diff --git a/src/complex/__test__/arithmetic.spec.js b/src/complex/__test__/arithmetic.spec.js index 536277a..4a54f42 100644 --- a/src/complex/__test__/arithmetic.spec.js +++ b/src/complex/__test__/arithmetic.spec.js @@ -32,6 +32,15 @@ describe('complex arithmetic operations', () => { const addFull = math.add.resolve([CplxNum, CplxNum], full) assert.deepStrictEqual(addFull(z, math.conj(z)), cplx(6, 0)) }) + it('adds real and complex numbers', () => { + const z = cplx(3, 4) + const zp3 = cplx(6,4) + const add = math.add + assert.deepStrictEqual(add(z, 3), zp3) + assert.deepStrictEqual(add(3, z), zp3) + assert.deepStrictEqual(add(3, cplx(z, z)), cplx(zp3, z)) + assert.deepStrictEqual(add(cplx(z, z), 3), cplx(zp3, z)) + }) it('conjugates complex numbers', () => { const conj = math.conj const z = cplx(3, 4) @@ -75,6 +84,25 @@ describe('complex arithmetic operations', () => { mult(q0, cplx(cplx(2, 0.1), cplx(1, 0.1))), cplx(cplx(1.9, 2.1), cplx(1.1, -0.9))) }) + it('takes the square roots of complex numbers', () => { + const {sqrt, multiply} = math + const rhalf = Math.sqrt(1 / 2) + assert.deepStrictEqual(sqrt(cplx(1, 0)), 1) + assert.deepStrictEqual(sqrt(cplx(-1, 0)), cplx(0, 1)) + assert(math.equal(sqrt(cplx(0, 1)), cplx(rhalf, rhalf))) + assert(math.equal(sqrt(cplx(0, -1)), cplx(rhalf, -rhalf))) + assert.deepStrictEqual(sqrt(cplx(5, 12)), cplx(3, 2)) + const z = cplx(3, 4) + const rz = sqrt(z) + assert.deepStrictEqual(multiply(rz, rz), z) + // quaternions, too: + assert.deepStrictEqual( + sqrt(cplx(cplx(-46, 20), cplx(12, 16))), + cplx(cplx(2, 5), cplx(3, 4))) + const q = cplx(z, z) + const rq = sqrt(q) + assert(math.equal(multiply(rq, rq), q)) + }) it('subtracts complex numbers', () => { const z = cplx(3, 4) const sub = math.subtract diff --git a/src/complex/__test__/utils.spec.js b/src/complex/__test__/utils.spec.js index 9b88c75..d7f2469 100644 --- a/src/complex/__test__/utils.spec.js +++ b/src/complex/__test__/utils.spec.js @@ -36,5 +36,16 @@ describe('complex utilities', () => { assert(isReal(cplx(5, 0))) assert(isReal(cplx(5, -1e-17))) assert(!isReal(cplx(5, 0.000001))) + assert(isReal(cplx(cplx(5, 1e-16), cplx(-0, 0)))) + assert(!isReal(cplx(cplx(5, 2), cplx(0, 0)))) + assert(!isReal(cplx(cplx(5, 0), cplx(0, 0.00002)))) + }) + it('identifies complex numbers that only have a real part', () => { + const noImag = math.nonImaginary + assert(noImag(cplx(5, 0))) + assert(noImag(cplx(5, -1e-17))) + assert(!noImag(cplx(5, 0.000001))) + assert(noImag(cplx(cplx(5, 2), cplx(0, 0)))) + assert(!noImag(cplx(cplx(5, 0), cplx(0, 0.00002)))) }) }) diff --git a/src/complex/arithmetic.js b/src/complex/arithmetic.js index 9b72697..f839c5b 100644 --- a/src/complex/arithmetic.js +++ b/src/complex/arithmetic.js @@ -1,11 +1,11 @@ import {Complex} from './Complex.js' import {maybeComplex, promoteBinary, promoteUnary} from './helpers.js' import {ResolutionError} from '#core/helpers.js' -import {ReturnTyping} from '#core/Type.js' +import {Returns, ReturnTyping} from '#core/Type.js' import {match} from '#core/TypePatterns.js' import {ReturnsAs} from '#generic/helpers.js' -const {full} = ReturnTyping +const {conservative, full, free} = ReturnTyping export const absquare = match(Complex, (math, C, strategy) => { const compAbsq = math.absquare.resolve(C.Component, full) @@ -50,23 +50,87 @@ export const invert = match(Complex, (math, C, strategy) => { // We want this to work for complex numbers, quaternions, octonions, etc // See https://math.ucr.edu/home/baez/octonions/node5.html -export const multiply = match([Complex, Complex], (math, [W, Z], strategy) => { - const conj = math.conj.resolve(W.Component, full) - if (conj.returns !== W.Component) { - throw new ResolutionError( - `conjugation on ${W.Component} returns other type (${conj.returns})`) - } - const mWZ = math.multiply.resolve([W.Component, Z.Component], full) - const mZW = math.multiply.resolve([Z.Component, W.Component], full) - const sub = math.subtract.resolve([mWZ.returns, mZW.returns], full) - const add = math.add.resolve([mWZ.returns, mZW.returns], full) - const cplx = maybeComplex(math, strategy, sub.returns, add.returns) - return ReturnsAs(cplx, (w, z) => { - const real = sub(mWZ( w.re, z.re), mZW(z.im, conj(w.im))) - const imag = add(mWZ(conj(w.re), z.im), mZW(z.re, w.im)) - return cplx(real, imag) +export const multiply = [ + match([T => !T.complex, Complex], (math, [R, C], strategy) => { + const mult = math.multiply.resolve([R, C.Component], full) + const cplx = maybeComplex(math, strategy, mult.returns, mult.returns) + return ReturnsAs(cplx, (r, z) => cplx(mult(r, z.re), mult(r, z.im))) + }), + match([Complex, T => !T.complex], (math, [C, R], strategy) => { + const mult = math.multiply.resolve([R, C], strategy) + return ReturnsAs(mult, (z, r) => mult(r, z)) + }), + match([Complex, Complex], (math, [W, Z], strategy) => { + const conj = math.conj.resolve(W.Component, full) + if (conj.returns !== W.Component) { + throw new ResolutionError( + `conjugation on ${W.Component} returns type (${conj.returns})`) + } + const mWZ = math.multiply.resolve([W.Component, Z.Component], full) + const mZW = math.multiply.resolve([Z.Component, W.Component], full) + const sub = math.subtract.resolve([mWZ.returns, mZW.returns], full) + const add = math.add.resolve([mWZ.returns, mZW.returns], full) + const cplx = maybeComplex(math, strategy, sub.returns, add.returns) + return ReturnsAs(cplx, (w, z) => { + const real = sub(mWZ( w.re, z.re), mZW(z.im, conj(w.im))) + const imag = add(mWZ(conj(w.re), z.im), mZW(z.re, w.im)) + return cplx(real, imag) + }) }) -}) +] export const negate = promoteUnary('negate') + +// Should work for complex, quaternions, octonions, etc., even with +// integer coordinates. +export const sqrt = match(Complex, (math, C, strategy) => { + const re = math.re.resolve(C) + const R = re.returns + const isReal = math.isReal.resolve(C) + // dependencies for the real case: + const zComp = math.zero(C.Component) + const sign = math.sign.resolve(R, conservative) + const oneR = math.one(R) + const zR = math.zero(R) + const addRComp = math.add.resolve([R, C.Component], full) + const sqrtR = math.sqrt.resolve(R, conservative) + const neg = math.negate.resolve(R, conservative) + const cplx = math.complex.resolve([C.Component, C.Component], full) + // additional dependencies for the complex case + const abs = math.abs.resolve(C, full) + if (abs.returns !== R) { + throw new TypeError(`abs on ${C} returns ${abs.returns}, not ${R}`) + } + const addRR = math.add.resolve([R, R], conservative) + const twoR = addRR(oneR, oneR) + const multRR = math.multiply.resolve([R, R], conservative) + const im = math.im.resolve(C, full) + const addRC = math.add.resolve([R, C], full) + const divRR = math.divide.resolve([R, R], conservative) + const divCR = math.divide.resolve([C, R], full) + + // The guts of the computation: + const sqrtImp = Returns(C, z => { + const a = re(z) + if (isReal(z)) { // always a special case + let rp = zComp + let ip = zComp + const sgn = sign(a) + if (sgn === oneR) rp = addRComp(sqrtR(a), rp) + else if (sgn !== zR) ip = addRComp(sqrtR(neg(a)), ip) + return cplx(rp, ip) + } + // Complex case: + // We can write z = a + q where q is pure imaginary. + // Let s = sqrt(2(|z|+a)). Then sqrt(z) = (s/2) + (q/s). + const s = sqrtR(multRR(twoR, addRR(abs(z), a))) + const q = im(z) + return addRC(divRR(s, twoR), divCR(q, s)) + }) + + if (strategy != free) return sqrtImp + const prune = math.pruneImaginary.resolve(C, free) + return ReturnsAs(prune, z => prune(sqrtImp(z))) +}) + export const subtract = promoteBinary('subtract') diff --git a/src/complex/type.js b/src/complex/type.js index cc3ee58..812ee06 100644 --- a/src/complex/type.js +++ b/src/complex/type.js @@ -1,5 +1,5 @@ import {Complex} from './Complex.js' -import {OneOf, Returns, ReturnTyping} from "#core/Type.js" +import {OneOf, Returns, ReturnTyping, TypeOfTypes} from "#core/Type.js" import {Any, match} from "#core/TypePatterns.js" import {BooleanT} from '#boolean/BooleanT.js' import {NumberT} from '#number/NumberT.js' @@ -29,8 +29,20 @@ export const complex = [ }) ] -export const arg = match( - Complex(NumberT), Returns(NumberT, z => Math.atan2(z.im, z.re))) +export const arg = // [ // enable when we have atan2 in mathjs +// match(Complex, (math, C) => { +// const re = math.re.resolve(C) +// const R = re.returns +// const im = math.im.resolve(C) +// const abs = math.abs.resolve(C) +// const atan2 = math.atan2.resolve([R, R], conservative) +// return Returns(R, z => atan2(abs(im(z)), re(z))) +// }), // note always between 0 and tau/2; need to use in conjunction +// // with a complex unit function that gives you the proper +// // imaginary unit, ±i in the simple complex case, to restore the +// // full circle of values for the direction of a complex number + match(Complex(NumberT), Returns(NumberT, z => Math.atan2(z.im, z.re))) +//] /* Returns true if w is z multiplied by a complex unit */ export const associate = match( @@ -82,10 +94,16 @@ export const pruneImaginary = match(Complex, (math, C) => { T = T.Component outcomes.push(T) } - const real = math.isReal.resolve(C, full) + const noImag = math.nonImaginary.resolve(C, full) if (C.Component.complex) { const compPrune = math.pruneImaginary.resolve(C.Component) - return Returns(OneOf(...outcomes), z => real(z) ? compPrune(z.re) : z) + return Returns(OneOf(...outcomes), z => noImag(z) ? compPrune(z.re) : z) } - return Returns(OneOf(...outcomes), z => real(z) ? z.re : z) + return Returns(OneOf(...outcomes), z => noImag(z) ? z.re : z) }) + +// Returns the type of re(z) +export const Real = match(TypeOfTypes, Returns(TypeOfTypes, t => { + while (t.complex) t = t.Component + return t +})) diff --git a/src/complex/utils.js b/src/complex/utils.js index 0aae100..04cdd45 100644 --- a/src/complex/utils.js +++ b/src/complex/utils.js @@ -1,7 +1,7 @@ import {Complex} from './Complex.js' import {promotePredicateAnd, promoteUnary} from './helpers.js' -import {ReturnTyping} from '#core/Type.js' +import {Returns, ReturnTyping} from '#core/Type.js' import {match} from '#core/TypePatterns.js' import {ReturnsAs} from '#generic/helpers.js' @@ -14,8 +14,41 @@ export const isfinite = promotePredicateAnd('isfinite') // just so-called rational integers. export const isInteger = promotePredicateAnd('isInteger') +// true if the Complex is truly a real number (nonImaginary and its +// real part is a real number) export const isReal = match(Complex, (math, C) => { + const nonImag = math.nonImaginary.resolve(C, full) + const realComp = math.isReal.resolve(C.Component) + return ReturnsAs(realComp, z => nonImag(z) && realComp(z.re)) +}) + +// true if the imaginary part of a Complex is negligible compared to the real +export const nonImaginary = match(Complex, (math, C) => { const eq = math.equal.resolve([C.Component, C.Component]) const add = math.add.resolve([C.Component, C.Component], full) return ReturnsAs(eq, z => eq(z.re, add(z.re, z.im))) }) + +// Always returns the "true" real part of a complex number z (i.e., the part +// that is actually in the real numbers mathematically). In other words, +// performs z.re recursively until it gets something not Complex. +export const re = match(Complex, (math, C) => { + return Returns(math.Real(C), z => { + let T = C + while (T.complex) { + T = T.Component + z = z.re + } + return z + }) +}) + +// Returns everything but the real part of z. NOTE this is a complex +// number with zero real part, _not_ the coefficient of the imaginary +// component of z. If you want the latter, just use `z.im` + +export const im = match(Complex, (math, C) => { + const imComp = math.im.resolve(C.Component, full) + const cplx = math.complex.resolve([C.Component, C.Component], full) + return ReturnsAs(cplx, z => cplx(imComp(z.re), z.im)) +}) diff --git a/src/generic/__test__/utils.spec.js b/src/generic/__test__/utils.spec.js index f8b6d53..41efea2 100644 --- a/src/generic/__test__/utils.spec.js +++ b/src/generic/__test__/utils.spec.js @@ -22,4 +22,7 @@ describe('generic utility functions', () => { assert(isReal(math.complex(-3.25, 4e-16))) assert(!isReal(math.complex(3, 4))) }) + it('tests for no imaginary part', () => { + assert(math.nonImaginary(true)) + }) }) diff --git a/src/generic/arithmetic.js b/src/generic/arithmetic.js index 5168111..d561912 100644 --- a/src/generic/arithmetic.js +++ b/src/generic/arithmetic.js @@ -1,8 +1,15 @@ -import {Returns} from '#core/Type.js' +import {ReturnsAs} from './helpers.js' + +import {Returns, ReturnTyping} from '#core/Type.js' import {match, Any} from '#core/TypePatterns.js' -export const conj = match(Any, (_math, T) => Returns(T, a => a)) +export const abs = match(Any, (math, T) => { + const absq = math.absquare.resolve(T) + const sqrt = math.sqrt.resolve(absq.returns, ReturnTyping.conservative) + return ReturnsAs(sqrt, t => sqrt(absq(t))) +}) +export const conj = match(Any, (_math, T) => Returns(T, t => t)) export const square = match(Any, (math, T, strategy) => { const mult = math.multiply.resolve([T, T], strategy) - return Returns(mult.returns, a => mult(a, a)) + return Returns(mult.returns, t => mult(t, t)) }) diff --git a/src/generic/relational.js b/src/generic/relational.js index e128cc0..2d70a1e 100644 --- a/src/generic/relational.js +++ b/src/generic/relational.js @@ -88,6 +88,12 @@ export const larger = match([Any, Any], (math, [T, U]) => { return boolnum((t, u) => !eq(t, u) && bigger(t, u))(math) }) +export const isPositive = match(Any, (math, T) => { + const zero = math.zero(T) + const larger = math.larger.resolve([T, T]) + return boolnum(t => larger(t, zero)) +}) + export const largerEq = match([Any, Any], (math, [T, U]) => { const eq = math.equal.resolve([T, U]) const bigger = math.exceeds.resolve([T, U]) diff --git a/src/generic/utils.js b/src/generic/utils.js index d696d6c..716609e 100644 --- a/src/generic/utils.js +++ b/src/generic/utils.js @@ -1,10 +1,15 @@ import {ReturnsAs} from './helpers.js' import {ResolutionError} from '#core/helpers.js' -import {Passthru, match} from '#core/TypePatterns.js' +import {Returns} from '#core/Type.js' +import {Any, Passthru, match} from '#core/TypePatterns.js' import {boolnum} from '#number/helpers.js' -// Most types are real. Have to make sure to redefine on all non-real types -export const isReal = match(Passthru, boolnum(() => true)) +// Most types are real, so we just define them that way generically. +// We have to make sure to redefine isReal on all non-real types. +// We use Any here so that it will match before a specific type matches +// with conversion; such a match runs the risk of not producing the correct +// result, or worse, leading to a resolution loop. +export const isReal = match(Any, boolnum(() => true)) export const isZero = match(Passthru, (math, [T]) => { if (!T) { // called with no arguments throw new ResolutionError('isZero() requires one argument') @@ -13,3 +18,9 @@ export const isZero = match(Passthru, (math, [T]) => { const eq = math.equal.resolve([T, T]) return ReturnsAs(eq, x => eq(z, x)) }) +export const nonImaginary = match(Passthru, boolnum(() => true)) +export const re = match(Any, (_math, T) => Returns(T, t => t)) +export const im = match(Any, (math, T) => { + const z = math.zero(T) + return Returns(T, () => z) +})