feat(polynomialRoot) (#57)
Implements a simply polynomial root finder function polynomialRoot, intended to be used for benchmarking against mathjs. For this purpose, adds numerous other functions (e.g. cbrt, arg, cis), refactors sqrt (so that you can definitely get the complex square root when you want it), and makes numerous enhancements to the core so that a template can match after conversions. Co-authored-by: Glen Whitney <glen@studioinfinity.org> Reviewed-on: #57
This commit is contained in:
parent
31add66f4c
commit
0dbb95bbbe
18 changed files with 634 additions and 264 deletions
|
|
@ -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
|
||||
|
|
|
|||
7
src/complex/arg.mjs
Normal file
7
src/complex/arg.mjs
Normal file
|
|
@ -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))
|
||||
}
|
||||
28
src/complex/cbrtc.mjs
Normal file
28
src/complex/cbrtc.mjs
Normal file
|
|
@ -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))
|
||||
)
|
||||
})
|
||||
}
|
||||
9
src/complex/cis.mjs
Normal file
9
src/complex/cis.mjs
Normal file
|
|
@ -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))
|
||||
)
|
||||
}
|
||||
7
src/complex/isReal.mjs
Normal file
7
src/complex/isReal.mjs
Normal file
|
|
@ -0,0 +1,7 @@
|
|||
import Returns from '../core/Returns.mjs'
|
||||
export * from './Types/Complex.mjs'
|
||||
|
||||
export const isReal = {
|
||||
'Complex<T>': ({'equal(T,T)': eq, 'add(T,T)': plus}) => Returns(
|
||||
'boolean', z => eq(z.re, plus(z.re, z.im)))
|
||||
}
|
||||
|
|
@ -3,18 +3,24 @@ 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'
|
||||
export {gcd} from './gcd.mjs'
|
||||
export {invert} from './invert.mjs'
|
||||
export {isReal} from './isReal.mjs'
|
||||
export {isZero} from './isZero.mjs'
|
||||
export {multiply} from './multiply.mjs'
|
||||
export {negate} from './negate.mjs'
|
||||
export {polynomialRoot} from './polynomialRoot.mjs'
|
||||
export {quaternion} from './quaternion.mjs'
|
||||
export {quotient} from './quotient.mjs'
|
||||
export {roundquotient} from './roundquotient.mjs'
|
||||
export {sqrt} from './sqrt.mjs'
|
||||
export {sqrtc} from './sqrtc.mjs'
|
||||
export {zero} from './zero.mjs'
|
||||
|
||||
|
|
|
|||
118
src/complex/polynomialRoot.mjs
Normal file
118
src/complex/polynomialRoot.mjs
Normal file
|
|
@ -0,0 +1,118 @@
|
|||
import Returns from '../core/Returns.mjs'
|
||||
export * from './Types/Complex.mjs'
|
||||
|
||||
export const polynomialRoot = {
|
||||
'Complex<T>,...Complex<T>': ({
|
||||
T,
|
||||
'tuple(...Complex<T>)': tupCplx,
|
||||
'tuple(...T)': tupReal,
|
||||
'isZero(Complex<T>)': zero,
|
||||
'complex(T)': C,
|
||||
'multiply(Complex<T>,Complex<T>)': mul,
|
||||
'divide(Complex<T>,Complex<T>)': div,
|
||||
'negate(Complex<T>)': neg,
|
||||
'isReal(Complex<T>)': real,
|
||||
'equalTT(Complex<T>,Complex<T>)': eq,
|
||||
'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 => {
|
||||
if (arr.every(real)) {
|
||||
return tupReal.apply(tupReal, arr.map(z => z.re))
|
||||
}
|
||||
return tupCplx.apply(tupCplx, arr)
|
||||
}
|
||||
|
||||
const coeffs = [constant, ...rest]
|
||||
while (coeffs.length > 0 && zero(coeffs[coeffs.length - 1])) {
|
||||
coeffs.pop()
|
||||
}
|
||||
if (coeffs.length < 2) {
|
||||
}
|
||||
switch (coeffs.length) {
|
||||
case 0: case 1:
|
||||
throw new RangeError(
|
||||
`Polynomial [${constant}, ${rest}] must have at least one`
|
||||
+ 'non-zero non-constant coefficient')
|
||||
case 2: // linear
|
||||
return typedTup([neg(div(coeffs[0], coeffs[1]))])
|
||||
case 3: { // quadratic
|
||||
const [c, b, a] = coeffs
|
||||
const denom = mul(C(2), a)
|
||||
const d1 = mul(b, b)
|
||||
const d2 = mul(C(4), mul(a, c))
|
||||
if (eq(d1, d2)) {
|
||||
return typedTup([div(neg(b), denom)])
|
||||
}
|
||||
let discriminant = sqt(sub(d1, d2))
|
||||
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, '
|
||||
+ `not ${JSON.stringify(coeffs)}`)
|
||||
}
|
||||
})
|
||||
}
|
||||
|
|
@ -4,49 +4,30 @@ export * from './Types/Complex.mjs'
|
|||
export const sqrt = {
|
||||
'Complex<T>': ({
|
||||
config,
|
||||
'sqrtc(Complex<T>)': predictableSqrt,
|
||||
'isZero(T)': isZ,
|
||||
'sign(T)': sgn,
|
||||
'one(T)': uno,
|
||||
'add(T,T)': plus,
|
||||
'complex(T)': cplxU,
|
||||
'complex(T,T)': cplxB,
|
||||
'multiply(T,T)': mult,
|
||||
'self(T)': me,
|
||||
'divide(T,T)': div,
|
||||
'absquare(Complex<T>)': absqC,
|
||||
'subtract(T,T)': sub
|
||||
}) => {
|
||||
let baseReturns = returnTypeOf(me)
|
||||
if (baseReturns.includes('|')) {
|
||||
// Bit of a hack, because it is relying on other implementations
|
||||
// to list the "typical" value of sqrt first
|
||||
baseReturns = baseReturns.split('|', 1)[0]
|
||||
}
|
||||
|
||||
if (config.checkingDependency) return undefined
|
||||
const complexReturns = returnTypeOf(predictableSqrt)
|
||||
const baseReturns = complexReturns.slice(8, -1); // Complex<WhatWeWant>
|
||||
if (config.predictable) {
|
||||
return Returns(`Complex<${baseReturns}>`, z => {
|
||||
const reOne = uno(z.re)
|
||||
if (isZ(z.im) && sgn(z.re) === reOne) return cplxU(me(z.re))
|
||||
const reTwo = plus(reOne, reOne)
|
||||
const myabs = me(absqC(z))
|
||||
return cplxB(
|
||||
mult(sgn(z.im), me(div(plus(myabs, z.re), reTwo))),
|
||||
me(div(sub(myabs, z.re), reTwo))
|
||||
)
|
||||
})
|
||||
return Returns(complexReturns, z => predictableSqrt(z))
|
||||
}
|
||||
|
||||
return Returns(
|
||||
`Complex<${baseReturns}>|${baseReturns}|undefined`,
|
||||
z => {
|
||||
const reOne = uno(z.re)
|
||||
if (isZ(z.im) && sgn(z.re) === reOne) return me(z.re)
|
||||
const reTwo = plus(reOne, reOne)
|
||||
const myabs = me(absqC(z))
|
||||
const reSqrt = me(div(plus(myabs, z.re), reTwo))
|
||||
const imSqrt = me(div(sub(myabs, z.re), reTwo))
|
||||
if (reSqrt === undefined || imSqrt === undefined) return undefined
|
||||
return cplxB(mult(sgn(z.im), reSqrt), imSqrt)
|
||||
let complexSqrt
|
||||
try {
|
||||
complexSqrt = predictableSqrt(z)
|
||||
} catch (e) {
|
||||
return undefined
|
||||
}
|
||||
if (complexSqrt.re === undefined || complexSqrt.im === undefined) {
|
||||
return undefined
|
||||
}
|
||||
if (isZ(complexSqrt.im)) return complexSqrt.re
|
||||
return complexSqrt
|
||||
}
|
||||
)
|
||||
}
|
||||
|
|
|
|||
41
src/complex/sqrtc.mjs
Normal file
41
src/complex/sqrtc.mjs
Normal file
|
|
@ -0,0 +1,41 @@
|
|||
import {Returns, returnTypeOf} from '../core/Returns.mjs'
|
||||
export * from './Types/Complex.mjs'
|
||||
|
||||
export const sqrtc = {
|
||||
'Complex<T>': ({
|
||||
'isZero(T)': isZ,
|
||||
'sign(T)': sgn,
|
||||
'one(T)': uno,
|
||||
'add(T,T)': plus,
|
||||
'complex(T)': cplxU,
|
||||
'complex(T,T)': cplxB,
|
||||
'multiply(T,T)': mult,
|
||||
'sqrt(T)': sqt,
|
||||
'divide(T,T)': div,
|
||||
'absquare(Complex<T>)': absqC,
|
||||
'subtract(T,T)': sub
|
||||
}) => {
|
||||
if (isZ.checkingDependency) return undefined
|
||||
let baseReturns = returnTypeOf(sqt)
|
||||
if (baseReturns.includes('|')) {
|
||||
// Bit of a hack, because it is relying on other implementations
|
||||
// to list the "typical" value of sqrt first
|
||||
baseReturns = baseReturns.split('|', 1)[0]
|
||||
}
|
||||
return Returns(`Complex<${baseReturns}>`, z => {
|
||||
const reOne = uno(z.re)
|
||||
if (isZ(z.im) && sgn(z.re) === reOne) return cplxU(sqt(z.re))
|
||||
const myabs = sqt(absqC(z))
|
||||
const reTwo = plus(reOne, reOne)
|
||||
const reQuot = div(plus(myabs, z.re), reTwo)
|
||||
const imQuot = div(sub(myabs, z.re), reTwo)
|
||||
if (reQuot === undefined || imQuot === undefined) {
|
||||
throw new TypeError(`Cannot compute sqrt of ${z.re} + {z.im}i`)
|
||||
}
|
||||
return cplxB(
|
||||
mult(sgn(z.im), sqt(div(plus(myabs, z.re), reTwo))),
|
||||
sqt(div(sub(myabs, z.re), reTwo))
|
||||
)
|
||||
})
|
||||
}
|
||||
}
|
||||
Loading…
Add table
Add a link
Reference in a new issue