summaryrefslogtreecommitdiff
path: root/lib/math/poly-impl.myr
blob: e7d645b515c6f74b07906f2b68d6a3f2440fe7f8 (plain)
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
use std

/* See [Mul16], section 5.1 */
pkg math =
        pkglocal const horner_poly32 : (f : flt32, a : flt32[:] -> flt32)
        pkglocal const horner_poly64 : (f : flt64, a : flt64[:] -> flt64)

        pkglocal const horner_polyu32 : (f : flt32, a : uint32[:] -> flt32)
        pkglocal const horner_polyu64 : (f : flt64, a : uint64[:] -> flt64)
;;

extern const fma32 : (x : flt32, y : flt32, z : flt32 -> flt32)
extern const fma64 : (x : flt64, y : flt64, z : flt64 -> flt64)

const horner_poly32 = {f : flt32, a : flt32[:]
        var r : flt32 = 0.0
        for var j = a.len - 1; j >= 0; j--
                r = fma32(r, f, a[j])
        ;;
        -> r
}

const horner_poly64 = {f : flt64, a : flt64[:]
        var r : flt64 = 0.0
        for var j = a.len - 1; j >= 0; j--
                r = fma64(r, f, a[j])
        ;;
        -> r
}

const horner_polyu32 = {f : flt32, a : uint32[:]
        var r : flt32 = 0.0
        for var j = a.len - 1; j >= 0; j--
                r = fma32(r, f, std.flt32frombits(a[j]))
        ;;
        -> r
}

const horner_polyu64 = {f : flt64, a : uint64[:]
        var r : flt64 = 0.0
        for var j = a.len - 1; j >= 0; j--
                r = fma64(r, f, std.flt64frombits(a[j]))
        ;;
        -> r
}

/* TODO: add Knuth's method */