From bc72c77a8a3d6a82cd234d40308dd8c47aa0e742 Mon Sep 17 00:00:00 2001 From: tomsmeding 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.. 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.. + 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.. 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.. 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.. - 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.. 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-70-g09d2