Crafting Formulas

Lambdas All the Way Down

Marvin Borner

2024-08-06

This project is based on the work of u/DaVinci103

Since its inception, the bruijn programming language had support for positive and negative integers.

One thing I love about lambda calculus encodings for numbers is that they are usually infinitely precise – only limited by your computer’s memory. Still, lazy (call-by-need) or optimal reducers can sometimes do magical optimizations automatically that would have to be hard-coded in other languages.

It’s natural to want to extend the integers to arbitrary-precision fractions, and later real and complex numbers. Encoding rational numbers, it turns out, is as easy as writing two integers as a (Church) pair. The relevant arithmetic operations are fairly trivial to implement as well.

However, for the real numbers I couldn’t come up with any elegant solution:

meme

This is, until I saw a great post on reddit about this exact topic.

So, in this article I want to show how the new arbitrary-precision arithmetic of bruijn is implemented, including some “more mathematical” real/complex formulas and ways to approximate their result. I also want to use this opportunity to introduce you to coding in bruijn and lambda calculus in general.

Note that all of the definitions (including library imports) have a 1:1 correspondence to pure lambda calculus – bruijn is just syntactic sugar (syntax explanation). We’re literally crafting complex formulas from the bottom up. You can click on any identifier to jump to its definition.

Integer

Natural numbers in lambda calculus are most commonly encoded as Church numerals. As in the meme above, they can be extended to integers by writing them as a pair.

Z={aba,bN}\mathbb{Z}=\{a-b\mid a,b\in\mathbb{N}\}

While the implementations are quite small and elegant, any number now takes O(2n)\mathcal{O}(2n) space:

:import std/Pair .
 

35=23-5=-2

# syntactic sugar with 'u' suffix for unary/Church encoding
:test ((+3u) : (+5u)) ([[1 (1 (1 0))]] : [[1 (1 (1 (1 (1 0))))]])

Instead, I use balanced ternary numbers. Theoretically one could use any base, but an investigation by Torben Mogensen showed that bigger bases result in much more complicated implementations without giving significant performance benefits.

The encoding of a number nn in arbitrary base b>1b>1:

nb=λS(b)(d0dlogb(n) b) \langle n\rangle_b=\lambda^{S(b)}(d_0\circ\dots\circ d_{\lfloor\log_b(n)\rfloor}\ b) where de Bruijn index  dinbi(modb) \mathrm{where\ de\ Bruijn\ index\ }\ d_i\equiv\frac{n}{b^i}\pmod{b}

Balanced ternary numbers provide a very elegant way to denote negative numbers by assigning one of its three digits to a negative multiplier. I use the second de Bruijn index as the negative trit, the zeroth index as the zero trit, and the first index as the positive trit:

# from left to right:
 

2=130+1312=-1\cdot3^0 + 1\cdot3^1

:test ((+2)) ([[[[2 (1 3)]]]])
 

2=130+131-2=1\cdot3^0 + -1\cdot3^1

:test ((-2)) ([[[[1 (2 3)]]]])

Now we only need O(logn)\mathcal{O}(\log n) space!

You can find details about definitions of the arithmetic operations in std/Number/Ternary and my Rosetta Code entry. All of the implementations below import integer maths into the N namespace.

Rational

The rational numbers represent the ratio between two integers:

Q={aba,bZ, b0}\mathbb{Q}=\left\{\frac ab\mid a,b\in\mathbb{Z},\ b\ne0\right\}

The obvious implementation is therefore to use a pair of two balanced ternary numbers! Since we have the restriction of a non-zero bb, we subtract 11 from the denominator in the encoding:

:import std/Pair .
:import std/Math N
 

5.0=515.0=\frac51

# new syntactic sugar with 'q' suffix!
:test ((+5.0q)) ((+5) : (+0))
 

3.14=157503.14=\frac{157}{50}

:test ((+3.14q)) ((+157) : (+49))

# converts a balanced ternary number to a rational number
number→rational [0 : (+0)] ⧗ Number → Rational

:test (number→rational (+5)) ((+5.0q))

Operators

Now, to check equality of two rational numbers pp and qq, we only need to compare the product of pp’s denominator and qq’s numerator and vice versa. In the following snippets I uncurry the pairs using the thrush combinator &.

:import std/List L
:import std/Combinator .
:import std/Logic/Binary .

# returns true if two rational numbers are equal

p=q    npdp+1=nqdq+1    np(dq+1)=(dp+1)nqp=q\iff\frac{n_p}{d_p+1}=\frac{n_q}{d_q+1}\iff n_p(d_q+1)=(d_p+1)n_q

eq? &[[&[[N.eq? (N.mul 3 N.++0) (N.mul N.++2 1)]]]] ⧗ Rational → Rational → Boolean

…=?… eq?

:test (((+1) : (+3)) =? ((+2) : (+7))) (true)
:test ((+0.25q) =? (+0.25q)) (true)

# returns true if a rational number is greater than another

p>q    npdp+1>nqdq+1    np(dq+1)>(dp+1)nqp>q\iff\frac{n_p}{d_p+1}>\frac{n_q}{d_q+1}\iff n_p(d_q+1)>(d_p+1)n_q

gt? &[[&[[N.gt? (N.mul 3 N.++0) (N.mul N.++2 1)]]]] ⧗ Rational → Rational → Boolean

Other basic operators can be implemented just as easily:

# negates a rational number
# this could also be done in the denominator
negate &[[N.-1 : 0]] ⧗ Rational → Rational

-‣ negate

# adds two rational numbers

p+q=npdp+1+nqdq+1=np(dq+1)+(dp+1)nqdpdq+dp+dq+1p+q=\frac{n_p}{d_p+1}+\frac{n_q}{d_q+1}=\frac{n_p(d_q+1) + (d_p+1)n_q}{d_pd_q+d_p+d_q+1}

