feat: add cbrt for complex numbers and quaternions, returning three values

src/complex/__test__/arithmetic.spec.js

@@ -132,4 +132,15 @@ describe('complex arithmetic operations', () => {
       assert.deepStrictEqual(
          sub(cplx(z,z), cplx(10,20)), cplx(cplx(-7, 4), cplx(-17, 4)))
    })
+   it('cube roots complex numbers and quaternions', () => {
+      assert(math.equal(math.cbrt(cplx(0, 8))[0], cplx(Math.sqrt(3), 1)))
+      assert.deepStrictEqual(
+         math.cbrt(cplx(2, 3)),
+         [cplx(1.4518566183526649, 0.49340353410400467),
+            cplx(-1.1532283040274218, 1.0106429470939742),
+            cplx(-0.29862831432524256, -1.5040464811979786)])
+      const quat = cplx(cplx(1, 2), cplx(2, 1))
+      const root = math.cbrt(quat)[0]
+      assert(math.equal(math.multiply(root, root, root), quat))
+   })
 })

src/complex/__test__/type.spec.js

@@ -19,6 +19,13 @@ describe('complex type operations', () => {
       assert.strictEqual(math.arg(cplx(1, Math.sqrt(3))), Math.PI/3)
       assert.strictEqual(math.arg(cplx(true, true)), Math.PI/4)
    })
+   it('calculates the real parts of complex numbers and quaternions', () => {
+      assert.strictEqual(math.re(cplx(2, -3)), 2)
+      assert.strictEqual(math.re(cplx(cplx(0.5), cplx(8, -7))), 0.5)
+   })
+   it('calculates the argument of a quaternion', () => {
+      assert(math.equal(math.arg(cplx(cplx(1, 1), cplx(1, 1))), Math.acos(0.5)))
+   })
    it('detects associates of a complex number', () => {
       const z = cplx(3, 4)
       const assoc = math.associate

src/complex/arithmetic.js

@@ -4,6 +4,7 @@ import {ResolutionError} from '#core/helpers.js'
 import {Returns, ReturnType, ReturnTyping} from '#core/Type.js'
 import {match} from '#core/TypePatterns.js'
 import {ReturnsAs} from '#generic/helpers.js'
+import {NumberT} from '#number/NumberT.js'
 
 const {conservative, full, free} = ReturnTyping
 
@@ -141,4 +142,44 @@ export const sqrt = match(Complex, (math, C, strategy) => {
    return ReturnsAs(prune, z => prune(sqrtImp(z)))
 })
 
+const TAU3 = 2 * Math.PI / 3
+
+export const cbrt = match(Complex, (math, C) => {
+   const arg = math.arg.resolve(C)
+   const divArg = math.divide.resolve([ReturnType(arg), NumberT])
+   const absC = math.abs.resolve(C)
+   const cbrtR = math.cbrt.resolve(ReturnType(absC), conservative)
+   const im = math.im.resolve(C)
+   const absIm = math.abs.resolve(ReturnType(im))
+   const divIm = math.divide.resolve([ReturnType(im), ReturnType(absIm)])
+   // TODO: replace with nanomath cos and sin when available
+   // const cos = math.cos.resolve(ReturnType(divArg))
+   // const CosType = ReturnType(cos)  // and similarly for sin
+   const cos = Math.cos
+   const CosType = NumberT
+   const sin = Math.sin
+   const SinType = NumberT
+   const mulIm = math.multiply.resolve([ReturnType(divIm), SinType])
+   const addUp = math.add.resolve([CosType, ReturnType(mulIm)])
+   const addAngle = math.add.resolve([ReturnType(divArg), math.typeOf(TAU3)])
+   const subAngle = math.subtract.resolve(
+      [ReturnType(divArg), math.typeOf(TAU3)])
+   const scale = math.multiply.resolve([ReturnType(cbrtR), ReturnType(addUp)])
+   const Cout = ReturnType(scale)
+   const vec = math.vector.resolve([Cout, Cout, Cout])
+   return ReturnsAs(vec, z => {
+      const arg3 = divArg(arg(z), 3)
+      const mag = absC(z)
+      const r = cbrtR(mag)
+      const imz = im(z)
+      const unit = divIm(imz, absIm(imz))
+      // At this point, z = mag*(cos(arg) + unit * sin(arg))
+      // so one cube root is r*(cos(arg3) + unit * sin(arg3))
+      // and we get the other two by adjusting arg3 by +- TAU3
+      const cus = theta => scale(r, addUp(cos(theta), mulIm(unit, sin(theta))))
+      return vec(
+         cus(arg3), cus(addAngle(arg3, TAU3)), cus(subAngle(arg3, TAU3)))
+   })
+})
+
 export const subtract = promoteBinary('subtract')

src/complex/type.js

@@ -31,20 +31,26 @@ export const complex = [
    })
 ]
 
-export const arg = // [  // enable when we have atan2 in mathjs
-//   match(Complex, (math, C) => {
-//      const re = math.re.resolve(C)
-//      const R = ReturnType(re)
-//      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)))
-//]
+export const arg = match(Complex, (math, C) => {
+   const re = math.re.resolve(C)
+   const atan2 = Math.atan2
+   const ArgType = NumberT
+   if (ReturnType(re) === C.Component) {
+      // "plain" Complex (not nested, i.e. not quaternions)
+      // FIXME: use math.atan2 once available
+      // const atan2 = math.atan2.resolve(
+      //   [C.Component, C.Component], conservative)
+      // const argType = ReturnType(atan2)
+      return Returns(ArgType, z => atan2(z.im, z.re))
+   }
+   const im = math.im.resolve(C)
+   const abs = math.abs.resolve(ReturnType(im))
+   // FIXME: use math.atan2 once available
+   // const atan2 = math.atan2.resolve(
+   //    [ReturnType(abs), ReturnType(re)], conservative)
+   // const ArgType = ReturnType(atan2)
+   return Returns(ArgType, z => atan2(abs(im(z)), re(z)))
+})
 
 /* Returns true if w is z multiplied by a complex unit */
 export const associate = match(