Numeric algorithms, as a sanity check on the calculus

Λs is pure and total: no recursion, no conditionals, no mutable state. The question this file answers is whether that leaves enough to write recognizable numerical methods, or only formulas.

The answer is that straight-line numerical methods transcribe directly, and their types are their specifications. A central difference has type (Q v → Q u) → Q v → Q v → Q (u/v): differentiate a position with respect to a time and you get a velocity, by derivation rather than by comment. What is missing is iteration: one Runge–Kutta step is expressible, a solver loop is not.

What the checker is doing that a syntactic system could not

Unit equality here is decided by arithmetic on exponent vectors, not by normalizing syntax. In the Runge–Kutta step below, the intermediate (h/2)·k₁ has unit t · (y/t), and add demands it equal y. As exponent vectors those are t + (y − t) and y, and the checker settles it by adding rationals. There is no normalization pass, so there is nothing to get wrong.

Scopes

Two nested dimension binders and two nested unit binders. Inside the body the innermost unit variable is de Bruijn 0 and the outer one is 1; likewise for dimensions.

@[reducible, inline]

abbrev LambdaS.Algorithms.uU : UExp Examples.Base 2

Equations

• LambdaS.Algorithms.uU = LambdaS.Term.ofVar 1

@[reducible, inline]

abbrev LambdaS.Algorithms.vU : UExp Examples.Base 2

Equations

• LambdaS.Algorithms.vU = LambdaS.Term.ofVar 0

@[reducible, inline]

abbrev LambdaS.Algorithms.dU : DExp Examples.Dim 1

Equations

• LambdaS.Algorithms.dU = LambdaS.Term.ofVar 0

@[reducible, inline]

abbrev LambdaS.Algorithms.dV : DExp Examples.Dim 2

Equations

• LambdaS.Algorithms.dV = LambdaS.Term.ofVar 0

@[reducible, inline]

abbrev LambdaS.Algorithms.poly (body : Tm Examples.Base Examples.Dim 2 2) : Examples.Term₀

Equations

• LambdaS.Algorithms.poly body = (LambdaS.Tm.ulam LambdaS.Algorithms.dU (LambdaS.Tm.ulam LambdaS.Algorithms.dV body).dlam).dlam

@[reducible, inline]

abbrev LambdaS.Algorithms.at₂ (du : DExp Examples.Dim 0) (u : UExp Examples.Base 0) (dv : DExp Examples.Dim 0) (v : UExp Examples.Base 0) (e : Examples.Term₀) : Examples.Term₀

Equations

• LambdaS.Algorithms.at₂ du u dv v e = (((LambdaS.Tm.dapp e du).uapp u).dapp dv).uapp v

@[reducible, inline]

abbrev LambdaS.Algorithms.dLen : DExp Examples.Dim 0

Equations

• LambdaS.Algorithms.dLen = LambdaS.Term.ofBase LambdaS.Examples.Dim.length

@[reducible, inline]

abbrev LambdaS.Algorithms.dTime : DExp Examples.Dim 0

Equations

• LambdaS.Algorithms.dTime = LambdaS.Term.ofBase LambdaS.Examples.Dim.time

Central difference

f'(x) ≈ (f(x+h) − f(x−h)) / 2h. Inside the three value binders, %0 is h, %1 is x, %2 is f.

def LambdaS.Algorithms.deriv : Examples.Term₀

Equations

• One or more equations did not get rendered due to their size.

Simpson's rule

∫ₐᵇ f ≈ ((b−a)/6)·(f(a) + 4·f((a+b)/2) + f(b)). %0 is b, %1 is a, %2 is f.

def LambdaS.Algorithms.simpson : Examples.Term₀

Equations

• One or more equations did not get rendered due to their size.

One Runge–Kutta step

y' = y + (h/6)(k₁ + 2k₂ + 2k₃ + k₄) for dy/dt = f(t,y). Binders, innermost first: %0 is y, %1 is t, %2 is h, %3 is f.

This is where unit equality does real work. k₁ : Q (y/t), so (h/2)·k₁ has unit t·(y/t), and adding it to y requires the checker to see those as the same exponent vector.

def LambdaS.Algorithms.rk4 : Examples.Term₀

Equations

• One or more equations did not get rendered due to their size.

Tolerances, where dimensioned types catch the common bug

An absolute tolerance has the unit of the quantity; a relative one is dimensionless. Writing |x − y| < 1e-6 with a bare literal is a type error, and the fix the checker forces is the correct one.