add &[[&[[p : q]]]] ⧗ Rational → Rational → Rational
    p N.add (N.mul 3 N.++0) (N.mul N.++2 1)
    q N.add (N.mul 2 0) (N.add 2 0)

…+… add

# subtracts two rational numbers

pq=npdp+1nqdq+1=np(dq+1)(dp+1)nqdpdq+dp+dq+1p-q=\frac{n_p}{d_p+1}-\frac{n_q}{d_q+1}=\frac{n_p(d_q+1) - (d_p+1)n_q}{d_pd_q+d_p+d_q+1}

sub &[[&[[p : q]]]] ⧗ Rational → Rational → Rational
    p N.sub (N.mul 3 N.++0) (N.mul N.++2 1)
    q N.add (N.mul 2 0) (N.add 2 0)

…-… sub

# multiplies two rational numbers

pq=npdp+1nqdq+1=npnqdpdq+dp+dq+1p\cdot q=\frac{n_p}{d_p+1}\cdot\frac{n_q}{d_q+1}=\frac{n_pn_q}{d_pd_q+d_p+d_q+1}

mul &[[&[[p : q]]]] ⧗ Rational → Rational → Rational
    p N.mul 3 1
    q N.add (N.mul 2 0) (N.add 2 0)

…⋅… mul

# calculates the multiplicative inverse of a rational number

p1=(npdp+1)1={undefinednp=0dp+1npnp>0(dp+1)npnp<0p^{-1}=\left(\frac{n_p}{d_p+1}\right)^{-1}=\begin{cases}\mathrm{undefined}&n_p=0\\\frac{d_p+1}{n_p}&n_p>0\\\frac{-(d_p+1)}{-n_p}&n_p<0\end{cases}

invert &[[N.compare-case eq gt lt 1 (+0)]] ⧗ Rational → Rational
    eq Ω
    gt N.++0 : N.--1
    lt N.-(N.++0) : N.--(N.-1)

~‣ invert

# divides two rational numbers

pq=pq1\frac{p}{q}=p\cdot q^{-1}

div [[1  ~0]] ⧗ Rational → Rational → Rational

…/… div

# approximates a rational number to n digits
approximate &[[[&[[int : float]] (N.quot-rem 2 N.++1)]]] ⧗ Rational → Number → (List Number)
    int 1
    float N.number→list (N.div (N.mul 0 (N.pow (+10) 2)) N.++3)

This is nothing special really and appears similar in various libraries like Haskell’s Rational (although somewhat obfuscated).

Reduction

Such libraries often do a reduction step that shortens fractions similar to how we do intuitively. These steps require a gcd\mathrm{gcd} (greatest common divisor) operation, which is currently quite expensive in bruijn1. Still, one could do something like this:

# finds the smallest equivalent representation of a rational number

p=npdp+1=np/gcd(np,dp+1)(dp+1)/gcd(np,dp+1)p=\frac{n_p}{d_p+1}=\frac{n_p/\mathrm{gcd}(n_p,d_p+1)}{(d_p+1)/\mathrm{gcd}(n_p,d_p+1)}

compress &[[[(N.div 2 0) : N.--(N.div N.++1 0)] (N.gcd 1 N.++0)]] ⧗ Rational → Rational

%‣ compress

# equivalent of Haskell's Fractional reduce function
# similar to compress, but without uncurrying
reduce [[[(N.div 2 0) : N.--(N.div 1 0)] (N.gcd 1 0)]] ⧗ Number → Number → Rational

# shortened (and much slower) add operation
add' &[[&[[reduce p q]]]] ⧗ Rational → Rational → Rational
    p N.add (N.mul 3 N.++0) (N.mul 1 N.++2)
    q N.mul N.++0 N.++2

All of this could already be used to do some pretty advanced calculations.

Example

a0=4.2an=(an1+0.1)2a4= ??? \begin{align*} a_0&=4.2\\ a_n&=(a_{n-1} + 0.1)^2\\ a_{4}&=\ ??? \end{align*}

λ pow-n = [L.nth-iterate (mul 0) (+1.0q)]
λ n = L.nth-iterate [pow-n (0 + (+0.1q)) (+2)] (+4.2q) (+4)
  2179006444280423418821440444946289062500000000000000/152587890625000000000000000000000000000000
  (~14280336633.23618300, time: 5.7e-2s, length: 1018)
λ approximate n (+5)
  {(+14280336633), (+2), (+3), (+6), (+1), (+8)}

Note: All REPL examples ran on my ThinkPad T14G1 AMD using bruijn’s HigherOrder reduction mode (ergo Haskell).

The real magic, of course, happens with real numbers!

Real

The missing insight I got from reddit is the fact that approximations of real numbers can be seen as the infinite limit of a function that maps natural to rational numbers. For example, for specific xRx\in\mathbb{R} there exists fx:NQf_x:\mathbb{N}\to\mathbb{Q}, such that

x=limnfx(n). x=\lim_{n\to\infty}f_x(n).

This fxf_x would typically be defined as a sequence that converges to the desired real number2.

Therefore we represent a real number as an abstraction over the rational numbers:

:import std/Pair .
:import std/Math/Rational Q

# if n is not irrational, just abstract the rational directly

fx:nxf_x:n\mapsto x


# new syntactic sugar with 'r' suffix!
:test ((+5.0r)) ([(+5) : (+0)])

# converts a balanced ternary number to a real number
number→real [[Q.number→rational 1]] ⧗ Number → Real

:test (number→real (+5)) ((+5.0r))

Approximating a real number to an arbitrarily precise rational number can then be done by applying some natural number. In theory we could use any base here, yet I still use balanced ternary. This wastes some space for negative numbers, but is best supported in bruijn’s standard library.

Operators

fx(42)=x+δf_x(42) = x+\delta

:test ((+4.2r) (+42)) ((+4.2q))

# returns true if two real numbers are approximately equal

x1=x2    limnfx1(n)=limnfx2(n)x_1=x_2\iff\lim_{n\to\infty}f_{x_1}(n)=\lim_{n\to\infty}f_{x_2}(n)

# one could further approximate by comparing the difference with some threshold
approx-eq? [[[Q.eq? (1 2) (0 2)]]] ⧗ Number → Real → Real → Boolean

# initialized with some large enough value
…≈?… approx-eq? (+42) ⧗ Real → Real → Boolean

:test ((+4.2r) ≈? (+4.2r)) (true)

This is one of the drawbacks of the implementation: It’s not possible to be 100% certain that two real numbers are equivalent! This is okay though, usecases for this are rare anyway.

For other binary operators I use a small combinator trick that “joins”/“lifts” the argument nn of both sides such that, inductively, we only ever need to apply a single number nn to any real number to get its approximation:

:import std/Combinator .

# to remind you:
# phoenix combinator / liftM2
φ [[[[3 (2 0) (1 0)]]]] ⧗ (b → c → d) → (a → b) → (a → c) → (a → d)

#    φ (Rational → Rational → Rational) (Number → Rational) (Number → Rational)
# ~> (Number → Rational) ~> Real

# adds two real numbers
add φ Q.add ⧗ Real → Real → Real

…+… add

# subtracts two real numbers
sub φ Q.sub ⧗ Real → Real → Real

…-… sub

# multiplies two real numbers
mul φ Q.mul ⧗ Real → Real → Real

…⋅… mul

# divides two real numbers
div φ Q.div ⧗ Real → Real → Real

…/… div

We can do something similar for the missing unary operators:

# to remind you:
# bluebird combinator / composition
b [[[2 (1 0)]]] ⧗ (b → c) → (a → b) → (a → c)

#    b (Rational → Rational) (Number → Rational)
# ~> (Number → Rational) ~> Real

# negates a real number
negate b Q.negate ⧗ Real → Real

-‣ negate

# inverts a real number
invert b Q.invert ⧗ Real → Real

~‣ invert

# finds smallest equivalent representation of a real number
compress b Q.compress ⧗ Real → Real

%‣ compress

Like with the rationals before, we can now go calculate some things!

:import std/List L

# power function: real^integer
pow-n [L.nth-iterate (mul 0) (+1.0r)] ⧗ Real → Number → Real

...

This is boring though. How can we calculate pow:RR\mathrm{pow}:\mathbb{R}\to\mathbb{R}? The trick is that for positive aa, ab=eblna.a^b=e^{b\ln a}.

Power

Therefore we only need to implement exp and ln.

# simple exp using infinite limit

ex=limn(1+xn)ne^x=\lim_{n\to\infty}(1+\frac{x}{n})^n

error:O(x2ex2n)\mathrm{error}: \mathcal{O}\left(\frac{x^2e^x}{2n}\right)

exp-lim [[pow-n ++(1 / (number→real 0)) 0]] ⧗ Real → Real

exp-lim converges very slowly. Since we want to combine exp with other (faster converging) formulas, we need to use a different technique using Taylor series. This is not as trivial:

# exp using Taylor series

ex=n=0xnn!e^x=\sum_{n=0}^\infty \frac{x^n}{n!}

error:O(xn+1(n+1)!)\mathrm{error}: \mathcal{O}\left(\frac{x^{n+1}}{(n+1)!}\right)

# the factorial and power is embedded into iteration
#   such that we don't always calculate them redundantly
exp [[L.nth-iterate &[[[[op]]]] start 0 [[[[1]]]] 0]] ⧗ Real → Real
    start [0 (+1.0r) (+1.0r) (+1.0r) (+1)]
    op [[[2 1 0 (4 + (1 / 0)) N.++3]] pow fac]
        pow 6  4
        fac 3  (number→real 1)

The iteration can be described by a recursive sequence:

numerator0=1.0denominator0=1.0result0=1.0iteration0=1numeratorn+1=numeratornxdenominatorn+1=denominatorniterationnresultn+1=resultn+numeratorn+1denominatorn+1iterationn+1=iterationn+1 \begin{align*} \mathrm{numerator}_0 &= 1.0\\ \mathrm{denominator}_0 &= 1.0\\ \mathrm{result}_0 &= 1.0\\ \mathrm{iteration}_0 &= 1\\ &\\ \mathrm{numerator}_{n+1} &= \mathrm{numerator}_n ⋅ x\\ \mathrm{denominator}_{n+1} &= \mathrm{denominator}_n ⋅ \mathrm{iteration}_n\\ \mathrm{result}_{n+1} &= \mathrm{result}_n + \frac{\mathrm{numerator}_{n+1}}{\mathrm{denominator}_{n+1}}\\ \mathrm{iteration}_{n+1} &= \mathrm{iteration}_n + 1 \end{align*}

The [[[[1]]]] then extracts the third value, resultn\mathrm{result}_n, of the 4-tuple start.

  *\ *\ *

ln is similarly complicated:

# natural logarithm using Taylor series

ln(x)=n=022n+1(x1x+1)2n+1\ln(x)=\sum_{n=0}^\infty\frac{2}{2n+1}(\frac{x-1}{x+1})^{2n+1}

error:O((x12)2n+1)\mathrm{error}: \mathcal{O}((\frac{x-1}{2})^{2n+1})

# again, the power is embedded into iteration (see start)
ln [[[L.nth-iterate &[[[op]]] start 1] (--1 / ++1 0) [[[1]]]]] ⧗ Real → Real
    start [0 1 (+0.0q) (+0)]
    op [0 pow (Q.add 2 go) N.++1]
        pow Q.mul 3 (Q.mul 4 4)
        go Q.mul ((+2) : (N.mul (+2) 1)) 3

# real power function: real^real
pow [[exp (0  (ln 1))]] ⧗ Real → Real → Real

…**… pow

You can already see a common problem: Different approximations converge with different speed (here notated as error in iteration step nn). Combining them into different terms can yield unbalanced ratios of convergence, in turn impacting performance of more precise calculations. I could not find a solution yet, so I just try to make most formulas converge roughly with the same speed.

Some example REPL usage, with precision n=10n=10:

λ exp (+1.0r) (+10)
  18099969098565397826764800000/6658606584104736522240000000
  (~2.7182818011463845, error of ~1.73e-8, time: ~0.03s, 632bit)
λ ln (+2.0r) (+10)
  233890488187863915541701543725460185763128364510093286440/337432647424663101659020986764900182676573608430256879075
  (~0.6931471805498117, error of ~1.01e-11, time: ~0.03s, length: 1051bit)
λ ((+2.0r) ** (+5.0r)) (+10)
  598079134110276147827564883538659081286604734987944789387834372340947075110943453837234618373552693185584716635684189972636644241652171376916943525486461022728058976842778463558150783565917373041554240375907383228706705686070473570707758765671846250321276278714131843233090967960950100066733563643600273605634465180162835902914836231104763335304562987169953913286974135123478137546638315618992875001966578776051314748446370950239017032303780045954608544272508524030682516609446564296943397047151770455802591948913236771550241140628206362089337306080729456884605657473736006900750062127183174575149591063756525402897186089732707298165462679140041446762061630523298302176019472132682528624760986533490881586211470148319049682387677590759715636301654414721058038524158625347730839117574417904657884025887070370668002870317851865610637261303431591432712751211868675681918441438181009236310163455557218485924474140212254727259954387270719434271208143299910616186908085375881165357827436854844415894159309606791014794031897082269022230707092560145786402326155316644252639477330086223550047352406410530249264625980127611616014612744110563965082672579163660755936886822130116400763646980101508854615001526358232971088757782801868634665475194968733501724384816076327046316070654670769403132424875856535131094386148667653955301369622336790381866119818316339256930323110721023022331415394093223237411495404394416682243161354096950645896692264278922066316269661717458816226041291616605553523205391825489419079585537561636737655730488488440777060126875525392906188334680281882510800224189036509975301509573328453196523002298919394859926114954914787468282500028398858785002799476737665131442186918408368665754947375031323023994205757187950886526437405495638105657127387780793248754925231819696930561066681593613019051902903759963710577600178126551532374230654018917071744634056944518994438021779780944427419530474781506626479456134411453643022372837732547369414691303401174003771686180236856810234420872033503572191499213696774228469644141715532392472118731282252789583127278677420087233498337120738433880137357508449801704906102451094822409824236222734766813101302782137436899378522610042953094252716889434203754695406049119679351289357796953897381003734393866725685601294639966545608927832846146076337329964777111812830183537916681463188543158146621245507172569192113334571731610610293877626003372555010124576440160375294328387964074220385933549089782641666077677368723655537628973398985260911246931635785708900043886771153987974930564246498556798449367380515140417080137186355779417955734278686914202430144848198772372684160394117644408417711866459815541733128505818698901992840753593067270178536796957546332531748299524902599475784658804350813184544953176440746444798346737996269153523375711529045417264644913503743443882289016496822187252992674242720739998853511069087517775444431645507009913242118489351376976411894580661365924273926755289179739713122277279765888522997890201654099882914400318594082414541185666356171548192583342325223014342107882723075838412735160971173131727596228301990777254104614257812500000000000000000000000000000000000000/74760779145070733867908932761413155762971594365211605875643688715795658071296489021978852498729185241962654371126198985856785158579859308675958564026242002487212659492193987786684537196431214280757975925264758942883634825715655210754259881854161850075555156727957799415478911599087787338759644279029501022468741585284712250870563608807772216626319992978677936038223391819445771394867916047740401921773908140195393144179369141320178364174064953698343981949644450414632898184304796423731597697465375566339947566525590816014253540526958504467447037286922629637764004487538446911081259943320301928694733654797538871623891971353728913964342443511026866560197170972549628045705440407603312893994442950853177415705720888410553006515572923881893088120078667252896080200863364262936348551786574482904929604768379106016250412338548222219088828801332785047621040579051084935367004882760987761820574788206148525328752068539097310214079835762720917029086377377340872933055196076185798222540307196361313399629934793229668735088328473360784158080026141115666746992410622458258295577175907820466408926620132964962486225515990672744837581085043326597011264706252949090969847322322495952154924286412491550689193904192548701376743857571103786655274607787837390201233431227051287558269947139643643354748534743176311469903335733573895699231869252769229007514872672526870460166704215242799424899976518534714924877197462712550853495447879709422302180307953400983205499497427437609334221946874236442601011590334150274580045185317079861788907309774171436652200887227409217065183012520489092932128741216386785887485237903535011967264310628081772805865418646841771828566494422331918960006779437847077838135379228767020888503346221487625606989561464766098979814616703315507210693772610145034125230859983749588645164959428343521078143981388069312327590834893754353340597878304329885934981702091972861233235112086702432657823805894511173586979592324513656316233277138188041726134730813327170727412918157800292547137553599968223539455477182881451302345840397262591228570188639075069091419714621221515228131751997428030553750909780200688042650201197354138784925439398522202355311797657482825417733561554769273059694184443577163344456922849979375864932918277480426922461006181579979595175551025678507687767005418040140775325648296021675347480784155313735196658541390475382786035917854927202267064655575298997606246563002418114754520547111452824170514469941660203289142179436597789664589790956742643067949813175447395938470322035555857103944544604786924840964930980034735034061609181845620607831460713770653133250923402551987297201485964824552332216241734485472613124394181347966082706956069628488992962113147532752637402860615037429791398408738208746234630119201878416063821300338439170456002649326583102180899337539911353548588857412013494440755404189179361325795357642589917024984308672607471573171033808418901434055060528412896195791352491112725007239345428334121420462709200860489873502414825861884547983847540875126400700037685065643942537861409047019979211343772995432655685736347752224500995410494397219736129045486450195312500000000000000000000000000000000000000
  (~7.9999050431205925, error of ~9.5e-5, time: ~49.2s, length: 54kbit)

The numbers getting that large is the implication of missing reduction of rationals as explained in the previous section, so yes, this definitely needs some improvements. The calculations then also get incredibly slow once multiple formulas are involved. This could potentially also be improved with better reducers (currently almost nothing is shared/parallelized).

Roots

Or, one could go a different route and just find a better formula that converges faster. For example, for sqrt and cbrt:

# square root via x^0.5
pow-sqrt \pow (+0.5r) ⧗ Real → Real

# real square root
# Newton's/Heron's method, quadratic convergence

xn+1=xn+axn2x_{n+1}=\frac{x_n+\frac{a}{x_n}}{2}

sqrt [[y [[[N.=?0 guess go]]] (1 0) 0]] ⧗ Real → Real
    guess (+1.0q)
    go [Q.div (Q.add 0 (Q.div 2 0)) (+2.0q)] (2 1 N.--0)

# real cube root
# Newton's/Heron's method, quadratic convergence

xn+1=3xn+axn24x_{n+1}=\frac{3x_n+\frac{a}{x_n^2}}{4}

cbrt [[y [[[N.=?0 guess go]]] (1 0) 0]] ⧗ Real → Real
    guess (+1.0q)
    go [Q.div (Q.add (Q.mul (+3.0q) 0) (Q.div 2 (Q.mul 0 0))) (+4.0q)] (2 1 N.--0)

Comparison:

λ pow-sqrt (+2.0r) (+7)
  89033367892848420325715129500331038667412886056368323573159431977310238895439218094980112725702840342993272134600096124077626686964733473513621053153947308252557325158480398438676319722582712899722911291367231026263483599926951994934434361052597331385548045264275022690705284987673823814184868974677551776642419208825170717874123852968703476348321320937471629780426469759321294332061735019180158114422592282352023648657975088457689473927068241499703325954672642420561944499393423107013113058488991800090967348390600054350006070192488093414571502040824332030927892137510876992387627104157234890750860521744152469018399002283691258127609527505782738117052145118472073066393963385416112881333398245475247711810954299587984415353368400124041917003009441829548789670468653120761682572674942660955340800000000000000000000000/62956098752232057387808350112841270755951480654727337171842489892946126359947885037380093400728241011317082442072604878261728100965154427911444757570897439621083387621109809746375368978651533367238019668081714567583330572204727730673367194112699820403536457063036501180427572233579383441923086911895434098875666345475230727812250574821804684019255992429473998197711015819847678940543489342616860950975760163509063079641664150888593726203230408507924318543814430554711334286943508216430250098869852834437753281169571221961391204635543888995125143434949553386638535183467482640471944118689238928081592507454317778459894864002701265511068264849228331714462408887817411373838522061660807386673012997199144341243681258637527781502269807089280818973625945249329852752982953473633532574022778880000000000000000000000000000000
  (~1.4142135497189112, error: ~1.26e-8, time: ~3.6s, length: 14kbit)
λ sqrt (+2.0r) (+7)
  4946041176255201878775086487573351061418968498177/3497379255757941172020851852070562919437964212608
  (~1.4142135623730951, error: ~2.89e-98, time: ~0.069s, length: 1071bit)

For x=2x=2 and n=7n=7 the previous method takes 50×50{\times} longer and is 12×12{\times} less precise! You can see how different formulas result in vastly different performance. This way, I believe, one could probably craft reasonably efficient implementations of many other mathematical functions.

Differentiation

Also fun:

# returns the derivative of a function

dydx=limh0f(x+h)f(x)h\frac{\mathrm dy}{\mathrm dx}=\lim_{h\to0}\frac{f(x+h)-f(x)}{h}

# (so much in that excellent formula!!1)
derive [[[[((3 (2 + 0)) - (3 2)) / 0] ((+1.0r) / conv) 0]]] ⧗ (Real → Real) → (Real → Real)
    conv number→real (N.pow (+2) 0) # custom convergence ratio

# e.g. derive sqrt (+3.0r) n, n→∞

     123120.2287\ \ \ \ \ \leadsto\frac12\cdot3^{-\frac12}\approx0.2287

In REPL:

λ derive sqrt (+3.0r) (+5)
  4548499328554333130253458031438461086963398946474394806600894005568776844296474084125321695485140926464/15756464483850753172492714095995675761947988863890026857797915483539179836951188627907703817126616760320
  (~0.2886751233575412, error: ~1.12e-8, time: ~0.6s, length: 2201bit)

The results could be improved by using a better adapted convergence ratio in derive, although this would also imply much longer results.

Constants

I’ve already shown you how to calculate ee using the Taylor series of exe^x. Other interesting constants:

π2=n=02nn!2(2n+1)!\frac{\pi}2=\sum_{n=0}^\infty\frac{2^n n!^2}{(2n+1)!}

π/2 [L.nth-iterate &[[[[[op]]]]] start 0 [[[[[3]]]]]] ⧗ Real
    start [0 (+1) (+0.0q) (+1) (+1) (+1)]
    op [0 N.++5 (Q.add 4 ((N.mul 3 2) : N.--1)) num-pow num-fac denom]
        num-pow N.mul 3 (+2)
        num-fac N.mul 2 (N.mul 5 5)
        denom [N.mul 2 (N.mul 0 N.++0)] (N.mul (+2) 5)

# ratio of circle's circumference to its diameter
π π/2  (+2.0r) ⧗ Real

# golden ratio from direct formula
φ ++(sqrt (+5.0r)) / (+2.0r) ⧗ Real

# conjugate golden ratio
ψ -(~φ) ⧗ Real

# golden ratio from fibonacci convergence

φ=limnFn+1Fn\varphi=\lim_{n\to\infty}\frac{F_{n+1}}{F_n}

φ* ++[(L.nth-iterate &[[0 : (N.add 1 0)]] ((+0) : (+1)) 0) [[1 : N.--0]]] ⧗ Real

# real fibonacci becomes trivial

Fn=φn(1φ)n5F_n=\frac{\varphi^n-(1-\varphi)^n}{\sqrt5}

fib [((pow φ 0) - (pow ψ 0)) / (sqrt (+5.0r))] ⧗ Real → Real

Results:

λ π (+10)
  58161346883519614707939305922472528176805953649166549778432000000000000000/18519311967783166520148611995625183290053887455979435458560000000000000000
  (~3.140578169680337, error: ~0.001, time: ~0.141s, length: 1511bit)
λ π (+20)
  23770718233993033585696255811682815853064495117892501617611447251371134322557383888051869648991848387499657110718132777488467529319755142727319452657262872133986736812397012330767068148003469821569121725937860722210965191933852585434483132680685165347732376029802802822656597172962731748393261075522609795620777191496744960000000000000000000000000000000000000000000000000000000000000000000000000000000/7566456363918636918448056774706264206207563561731028263307703141167935514660420425628485969667487374826285184836250356538021409394395657113901240231996338685382563709519924768061981073010130735334419464639463771932114561342847927309380840121890195187487927033515312915490616963923270255231275266462271264256491368166195200000000000000000000000000000000000000000000000000000000000000000000000000000000
  (~3.141591927675147, error: ~7.26e-7, time: ~8.3s, length: 8089bit)
λ φ (+7)
  138598980290857213752199295710131342640539745689910576085702017024/85658880625826545040363067218944712997310200964513375659560534016
  (~1.618033988749895, error: 7.06e-54, time: ~0.17s, length: 1418bit)
λ # another example of the convergence rate problem:
λ φ* (+500)
  1206484255615496768210420703829205488386909032955899056732883572731058504300529011051/745648276861993189984204820755333242001144560869101973859680384188513887852280667476
  (~1.618033988749895, error: 5.12e-85, time: ~0.17s, length: 1757bit)

Also see std/Math for a list-based approximation of π, it’s actually more performant.

Trigonometry

Some trigonometric functions (kinda arbitrarily chosen):

# hypotenuse
hypot [[sqrt ((0  0) + (1  1))]] ⧗ Real → Real → Real

# arctan by Taylor series

arctan(x)=n=0(1)nx2n+12n+1,x1\arctan(x)=\sum_{n=0}^\infty(-1)^n \frac{x^{2n+1}}{2n+1},\quad |x|\le1

atan* [[[L.nth-iterate &[[[[op]]]] start 1] (1 0) [[[[3]]]]]] ⧗ Real → Real
    start [0 1 [[0]] (Q.pow-n 1 (+3)) (+3.0q)]
    op [0 ((3 Q.add Q.sub) 4 (Q.div 2 1)) \3 num denom]
        num Q.mul 2 (Q.mul 5 5)
        denom Q.add 1 (+2.0q)

# actual arctan for arbitrary x
atan [[[Q.lt? Q.|0 (+1.0) (atan* [1] 1) go] (1 0)]] ⧗ Real → Real
    go Q.sub π/2-pm conj-atan
        π/2-pm Q.mul (Q.>?0 (+1.0q) (-1.0q)) (π/2 1)
        conj-atan atan* [Q.div (+1.0q) 1] 1

# 2-argument arctan
# this one's really weird
atan2 [[[[[Q.add a (Q.mul b (Q.add c d))]] (2 0) (1 0)]]] ⧗ Real → Real → Real
    z (+0.0q)
    a Q.=?0 z (atan [Q.div 2 1] 2)
    b Q.sub (+1.0q) (Q.<?1 (+2.0q) z)
    c Q.<?0 (π 2) z
    d Q.=?0 (π/2 2) z

Formulas like sin and cos are similar to atan and are left as an exercise to the reader.

Methods

There are many ways to translate real approximations like Taylor series to lambda calculus. Depending on your needs, my method using L.nth-iterate may not be the optimal approach.

Unary

If you use unary/Church numerals instead, you actually get repeated iteration for free! A Church numeral is [[1 (1 (1 (... (1 0))))]], so therefore we can just apply the iteration function as first argument and the starting value as second. The previous exp implementation could easily be translated:

# unary exp using Taylor series

n=0xnn!\sum_{n=0}^\infty \frac{x^n}{n!}

unary-exp [[0 &[[[[[op]]]]] start [[[[1]]]]]] ⧗ Number → (Unary → Rational)
    start [0 (+1) (+1) (+1.0q) (+1)]
    op [[2 1 0 (Q.add 4 (1 : N.--0)) N.++3]] pow fac
        pow N.mul 6 4
        fac N.mul 3 1

Things like list indexing are also easier with Church numerals. However, they also take more space, since they’re not only O(n)\mathcal{O}(n), but also because they literally duplicate the iteration function nn times! They could still be a good option for rapidly converging series. Also, if for whatever reason you’d want to golf an exp function or similar, you could use unary numerals at every step starting with defining integers as a pair of two unary numerals. unary-exp could then get truly tiny.

Recurrence

I’ve demonstrated some conversions from Taylor series to recursive sequences. In theory, such recurrence relations could be resolved using a variadic fixed-point combinator.

In my experience, the use of this combinator doesn’t give performance improvements over normal iteration using y-recursion (like I’ve used here) and can be quite confusing. Still, L.y* can be fairly small, so this could also be a good opportunity for golfing.

Complex

Complex numbers are, contrary to their name, quite simple. They are just a pair of real numbers!

C={a+bia,bR}\mathbb{C}=\{a+bi\mid a,b\in\mathbb{R}\}

Still, we need to lift the approximation arguments of the two real numbers to the outside, such that we only ever need to apply a single natural number:

:import std/Pair .
:import std/Combinator .
:import std/Math/Real R
:import std/Math/Rational Q

# new syntactic sugar with complex part as suffix!
:test ((+0.0+1.0i)) ([((+0.0r) 0) : ((+1.0r) 0)])

# imaginary unit
i (+0.0+1.0i)

# converts a balanced ternary number to a complex number
number→complex [[(R.number→real 1 0) : ((+0.0r) 0)]] ⧗ Number → Complex

:test (number→complex (+5)) ((+5.0+0.0i))

We can then extract the real/imaginary part using traditional pair selectors mixed with a lifted argument, therefore resembling a real number:

# returns real part of a complex number
real [[1 0 [[1]]]] ⧗ Complex → Real

:test (real (+5.0+2.0i)) ((+5.0r))

# returns imaginary part of a complex number
imag [[1 0 [[0]]]] ⧗ Complex → Real

:test (imag (+5.0+2.0i)) ((+2.0r))

Operators

Here we use φ and b again to lift the approximator argument. After their application, the real numbers inside the complex number will actually appear rational, such that we can use rational operations directly.

# adds two complex numbers

p+q=(ap+bpi)+(aq+bqi)=(ap+aq)+(bp+bq)ip+q=(a_p+b_pi)+(a_q+b_qi)=(a_p+a_q)+(b_p+b_q)i

add φ &[[&[[(Q.add 3 1) : (Q.add 2 0)]]]] ⧗ Complex → Complex → Complex

…+… add

# subtracts two complex numbers

pq=(ap+bpi)(aq+bqi)=(apaq)+(bpbq)ip-q=(a_p+b_pi)-(a_q+b_qi)=(a_p-a_q)+(b_p-b_q)i

sub φ &[[&[[(Q.sub 3 1) : (Q.sub 2 0)]]]] ⧗ Complex → Complex → Complex

…-… sub

# multiplies two complex numbers

pq=(ap+bpi)(aqbqi)=(apaqbpbq)+(apbq+bpaq)ipq=(a_p+b_pi)(a_q-b_qi)=(a_pa_q-b_pb_q)+(a_pb_q+b_pa_q)i

mul φ &[[&[[p : q]]]] ⧗ Complex → Complex → Complex
    p Q.sub (Q.mul 3 1) (Q.mul 2 0)
    q Q.add (Q.mul 3 0) (Q.mul 2 1)

…⋅… mul

# divides two complex numbers

pq=ap+bpiaqbqi=apaq+bpbqaq2+bq2+bpaqapbqaq2+bq2i\frac{p}{q}=\frac{a_p+b_pi}{a_q-b_qi}=\frac{a_pa_q+b_pb_q}{a_q^2+b_q^2}+\frac{b_pa_q-a_pb_q}{a_q^2+b_q^2}i

div φ &[[&[[p : q]]]] ⧗ Complex → Complex → Complex
    p Q.div (Q.add (Q.mul 3 1) (Q.mul 2 0)) (Q.add (Q.mul 1 1) (Q.mul 0 0))
    q Q.div (Q.sub (Q.mul 2 1) (Q.mul 3 0)) (Q.add (Q.mul 1 1) (Q.mul 0 0))

…/… div

# negates a complex number

p=(ap+bpi)=ap+bpi-p=-(a_p+b_pi)=-a_p+-b_pi

negate b &[[(Q.negate 1) : (Q.negate 0)]] ⧗ Complex → Complex

-‣ negate

# inverts a complex number

1p=1ap+bpi=apap2+bp2+bpap2+bp2i\frac{1}{p}=\frac{1}{a_p+b_pi}=\frac{a_p}{a_p^2+b_p^2}+\frac{b_p}{a_p^2+b_p^2}i

invert b &[[p : q]] ⧗ Complex → Complex
    p Q.div 1 (Q.add (Q.mul 1 1) (Q.mul 0 0))
    q Q.div 0 (Q.add (Q.mul 1 1) (Q.mul 0 0))

~‣ invert

Pie

Several of the defined functions can now be combined to (attempt to) calculate eπie^{\pi i}:

# very similar to the real one!

ex=n=0xnn!e^x=\sum_{n=0}^\infty \frac{x^n}{n!}

exp [[L.nth-iterate &[[[[op]]]] start 0 [[[[1]]]] 0]] ⧗ Complex → Complex
    start [0 (+1.0+0.0i) (+1.0+0.0i) (+1.0+0.0i) (+1)]
    op [[[2 1 0 (4 + (1 / 0)) N.++3]] pow fac]
        pow 6  4
        fac 3  (number→complex 1)

# complex natural logarithm
ln [[[p : q] (1 0)]] ⧗ Complex → Complex
    p R.ln (&[[R.hypot [2] [1]]] 0) 1
    q &[[\R.atan2 [2] [1]]] 0 1

# complex power function: complex^complex
pow [[exp (0  (ln 1))]] ⧗ Complex → Complex → Complex

…**… pow

In REPL:

λ ln (+2.0+4.2i) (+3)
  131333047804453819160134580326741441990719587929624714214939549619500790973864294006194984228668428249332105635341198746982341517036711701337460293632000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000/83484794802006666694721421339938768091461941633007787921780106775988824518415759807105367379415841350472645076302765692152680960498897197586158714880000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
  : 110569850108556581649876480/108150992982294090821067600
  (~1.57313734+1.02236556i, error: ~0.03+0.1i, time: ~1.5s, length: 6376bit)
λ # ugh, numerator converges a lot faster than denominator, anyway...

λ (+2.0+0.0i) ** (+3.0+0.0i) (+2)
  46279959818696076004546197821636492728892381317705253201507013980345318470977416794052934030379365364822245947734236637528515350739555014550295902286727588227228815359769649963518843577958400000000000000000000000000000000000000000000000000000000000000000000000000000000/8488066053046236477811349396282461894378738844356320153510331083656243451009448569006437421681264652074166946504077170970137877639991273678242184865776333754718651791926832629500105614950400000000000000000000000000000000000000000000000000000000000000000000000000000000
  : 0/8488066053046236477811349396282461894378738844356320153510331083656243451009448569006437421681264652074166946504077170970137877639991273678242184865776333754718651791926832629500105614950400000000000000000000000000000000000000000000000000000000000000000000000000000000
  (~5.45235623+0i, error: ~2.58+0i, time: ~2.9s, length: 10740bit)
λ # yes, many abstractions make everything worse. Looking at you, Electron!

λ exp ([(R.π 0) : (+0.0)]  i) (+4)
  -67479993202306472714667916921031425556099775398952539611521667570250447344805122008169714475335413570680821954970430751292829585907152642903675057639798448289757605848159425094599290482701423963369343202081203932066765475487409153721693127606181916420307244302193597548169505958661717082177471539283605044618375306864059617331387590135051834499754896430296524931141329329805419250287631708024288329733458404003715292655399963717259776754256684136966566721257811784729486080081920000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000/1361884868862060097549690181333056765046439560624409926616798792493328435787994562240567391560867539517753322628725069314972581967222209672892486307189350239216889094360499479590034032207920263636241298759966234243205570905766386207577079921475443490332191343090551009522870423225966591704205283036225548377573722514188921626990790748596113579717060293956259137133760787580672566454929838687772488012891423833655251357457695563011145177240090717870704091720599864440759816908963840000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
  -2274464885548702030965181024866063200451887687181824158478906668527172710670660032680499836893937077839821263721442724551514516189004446201713686954381666886213283445931943125399617719268887457535960704836138403659962512688606713459006461086148384909887833506799178160812769954205264005149736354349093526426930883081681818503846693688041596303454328101341483063852229800372285465338831102296823361184831086127545753265883348598917220466614124079804220921276889958249694255644672000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000/1361884868862060097549690181333056765046439560624409926616798792493328435787994562240567391560867539517753322628725069314972581967222209672892486307189350239216889094360499479590034032207920263636241298759966234243205570905766386207577079921475443490332191343090551009522870423225966591704205283036225548377573722514188921626990790748596113579717060293956259137133760787580672566454929838687772488012891423833655251357457695563011145177240090717870704091720599864440759816908963840000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
  (~-0.05-1.67i, error: 1.05+1.57i, time: ~19s, length: 27219bit) 
λ exp ([(+3.14159265358) : (+0.0)]  i) (+4)
  5067240574827142493824884283487389404038327255808591266950606056410685157365136888074630405753850936889648437500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000/40894549326657251992544008824259638146259021596051752567291259765625000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
  : -82857269451682879813602954637549515087487631913625936626050655364750241460569668561220169067382812500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000/40894549326657251992544008824259638146259021596051752567291259765625000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
  (~0.12-2.02i, error: 0.88+1.57i, time: ~80s, length: 46851bit)
λ # oof

Sadly, these results are really disappointing. Unfortunately the complex functions converge rather slowly and became a victim of many layered abstractions and recursions.

Also, if you’re interested: The unevaluated REPL term with R.π (after complete substitution of definitions) has a BLC length of 1,106,776bit – that’s around 138kB! Compressing the term with my BLoC tool results in a size of 704B. Many duplicated terms could obviously be shared by abstraction directly, the large size doesn’t seem to impact performance significantly, though.

  *\ *\ *

Otherwise the techniques shown for real numbers apply here, too, so I won’t demonstrate them again. However, feel free to apply your gained knowledge by PRing new functions to bruijn’s standard library!

Lessons

In my experience many readers won’t immediately understand, but working with bruijn and pure lambda calculus in general is really satisfying and enjoyable. I have yet to find another language that allows the same freedom and elegancy. I can only recommend digging into all the definitions and trying to understand everything – it’s definitely possible and rewarding.

Especially in the end, we’ve seen that this technique comes with large performance disadvantages. If someone wanted to use this for actual programs, they would need to invest a huge amount of work to make all of these functions even remotely efficient and accurate. This article should mostly serve as a fun proof of concept and notebook. I’ll probably still work on some obvious performance and accuracy improvements – you will find them on GitHub.

Aside from showing bruijn’s weaknesses, this project also showed me how fast bruijn can actually be. If you think about it, there are literally millions of reductions happening for the last formulas and the result still appears within a few seconds! Still, having to optimize even the smallest operations is something I don’t experience in normal programming languages. If I were to implement something similar in a language like Python or Java, decomposing powers or sqrt into the iteration could even hurt performance, since the traditional functions are already hyper-optimized.

Here, we can’t trust any optimizations and need to consider the usefulness of every additional abstraction and, especially, layered y-recursion – it seems like I overestimated the performance of fixed-point recursion a lot. While some of this could be improved with specially optimized reducers, I think this goes against the spirit of the project and language. For me, it’s more interesting trying to implement procedures that run efficiently on off the shelf reducers. I’m obviously not there yet, though.

Thanks for reading. Contact me via email. Support on Ko-fi. Subscribe on RSS. Follow on Mastodon. Discuss on hacker news. Program in bruijn.

You can find the source of every described function (and a lot more) here or on GitHub.


Imprint · AI Statement · RSS