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
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
|
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
-- First element of pair steps fastest
let iota2 (n: i32) (m: i32): [](i32, i32) =
flatten (map (\y -> map (\x -> (x, y)) (iota n)) (iota m))
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) (take (PolyN - 1) p)) ++ [0]
-- Cauchy's bound: https://en.wikipedia.org/wiki/Geometrical_properties_of_polynomial_roots#Lagrange's_and_Cauchy's_bounds
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
-- bound is 's' in the stop condition formulated at p.189-190 of
-- https://link.springer.com/article/10.1007%2FBF02207694
type context = {p: poly, deriv: poly, bound: poly, 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 generate (ctx: context) (rng: *rand_engine.rng): *(rand_engine.rng, approx) =
let rngs = rand_engine.split_rng N rng
let (rngs, approx) = unzip (map (gen_coord_c ctx.radius) rngs)
let rng = rand_engine.join_rng rngs
in (rng, approx)
let compute_bound_poly (p: poly): *poly =
map2 (\coef i -> f64.abs coef * f64.i32 (4 * i + 1)) p (0..<PolyN)
let initialise (p: *poly): *context =
let deriv = derivative p
let bound = compute_bound_poly p
let radius = max_root_norm p
in {p, deriv, bound, radius}
-- Lagrange-style step where the new elements are computed in parallel from
-- the previous values
let step (ctx: context) (approx: *approx): *(bool, approx) =
let pvals = map (eval_c ctx.p) approx
let derivvals = map (evaln_c ctx.deriv (PolyN - 1)) approx
let quos = map2 (c64./) pvals derivvals
let sums = map (\i ->
reduce_comm (c64.+) (c64.mk_re 0.0)
(map (\j ->
if i == j then c64.mk_re 0.0
else c64.mk_re 1.0 c64./
(approx[i] c64.- approx[j]))
(0..<N)))
(0..<N)
let offsets = map2 (\quo sum -> quo c64./ (c64.mk_re 1.0 c64.- quo c64.* sum))
quos sums
let approx = map2 (c64.-) approx offsets
let svals = map (eval_d ctx.bound <-< c64.mag) approx
let conditions = map2 (\p s -> c64.mag p > 1e-9 * s) pvals svals
let all_converged = all id conditions
in (all_converged, approx)
let iterate (ctx: context) (rng: *rand_engine.rng): (rand_engine.rng, *approx) =
let (rng, approx) = generate ctx rng
let (init_conv, approx) = step ctx approx
let (rng, _, _, _, approx) =
loop (rng, conv, tries, step_idx, approx) =
(rng, init_conv, 1, 1: i32, approx)
while !conv
do if step_idx + 1 > tries * 100
then let (rng, approx) = generate ctx rng
let (conv, approx) = step ctx approx
in (rng, conv, tries + 1, 0, approx)
else let (conv, approx) = step ctx approx
in (rng, conv, tries, step_idx + 1, approx)
in (rng, approx)
let aberth (p: *poly) (rng: *rand_engine.rng): *(rand_engine.rng, approx) =
iterate (initialise p) rng
}
-- Set the constant coefficient to 1; nextDerbyshire will never change it
let init_derbyshire: poly =
[1] ++ replicate (PolyN - 1) (-1)
let next_derbyshire (p: *poly): *(bool, poly) =
let (_, p, looped) =
loop (i, p, cont) = (0, p, true)
while cont && i < length p
do if p[i] == -1
then (i, p with [i] = 1, false)
else (i + 1, p with [i] = -1, true)
in (looped, p)
let derbyshire_at_index (index: i32): *poly =
let bitfield = (index << 1) + 1
in tabulate PolyN (\i -> f64.i32 (i32.get_bit i bitfield * 2 - 1))
entry main: i32 = 42
|