@[reducible, inline]

abbrev LambdaS.Algorithms.tolCtx : Ctx Examples.Base Examples.Dim 0 0

Equations

• LambdaS.Algorithms.tolCtx = [LambdaS.Ty.Q LambdaS.Examples.m, LambdaS.Ty.Q LambdaS.Examples.m]

Running them

Free fall: s(t) = ½·g·t² with g = 9.81 m/s². A central difference is exact on a quadratic and Simpson's rule is exact on a cubic, so both answers are exact and the guards can be tight.

@[reducible, inline]

abbrev LambdaS.Algorithms.accel : UExp Examples.Base 0

Equations

• LambdaS.Algorithms.accel = LambdaS.Term.div LambdaS.Examples.m (LambdaS.Term.mul LambdaS.Examples.sec LambdaS.Examples.sec)

def LambdaS.Algorithms.gravity : Examples.Term₀

Equations

• LambdaS.Algorithms.gravity = (LambdaS.Tm.lit (981 / 100)).mul (LambdaS.Tm.ucon LambdaS.Algorithms.accel)

def LambdaS.Algorithms.freeFall : Examples.Term₀

Equations

• One or more equations did not get rendered due to their size.

@[reducible, inline]

abbrev LambdaS.Algorithms.derivMS : Examples.Term₀

Equations

• LambdaS.Algorithms.derivMS = LambdaS.Algorithms.at₂ LambdaS.Algorithms.dLen LambdaS.Examples.m LambdaS.Algorithms.dTime LambdaS.Examples.sec LambdaS.Algorithms.deriv

@[reducible, inline]

abbrev LambdaS.Algorithms.simpsonMS : Examples.Term₀

Equations

• LambdaS.Algorithms.simpsonMS = LambdaS.Algorithms.at₂ LambdaS.Algorithms.dLen LambdaS.Examples.m LambdaS.Algorithms.dTime LambdaS.Examples.sec LambdaS.Algorithms.simpson

def LambdaS.Algorithms.twoSec : Examples.Term₀

Equations

• LambdaS.Algorithms.twoSec = (LambdaS.Tm.lit 2).mul (LambdaS.Tm.ucon LambdaS.Examples.sec)

def LambdaS.Algorithms.zeroSec : Examples.Term₀

Equations

• LambdaS.Algorithms.zeroSec = (LambdaS.Tm.lit 0).mul (LambdaS.Tm.ucon LambdaS.Examples.sec)

def LambdaS.Algorithms.milliSec : Examples.Term₀

Equations

• LambdaS.Algorithms.milliSec = (LambdaS.Tm.lit (1 / 1000)).mul (LambdaS.Tm.ucon LambdaS.Examples.sec)

def LambdaS.Algorithms.speedAt2 : Examples.Term₀

Equations

• LambdaS.Algorithms.speedAt2 = ((LambdaS.Tm.app LambdaS.Algorithms.derivMS LambdaS.Algorithms.freeFall).app LambdaS.Algorithms.twoSec).app LambdaS.Algorithms.milliSec

def LambdaS.Algorithms.absementTo2 : Examples.Term₀

Equations

• LambdaS.Algorithms.absementTo2 = ((LambdaS.Tm.app LambdaS.Algorithms.simpsonMS LambdaS.Algorithms.freeFall).app LambdaS.Algorithms.zeroSec).app LambdaS.Algorithms.twoSec

def LambdaS.Algorithms.run (e : Examples.Term₀) : Option Float

Equations

• LambdaS.Algorithms.run e = match LambdaS.evalC (fun (x x_1 : LambdaS.UExp LambdaS.Examples.Base 0) => 1.0) 10000 [] e with | some (LambdaS.Val.scalar x) => some x.mag | x => none

def LambdaS.Algorithms.report : String

Declared conversions, compiled

The conversion oracle below is not fun _ _ => 1.0: it is computed from the declared magnitudes (meter 1, foot 0.3048, yard 0.9144). magB is a transcription by hand, at Float, of the valuation ψyd that LambdaS.Examples proves satisfies the declaration set; nothing derives it from ψyd, and the agreement between the two is checked by the numbers below rather than by a theorem. evalC_convert_declared is the theorem that the instrumented evaluator at carrier ℝ multiplies by exactly the declared factor; this is that theorem's number coming out of the compiled binary.

Two routes from yards to meters (direct, and through feet) print the same number, which is convChain_eq and yard_forced made observable: the factors are forced by the declarations, so there is no route to get wrong.

