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 = 20i32
let PolyN = N + 1

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) (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 =
    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))


let calc_index (value: f64) (left: f64) (right: f64) (steps: i32): i32 =
    t64 ((value - left) / (right - left) * (r64 steps - 1) + 0.5)

let point_index
        (width: i32) (height: i32)
        (bottom_left: complex) (top_right: complex)
        (pt: complex)
        : i32 =
    -- Range for 'yi' is reversed because image coordinates go down in the y
    -- direction, while complex coordinates go up in the y direction
    let xi = calc_index (c64.re pt) (c64.re bottom_left) (c64.re top_right) width
    let yi = calc_index (c64.im pt) (c64.im top_right) (c64.im bottom_left) height
    in if 0 <= xi && xi < width && 0 <= yi && yi < height
        then width * yi + xi
        else -1

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 bottom_left = c64.mk left bottom
    let top_right = c64.mk right top
    let indices = flatten
            (map (\idx ->
                    let p = derbyshire_at_index idx
                    let (_, pts) = aberth.aberth p 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)

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