I run into precision problems very easily

Thanks. Indeed with 24 bit float you can't do much. Another example:

But! You can use the derivative of the inflections (before Julia iterations) to calculate the required precision needed

Something like

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

Using this precision (different for each pixel) gives a clean image:

Here's a full working Haskell program using the 'rounded' package for MPFR bindings:

{-

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

]

You can download the program here:

https://mathr.co.uk/misc/2017-03-21_inflection-mapping-dynamic-precision.hsIt's not really practical for live exploring though, 24bit float on GPU is pretty much instant, while this variable precision software floating point on CPU takes several minutes at the same resolution - probably it would be beneficial to just use double precision, or double-double if required for some pixels...