path: root/std/Math/Real.bruijn

# some ideas by u/DaVinci103
# MIT License, Copyright (c) 2024 Marvin Borner

:import std/Logic .
:import std/Combinator .
:import std/Pair .
:import std/Math/Rational Q
:import std/Math N

# a Real is just a Number → Rational!

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

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

# returns true if two real numbers are equal approximately
approx-eq? [[[Q.eq? (1 2) (0 2)]]] ⧗ Number → Real → Real → Boolean

# TODO: bigger value??
…≈?… approx-eq? (+100)

# returns true if a real number is greater than another
gt? [[[Q.gt? (1 2) (0 2)]]] ⧗ Number → Real → Real → Boolean

…>?… gt?

# returns true if a real number is less than another
lt? [[[Q.lt? (1 2) (0 2)]]] ⧗ Number → Real → Real → Boolean

…<?… lt?

# returns true if a real number is positive
positive? [\(gt? 0) (+0.0r)] ⧗ Number → Real → Boolean

>?‣ positive?

# returns true if a real number is negative
negative? [\(lt? 0) (+0.0r)] ⧗ Number → Real → Boolean

<?‣ negative?

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

…+… add

:test ((+1.0r) + (+0.0r) ≈? (+1.0r)) (true)
:test ((+0.0r) + (-1.0r) ≈? (-1.0r)) (true)

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

…-… sub

:test ((+1.0r) - (+0.5r) ≈? (+0.5r)) (true)
:test ((+0.0r) - (-1.0r) ≈? (+1.0r)) (true)

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

…⋅… mul

:test ((+5.0r) ⋅ (+5.0r) ≈? (+25.0r)) (true)
:test ((+1.8r) ⋅ (+1.2r) ≈? (+2.16r)) (true)

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

…/… div

:test ((+8.0r) / (+4.0r) ≈? (+2.0r)) (true)
:test ((+18.0r) / (+12.0r) ≈? (+1.5r)) (true)

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

-‣ negate

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

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

~‣ invert

:test (~(+0.5r) ≈? (+2.0r)) (true)
:test (~(-0.5r) ≈? (-2.0r)) (true)

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

%‣ compress

:test (%[(+4) : (+1)] ≈? (+2.0r)) (true)
:test (%[(+4) : (+7)] ≈? (+0.5r)) (true)

# increments a real number
inc add (+1.0r) ⧗ Real → Real

++‣ inc

# decrements a real number
dec \sub (+1.0r) ⧗ Real → Real

--‣ dec

# ---

:import std/List L

# e^x using Taylor expansion
# tex: \sum_{n=0}^\infty \frac{x^n}{n!}
unary-exp [[0 &[[[[[[[2 1 0 (Q.add 4 (1 : N.--0)) N.++3]] pow fac]]]]] start [[[[1]]]]]] ⧗ Unary → Real
	pow N.mul 6 4
	fac N.mul 1 3
	start [0 (+1) (+1) (+1.0q) (+1)]

# equivalent to unary-exp but with ternary index using infinite list iteration
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 (number→real 1) ⋅ 3

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

# e^x using infinite limit
# tex: \lim_{n\to\infty}(1+x/n)^n
lim-exp [[pow-n [(N.add 2 1) : N.--1] 0]] ⧗ Number → Real

# natural logarithm using Taylor expansion
# tex: \sum_{n=0}^\infty\frac{2}{2n+1}(\frac{x-1}{x+1})^{2n+1}
# error: O(((x-1)/2)^{2n+1})
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

:test (Q.eq? (ln (+2.0r) (+2)) ((+168) : (+242))) (true)

derive [[[[((3 (2 + 0)) - (3 2)) / 0] ((+1.0r) / conv) 0]]] ⧗ (Real → Real) → (Real → Real)
	conv number→real (N.pow (+2) (N.mul (+2) 0))

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

…**… pow

# square root by x^{0.5}
sqrt* \pow (+0.5r) ⧗ Real → Real

# Newton's/Heron's method, quadratic convergence
# tex: x_{n+1}=\frac{x_n+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)

# Newton's/Heron's method, quadratic convergence
# tex: x_{n+1}=\frac{3x_n+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)

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

# tex: \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

# Gauss-Legendre, quadratic convergence
# Chudnovsky would be even faster but is ugly (i.e. magic numbers)
π-gauss [L.nth-iterate &[[[[op]]]] start 0 final]
	start [0 (+1.0q) Q.~(sqrt (+2.0r) 1) (+0.25q) (+1.0q)]
	op [[1 0 b t p] a]
		a Q.div (Q.add 4 3) (+2.0q)
		b sqrt [Q.mul 6 5] 6
		t Q.sub 3 (Q.mul 2 (Q.pow-n (Q.sub 5 0) (+2)))
		p Q.mul (+2.0q) 2
	final [[[[Q.div (Q.pow-n (Q.add 3 2) (+2)) (Q.mul (+4.0q) 1)]]]]

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

# conjugate golden ratio
ψ -(~φ)

# golden ratio from fibonacci convergence
φ* ++[(L.nth-iterate &[[0 : (N.add 1 0)]] ((+0) : (+1)) 0) [[1 : N.--0]]]

# real fibonacci
fib [((pow φ 0) - (pow ψ 0)) / (sqrt (+5.0r))]

# arctan by Taylor expansion, only for |x|<=1
# tex: \sum_{n=0}^\infty(-1)^n \frac{x^{2n+1}}{2n+1}
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 [[[fallback] (1 0)]] ⧗ Real → Real
	normal 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
	fallback Q.lt? Q.|0 (+1.0) (atan* [1] 1) normal

# 2-argument arctan
atan2 [[[[[go]] (2 0) (1 0)]]] ⧗ Real → Real → Real
	go Q.add a (Q.mul b (Q.add c d))
		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