diff options
Diffstat (limited to 'aberth/aberth_kernel.fut')
| -rw-r--r-- | aberth/aberth_kernel.fut | 45 |
1 files changed, 40 insertions, 5 deletions
diff --git a/aberth/aberth_kernel.fut b/aberth/aberth_kernel.fut index a52d8b6..815c254 100644 --- a/aberth/aberth_kernel.fut +++ b/aberth/aberth_kernel.fut @@ -29,7 +29,7 @@ let evaln_d (p: poly) (nterms: i32) (pt: f64): f64 = 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] + map (\(i, v) -> f64.i32 i * v) (zip (1..<PolyN) (drop 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 = @@ -37,8 +37,8 @@ let max_root_norm (p: poly): f64 = 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 + -- 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) = @@ -82,7 +82,7 @@ module aberth = { 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 conditions = map2 (\p s -> c64.mag p <= 1e-9 * s) pvals svals let all_converged = all id conditions in (all_converged, approx) @@ -123,4 +123,39 @@ let derbyshire_at_index (index: i32): *poly = in tabulate PolyN (\i -> f64.i32 (i32.get_bit i bitfield * 2 - 1)) -entry main: i32 = 42 +let point_index + (width: i32) (height: i32) + (top_left: complex) (bottom_right: complex) + (pt: complex) + : i32 = + let x = (c64.re pt - c64.re top_left) / (c64.re bottom_right - c64.re top_left) + let y = (c64.im pt - c64.im top_left) / (c64.im bottom_right - c64.im top_left) + let xi = t64 (x * r64 (width - 1)) + let yi = t64 (y * r64 (height - 1)) + in width * yi + xi + +entry main_job + (start_index: i32) (num_polys: i32) + (width: i32) (height: i32) + (left: f64) (top: f64) (right: f64) (bottom: f64) + (seed: i32) + : []i32 = + -- Unnecessary to give each polynomial a different seed + let rng = rand_engine.rng_from_seed [seed] + let top_left = c64.mk left top + let bottom_right = c64.mk right bottom + let indices = flatten + (map (\idx -> + let p = derbyshire_at_index idx + let (_, pts) = aberth.aberth p rng + let indices = map (point_index width height top_left bottom_right) pts + in indices) + (start_index ..< start_index + num_polys)) + in reduce_by_index (replicate (width * height) 0) (+) 0 indices (replicate (length indices) 1) + +entry main_all + (width: i32) (height: i32) + (left: f64) (top: f64) (right: f64) (bottom: f64) + (seed: i32) + : []i32 = + main_job 0 (1 << N) width height left top right bottom seed |