From bc72c77a8a3d6a82cd234d40308dd8c47aa0e742 Mon Sep 17 00:00:00 2001
From: tomsmeding <tom.smeding@gmail.com>
Date: Fri, 19 Apr 2019 11:15:28 +0200
Subject: Rename lib.fut to more reasonable name

---
 aberth/aberth_kernel.fut | 126 +++++++++++++++++++++++++++++++++++++++++++++++
 aberth/lib.fut           | 126 -----------------------------------------------
 2 files changed, 126 insertions(+), 126 deletions(-)
 create mode 100644 aberth/aberth_kernel.fut
 delete mode 100644 aberth/lib.fut

diff --git a/aberth/aberth_kernel.fut b/aberth/aberth_kernel.fut
new file mode 100644
index 0000000..a52d8b6
--- /dev/null
+++ b/aberth/aberth_kernel.fut
@@ -0,0 +1,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
diff --git a/aberth/lib.fut b/aberth/lib.fut
deleted file mode 100644
index a52d8b6..0000000
--- a/aberth/lib.fut
+++ /dev/null
@@ -1,126 +0,0 @@
-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
-- 
cgit v1.2.3-54-g00ecf