From d76eddc7a57e28ba98be427ef638e560e478d6e4 Mon Sep 17 00:00:00 2001 From: Glen Whitney Date: Thu, 1 Dec 2022 12:40:05 -0500 Subject: [PATCH] feat(polynomialRoot): Complete implementation --- src/complex/abs.mjs | 1 - src/complex/arg.mjs | 7 ++++ src/complex/cbrtc.mjs | 28 +++++++++++++ src/complex/cis.mjs | 9 +++++ src/complex/native.mjs | 3 ++ src/complex/polynomialRoot.mjs | 67 ++++++++++++++++++++++++++++---- src/complex/sqrtc.mjs | 1 - src/number/cbrt.mjs | 19 +++++++++ src/number/native.mjs | 1 + test/complex/_polynomialRoot.mjs | 33 +++++++++++++++- tools/approx.mjs | 46 ++++++++++++++++++++++ 11 files changed, 204 insertions(+), 11 deletions(-) create mode 100644 src/complex/arg.mjs create mode 100644 src/complex/cbrtc.mjs create mode 100644 src/complex/cis.mjs create mode 100644 src/number/cbrt.mjs create mode 100644 tools/approx.mjs diff --git a/src/complex/abs.mjs b/src/complex/abs.mjs index 53a2464..3ea0274 100644 --- a/src/complex/abs.mjs +++ b/src/complex/abs.mjs @@ -3,7 +3,6 @@ export * from './Types/Complex.mjs' export const abs = { 'Complex': ({ - T, sqrt, // Unfortunately no notation yet for the needed signature 'absquare(T)': baseabsq, 'absquare(Complex)': absq diff --git a/src/complex/arg.mjs b/src/complex/arg.mjs new file mode 100644 index 0000000..d654795 --- /dev/null +++ b/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': () => Returns('number', z => Math.atan2(z.im, z.re)) +} diff --git a/src/complex/cbrtc.mjs b/src/complex/cbrtc.mjs new file mode 100644 index 0000000..118da60 --- /dev/null +++ b/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': ({ + 'arg(T)': theta, + 'divide(T,T)': div, + 'abs(Complex)': absval, + 'complex(T)': cplx, + 'cbrt(T)': cbrtT, + 'multiply(Complex,Complex)': mult, + 'cis(T)': cisT, + 'tuple(...Complex)': tup + }) => Returns('Tuple>', 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)) + ) + }) +} diff --git a/src/complex/cis.mjs b/src/complex/cis.mjs new file mode 100644 index 0000000..fd541e7 --- /dev/null +++ b/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', t => cplx(Math.cos(t), Math.sin(t)) + ) +} diff --git a/src/complex/native.mjs b/src/complex/native.mjs index 1668743..eba3859 100644 --- a/src/complex/native.mjs +++ b/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' diff --git a/src/complex/polynomialRoot.mjs b/src/complex/polynomialRoot.mjs index fc7c255..7a6b9a3 100644 --- a/src/complex/polynomialRoot.mjs +++ b/src/complex/polynomialRoot.mjs @@ -12,9 +12,11 @@ export const polynomialRoot = { 'divide(Complex,Complex)': div, 'negate(Complex)': neg, 'isReal(Complex)': real, - 'equalTT(Complex, Complex)': eq, - 'subtract(Complex, Complex)': sub, + 'equalTT(Complex,Complex)': eq, + 'add(Complex,Complex)': plus, + 'subtract(Complex,Complex)': sub, 'sqrtc(Complex)': sqt, + 'cbrtc(Complex)': cbt }) => Returns(`Tuple<${T}>|Tuple>`, (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, ' diff --git a/src/complex/sqrtc.mjs b/src/complex/sqrtc.mjs index 3a7ea1d..e909ff7 100644 --- a/src/complex/sqrtc.mjs +++ b/src/complex/sqrtc.mjs @@ -3,7 +3,6 @@ export * from './Types/Complex.mjs' export const sqrtc = { 'Complex': ({ - T, 'isZero(T)': isZ, 'sign(T)': sgn, 'one(T)': uno, diff --git a/src/number/cbrt.mjs b/src/number/cbrt.mjs new file mode 100644 index 0000000..af7587b --- /dev/null +++ b/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 + }) +} + diff --git a/src/number/native.mjs b/src/number/native.mjs index fcebecc..ad2de12 100644 --- a/src/number/native.mjs +++ b/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') diff --git a/test/complex/_polynomialRoot.mjs b/test/complex/_polynomialRoot.mjs index 1224bcf..22ad90c 100644 --- a/test/complex/_polynomialRoot.mjs +++ b/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))) + }) }) diff --git a/tools/approx.mjs b/tools/approx.mjs new file mode 100644 index 0000000..cab5483 --- /dev/null +++ b/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) +}