Math.bruijn

# MIT License, Copyright (c) 2022 Marvin Borner
# experimental functions; sometimes list-based; could work on any base

:input std/Number

:import std/List L

# adds all values in list
sum L.foldl add (+0) ⧗ (List Number) → Number

∑‣ sum

:test (∑((+1) : ((+2) : L.{}(+3)))) ((+6))

# digit sum of all values
digit-sum sum ∘ number→list ⧗ Number → Number

:test ((digit-sum (+0)) =? (+0)) (true)
:test ((digit-sum (+10)) =? (+1)) (true)
:test ((digit-sum (+19)) =? (+10)) (true)

# returns max value of list
lmax L.foldl1 max ⧗ (List Number) → Number

:test (lmax ((+1) : ((+3) : L.{}(+2)))) ((+3))

# returns min value of list
lmin L.foldl1 min ⧗ (List Number) → Number

:test (lmin ((+2) : ((+1) : L.{}(+0)))) ((+0))

# list from num to num
{…→…} z [[[rec]]] ⧗ Number → Number → (List Number)
	rec 1 =? ++0 case-end case-list
		case-list 1 : (2 ++1 0)
		case-end L.empty

:test ({ (+0) → (+2) }) ((+0) : ((+1) : L.{}(+2)))

# equivalent of mathematical sum function
∑…→…|… z [[[[[rec]]]]] (+0) ⧗ Number → Number → (Number → Number) → Number
	rec 2 =? ++1 case-end case-sum
		case-sum 4 (3 + (0 2)) ++2 1 0
		case-end 3

:test (∑ (+1) → (+3) | ++‣) ((+9))

# multiplies all values in list
product L.foldl mul (+1) ⧗ (List Number) → Number

∏‣ product

:test (∏((+1) : ((+2) : L.{}(+3)))) ((+6))

# equivalent of mathematical product function
∏…→…|… z [[[[[rec]]]]] (+1) ⧗ Number → Number → (Number → Number) → Number
	rec 2 =? ++1 case-end case-sum
		case-sum 4 (3 ⋅ (0 2)) ++2 1 0
		case-end 3

:test (∏ (+1) → (+3) | ++‣) ((+24))

# greatest common divisor using repeated subtraction
gcd* z [[[1 =? 0 case-eq (1 >? 0 case-gre case-les)]]] ⧗ Number → Number → Number
	case-eq 1
	case-gre 2 (1 - 0) 0
	case-les 2 1 (0 - 1)

# greatest common divisor using modulo (mostly faster than gcd*)
gcd z [[[=?0 1 (2 0 (1 % 0))]]] ⧗ Number → Number → Number

:test ((gcd (+2) (+4)) =? (+2)) (true)
:test ((gcd (+10) (+5)) =? (+5)) (true)
:test ((gcd (+3) (+8)) =? (+1)) (true)

# least common multiple using gcd
lcm [[=?1 1 (=?0 0 |(1 / (gcd 1 0) ⋅ 0))]] ⧗ Number → Number → Number

:test ((lcm (+12) (+18)) =? (+36)) (true)
:test ((lcm (+42) (+25)) =? (+1050)) (true)

# power function
pow […!!… (iterate (…⋅… 0) (+1))] ⧗ Number → Number → Number

…**… pow

:test (((+2) ** (+3)) =? (+8)) (true)

# modulo exponentiation
pow-mod [[[(f (2 % 0) 1 (+1)) % 0]]] ⧗ Number → Number → Number → Number
	f y [[[[=?1 0 rec]]]]
		rec 3 (2 ⋅ 2 % 4) /²1 (=²?1 0 (2 ⋅ 0 % 4))

:test ((pow-mod (+2) (+3) (+5)) =? (+3)) (true)

# power function using ternary exponentiation (TODO: fix, wrong..)
pow* z [[[rec]]] ⧗ Number → Number → Number
	rec =?0 case-end case-pow
		case-pow =?(lst 0) ³(2 1 /³0) (³(2 1 /³0) ⋅ 1)
			³‣ [0 ⋅ 0 ⋅ 0]
		case-end (+1)

# factorial function
fac [∏ (+1) → 0 | i] ⧗ Number → Number

:test ((fac (+3)) =? (+6)) (true)

# super factorial function
superfac [∏ (+1) → 0 | fac] ⧗ Number → Number

:test ((superfac (+4)) =? (+288)) ([[1]])

# hyper factorial function
hyperfac [∏ (+1) → 0 | [0 ** 0]] ⧗ Number → Number

:test ((hyperfac (+2)) =? (+4)) (true)
:test ((hyperfac (+3)) =? (+108)) (true)
:test ((hyperfac (+4)) =? (+27648)) ([[1]])

# alternate factorial function
altfac y [[=?0 0 ((fac 0) - (1 --0))]]

:test ((altfac (+3)) =? (+5)) ([[1]])

# exponential factorial function
expfac y [[(0 =? (+1)) 0 (0 ** (1 --0))]]

:test ((expfac (+4)) =? (+262144)) ([[1]])

# inverse factorial function
invfac y [[[compare-case 1 (2 ++1 0) (-1) 0 (∏ (+0) → --1 | ++‣)]]] (+0)

:test ((invfac (+1)) =? (+0)) ([[1]])
:test ((invfac (+2)) =? (+2)) ([[1]])
:test ((invfac (+120)) =? (+5)) ([[1]])
:test ((invfac (+119)) =? (-1)) ([[1]])