def LambdaS.Algorithms.magB : Examples.Base → Float

The declared magnitude of each base unit, in meters (and SI mates): ψyd of LambdaS.Examples, transcribed at Float.

Equations

• LambdaS.Algorithms.magB LambdaS.Examples.Base.meter = 1
• LambdaS.Algorithms.magB LambdaS.Examples.Base.foot = 0.3048
• LambdaS.Algorithms.magB LambdaS.Examples.Base.yard = 0.9144
• LambdaS.Algorithms.magB LambdaS.Examples.Base.kilogram = 1
• LambdaS.Algorithms.magB LambdaS.Examples.Base.second = 1

def LambdaS.Algorithms.magF (w : UExp Examples.Base 0) : Float

The magnitude of a compound unit: the product of its bases' declared magnitudes, at their exponents.

Equations

• One or more equations did not get rendered due to their size.

def LambdaS.Algorithms.cfDecl (u v : UExp Examples.Base 0) : Float

The conversion oracle the declarations determine: a ratio of magnitudes, which is conv at the Float shadow of ψyd.

Equations

• LambdaS.Algorithms.cfDecl u v = LambdaS.Algorithms.magF u / LambdaS.Algorithms.magF v

def LambdaS.Algorithms.hundredYards : Examples.Term₀

Equations

• LambdaS.Algorithms.hundredYards = (LambdaS.Tm.lit 100).mul (LambdaS.Tm.ucon LambdaS.Examples.yd)

def LambdaS.Algorithms.inFeet : Examples.Term₀

Equations

• LambdaS.Algorithms.inFeet = LambdaS.Tm.convert LambdaS.Algorithms.hundredYards LambdaS.Examples.yd LambdaS.Examples.ft

def LambdaS.Algorithms.inMeters : Examples.Term₀

Equations

• LambdaS.Algorithms.inMeters = LambdaS.Tm.convert LambdaS.Algorithms.hundredYards LambdaS.Examples.yd LambdaS.Examples.m

def LambdaS.Algorithms.viaFeet : Examples.Term₀

Equations

• LambdaS.Algorithms.viaFeet = (LambdaS.Tm.convert LambdaS.Algorithms.hundredYards LambdaS.Examples.yd LambdaS.Examples.ft).convert LambdaS.Examples.ft LambdaS.Examples.m

def LambdaS.Algorithms.runDecl (e : Examples.Term₀) : Option Float

Run a closed term under the declared-conversion oracle cfDecl, returning its scalar magnitude; the declared-factor counterpart of run.

Equations

• LambdaS.Algorithms.runDecl e = match LambdaS.evalC LambdaS.Algorithms.cfDecl 10000 [] e with | some (LambdaS.Val.scalar x) => some x.mag | x => none

One yard per foot

1_yd / 1_ft denotes the magnitude 1, at unit yd/ft: unit constants denote the number one, so the quotient is convert-free and, by den_indep, cannot see the declared magnitudes. The number 3 is one conversion away, (1_yd / 1_ft) in 1, and it must be: the 3 is data from the declaration table, and conversion is the only construct that reads it.

def LambdaS.Algorithms.ydPerFt : Examples.Term₀

Equations

• LambdaS.Algorithms.ydPerFt = (LambdaS.Tm.ucon (LambdaS.Term.ofBase LambdaS.Examples.Base.yard)).div (LambdaS.Tm.ucon (LambdaS.Term.ofBase LambdaS.Examples.Base.foot))

@[reducible, inline]

abbrev LambdaS.Algorithms.one : UExp Examples.Base 0

The unit one: the SI's name for the unit of dimension-one quantities (9th SI Brochure, 5.4.7), which the brochure calls the neutral element of any system of units. Here it is the empty exponent vector. The surface reading of the conversion below is (1 yd / 1 ft) in one.

Equations

• LambdaS.Algorithms.one = 1

def LambdaS.Algorithms.ydPerFtIn1 : Examples.Term₀

Equations

• One or more equations did not get rendered due to their size.

def LambdaS.Algorithms.ydPerFtViaFt : Examples.Term₀

The other route to the same number: convert one operand before dividing, (1_yd in ft) / 1_ft. The quotient is at unit ft/ft = 1 with no conversion of its own, and the 3 again enters through the one conversion in the term.

Equations

• One or more equations did not get rendered due to their size.

def LambdaS.Algorithms.reportDecl : String