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
|
import "lib/github.com/diku-dk/complex/complex"
import "lib/github.com/diku-dk/cpprandom/random"
module uniform_real = uniform_real_distribution f64 minstd_rand
module rand_engine = minstd_rand
module c64 = mk_complex f64
type complex = c64.complex
let N = 18i32
let PolyN = 19i32
type poly = [PolyN]f64
let evaln_c (p: poly) (nterms: i32) (pt: complex): complex =
foldr (\coef accum -> c64.mk_re coef c64.+ pt c64.* accum)
(c64.mk_re p[nterms-1]) (take (nterms - 1) p)
let eval_c (p: poly) (pt: complex): complex = evaln_c p (length p) pt
let evaln_d (p: poly) (nterms: i32) (pt: f64): f64 =
foldr (\coef accum -> coef + pt * accum)
p[nterms-1] (take (nterms - 1) p)
let eval_d (p: poly) (pt: f64): f64 = evaln_d p (length p) pt
let derivative (p: poly): *poly =
map (\(i, v) -> f64.i32 i * v) (zip (1...PolyN-1) (take (PolyN - 1) p)) ++ [0]
let max_root_norm (p: poly): f64 =
1 + f64.maximum (map (\coef -> f64.abs (coef / p[PolyN-1])) p)
module aberth = {
type approx = [N]complex
type state = {p: poly, deriv: poly, bound: poly, approx: approx, radius: f64}
let gen_coord (r: f64) (rng: rand_engine.rng): (rand_engine.rng, f64) =
uniform_real.rand (-r, r) rng
let gen_coord_c (r: f64) (rng: rand_engine.rng): (rand_engine.rng, complex) =
let (rng, x) = gen_coord r rng
let (rng, y) = gen_coord r rng
in (rng, c64.mk x y)
let regenerate (rng: rand_engine.rng) (s: state): (rand_engine.rng, *approx) =
let rngs = rand_engine.split_rng N rng
let (rngs, approx) = unzip (map (gen_coord_c s.radius) rngs)
let rng = rand_engine.join_rng rngs
in (rng, approx)
}
entry main: i32 = 42
|
