Welcome to Fractal Forums

max(16, ceil(16 - log2(abs(derivative))))

{-

Simple inflection mapping with dynamic precision.  First compute the inflections
with low constant precision, then use the computed derivative to decide how much
precision to use for each pixel in the final render.

Usage:

    ghc -O2 Main.hs
    ./Main > output.pgm

Runtime is about 5m30s on one core of a 3GHz amd64 desktop running Debian.

-}

{-# LANGUAGE DataKinds #-}
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE ScopedTypeVariables #-}

module Main (main) where

-- MPFR binding, not on Hackage yet, get it from either of these:
-- https://code.mathr.co.uk/rounded/shortlog/refs/heads/claude5
-- https://github.com/ekmett/rounded/pull/27
import Numeric.Rounded

import Data.Complex (Complex((:+)), magnitude)
import Data.List (foldl')
import Data.Monoid ((<>))
import Data.Proxy (Proxy)
import Data.Word (Word8)

import qualified Data.ByteString as BS
import qualified Data.ByteString.Char8 as BSC

-- arbitrary precision types (p bits of precision)

type R (p :: k) = Rounded 'TowardNearest p
type C (p :: k) = Complex (R p)

-- dual numbers for automatic differentiation

data AD p = AD !(C p) !(C p)
  deriving (Eq, Read, Show)

instance Precision p => Num (AD p) where
  negate (AD a b) = AD (negate a) (negate b)
  AD a b + AD c d = AD (a + c) (b + d)
  AD a b - AD c d = AD (a - c) (b - d)
  AD a b * AD c d = AD (a * c) (a * d + b * c)
  fromInteger = constant . fromInteger

constant :: Precision p => C p -> AD p
constant p = AD p 0

variable :: Precision p => C p -> C p -> AD p
variable d p = AD p d

getVariable :: Precision p => AD p -> C p
getVariable (AD x _) = x

getDerivative :: Precision p => AD p -> C p
getDerivative (AD _ y) = y

-- inflection mapping

inflect :: Precision p => C p -> AD p -> AD p
inflect p c =
  let f = constant p
      d = c - f
  in  d * d + f

inflects :: Precision p => [C p] -> AD p -> AD p
inflects ps z = foldl' (flip inflect) z ps

-- estimate required precision from derivative

precisionRequired :: Precision p => AD p -> Int
precisionRequired (AD _ d) = max 16 . ceiling $ 16 - logBase 2 (magnitude d)

-- generate a grid of coordinates

grid :: Precision p => Int -> Int -> C p -> R p -> [[C p]]
grid w h c r =
    [ [ c + d * (x i :+ y j) | i <- take w [0..] ] | j <- take h [0..] ]
  where
    d = 2 * r :+ 0
    aspect = fromInt w / fromInt h
    x i = ((fromInt i + 0.5) / fromInt w - 0.5) * aspect
    y j = (fromInt h - fromInt j + 0.5) / fromInt h - 0.5

clamp :: Ord a => a -> a -> a -> a
clamp mi ma = min ma . max mi

-- simple greyscale colouring

grey :: Double -> Word8
grey = fromIntegral . clamp 0 (255 :: Int) . floor . (255 *) . tanh . clamp 0 8

-- distance estimator

de :: Precision p => AD p -> Double
de (AD z' dz')
  | isInfinite mdz || isNaN mdz || isInfinite mz || isNaN mz = 1e10
  | otherwise = 2 * mz * log mz / mdz
  where
    z = fmap toDouble z'
    dz = fmap toDouble dz'
    mz = magnitude z
    mdz = magnitude dz

-- Julia iteration loop

calculate :: Precision p => AD p -> AD p -> AD p
calculate c
    = head . (++ [AD 0 0])
    . dropWhile ((< escapeRadius) . magnitude . getVariable)
    . take maxIters
    . iterate (z -> z * z + c)

-- calculate and colour a point

render :: Precision p => [C p] -> AD p -> Proxy p -> Word8
render ps z0 _proxy = grey . de . calculate c . inflects ps $ z0
  where
    c = constant . fmap precRound . last $ ps

-- compute a pixel

pixel :: Precision p => [C p] -> C p -> Word8
pixel ps z0 = reifyPrecision prec render'
  where
    -- compute precision at base precision of input
    -- prec = 24 -- uncomment this and comment the next line to see the effect
    prec = precisionRequired . inflects ps . variable delta0 $ z0
    delta0 = fmap precRound delta
    -- round to computed necessary precision for rendering
    render' :: forall p . Precision p => Proxy p -> Word8
    render' = render ps' (variable delta' z0')
      where
        z0' = fmap precRound z0
        delta' = fmap precRound delta
        ps' = map (fmap precRound) ps

-- pixel size
delta :: C 24
delta = radius * 2 / fromInt height :+ 0

main' :: [C 24] -> BS.ByteString
main' ps
  =  BSC.pack ("P5
" ++ show width ++ " " ++ show height ++ "
255
")
  <> BS.pack (concatMap (map (pixel ps)) (grid width height (head ps) radius))

main :: IO ()
main = BS.putStr (main' example)

-- image parameters

width :: Int
width = 1280

height :: Int
height = 702

-- iteration parameters

escapeRadius :: Precision p => R p
escapeRadius = 65536

maxIters :: Int
maxIters = 1000

-- view size

radius :: R 24
radius = 1.1

-- example inflections

example :: [C 24]
example =
  [   0.29166666  :+ 1.20277774
  ,   0.15833333  :+ 1.21944439
  ,   0.06388889  :+ 1.20277774
  , (-0.05833333) :+ 1.21388888
  , (-0.22500001) :+ 1.17499995
  , (-0.34166667) :+ 1.16388893
  , (-0.48055553) :+ 1.19166672
  , (-0.59722221) :+ 1.13611114
  , (-0.69166666) :+ 1.06388891
  , (-0.80833334) :+ 0.96388888
  , (-0.76944447) :+ 0.13055556
  ]