feat(polynomialRoot): Complete implementation

src/complex/abs.mjs

@@ -3,7 +3,6 @@ export * from './Types/Complex.mjs'
 
 export const abs = {
    'Complex<T>': ({
-      T,
       sqrt, // Unfortunately no notation yet for the needed signature
       'absquare(T)': baseabsq,
       'absquare(Complex<T>)': absq

src/complex/arg.mjs

@@ -0,0 +1,7 @@
+import {Returns, returnTypeOf} from '../core/Returns.mjs'
+export * from './Types/Complex.mjs'
+
+/* arg is the "argument" or angle theta of z in its form r cis theta */
+export const arg = {
+   'Complex<number>': () => Returns('number', z => Math.atan2(z.im, z.re))
+}

src/complex/cbrtc.mjs

@@ -0,0 +1,28 @@
+import {Returns, returnTypeOf} from '../core/Returns.mjs'
+export * from './Types/Complex.mjs'
+
+const TAU3 = 2 * Math.PI / 3
+
+/* Complex cube root that returns all three roots as a tuple of complex. */
+/* follows the implementation in mathjs */
+/* Really only works for T = number at the moment because of arg and cbrt */
+export const cbrtc = {
+   'Complex<T>': ({
+      'arg(T)': theta,
+      'divide(T,T)': div,
+      'abs(Complex<T>)': absval,
+      'complex(T)': cplx,
+      'cbrt(T)': cbrtT,
+      'multiply(Complex<T>,Complex<T>)': mult,
+      'cis(T)': cisT,
+      'tuple(...Complex<T>)': tup
+   }) => Returns('Tuple<Complex<T>>', z => {
+      const arg3 = div(theta(z), 3)
+      const r = cplx(cbrtT(absval(z)))
+      return tup(
+         mult(r, cisT(arg3)),
+         mult(r, cisT(arg3 + TAU3)),
+         mult(r, cisT(arg3 - TAU3))
+      )
+   })
+}

src/complex/cis.mjs

@@ -0,0 +1,9 @@
+import {Returns, returnTypeOf} from '../core/Returns.mjs'
+export * from './Types/Complex.mjs'
+
+/* Returns cosine plus i sin theta */
+export const cis = {
+   'number': ({'complex(number,number)': cplx}) => Returns(
+      'Complex<number>', t => cplx(Math.cos(t), Math.sin(t))
+   )
+}

src/complex/native.mjs

@@ -3,7 +3,10 @@ export * from './Types/Complex.mjs'
 export {abs} from './abs.mjs'
 export {absquare} from './absquare.mjs'
 export {add} from './add.mjs'
+export {arg} from './arg.mjs'
 export {associate} from './associate.mjs'
+export {cbrtc} from './cbrtc.mjs'
+export {cis} from './cis.mjs'
 export {complex} from './complex.mjs'
 export {conjugate} from './conjugate.mjs'
 export {equalTT} from './equalTT.mjs'

src/complex/polynomialRoot.mjs

@@ -12,9 +12,11 @@ export const polynomialRoot = {
       'divide(Complex<T>,Complex<T>)': div,
       'negate(Complex<T>)': neg,
       'isReal(Complex<T>)': real,
-      'equalTT(Complex<T>, Complex<T>)': eq,
+      'equalTT(Complex<T>,Complex<T>)': eq,
-      'subtract(Complex<T>, Complex<T>)': sub,
+      'add(Complex<T>,Complex<T>)': plus,
+      'subtract(Complex<T>,Complex<T>)': sub,
       'sqrtc(Complex<T>)': sqt,
+      'cbrtc(Complex<T>)': cbt
    }) => Returns(`Tuple<${T}>|Tuple<Complex<${T}>>`, (constant, rest) => {
       // helper to convert results to appropriate tuple type
       const typedTup = arr => {
@@ -39,23 +41,74 @@ export const polynomialRoot = {
             return typedTup([neg(div(coeffs[0], coeffs[1]))])
          case 3: { // quadratic
             const [c, b, a] = coeffs
-            console.log('solving', a, b, c)
             const denom = mul(C(2), a)
             const d1 = mul(b, b)
             const d2 = mul(C(4), mul(a, c))
-            console.log('Whoa', denom, d1, d2)
             if (eq(d1, d2)) {
-               console.log('Hello', b, denom, div(neg(b), denom))
                return typedTup([div(neg(b), denom)])
             }
             let discriminant = sqt(sub(d1, d2))
-            console.log('Uhoh', discriminant)
-            console.log('Roots', div(sub(discriminant, b), denom), div(sub(neg(discriminant), b), denom))
             return typedTup([
                div(sub(discriminant, b), denom),
                div(sub(neg(discriminant), b), denom)
             ])
          }
+         case 4: { // cubic, cf. https://en.wikipedia.org/wiki/Cubic_equation
+            const [d, c, b, a] = coeffs
+            const denom = neg(mul(C(3), a))
+            const asqrd = mul(a, a)
+            const D0_1 = mul(b, b)
+            const bcubed = mul(D0_1, b)
+            const D0_2 = mul(C(3), mul(a, c))
+            const D1_1 = plus(
+               mul(C(2), bcubed), mul(C(27), mul(asqrd, d)))
+            const abc = mul(a, mul(b, c))
+            const D1_2 = mul(C(9), abc)
+            // Check for a triple root
+            if (eq(D0_1, D0_2) && eq(D1_1, D1_2)) {
+               return typedTup([div(b, denom)])
+            }
+            const Delta0 = sub(D0_1, D0_2)
+            const Delta1 = sub(D1_1, D1_2)
+            const csqrd = mul(c, c)
+            const discriminant1 = plus(
+               mul(C(18), mul(abc, d)), mul(D0_1, csqrd))
+            const discriminant2 = plus(
+               mul(C(4), mul(bcubed, d)),
+               plus(
+                  mul(C(4), mul(a, mul(csqrd, c))),
+                  mul(C(27), mul(asqrd, mul(d, d)))))
+            // See if we have a double root
+            if (eq(discriminant1, discriminant2)) {
+               return typedTup([
+                  div(
+                     sub(
+                        mul(C(4), abc),
+                        plus(mul(C(9), mul(asqrd, d)), bcubed)),
+                     mul(a, Delta0)), // simple root
+                  div(
+                     sub(mul(C(9), mul(a, d)), mul(b, c)),
+                     mul(C(2), Delta0)) // double root
+               ])
+            }
+            // OK, we have three distinct roots
+            let Ccubed
+            if (eq(D0_1, D0_2)) {
+               Ccubed = Delta1
+            } else {
+               Ccubed = div(
+                  plus(
+                     Delta1,
+                     sqt(sub(
+                        mul(Delta1, Delta1),
+                        mul(C(4), mul(Delta0, mul(Delta0, Delta0)))))
+                  ),
+                  C(2))
+            }
+            const croots = cbt(Ccubed)
+            return typedTup(cbt(Ccubed).elts.map(
+               C => div(plus(b, plus(C, div(Delta0, C))), denom)))
+         }
          default:
             throw new RangeError(
                'only implemented for cubic or lower-order polynomials, '

src/complex/sqrtc.mjs

@@ -3,7 +3,6 @@ export * from './Types/Complex.mjs'
 
 export const sqrtc = {
    'Complex<T>': ({
-      T,
       'isZero(T)': isZ,
       'sign(T)': sgn,
       'one(T)': uno,

src/number/cbrt.mjs

@@ -0,0 +1,19 @@
+import Returns from '../core/Returns.mjs'
+export * from './Types/number.mjs'
+
+/* Returns just the real cube root, following mathjs implementation */
+export const cbrt = {
+   number: ({'negate(number)': neg}) => Returns('number', x => {
+      if (x === 0) return x
+      const negate = x < 0
+      if (negate) x = neg(x)
+      let result = x
+      if (isFinite(x)) {
+         result = Math.exp(Math.log(x) / 3)
+         result = (x / (result * result) + (2 * result)) / 3
+      }
+      if (negate) return neg(result)
+      return result
+   })
+}
+                                                

src/number/native.mjs

@@ -6,6 +6,7 @@ export * from './Types/number.mjs'
 export {abs} from './abs.mjs'
 export {absquare} from './absquare.mjs'
 export {add} from './add.mjs'
+export {cbrt} from './cbrt.mjs'
 export {compare} from './compare.mjs'
 export const conjugate = {'T:number': identitySubTypes('number')}
 export const gcd = gcdType('NumInt')

test/complex/_polynomialRoot.mjs

@@ -1,4 +1,5 @@
 import assert from 'assert'
+import * as approx from '../../tools/approx.mjs'
 import math from '../../src/pocomath.mjs'
 
 describe('polynomialRoot', () => {
@@ -29,6 +30,34 @@ describe('polynomialRoot', () => {
       assert.deepEqual(
          pRoot(complex(3, 1), -3, 1), tup(complex(1, 1), complex(2, -1)))
    })
-
+   it('should solve a cubic with a triple root', function () {
-
+      assert.deepEqual(pRoot(8, 12, 6, 1), tup(-2))
+      assert.deepEqual(
+         pRoot(complex(-2, 11), complex(9, -12), complex(-6, 3), 1),
+         tup(complex(2, -1)))
+   })
+   it('should solve a cubic with one simple and one double root', function () {
+      assert.deepEqual(pRoot(4, 0, -3, 1), tup(-1, 2))
+      assert.deepEqual(
+         pRoot(complex(9, 9), complex(15, 6), complex(7, 1), 1),
+         tup(complex(-1, -1), -3))
+      assert.deepEqual(
+         pRoot(complex(0, 6), complex(6, 8), complex(5, 2), 1),
+         tup(-3, complex(-1, -1)))
+      assert.deepEqual(
+         pRoot(complex(2, 6), complex(8, 6), complex(5, 1), 1),
+         tup(complex(-3, 1), complex(-1, -1)))
+   })
+   it('should solve a cubic with three distinct roots', function () {
+      approx.deepEqual(pRoot(6, 11, 6, 1), tup(-3, -1, -2))
+      approx.deepEqual(
+         pRoot(-1, -2, 0, 1),
+         tup(-1, (1 + math.sqrt(5)) / 2, (1 - math.sqrt(5)) / 2))
+      approx.deepEqual(
+         pRoot(1, 1, 1, 1),
+         tup(-1, complex(0, -1), complex(0, 1)))
+      approx.deepEqual(
+         pRoot(complex(0, -10), complex(8, 12), complex(-6, -3), 1),
+         tup(complex(1, 1), complex(3, 1), complex(2, 1)))
+   })
 })

tools/approx.mjs

@@ -0,0 +1,46 @@
+import assert from 'assert'
+
+export const epsilon = 1e-12
+
+const isNumber = entity => (typeof entity === 'number')
+
+export function equal(a, b) {
+   if (isNumber(a) && isNumber(b)) {
+      if (a === b) return true
+      if (isNaN(a)) return assert.strictEqual(a.toString(), b.toString())
+      const message = `${a} ~= ${b} (to ${epsilon})`
+      if (a === 0) return assert.ok(Math.abs(b) < epsilon, message)
+      if (b === 0) return assert.ok(Math.abs(a) < epsilon, message)
+      const diff = Math.abs(a - b)
+      const maxDiff = Math.abs(epsilon * Math.max(Math.abs(a), Math.abs(b)))
+      return assert.ok(diff <= maxDiff, message)
+   }
+   return assert.strictEqual(a, b)
+}
+
+export function deepEqual(a, b) {
+   if (Array.isArray(a) && Array.isArray(b)) {
+      const alen = a.length
+      assert.strictEqual(alen, b.length, `${a} ~= ${b}`)
+      for (let i = 0; i < alen; ++i) deepEqual(a[i], b[i])
+      return true
+   }
+   if (typeof a === 'object' && typeof b === 'object') {
+      for (const prop in a) {
+         if (a.hasOwnProperty(prop)) {
+            assert.ok(
+               b.hasOwnProperty(prop), `a[${prop}] = ${a[prop]} ~= ${b[prop]}`)
+            deepEqual(a[prop], b[prop])
+         }
+      }
+
+      for (const prop in b) {
+         if (b.hasOwnProperty(prop)) {
+            assert.ok(
+               a.hasOwnProperty(prop), `${a[prop]} ~= ${b[prop]} = b[${prop}]`)
+         }
+      }
+      return true
+   }
+   return equal(a, b)
+}