feat: implement complex sqrt

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

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))))
    })
 })

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')

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
+}))

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))
+})

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))
+   })
 })

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))
 })

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])

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)
+})