sounio-examples/newton_root.sio at main · Sounio-lang/sounio-examples · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
// Sprint 230: Functional Calculus - Native ELF Demo
//
// Real f64 numerical computing: differentiate, integrate, find roots, optimize.
// Every algorithm takes fn(f64) -> f64 as argument.
// Pure Sounio - no libc, no runtime, no extern C.
//
// Compile:
//   $SOUC run self-hosted/compiler/render_native_compile_driver_lean.sio \
//     -- examples/native_calculus.sio /tmp/calculus.elf
// Run:
//   /tmp/calculus.elf; echo $?  # -> 50

fn fc_abs(x: f64) -> f64 { if x < 0.0 { 0.0 - x } else { x } }

fn my_sin(x: f64) -> f64 with Mut, Panic, Div {
    let pi = 3.14159265358979
    let two_pi = 6.28318530717959
    var r = x
    while r > pi { r = r - two_pi }
    while r < 0.0 - pi { r = r + two_pi }
    let x2 = r * r
    let t1 = r
    let t3 = 0.0 - r * x2 / 6.0
    let t5 = 0.0 - t3 * x2 / 20.0
    let t7 = 0.0 - t5 * x2 / 42.0
    let t9 = 0.0 - t7 * x2 / 72.0
    let t11 = 0.0 - t9 * x2 / 110.0
    t1 + t3 + t5 + t7 + t9 + t11
}

fn parabola(x: f64) -> f64 { x * x - 4.0 * x + 3.0 }

fn deriv(f: fn(f64) -> f64, x: f64) -> f64 with Mut, Panic, Div {
    let h = 0.00000001
    (f(x + h) - f(x - h)) / (2.0 * h)
}

fn integrate(f: fn(f64) -> f64, a: f64, b: f64) -> f64 with Mut, Panic, Div {
    let n: i64 = 100
    let h = (b - a) / 100.0
    var sum = f(a) + f(b)
    var i: i64 = 1
    var x = a + h
    while i < n {
        if i % 2 == 0 { sum = sum + 2.0 * f(x) }
        else { sum = sum + 4.0 * f(x) }
        i = i + 1
        x = x + h
    }
    sum * h / 3.0
}

fn minimize(f: fn(f64) -> f64, a: f64, b: f64) -> f64 with Mut, Panic, Div {
    let resphi = 0.3819660113
    var lo = a
    var hi = b
    var x1 = lo + resphi * (hi - lo)
    var x2 = hi - resphi * (hi - lo)
    var f1 = f(x1)
    var f2 = f(x2)
    var i: i64 = 0
    while i < 100 {
        if fc_abs(hi - lo) < 0.0000000001 { return 0.5 * (lo + hi) }
        if f1 < f2 {
            hi = x2
            x2 = x1
            f2 = f1
            x1 = lo + resphi * (hi - lo)
            f1 = f(x1)
        } else {
            lo = x1
            x1 = x2
            f1 = f2
            x2 = hi - resphi * (hi - lo)
            f2 = f(x2)
        }
        i = i + 1
    }
    0.5 * (lo + hi)
}

fn f2exit(x: f64) -> i64 with Mut, Panic, Div {
    var result: i64 = 0
    var rem = if x < 0.0 { 0.0 - x } else { x }
    while rem >= 1.0 {
        rem = rem - 1.0
        result = result + 1
    }
    if x < 0.0 { 0 - result } else { result }
}

fn main() -> i64 with Mut, Panic, Div {
    let pi = 3.14159265358979

    // d/dx sin(0) = cos(0) = 1.0 -> f2exit(1.0 * 10) = 10
    let d = deriv(my_sin, 0.0)

    // integral of sin from 0 to pi = 2.0 -> f2exit(2.0 * 10) = 20
    let area = integrate(my_sin, 0.0, pi)

    // minimize x^2-4x+3 on [0,5] = 2.0 -> f2exit(2.0 * 10) = 20
    let min_x = minimize(parabola, 0.0, 5.0)

    // Total: 10 + 20 + 20 = 50
    f2exit(d * 10.0) + f2exit(area * 10.0) + f2exit(min_x * 10.0)
}