# calculates a powertower
# also: [[foldr pow (+1) (replicate 0 1)]]
powertower z [[[rec]]] ⧗ Number → Number → Number
	rec =?0 case-end case-rec
		case-end (+1)
		case-rec 1 ** (2 1 --0)

:test ((powertower (+2) (+1)) =? (+2)) (true)
:test ((powertower (+2) (+2)) =? (+4)) (true)
:test ((powertower (+2) (+3)) =? (+16)) (true)
:test ((powertower (+2) (+4)) =? (+65536)) (true)

# knuth's up-arrow notation
# arrow count → base → exponent
arrow z [[[[rec]]]] ⧗ Number → Number → Number → Number
	rec =?2 case-end case-rec
		case-end 1 ⋅ 0
		case-rec L.foldr (3 --2) 1 (L.replicate --0 1)

:test ((arrow (+1) (+1) (+1)) =? (+1)) (true)
:test ((arrow (+1) (+2) (+4)) =? (+16)) (true)
:test ((arrow (+2) (+2) (+4)) =? (+65536)) (true)

# fibonacci sequence
# TODO: faster fib?
fibs L.map L.head (L.iterate &[[0 : (1 + 0)]] ((+0) : (+1))) ⧗ (List Number)

fib [L.index fibs ++0] ⧗ Number

:test (fib (+5)) ((+8))

# floored integer square root using Babylonian method
sqrt [z [[[[rec]]]] (+1) 0 0] ⧗ Number → Number
	rec (1 >? 2) case-rec case-end
		case-rec [4 (1 / 0) 0 1] /²(2 + 1)
		case-end 1

:test ((sqrt (+0)) =? (+0)) (true)
:test ((sqrt (+1)) =? (+1)) (true)
:test ((sqrt (+2)) =? (+1)) (true)
:test ((sqrt (+5)) =? (+2)) (true)
:test ((sqrt (+9)) =? (+3)) (true)

# integer logarithm
# TODO: could we somehow use the change-of-base rule and efficient log3?
log z [[[[rec]]]] (+1) ⧗ Number → Number → Number
	rec [((3 ≤? 1) ⋀? (1 <? 0)) case-end case-rec] (2 ⋅ 1)
		case-end (+0)
		case-rec ++(4 0 2 1)

:test ((log (+2) (+1)) =? (+0)) (true)
:test ((log (+2) (+2)) =? (+1)) (true)
:test ((log (+2) (+3)) =? (+1)) (true)
:test ((log (+2) (+4)) =? (+2)) (true)
:test ((log (+2) (+32)) =? (+5)) (true)
:test ((log (+2) (+48)) =? (+5)) (true)

# iterated logarithm
# note that log! 1 is defined as 1
log! [z [[rec]] --0] ⧗ Number → Number
	rec (0 ≤? (+1)) case-end case-rec
		case-end (+1)
		case-rec ++(1 (log (+2) 0))

:test ((log! (+1)) =? (+1)) (true)
:test ((log! (+2)) =? (+1)) (true)
:test ((log! (+3)) =? (+2)) (true)
:test ((log! (+4)) =? (+2)) (true)
:test ((log! (+5)) =? (+3)) (true)
:test ((log! (+16)) =? (+3)) (true)
:test ((log! (+17)) =? (+4)) (true)
:test ((log! (+65536)) =? (+4)) (true)
:test ((log! (+65537)) =? (+5)) (true)

# pascal triangle
# TODO: something is wrong in here
pascal L.iterate [L.zip-with …+… (L.{}(+0) ++ 0) (0 ; (+0))] (L.{}(+1))

# characteristic prime sequence by Tromp
characteristic-primes ki : (ki : (sieve s0)) ⧗ (List Bool)
	sieve y [[k : ([(2 0) (y' 0)] (ssucc 0))]]
		y' [[0 0] [1 (0 0)]]
		ssucc [[[[1 : (0 (3 2))]]]]
	s0 [[[ki : (0 2)]]]

# prime number sequence
primes L.map fst (L.filter snd (enumerate characteristic-primes)) ⧗ (List Number)

# slower but cooler prime number sequence
primes* L.nub ((…≠?… (+1)) ∘∘ gcd) (L.iterate ++‣ (+2)) ⧗ (List Number)

# prime factors
factors \divs primes ⧗ Number → (List Number)
	divs y [[&[[&[[3 ⋅ 3 >? 4 case-1 (=?0 case-2 case-3)]] (quot-rem 2 1)]]]]
		case-1 4 >? (+1) {}4 empty
		case-2 3 : (5 1 (3 : 2))
		case-3 5 4 2

# π as a list of decimal digits
# translation of unbounded spigot algorithm by Jeremy Gibbons
# TODO: faster!
#     → BBP/Bellard's formula with ternary base?
#       TODO: |log|, better primes/mod/div
π y [[[[[calc]]]]] (+1) (+180) (+60) (+2) ⧗ (List Number)
	calc [[0 : (6 q r t ++2)]] a b
		a ↑⁰(↑⁺0 ⋅ (↑⁰0 + (+2)))
		b (3 ⋅ ↑⁰(↑⁻(↑⁻0)) + ((+5) ⋅ 2)) / ((+5) ⋅ 1)
		q (+10) ⋅ 5 ⋅ 2 ⋅ --((+2) ⋅ 2)
		r (+10) ⋅ 1 ⋅ (5 ⋅ ((+5) ⋅ 2 - (+2)) + 4 - (0 ⋅ 3))
		t 3 ⋅ 1