diff options
Diffstat (limited to 'aberth/aberth_kernel.fut')
| -rw-r--r-- | aberth/aberth_kernel.fut | 14 |
1 files changed, 8 insertions, 6 deletions
diff --git a/aberth/aberth_kernel.fut b/aberth/aberth_kernel.fut index fd4d943..5d851f0 100644 --- a/aberth/aberth_kernel.fut +++ b/aberth/aberth_kernel.fut @@ -29,7 +29,8 @@ let evaln_d (p: poly) (nterms: i32) (pt: f32): f32 = let eval_d (p: poly) (pt: f32): f32 = evaln_d p (length p) pt let derivative (p: poly): *poly = - map (\(i, v) -> f32.i32 i * v) (zip (1..<PolyN) (tail p)) ++ [0] + let res = map (\(i, v) -> f32.i32 i * v) (tail (zip (0..<PolyN) p)) ++ [0] + in res :> poly -- 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): f32 = @@ -44,14 +45,14 @@ module aberth = { let gen_coord (r: f32) (rng: *rand_engine.rng): *(rand_engine.rng, f32) = uniform_real.rand (-r, r) rng - let gen_coord_c (r: f32) (rng: rand_engine.rng): (rand_engine.rng, complex) = + let gen_coord_c (r: f32) (rng: *rand_engine.rng): (rand_engine.rng, complex) = let (rng, x) = gen_coord r rng let (rng, y) = gen_coord r rng in (rng, cplx.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 (rngs, approx) = unzip (map (\rng -> gen_coord_c ctx.radius (copy rng)) rngs) let rng = rand_engine.join_rng rngs in (rng, approx) @@ -107,7 +108,8 @@ module aberth = { -- Set the constant coefficient to 1; nextDerbyshire will never change it let init_derbyshire: poly = - [1] ++ replicate (PolyN - 1) (-1) + let res = [1] ++ replicate (PolyN - 1) (-1) + in res :> poly let next_derbyshire (p: *poly): *(bool, poly) = let (_, p, looped) = @@ -152,14 +154,14 @@ entry main_job let indices = flatten (map (\idx -> let p = derbyshire_at_index idx - let (_, pts) = aberth.aberth p rng + let (_, pts) = aberth.aberth p (copy rng) in map (point_index width height bottom_left top_right) pts) (start_index ..< start_index + num_polys)) let filtered = filter (\i -> i != -1) indices in reduce_by_index (replicate (width * height) 0) (+) 0 filtered - (replicate (length filtered) 1) + (map (const 1) filtered) entry main_all (width: i32) (height: i32) |