open System.Numerics

open System.Diagnostics

​

[<EntryPoint>]

let main argv =

    let maxIteration = 35000

​

    let time f x =

        let s = new Stopwatch()

        s.Start()

        let r = f x

        s.Stop()

        (s.Elapsed.TotalSeconds, r)

​

    /// Complex math

    let mandelbrotc c =

        let rec loop (z:Complex) iterations =

            if Complex.Abs(Complex.Pow(z, 2.0)) >= 4.0 then iterations

            elif iterations = maxIteration then 0

            else loop (z * z + c) (iterations + 1)

        loop c 0

​

    /// Optimized the extremely slow Pow(z,2.0)

    let mandelbrotco c =

        let rec loop (z:Complex) iterations =

            if Complex.Abs(z) >= 2.0 then iterations

            elif iterations = maxIteration then 0

            else loop ((z * z) + c) (iterations + 1)

        loop c 0

​

    /// Optimized the slow Abs(z)

    let mandelbrotco2 c =

        let rec loop (z:Complex) iterations =

            if z.Real * z.Real + z.Imaginary * z.Imaginary >= 4.0 then iterations

            elif iterations = maxIteration then 0

            else loop ((z * z) + c) (iterations + 1)

        loop c 0

​

    /// Floating-point version

    let mandelbrotf (c:Complex) =

        let cr = c.Real

        let ci = c.Imaginary

        let rec loop zr zi iterations =

            if zr * zr + zi * zi  >= 4.0 then iterations

            elif iterations = maxIteration then 0

            else loop (zr * zr - zi * zi + cr) (2.0 * zr * zi + ci) (iterations + 1)

        loop cr ci 0

​

    /// Floating-point version with common subexpression elimination

    let mandelbrotfs (c:Complex) =

        let cr = c.Real

        let ci = c.Imaginary

        let rec loop zr zi iterations =

            let zr2 = zr * zr

            let zi2 = zi * zi

            let zri = zr * zi

            if zr2 + zi2  >= 4.0 then iterations

            elif iterations = maxIteration then 0

            else loop (zr2 - zi2 + cr) (zri + zri + ci ) (iterations + 1)

        loop cr ci 0

​

​

    let coord = [|

        for y in {-1.2..0.15..1.2} do

            for x in {-2.0..0.07..0.9} do

                yield Complex(x, y)  |]

​

​

    let mandelp f = Array.Parallel.map f coord

    let mandel f = Array.map f coord

​

    let pixelfunc =

        [(mandelbrotc, "Complex")

         (mandelbrotco, "Complex, optimized Pow(z,2.0)")

         (mandelbrotco2, "Complex, optimized Abs(z)")

         (mandelbrotf, "Floating-point")

         (mandelbrotfs, "Floating-point, common subexpression elimination")]

​

    let result f t = pixelfunc |> List.map (fun (f2,txt) -> ((time f f2), txt + t))

​

    let isOk xs = let ((_,r),_) = List.head xs

                  xs |> List.forall (fun ((_,r2),_) -> r = r2)

​

    let r = result mandel " (serial)" @ result mandelp " (parallel)"

    match isOk r with

    | false -> failwith "Results are not identical"

    | _     -> for ((t,_),txt) in r do printfn "%6.3fs  %s" t txt

​

    0 // return an integer exit code

​