3 files changed, 59 insertions, 0 deletions
diff --git a/aberth/.gitignore b/aberth/.gitignore
index 60d8b01..66c7e62 100644
--- a/aberth/.gitignore
+++ b/aberth/.gitignore
@@ -1 +1,2 @@
 aberth
+lib/
diff --git a/aberth/futhark.pkg b/aberth/futhark.pkg
new file mode 100644
index 0000000..963e457
--- /dev/null
+++ b/aberth/futhark.pkg
@@ -0,0 +1,4 @@
+require {
+  github.com/diku-dk/complex 0.1.2 #790bf261a6a36a9454b55f3f8109ac5cb1f4bf6e
+  github.com/diku-dk/cpprandom 1.1.4 #7dd47f3f2f5a8fa24d6847b34d9f8a9091b5d94a
+}
diff --git a/aberth/lib.fut b/aberth/lib.fut
new file mode 100644
index 0000000..9e79197
--- /dev/null
+++ b/aberth/lib.fut
@@ -0,0 +1,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