Check Nim's implementation of the Complex64 types by computing Julia sets

julia_set_benchmark.nim

# -*- encoding: utf-8 -*-

# Check different implementations of the algorithm to determine if a
# point z ∈ ℂ belongs to Julia set J_c

import complex
import math
import strformat
import times

type
  JuliaProc = proc(z: Complex64, c: Complex64, maxiter: int): int

# First implementation, the simplest
func julia1(z: Complex64, c: Complex64, maxiter: int = 256): int =
  var iteridx = 0
  var cur_z = z
  while (abs2(cur_z) < 4) and (iteridx < maxiter):
    cur_z = cur_z * cur_z + c
    iteridx += 1

  result = if iteridx == maxiter:
             -1
           else:
             iteridx

# Second implementation, use incremental operators
func julia2(z: Complex64, c: Complex64, maxiter: int = 256): int =
  var iteridx = 0
  var cur_z = z
  while (abs2(cur_z) < 4) and (iteridx < maxiter):
    cur_z *= cur_z
    cur_z += c
    iteridx += 1

  result = if iteridx == maxiter:
             -1
           else:
             iteridx

# Third implementation, spell out the calculations for the real and
# imaginary parts
func julia3(z: Complex64, c: Complex64, maxiter: int = 256): int =
  var iteridx = 0
  var cur_z = z
  while (cur_z.re * cur_z.re + cur_z.im * cur_z.im < 4) and (iteridx < maxiter):
    let tmp = cur_z.re * cur_z.re - cur_z.im * cur_z.im
    cur_z.im = 2 * cur_z.re * cur_z.im + c.im
    cur_z.re = tmp + c.re
    
    iteridx += 1

  result = if iteridx == maxiter:
             -1
           else:
             iteridx

# Basic interpolation: return out_min if val = 0, out_max if val =
# (max-1), and apply linear interpolation if 0 < val < max - 1
func linspace(val: int, max: int, out_min: float, out_max: float): float {.inline.} =
    result = out_min + (float(val) * (out_max - out_min)) / float(max - 1)

# Run either "julia1", "julia2", or "julia3" on a square grid on the
# complex plane
func juliaset(fn: JuliaProc): seq[int] =
  const num_of_rows = 800
  const num_of_cols = 800
  const c = complex64(re = 0.2, im = 0.55)
  
  result = newSeq[int](num_of_rows * num_of_cols)
  
  var idx = 0
  for row in 1..num_of_rows:
    for col in 1..num_of_cols:
      let z = complex64(re = linspace(col, num_of_cols, -1.5, 1.5),
                        im = linspace(row, num_of_rows, -1.5, 1.5))

      result[idx] = fn(z, c, 10000)
      idx += 1

# Call "juliaset", passing an arbitrary "julia?" function and
# measuring the elapsed time
proc benchmark(fn: JuliaProc, title: string) =
  let start = now()
  let image = juliaset(fn)
  let stop = now()
  
  let duration = stop - start

  # To be sure that all the functions "julia?" are the same, compute
  # the sum of the pixels
  let pixel_sum = sum(image)
  
  echo fmt"{title}: {duration.inMilliseconds} ms (sum of pixels: {pixel_sum})"

benchmark(julia1, "julia1")
benchmark(julia2, "julia2")
benchmark(julia3, "julia3")