Skip to content

Instantly share code, notes, and snippets.

@siraben
Last active April 7, 2023 21:11
Show Gist options
  • Save siraben/e1b61057bb0074345192e61635183b10 to your computer and use it in GitHub Desktop.
Save siraben/e1b61057bb0074345192e61635183b10 to your computer and use it in GitHub Desktop.
Implementation of sine function in Nix using fourth-order Runge-Kutta
{ lib ? (import <nixpkgs> {}).lib }:
with lib;
let
# Fourth-order Runge-Kutta method
rk4 = f: x0: y0: z0: h:
let
k1_y = z0;
k1_z = f x0 y0 z0;
k2_y = z0 + h / 2 * k1_z;
k2_z = f (x0 + h / 2) (y0 + h / 2 * k1_y) (z0 + h / 2 * k1_z);
k3_y = z0 + h / 2 * k2_z;
k3_z = f (x0 + h / 2) (y0 + h / 2 * k2_y) (z0 + h / 2 * k2_z);
k4_y = z0 + h * k3_z;
k4_z = f (x0 + h) (y0 + h * k3_y) (z0 + h * k3_z);
in y0 + h / 6.0 * (k1_y + 2 * k2_y + 2 * k3_y + k4_y);
solve = f: x0: y0: z0: h: x:
let
iter = xi: yi: zi:
if xi >= x then yi
else iter (xi + h) (rk4 f xi yi zi h) (zi + h * f xi yi zi);
in iter x0 y0 z0;
in
{
rk4 = rk4;
solve = solve;
pi = 3.141592653589;
}
{ lib ? (import <nixpkgs> {}).lib, rk4 ? (import ./rk4.nix { inherit lib; }) }:
let
# Define the system of ODEs for the sine function.
f = x: y: z: -y;
# Initial conditions.
x0 = 0.0;
y0 = 0.0; # y(0) = 0
z0 = 1.0; # y'(0) = 1
h = 0.002; # Step size
fmod = x: y: (x - y * builtins.floor (x / y));
# Compute the sine value at a specific point.
sineValueAt = x:
let
pi = rk4.pi;
pi_2 = pi / 2;
adjustedX = let x' = fmod x (2 * pi); in
if x' >= 0 && x' <= pi_2 then x'
else if x' > pi_2 && x' <= pi then pi - x'
else if x' < 0 && x' >= -pi_2 then -x'
else if x' < -pi_2 && x' >= -pi then (pi + x')
else x'; # Should not reach this point
isNegMultiplier = if adjustedX <= 0 then - 1 else 1;
in isNegMultiplier * rk4.solve f x0 y0 z0 h adjustedX;
in
{
sineValueAt = sineValueAt;
}
@siraben
Copy link
Author

siraben commented Apr 7, 2023

Performance

$ nix run nixpkgs#hyperfine -- "nix-instantiate --eval -E 'let v = (import ./root.nix {}); f = n: x: if n == 0 then v.sineValueAt x else v.sineValueAt x + f (n - 1) x; in 3'"             
Benchmark 1: nix-instantiate --eval -E 'let v = (import ./root.nix {}); in 3'
  Time (mean ± σ):      30.4 ms ±   0.9 ms    [User: 20.0 ms, System: 8.4 ms]
  Range (min … max):    28.7 ms …  32.9 ms    89 runs
 

$ nix run nixpkgs#hyperfine -- "nix-instantiate --eval -E 'let v = (import ./root.nix {}); f = n: x: if n == 0 then v.sineValueAt x else v.sineValueAt x + f (n - 1) x; in f 10000 3.1415'"
Benchmark 1: nix-instantiate --eval -E 'let v = (import ./root.nix {}); f = n: x: if n == 0 then v.sineValueAt x else v.sineValueAt x + f (n - 1) x; in f 10000 3.1415'
  Time (mean ± σ):      53.4 ms ±   0.9 ms    [User: 39.8 ms, System: 10.9 ms]
  Range (min … max):    51.5 ms …  56.3 ms    55 runs

So it takes (53.4-30.4 ms)/10000 = 2.3 μs per function call to sineValueAt, not bad!

Error

Screen Shot 2023-04-07 at 16 11 11

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment