F# Mandelbrot

Mandelbrot.fsproj

<Project Sdk="Microsoft.NET.Sdk">

  <PropertyGroup>
    <OutputType>Exe</OutputType>
    <TargetFramework>net6.0</TargetFramework>
  </PropertyGroup>

  <ItemGroup>
    <Compile Include="Program.fs" />
  </ItemGroup>

  <ItemGroup>
    <PackageReference Include="SixLabors.ImageSharp" Version="2.1.1" />
  </ItemGroup>

</Project>

Program.fs

module Measure =
  open System
  open System.Diagnostics
  let private sw = Stopwatch.StartNew ()

  let now () = sw.ElapsedMilliseconds

  let time a = 
    let inline cc i = GC.CollectionCount i

    GC.Collect (2, GCCollectionMode.Forced)
    GC.WaitForFullGCComplete () |> ignore

    let bcc0, bcc1, bcc2 = cc 0, cc 1, cc 2
    let before = now ()

    let v = a ()

    let after = now ()
    let acc0, acc1, acc2 = cc 0, cc 1, cc 2

    v, (after-before), (acc0 - bcc0, acc1 - bcc1, acc2 - bcc2)

module Mandelbrot =
  open System
  open System.IO

  open SixLabors.ImageSharp
  open SixLabors.ImageSharp.PixelFormats

  let width   = 2000
  let height  = width

  let maxIter     = 1000        
  let cutOff      = 100.0       
  let centerX     = -0.863493827
  let centerY     = -0.23191307 
  let zoom        = 0.0004      
  let colorPeriod = 0.05        

  let drawPixel (img : Image<Bgr24>) x y iter l2 =
    if iter < maxIter then
      let d   = float iter - Math.Log2 (Math.Log2 (l2))
      let a   = 3.0 + d * colorPeriod

      let half = 255.0*0.5

      let b   = half + half * cos (a + 0.9)
      let g   = half + half * cos (a + 0.5)
      let r   = half + half * cos (a + 0.1)

      let b   = byte b;
      let g   = byte g;
      let r   = byte r;

      img[x,y] <- Bgr24 (r, g, b)
    else
      img[x,y] <- Bgr24 (0uy, 0uy, 0uy)

  let rec mandelbrotLoop zx zy cx cy rem =
    if rem > 0 then
      let re2   = zx*zx
      let im2   = zy*zy
      let reim  = zx*zy
      let l2    = re2 + im2

      if l2 > cutOff then 
        struct (maxIter - rem, l2)
      else
        let zx = re2 - im2 + cx
        let zy = 2.0*reim + cy
        mandelbrotLoop zx zy cx cy (rem - 1)

    else
      struct (maxIter, zx*zx + zy*zy)

  // Computes a mandelbrot set 1 pixel at a time
  let mandelbrot (img : Image<Bgr24>) =
    let w = float width
    let h = float height
    let r = w/h
    for y = 0 to height - 1 do
      let cy = centerY + zoom*(float y - 0.5*h) / h
      for x = 0 to width - 1 do
        let cx = centerX + r * zoom * (float x - 0.5 * w) / w
        let struct (iter, l2) = mandelbrotLoop cx cy cx cy maxIter
        drawPixel img x y iter l2

  let rec doubleMandelbrotLoop zx_0 zy_0 cx_0 cy_0 zx_1 zy_1 cx_1 cy_1 iter_0 iter_1 rem =
    if rem > 0 then
      let re2_0   = zx_0*zx_0
      let re2_1   = zx_1*zx_1

      let im2_0   = zy_0*zy_0
      let im2_1   = zy_1*zy_1

      let reim_0  = zx_0*zy_0
      let reim_1  = zx_1*zy_1

      let iter_0  =
        if iter_0 >= 0 then
          iter_0
        else
          if re2_0+im2_0 > cutOff then
            maxIter - rem
          else
            -1

      let iter_1  =
        if iter_1 >= 0 then
          iter_1
        else
          if re2_1+im2_1 > cutOff then
            maxIter - rem
          else
            -1

      if iter_0 >= 0 && iter_1 >= 0 then 
        struct (iter_0, re2_0+im2_0, iter_1, re2_1+im2_1)
      else
        let zx_0, zy_0 =
          if iter_0 >= 0 then
            zx_0, zy_0
          else
            re2_0 - im2_0 + cx_0, 2.0*reim_0 + cy_0

        let zx_1, zy_1 =
          if iter_1 >= 0 then
            zx_1, zy_1
          else
            re2_1 - im2_1 + cx_1, 2.0*reim_1 + cy_1

        doubleMandelbrotLoop zx_0 zy_0 cx_0 cy_0 zx_1 zy_1 cx_1 cy_1 iter_0 iter_1 (rem - 1)
    else
      let iter_0  =
        if iter_0 >= 0 then
          iter_0
        else
          maxIter

      let iter_1  =
        if iter_1 >= 0 then
          iter_1
        else
          maxIter


      struct (iter_0, zx_0*zx_0 + zy_0*zy_0, iter_1, zx_1*zx_1 + zy_1*zy_1)

  // Computes a mandelbrot set 2 pixel at a time
  //  Even though it is more complex than 1 pixel at a time it performs
  //  50% better on my hardware, why?
  let doubleMandelbrot (img : Image<Bgr24>) =
    let w = float width
    let h = float height
    let r = w/h
    for y = 0 to (height >>> 1) - 1 do
      let y_0 = y <<< 1
      let y_1 = y_0 + 1
      let cy_0 = centerY + zoom*(float y_0 - 0.5*h) / h
      let cy_1 = centerY + zoom*(float y_1 - 0.5*h) / h
      for x = 0 to width - 1 do
        let cx_0 = centerX + r * zoom * (float x - 0.5 * w) / w
        let cx_1 = cx_0
        let struct (iter_0, l2_0, iter_1, l2_1) = doubleMandelbrotLoop cx_0 cy_0 cx_0 cy_0 cx_1 cy_1 cx_1 cy_1 -1 -1 maxIter
        drawPixel img x y_0 iter_0 l2_0
        drawPixel img x y_1 iter_1 l2_1

  let saveImg alg name =
    use img = new Image<Bgr24> (width, height)

    printfn "Starting mandelbrot rendering for: %s" name
    let _, ms, cc = Measure.time (fun () -> alg img)
    let ms = Math.Round (float ms/1000., 2)
    printfn "  ... Rendering mandelbrot %dx%d took %f ms with %A CC" width height ms cc

    let home = Path.GetDirectoryName Environment.ProcessPath
    let path = Path.Combine (home, sprintf "%s.png" name)

    printfn "  Saving mandelbrot image to: %s" path
    img.SaveAsPng path

  let run () =
    saveImg mandelbrot          "simple"
    saveImg doubleMandelbrot    "double"

[<EntryPoint>]
let main argv =
  Mandelbrot.run ()
  0