starred/Foundry
1
0
Fork 0
mirror of https://github.com/Misterio77/Foundry.git synced 2026-08-24 10:04:09 -05:00
Code Issues Packages Projects Releases Wiki Activity

feat(lib): add pure-Nix float math library

Browse Source 
Hand-rolled transcendentals (sqrt, cbrt, exp, ln, pow, sin, cos, atan, atan2)
via range reduction + series/Newton, since Nix has no math builtins. Verified
against known values to ~1e-6. Foundation for a pure-Nix material-color-utilities
port (goal: drop the matugen IFD).

Assisted-by: pi (opus-4.8)
This commit is contained in:
Gabriel Fontes
2026-07-04 17:22:06 -03:00
parent 188c1837c1
commit b48623af3a
1 changed files with 170 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
lib/math.nix
+170
View File
@@ -0,0 +1,170 @@
# Pure-Nix floating-point math, enough to port material-color-utilities.
# Nix has floats + (+,-,*,/) + floor/ceil but no transcendentals, so these are
# hand-rolled via range reduction + series/Newton. Aiming for ~1e-9 accuracy.
{lib}: let
inherit (builtins) floor;
in rec {
pi = 3.141592653589793;
e = 2.718281828459045;
ln2 = 0.6931471805599453;
abs = x:
if x < 0.0
then 0.0 - x
else x;
# Float modulo
fmod = x: m: x - m * (1.0 * floor (x / m));
# sqrt via Newton-Raphson
sqrt = x:
if x <= 0.0
then 0.0
else let
step = g: (g + x / g) / 2.0;
iter = g: n:
if n <= 0
then g
else iter (step g) (n - 1);
in
iter (if x > 1.0 then x else 1.0) 40;
# cbrt via Newton, sign-aware
cbrt = x:
if x == 0.0
then 0.0
else let
s =
if x < 0.0
then -1.0
else 1.0;
a = abs x;
step = g: (2.0 * g + a / (g * g)) / 3.0;
iter = g: n:
if n <= 0
then g
else iter (step g) (n - 1);
in
s * iter (if a > 1.0 then a else 1.0) 60;
# exp via argument halving + Taylor, then repeated squaring
exp = x: let
# reduce until |r| < 1/1024 (halve 10 times max region)
k = 12;
scale = 4096.0; # 2^12
r = x / scale;
# Taylor around 0 for small r
term = i: acc: prod:
if i > 12
then acc
else let
p = prod * r / (1.0 * i);
in
term (i + 1) (acc + p) p;
small = term 1 1.0 1.0;
square = v: n:
if n <= 0
then v
else square (v * v) (n - 1);
in
square small k;
# ln via factoring powers of 2 into [2/3, 4/3], then artanh series
ln = x:
if x <= 0.0
then 0.0 # undefined; caller guards
else let
reduce = v: k:
if v > 1.3333333333
then reduce (v / 2.0) (k + 1)
else if v < 0.6666666667
then reduce (v * 2.0) (k - 1)
else {
inherit v k;
};
red = reduce x 0;
y = (red.v - 1.0) / (red.v + 1.0);
y2 = y * y;
term = i: acc: ypow:
if i > 25
then acc
else let
p = ypow / (1.0 * i);
in
term (i + 2) (acc + p) (ypow * y2);
in
2.0 * (term 1 0.0 y) + (1.0 * red.k) * ln2;
# pow for real exponents (base > 0); integer-exponent fast path for exactness
pow = b: p:
if b == 0.0
then 0.0
else if p == (1.0 * floor p)
then let
n = floor (abs p);
go = acc: i:
if i <= 0
then acc
else go (acc * b) (i - 1);
v = go 1.0 n;
in
if p < 0.0
then 1.0 / v
else v
else exp (p * ln b);
# sin/cos via range reduction to [-pi,pi] + Taylor
sin = x: let
r = fmod (x + pi) (2.0 * pi) - pi;
r2 = r * r;
term = i: acc: prod:
if i > 12
then acc
else let
p = prod * (0.0 - r2) / (1.0 * (2 * i) * (2 * i + 1));
in
term (i + 1) (acc + p) p;
in
term 1 r r;
cos = x: sin (pi / 2.0 - x);
# atan: half-angle reduce (atan z = 2*atan(z/(1+sqrt(1+z^2)))) until |z| small,
# so the alternating series converges fast, then scale back by 2^k.
atan = x: let
reduce = z: k:
if (abs z) < 0.2
then {inherit z k;}
else reduce (z / (1.0 + sqrt (1.0 + z * z))) (k + 1);
red = reduce x 0;
z = red.z;
z2 = z * z;
term = i: acc: zpow:
if i > 25
then acc
else let
s =
if (fmod (1.0 * i) 4.0) == 1.0
then 1.0
else -1.0;
in
term (i + 2) (acc + s * zpow / (1.0 * i)) (zpow * z2);
in
(pow 2.0 (1.0 * red.k)) * (term 1 0.0 z);
atan2 = y: x:
if x > 0.0
then atan (y / x)
else if x < 0.0 && y >= 0.0
then atan (y / x) + pi
else if x < 0.0 && y < 0.0
then atan (y / x) - pi
else if x == 0.0 && y > 0.0
then pi / 2.0
else if x == 0.0 && y < 0.0
then (0.0 - pi) / 2.0
else 0.0;
toDeg = r: r * 180.0 / pi;
toRad = d: d * pi / 180.0;
}