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/isReal.mjs

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

src/complex/native.mjs

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

src/complex/polynomialRoot.mjs

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

src/complex/sqrt.mjs

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

src/complex/sqrtc.mjs

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

src/core/PocomathInstance.mjs

@@ -7,6 +7,12 @@ import {typeListOfSignature, typesOfSignature, subsetOfKeys} from './utils.mjs'
 
 const anySpec = {} // fixed dummy specification of 'any' type
 
+/* Like `.some(predicate)` but for collections */
+function exists(collection, predicate) {
+   for (const item of collection) if (predicate(item)) return true;
+   return false;
+}
+
 /* Template/signature parsing stuff; should probably be moved to a
  * separate file, but it's a bit interleaved at the moment
  */
@@ -87,6 +93,7 @@ function substituteInSignature(signature, parameter, type) {
    return sig.replaceAll(pattern, type)
 }
 
+const UniversalType = 'ground' // name for a type that matches anything
 let lastWhatToDo = null // used in an infinite descent check
 
 export default class PocomathInstance {
@@ -97,6 +104,7 @@ export default class PocomathInstance {
    static reserved = new Set([
       'chain',
       'config',
+      'convert',
       'importDependencies',
       'install',
       'installType',
@@ -128,16 +136,30 @@ export default class PocomathInstance {
       // its onMismatch function, below:
       this._metaTyped = typed.create()
       this._metaTyped.clear()
+      this._metaTyped.addTypes([{name: UniversalType, test: () => true}])
+
       // And these are the meta bindings: (I think we don't need separate
       // invalidation for them as they are only accessed through a main call.)
       this._meta = {} // The resulting typed-functions
       this._metaTFimps = {} // and their implementations
       const me = this
-      const myTyped = this._typed
       this._typed.onMismatch = (name, args, sigs) => {
          if (me._invalid.has(name)) {
+            if (this._fixing === name) {
+               this._fixingCount += 1
+               if (this._fixingCount > this._maxDepthSeen + 2) {
+                  throw new ReferenceError(
+                     `Infinite descent rebuilding ${name} on ${args}`)
+               }
+            } else {
+               this._fixingCount = 0
+            }
             // rebuild implementation and try again
-            return me[name](...args)
+            const lastFixing = this._fixing
+            this._fixing = name
+            const value = me[name](...args)
+            this._fix = lastFixing
+            return value
          }
          const metaversion = me._meta[name]
          if (metaversion) {
@@ -183,6 +205,8 @@ export default class PocomathInstance {
       this._plainFunctions = new Set() // the names of the plain functions
       this._chainRepository = {} // place to store chainified functions
       this.joinTypes = this.joinTypes.bind(me)
+      // Provide access to typed function conversion:
+      this.convert = this._typed.convert.bind(this._typed)
    }
 
    /**
@@ -523,6 +547,10 @@ export default class PocomathInstance {
             }
          }
       }
+      // Need to metafy ground types
+      if (type === base) {
+         this._metafy(type)
+      }
       // update the typeOf function
       const imp = {}
       imp[type] = {uses: new Set(), does: () => Returns('string', () => type)}
@@ -640,7 +668,7 @@ export default class PocomathInstance {
       }
 
       // install the "base type" in the meta universe:
-      let beforeType = 'any'
+      let beforeType = UniversalType
       for (const other of spec.before || []) {
          if (other in this.templates) {
             beforeType = other
@@ -648,6 +676,13 @@ export default class PocomathInstance {
          }
       }
       this._metaTyped.addTypes([{name: base, test: spec.base}], beforeType)
+      // Add conversions to the base type:
+      if (spec.from && spec.from[theTemplateParam]) {
+         for (const ground of this._metafiedTypes) {
+            this._metaTyped.addConversion(
+               {from: ground, to: base, convert: spec.from[theTemplateParam]})
+         }
+      }
       this._instantiationsOf[base] = new Set()
 
       // update the typeOf function
@@ -750,45 +785,39 @@ export default class PocomathInstance {
    /**
     * Reset an operation to require creation of typed-function,
     * and if it has no implementations so far, set them up.
+    * name is the name of the operation, badType is a type that has been
+    * invalidated, and reasons is a set of specific operations/signatures
+    * that have been invalidated
     */
-   _invalidate(name, reason) {
+   _invalidate(name, badType = '', reasons = new Set()) {
       if (!(name in this._imps)) {
          this._imps[name] = {}
          this._TFimps[name] = {}
          this._metaTFimps[name] = {}
       }
-      if (reason) {
-         // Make sure no TF imps that depend on reason remain:
-         for (const [signature, behavior] of Object.entries(this._imps[name])) {
-            let invalidated = false
-            if (reason.charAt(0) === ':') {
-               const badType = reason.slice(1)
-               if (signature.includes(badType)) invalidated = true
+      // Go through each TF imp and invalidate it if need be
+      for (const [signature, imp] of Object.entries(this._TFimps[name])) {
+         if (imp.deferred
+             || (badType && signature.includes(badType))
+             || exists(imp.uses, dep => {
+                const [func, sig] = dep.split(/[()]/)
+                return reasons.has(dep)
+                   || (reasons.has(func) && !(sig in this._TFimps[func]))
+             })) {
+            // Invalidate this implementation:
+            delete this._TFimps[name][signature]
+            const behavior = imp.fromBehavior
+            if (behavior.explicit) {
+               behavior.resolved = false
             } else {
-               for (const dep of behavior.uses) {
-                  if (dep.includes(reason)) {
-                     invalidated = true
-                     break
-                  }
-               }
-            }
-            if (invalidated) {
-               if (behavior.explicit) {
-                  if (behavior.resolved) delete this._TFimps[signature]
-                  behavior.resolved = false
-               } else {
-                  for (const fullSig
-                       of Object.values(behavior.hasInstantiations)) {
-                     delete this._TFimps[fullSig]
-                  }
-                  behavior.hasInstantiations = {}
-               }
+               delete behavior.hasInstantiations[imp.instance]
             }
+            reasons.add(`${name}(${signature})`)
          }
       }
       if (this._invalid.has(name)) return
       this._invalid.add(name)
-      this._invalidateDependents(name)
+      this._invalidateDependents(name, badType, reasons)
       const self = this
       Object.defineProperty(this, name, {
          configurable: true,
@@ -802,11 +831,14 @@ export default class PocomathInstance {
 
    /**
     * Invalidate all the dependents of a given property of the instance
+    * reasons is a set of invalidated signatures
     */
-   _invalidateDependents(name) {
+   _invalidateDependents(name, badType, reasons = new Set()) {
+      if (name.charAt(0) === ':') badType = name.slice(1)
+      else reasons.add(name)
       if (name in this._affects) {
          for (const ancestor of this._affects[name]) {
-            this._invalidate(ancestor, name)
+            this._invalidate(ancestor, badType, reasons)
          }
       }
    }
@@ -847,7 +879,7 @@ export default class PocomathInstance {
       for (const [rawSignature, behavior] of usableEntries) {
          if (behavior.explicit) {
             if (!(behavior.resolved)) {
-               this._addTFimplementation(tf_imps, rawSignature, behavior)
+               this._addTFimplementation(name, tf_imps, rawSignature, behavior)
                tf_imps[rawSignature]._pocoSignature = rawSignature
                behavior.resolved = true
             }
@@ -867,11 +899,18 @@ export default class PocomathInstance {
          }
          /* First, add the known instantiations, gathering all types needed */
          if (ubType) behavior.needsInstantiations.add(ubType)
+         const nargs = typeListOfSignature(rawSignature).length
          let instantiationSet = new Set()
          const ubTypes = new Set()
          if (!ubType) {
             // Collect all upper-bound types for this signature
             for (const othersig in imps) {
+               const otherNargs = typeListOfSignature(othersig).length
+               if (nargs !== otherNargs) {
+                  // crude criterion that it won't match, that ignores
+                  // rest args, but hopefully OK for prototype
+                  continue
+               }
                const thisUB = upperBounds.exec(othersig)
                if (thisUB) ubTypes.add(thisUB[2])
                let basesig = othersig.replaceAll(templateCall, '')
@@ -881,7 +920,7 @@ export default class PocomathInstance {
                      basesig, theTemplateParam, '')
                   if (testsig === basesig) {
                      // that is not also top-level
-                     for (const templateType of typeListOfSignature(basesig)) {
+                     for (let templateType of typeListOfSignature(basesig)) {
                         if (templateType.slice(0,3) === '...') {
                            templateType = templateType.slice(3)
                         }
@@ -894,21 +933,20 @@ export default class PocomathInstance {
          for (const instType of behavior.needsInstantiations) {
             instantiationSet.add(instType)
             const otherTypes =
-                  ubType ? this.subtypesOf(instType) : this._priorTypes[instType]
+               ubType ? this.subtypesOf(instType) : this._priorTypes[instType]
             for (const other of otherTypes) {
                if (!(this._atOrBelowSomeType(other, ubTypes))) {
                   instantiationSet.add(other)
                }
             }
          }
-
          /* Prevent other existing signatures from blocking use of top-level
           * templates via conversions:
           */
          let baseSignature = rawSignature.replaceAll(templateCall, '')
          /* Any remaining template params are top-level */
          const signature = substituteInSignature(
-            baseSignature, theTemplateParam, 'any')
+            baseSignature, theTemplateParam, UniversalType)
          const hasTopLevel = (signature !== baseSignature)
          if (!ubType && hasTopLevel) {
             for (const othersig in imps) {
@@ -939,9 +977,9 @@ export default class PocomathInstance {
                }
             }
          }
-
          for (const instType of instantiationSet) {
-            this._instantiateTemplateImplementation(name, rawSignature, instType)
+            this._instantiateTemplateImplementation(
+               name, rawSignature, instType)
          }
          /* Now add the catchall signature */
          /* (Not needed if if it's a bounded template) */
@@ -996,6 +1034,7 @@ export default class PocomathInstance {
                            `Type inference failed for argument ${j} of ${name}`)
                      }
                      if (argType === 'any') {
+                        console.log('INCOMPATIBLE ARGUMENTS are', args)
                         throw TypeError(
                            `In call to ${name}, `
                               + 'incompatible template arguments:'
@@ -1018,9 +1057,10 @@ export default class PocomathInstance {
                   usedConversions = true
                   instantiateFor = self.joinTypes(argTypes, usedConversions)
                   if (instantiateFor === 'any') {
+                     let argDisplay = args.map(toString).join(', ')
                      throw TypeError(
                         `In call to ${name}, no type unifies arguments `
-                           + args.toString() + '; of types ' + argTypes.toString()
+                           + argDisplay + '; of types ' + argTypes.toString()
                            + '; note each consecutive pair must unify to a '
                            + 'supertype of at least one of them')
                   }
@@ -1042,7 +1082,9 @@ export default class PocomathInstance {
                for (j = 0; j < parTypes.length; ++j) {
                   if (wantTypes[j] !== parTypes[j] && parTypes[j].includes('<')) {
                      // actually used the param and is a template
-                     self._ensureTemplateTypes(parTypes[j], instantiateFor)
+                     const strippedType = parTypes[j].substr(
+                        parTypes[j].lastIndexOf('.') + 1)
+                     self._ensureTemplateTypes(strippedType, instantiateFor)
                   }
                }
 
@@ -1050,12 +1092,19 @@ export default class PocomathInstance {
                // But possibly since this resolution was grabbed, the proper
                // instantiation has been added (like if there are multiple
                // uses in the implementation of another method.
-               if (!(behavior.needsInstantiations.has(instantiateFor))) {
-                  behavior.needsInstantiations.add(instantiateFor)
+               let whatToDo
+               if (!(instantiateFor in behavior.hasInstantiations)) {
+                  const newImp = self._instantiateTemplateImplementation(
+                     name, rawSignature, instantiateFor)
+                  if (newImp) {
+                     whatToDo = {fn: newImp, implementation: newImp}
+                  }
                   self._invalidate(name)
                }
                const brandNewMe = self[name]
-               const whatToDo = self._typed.resolve(brandNewMe, args)
+               const betterToDo = self._typed.resolve(brandNewMe, args)
+               whatToDo = betterToDo || whatToDo
+
                // We can access return type information here
                // And in particular, if it might be a template, we should try to
                // instantiate it:
@@ -1069,8 +1118,13 @@ export default class PocomathInstance {
                }
                if (whatToDo === lastWhatToDo) {
                   throw new Error(
-                     `Infinite recursion in resolving $name called on`
-                        + args.map(x => x.toString()).join(','))
+                     `Infinite recursion in resolving ${name} called on `
+                        + args.map(x =>
+                           (typeof x === 'object'
+                            ? JSON.stringify(x)
+                            : x.toString())
+                        ).join(', ')
+                        + ` inferred to be ${wantSig}`)
                }
                lastWhatToDo = whatToDo
                const retval = whatToDo.implementation(...args)
@@ -1088,15 +1142,19 @@ export default class PocomathInstance {
          // correct return type a priori. Deferring because unclear what
          // aspects will be merged into typed-function.
          this._addTFimplementation(
-            meta_imps, signature, {uses: new Set(), does: patch})
+            name, meta_imps, signature,
+            {uses: new Set(), does: patch},
+            behavior)
          behavior.resolved = true
       }
-      this._correctPartialSelfRefs(name, tf_imps)
       // Make sure we have all of the needed (template) types; and if they
       // can't be added (because they have been instantiated too deep),
       // ditch the signature:
       const badSigs = new Set()
       for (const sig in tf_imps) {
+         if (!tf_imps[sig].uses) {
+            throw new ReferenceError(`MONKEY WRENCH: ${name} ${sig}`)
+         }
          for (const type of typeListOfSignature(sig)) {
             if (this._maybeInstantiate(type) === undefined) {
                badSigs.add(sig)
@@ -1117,11 +1175,13 @@ export default class PocomathInstance {
       if (Object.keys(tf_imps).length > 0) {
          tf = this._typed(name, tf_imps)
          tf.fromInstance = this
+         tf.isMeta = false
       }
       let metaTF
       if (Object.keys(meta_imps).length > 0) {
          metaTF = this._metaTyped(name, meta_imps)
          metaTF.fromInstance = this
+         metaTF.isMeta = true
       }
       this._meta[name] = metaTF
 
@@ -1215,10 +1275,12 @@ export default class PocomathInstance {
          return behavior.does(innerRefs)
       }
       const tf_imps = this._TFimps[name]
-      this._addTFimplementation(tf_imps, signature, {uses, does: patch})
+      this._addTFimplementation(
+         name, tf_imps, signature, {uses, does: patch}, behavior, instanceType)
       tf_imps[signature]._pocoSignature = templateSignature
       tf_imps[signature]._pocoInstance = instanceType
       behavior.hasInstantiations[instanceType] = signature
+      behavior.needsInstantiations.add(instanceType) // once we have it, keep it
       return tf_imps[signature]
    }
 
@@ -1226,17 +1288,23 @@ export default class PocomathInstance {
     * to typed-function implementations and inserts the result into plain
     * object imps
     */
-   _addTFimplementation(imps, signature, behavior) {
-      const {uses, does} = behavior
+   _addTFimplementation(
+      name, imps, signature, specificBehavior, fromImp, asInstance)
+   {
+      if (!fromImp) fromImp = specificBehavior
+      const {uses, does} = specificBehavior
       if (uses.length === 0) {
          const implementation = does()
+         implementation.uses = uses
+         implementation.fromInstance = this
+         implementation.fromBehavior = fromImp
+         implementation.instance = asInstance
          // could do something with return type information here
          imps[signature] = implementation
          return
       }
       const refs = {}
       let full_self_referential = false
-      let part_self_references = []
       for (const dep of uses) {
          let [func, needsig] = dep.split(/[()]/)
          /* Safety check that can perhaps be removed:
@@ -1252,82 +1320,71 @@ export default class PocomathInstance {
          }
          if (func === 'self') {
             if (needsig) {
-               /* Maybe we can resolve the self reference without troubling
-                * typed-function:
+               /* We now resolve all specific-signature self references
+                * here, without resorting to the facility in typed-function:
                 */
                if (needsig in imps && typeof imps[needsig] == 'function') {
                   refs[dep] = imps[needsig]
+                  continue
+               }
+               const needTypes = typesOfSignature(needsig)
+               const mergedTypes = Object.assign(
+                  {}, this.Types, this.Templates)
+               if (subsetOfKeys(needTypes, mergedTypes)) {
+                  func = name // just resolve it in limbo
                } else {
-                  if (full_self_referential) {
-                     throw new SyntaxError(
-                        'typed-function does not support mixed full and '
-                           + 'partial self-reference')
-                  }
-                  const needTypes = typesOfSignature(needsig)
-                  const mergedTypes = Object.assign(
-                     {}, this.Types, this.Templates)
-                  if (subsetOfKeys(needTypes, mergedTypes)) {
-                     part_self_references.push(needsig)
-                  }
+                  // uses an unknown type, so will get an undefined impl
+                  console.log(
+                     'WARNING: partial self-reference for', name, 'to',
+                     needsig, 'uses an unknown type')
+                  refs[dep] = undefined
+                  continue
                }
             } else {
-               if (part_self_references.length) {
-                  throw new SyntaxError(
-                     'typed-function does not support mixed full and '
-                        + 'partial self-reference')
-               }
                full_self_referential = true
-            }
-         } else {
-            if (this[func] === 'limbo') {
-               /* We are in the midst of bundling func */
-               let fallback =  true
-               /* So the first thing we can do is try the tf_imps we are
-                * accumulating:
-                */
-               if (needsig) {
-                  let typedUniverse
-                  let tempTF
-                  if (Object.keys(this._TFimps[func]).length > 0) {
-                     typedUniverse = this._typed
-                     tempTF = typedUniverse('dummy_' + func, this._TFimps[func])
-                  } else {
-                     typedUniverse = this._metaTyped
-                     tempTF = typedUniverse(
-                        'dummy_' + func, this._metaTFimps[func])
-                  }
-                  let result = undefined
-                  try {
-                     result = typedUniverse.find(tempTF, needsig, {exact: true})
-                  } catch {}
-                  if (result) {
-                     refs[dep] = result
-                     fallback = false
-                  }
-               }
-               if (fallback) {
-                  /* Either we need the whole function or the signature
-                   * we need is not available yet, so we have to use
-                   * an indirect reference to func. And given that, there's
-                   * really no helpful way to extract a specific signature
-                   */
-                  const self = this
-                  const redirect = function () { // is this the most efficient?
-                     return self[func].apply(this, arguments)
-                  }
-                  Object.defineProperty(redirect, 'name', {value: func})
-                  Object.defineProperty(redirect, 'fromInstance', {value: this})
-                  refs[dep] = redirect
-               }
-            } else {
-               // can bundle up func, and grab its signature if need be
-               let destination = this[func]
-               if (destination && needsig) {
-                  destination = this.resolve(func, needsig)
-               }
-               refs[dep] = destination
+               continue
             }
          }
+         if (this[func] === 'limbo') {
+            /* We are in the midst of bundling func (which may be ourself) */
+            /* So the first thing we can do is try the tf_imps we are
+             * accumulating:
+             */
+            if (needsig) {
+               const candidate = this.resolve(func, needsig)
+               if (typeof candidate === 'function') {
+                  refs[dep] = candidate
+                  continue
+               }
+            }
+            /* Either we need the whole function or the signature
+             * we need is not available yet, so we have to use
+             * an indirect reference to func. And given that, there's
+             * really no helpful way to extract a specific signature
+             */
+            const self = this
+            const redirect = function () { // is this the most efficient?
+               return self[func].apply(this, arguments)
+            }
+            Object.defineProperty(redirect, 'name', {value: func})
+            Object.defineProperty(redirect, 'fromInstance', {value: this})
+            refs[dep] = redirect
+            continue
+         }
+         // can bundle up func, and grab its signature if need be
+         let destination = this[func]
+         if (needsig) {
+            destination = this.resolve(func, needsig)
+         }
+         if (!destination) {
+            // Unresolved reference. This is allowed so that
+            // you can bundle up just some portions of the library,
+            // but let's warn.
+            console.log(
+               'WARNING: No definition found for dependency',
+               dep, 'needed by', name, '(', signature, ')')
+         }
+         refs[dep] = destination
       }
       if (full_self_referential) {
          imps[signature] = this._typed.referToSelf(self => {
@@ -1335,108 +1392,27 @@ export default class PocomathInstance {
             const implementation = does(refs)
             Object.defineProperty(implementation, 'name', {value: does.name})
             implementation.fromInstance = this
+            implementation.uses = uses
+            implementation.instance = asInstance
+            implementation.fromBehavior = fromImp
             // What are we going to do with the return type info in here?
             return implementation
          })
-         return
-      }
-      if (part_self_references.length) {
-         /* There is an obstruction here. The list part_self_references
-          * might contain a signature that requires conversion for self to
-          * handle. But I advocated this not be allowed in typed.referTo, which
-          * made sense for human-written functions, but is unfortunate now.
-          * So we have to defer creating these and correct them later, at
-          * least until we can add an option to typed-function.
-          */
-         imps[signature] = {
-            deferred: true,
-            builtRefs: refs,
-            sigDoes: does,
-            fromInstance: this,
-            psr: part_self_references
-         }
+         imps[signature].uses = uses
+         imps[signature].fromInstance = this
+         imps[signature].instance = asInstance
+         imps[signature].fromBehavior = fromImp
          return
       }
       const implementation = does(refs)
       implementation.fromInstance = this
+      implementation.fromBehavior = fromImp
+      implementation.instance = asInstance
+      implementation.uses = uses
       // could do something with return type information here?
       imps[signature] = implementation
    }
 
-   _correctPartialSelfRefs(name, imps) {
-      for (const aSignature in imps) {
-         if (!(imps[aSignature].deferred)) continue
-         const deferral = imps[aSignature]
-         const part_self_references = deferral.psr
-         const corrected_self_references = []
-         const remaining_self_references = []
-         const refs = deferral.builtRefs
-         for (const neededSig of part_self_references) {
-            // Have to find a match for neededSig among the other signatures
-            // of this function. That's a job for typed-function, but we will
-            // try here:
-            if (neededSig in imps) { // the easy case
-               corrected_self_references.push(neededSig)
-               remaining_self_references.push(neededSig)
-               continue
-            }
-            // No exact match, try to get one that matches with
-            // subtypes since the whole conversion thing in typed-function
-            // is too complicated to reproduce
-            let foundSig = this._findSubtypeImpl(name, imps, neededSig)
-            if (foundSig) {
-               corrected_self_references.push(foundSig)
-               remaining_self_references.push(neededSig)
-            } else {
-               // Maybe it's a template instance we don't yet have
-               foundSig = this._findSubtypeImpl(
-                  name, this._imps[name], neededSig)
-               if (foundSig) {
-                  const match = this._pocoFindSignature(name, neededSig)
-                  const neededTemplate = match.fn._pocoSignature
-                  const neededInstance = whichSigInstance(
-                     neededSig, neededTemplate)
-                  const neededImplementation =
-                     this._instantiateTemplateImplementation(
-                        name, neededTemplate, neededInstance)
-                  if (!neededImplementation) {
-                     refs[`self(${neededSig})`] = match.implementation
-                  } else {
-                     if (typeof neededImplementation === 'function') {
-                        refs[`self(${neededSig})`] = neededImplementation
-                     } else {
-                        corrected_self_references.push(neededSig)
-                        remaining_self_references.push(neededSig)
-                     }
-                  }
-               } else {
-                  throw new Error(
-                     'Implement inexact self-reference in typed-function for '
-                        + `${name}(${neededSig})`)
-               }
-            }
-         }
-         const does = deferral.sigDoes
-         if (remaining_self_references.length > 0) {
-            imps[aSignature] = this._typed.referTo(
-               ...corrected_self_references, (...impls) => {
-                  for (let i = 0; i < remaining_self_references.length; ++i) {
-                     refs[`self(${remaining_self_references[i]})`] = impls[i]
-                  }
-                  const implementation = does(refs)
-                  // What will we do with the return type info in here?
-                  return implementation
-               }
-            )
-         } else {
-            imps[aSignature] = does(refs)
-         }
-         imps[aSignature]._pocoSignature = deferral._pocoSignature
-         imps[aSignature]._pocoInstance = deferral._pocoInstance
-         imps[aSignature].fromInstance = deferral.fromInstance
-      }
-   }
-
    /* This function analyzes the template and makes sure the
     * instantiations of it for type and all prior types of type are present
     * in the instance.
@@ -1542,7 +1518,8 @@ export default class PocomathInstance {
       return wantsType
    })
 
-   _findSubtypeImpl(name, imps, neededSig) {
+   _findSubtypeImpl(name, imps, neededSig, raw = false) {
+      const detemplate = !raw
       if (neededSig in imps) return neededSig
       let foundSig = false
       const typeList = typeListOfSignature(neededSig)
@@ -1550,21 +1527,21 @@ export default class PocomathInstance {
          const otherTypeList = typeListOfSignature(otherSig)
          if (typeList.length !== otherTypeList.length) continue
          let allMatch = true
-         let paramBound = 'any'
+         let paramBound = UniversalType
          for (let k = 0; k < typeList.length; ++k) {
             let myType = typeList[k]
             let otherType = otherTypeList[k]
             if (otherType === theTemplateParam) {
-               otherTypeList[k] = paramBound
+               if (detemplate) otherTypeList[k] = paramBound
                otherType = paramBound
             }
             if (otherType === restTemplateParam) {
-               otherTypeList[k] = `...${paramBound}`
+               if (detemplate) otherTypeList[k] = `...${paramBound}`
                otherType = paramBound
             }
             const adjustedOtherType = otherType.replaceAll(templateCall, '')
             if (adjustedOtherType !== otherType) {
-               otherTypeList[k] = adjustedOtherType
+               if (detemplate) otherTypeList[k] = adjustedOtherType
                otherType = adjustedOtherType
             }
             if (myType.slice(0,3) === '...') myType = myType.slice(3)
@@ -1573,10 +1550,13 @@ export default class PocomathInstance {
             if (otherBound) {
                paramBound = otherBound[2]
                otherType = paramBound
-               otherTypeList[k] = otherBound[1].replaceAll(
-                  theTemplateParam, paramBound)
+               if (detemplate) {
+                  otherTypeList[k] = otherBound[1].replaceAll(
+                     theTemplateParam, paramBound)
+               }
             }
             if (otherType === 'any') continue
+            if (otherType === UniversalType) continue
             if (myType === otherType) continue
             if (otherType in this.Templates) {
                const [myBase] = splitTemplate(myType)
@@ -1584,7 +1564,7 @@ export default class PocomathInstance {
                if (this.instantiateTemplate(otherType, myType)) {
                   let dummy
                   dummy = this[name] // for side effects
-                  return this._findSubtypeImpl(name, this._imps[name], neededSig)
+                  return this._findSubtypeImpl(name, this._imps[name], neededSig, raw)
                }
             }
             if (!(otherType in this.Types)) {
@@ -1608,6 +1588,7 @@ export default class PocomathInstance {
          typedFunction = this[name]
       }
       const haveTF = this._typed.isTypedFunction(typedFunction)
+         && !(typedFunction.isMeta)
       if (haveTF) {
          // First try a direct match
          let result
@@ -1622,6 +1603,10 @@ export default class PocomathInstance {
               of typedFunction._typedFunctionData.signatureMap) {
             let allMatched = true
             const implTypes = typeListOfSignature(implSig)
+            if (implTypes.length > wantTypes.length) {
+               // Not enough arguments for that implementation
+               continue
+            }
             for (let i = 0; i < wantTypes.length; ++i) {
                const implIndex = Math.min(i, implTypes.length - 1)
                let implType = implTypes[implIndex]
@@ -1646,7 +1631,7 @@ export default class PocomathInstance {
          }
       }
       if (!(this._imps[name])) return undefined
-      const foundsig = this._findSubtypeImpl(name, this._imps[name], sig)
+      const foundsig = this._findSubtypeImpl(name, this._imps[name], sig, 'raw')
       if (foundsig) {
          if (haveTF) {
             try {
@@ -1654,19 +1639,74 @@ export default class PocomathInstance {
             } catch {
             }
          }
-         try {
-            return this._metaTyped.findSignature(this._meta[name], foundsig)
-         } catch {
+         const instantiationMatcher =
+            '^'
+            + substituteInSignature(foundsig, theTemplateParam, '(.*)')
+               .replaceAll(UniversalType, '(.*)')
+            + '$'
+         const instanceMatch = sig.match(instantiationMatcher)
+         let possibleInstantiator = false
+         if (instanceMatch) {
+            possibleInstantiator = instanceMatch[1]
+            for (let i = 2; i < instanceMatch.length; ++i) {
+               if (possibleInstantiator !== instanceMatch[i]) {
+                  possibleInstantiator = false
+                  break
+               }
+            }
+         }
+         if (possibleInstantiator) {
+            const behavior = this._imps[name][foundsig]
+            let newInstance
+            if (behavior) {
+               if (!(possibleInstantiator in behavior.hasInstantiations)) {
+                  newInstance = this._instantiateTemplateImplementation(
+                     name, foundsig, possibleInstantiator)
+               } else {
+                  // OK, so we actually have the instantiation. Let's get it
+                  newInstance = this._TFimps[name][sig]
+               }
+               // But we may not have taken advantage of conversions
+               this._invalidate(name)
+               const tryAgain = this[name]
+               let betterInstance
+               if (this._typed.isTypedFunction(tryAgain)) {
+                  betterInstance = this._typed.findSignature(tryAgain, sig)
+               }
+               if (betterInstance) {
+                  newInstance = betterInstance
+               } else {
+                  newInstance = {
+                     fn: newInstance,
+                     implementation: newInstance
+                  }
+               }
+               if (newInstance) return newInstance
+            }
+         }
+         const catchallSig = this._findSubtypeImpl(name, this._imps[name], sig)
+         if (catchallSig !== foundsig) {
+            try {
+               return this._metaTyped.findSignature(
+                  this._meta[name], catchallSig)
+            } catch {
+            }
          }
          // We have an implementation but not a typed function. Do the best
          // we can:
-         const foundImpl = this._imps[name][foundsig]
+         const restoredSig = foundsig.replaceAll('ground', theTemplateParam)
+         const foundImpl = this._imps[name][restoredSig]
          const needs = {}
          for (const dep of foundImpl.uses) {
-            const [base, sig] = dep.split('()')
-            needs[dep] = this.resolve(base, sig)
+            const [base, sig] = dep.split(/[()]/)
+            if (sig) {
+               needs[dep] = this.resolve(base, sig)
+            } else {
+               needs[dep] = this[dep]
+            }
          }
          const pseudoImpl = foundImpl.does(needs)
+         pseudoImpl.fromInstance = this
          return {fn: pseudoImpl, implementation: pseudoImpl}
       }
       // Hmm, no luck. Make sure bundle is up-to-date and retry:

src/core/extractors.mjs

@@ -6,7 +6,7 @@ export function dependencyExtractor(destinationSet) {
    return new Proxy({}, {
       get: (target, property) => {
          destinationSet.add(property)
-         return {}
+         return {checkingDependency: true}
       }
    })
 }

src/number/add.mjs

@@ -4,5 +4,5 @@ export * from './Types/number.mjs'
 export const add = {
     // Note the below assumes that all subtypes of number that will be defined
     // are closed under addition!
-    'T:number, T': ({T}) => Returns(T, (m,n) => m+n)
+    'T:number,T': ({T}) => Returns(T, (m,n) => m+n)
 }

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

src/tuple/tuple.mjs

@@ -6,5 +6,5 @@ export {Tuple} from './Types/Tuple.mjs'
  * are convertible to the same type.
  */
 export const tuple = {
-    '...T': ({T}) => Returns(`Tuple<${T}>`, args => ({elts: args}))
+   '...T': ({T}) => Returns(`Tuple<${T}>`, args => ({elts: args}))
 }

test/complex/_all.mjs

@@ -45,6 +45,11 @@ describe('complex', () => {
       assert.ok(!(math.equal(math.complex(45n, 3n), 45n)))
    })
 
+   it('tests for reality', () => {
+      assert.ok(math.isReal(math.complex(3, 0)))
+      assert.ok(!(math.isReal(math.complex(3, 2))))
+   })
+
    it('computes gcd', () => {
       assert.deepStrictEqual(
          math.gcd(math.complex(53n, 56n), math.complex(47n, -13n)),

test/complex/_polynomialRoot.mjs

@@ -0,0 +1,63 @@
+import assert from 'assert'
+import * as approx from '../../tools/approx.mjs'
+import math from '../../src/pocomath.mjs'
+
+describe('polynomialRoot', () => {
+   it('should solve a linear equation with real coefficients', function () {
+      assert.deepEqual(math.polynomialRoot(6, 3), math.tuple(-2))
+      assert.deepEqual(
+         math.polynomialRoot(math.complex(-3, 2), 2),
+         math.tuple(math.complex(1.5, -1)))
+      assert.deepEqual(
+         math.polynomialRoot(math.complex(3, 1), math.complex(-1, -1)),
+         math.tuple(math.complex(2, -1)))
+   })
+   // Should be safe now to capture the functions:
+   const complex = math.complex
+   const pRoot = math.polynomialRoot
+   const tup = math.tuple
+   it('should solve a quadratic equation with a double root', function () {
+      assert.deepEqual(pRoot(4, 4, 1), tup(-2))
+      assert.deepEqual(
+         pRoot(complex(0, 2), complex(2, 2), 1), tup(complex(-1, -1)))
+   })
+   it('should solve a quadratic with two distinct roots', function () {
+      assert.deepEqual(pRoot(-3, 2, 1), tup(1, -3))
+      assert.deepEqual(pRoot(-2, 0, 1), tup(math.sqrt(2), -math.sqrt(2)))
+      assert.deepEqual(
+         pRoot(4, 2, 1),
+         tup(complex(-1, math.sqrt(3)), complex(-1, -math.sqrt(3))))
+      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)
+}