re-organize into Numeric and Symbolic
[ruff:ruff.git] / Fractal / Mandelbrot / Numeric / RayTraceReverse.hs
1 {-# LANGUAGE BangPatterns, Rank2Types #-}
2 {- |
3 Module      :  Fractal.Mandelbrot.Numeric.RayTraceReverse
4 Copyright   :  (c) Claude Heiland-Allen 2011,2012
5 License     :  BSD3
6
7 Maintainer  :  claudiusmaximus@goto10.org
8 Stability   :  unstable
9 Portability :  BangPatterns, Rank2Types
10
11 Exterior parameters near the boundary can be traced outwards to compute
12 external angles.
13
14 Example usage:
15
16 > main :: IO ()
17 > main = do
18 >   let callback bs = RayTraceReverseCallback $ \m -> case m of
19 >         Nothing -> []
20 >         Just (continue, m') -> case m' of
21 >           Nothing -> let r = fromBits bs % (bit (length bs) - 1) in
22 >             "DONE" : ("\t" ++ show r) : ("\t" ++ show (fromRational r)) : rayTraceReverse continue (callback bs)
23 >           Just (c, m'') -> case m'' of
24 >             Nothing    -> (" \t" ++ show c) : rayTraceReverse continue (callback        bs )
25 >             Just False -> ("0\t" ++ show c) : rayTraceReverse continue (callback (False:bs))
26 >             Just True  -> ("1\t" ++ show c) : rayTraceReverse continue (callback (True :bs))
27 >       ray c = rayTraceReverse (rayTraceReverseStart c) (callback [])
28 >   putStr . unlines . ray . fromComplex fromDouble . read . head =<< getArgs
29
30 -}
31
32 module Fractal.Mandelbrot.Numeric.RayTraceReverse
33   ( RayTraceReverse(rayTraceReverse)
34   , RayTraceReverseCallback(RayTraceReverseCallback)
35   , rayTraceReverseStart
36   ) where
37
38 import qualified Data.Complex as X
39 import Numeric.VariablePrecision
40   ( NaturalNumber, adjustPrecision, VFloat
41   , VComplex, realPart, imagPart, magnitude2, sqr, magnitude2'
42   , scaleComplex, fromComplexDouble, fromDouble, toDouble
43   )
44
45 -- | For each point on the ray, a callback gets passed the current
46 --   coordinate and a bool whenever a dwell boundary is crossed
47 --   along with a continuation that can be used to get more points
48 --   if the ray hasn't passed the escape radius.
49 data RayTraceReverseCallback a = RayTraceReverseCallback{ rayTraceReverseCallback :: NaturalNumber p => Maybe (RayTraceReverse a, Maybe (VComplex p, Maybe Bool)) -> a }
50
51 -- | Step along the ray by one point, calling the callback.
52 data RayTraceReverse  a = RayTraceReverse { rayTraceReverse  :: RayTraceReverseCallback a -> a }
53
54 -- | Initialize ray tracing for a point.
55 rayTraceReverseStart :: NaturalNumber p => VComplex p -> RayTraceReverse a
56 rayTraceReverseStart c = RayTraceReverse{ rayTraceReverse = rayTraceReverseStep (RayTraceReverseContext 8 4 er eps2 c n z d) }
57   where
58     er = 1024
59     er2 = er * er
60     eps2 = encodeFloat 1 (2 * (4 - floatDigits (realPart c)))
61     (n, z) = iter er2 c 0 0
62     d = dwell er2 n z
63
64 data RayTraceReverseContext p =
65   RayTraceReverseContext
66   { rtrcSharpness :: !Int -- number of steps to take within each dwell band
67   --  , rtrcPrecision :: !Int -- enough bits must be available to represent the delta with this much effective precision  FIXME implement incremental precision reduction
68   , rtrcAccuracy  :: !Int -- scales epsilon relative to the length of the last step
69   , rtrcEscapeR :: !Double
70   , rtrcEpsilon2 :: !(VFloat p)
71   , rtrcLastC :: !(VComplex p)
72   , rtrcLastN :: !Integer
73   , rtrcLastZ :: !(X.Complex Double)
74   , rtrcLastD :: !Double
75   }
76
77 iter :: NaturalNumber p => Double -> VComplex p -> Integer -> VComplex p -> (Integer, X.Complex Double)
78 iter er2 c = go
79   where
80     er2' = fromDouble er2
81     go !n !z
82       | magnitude2 z' > er2' = (n, toDouble (realPart z') X.:+ toDouble (imagPart z'))
83       | otherwise = go (n + 1) (sqr z + c)
84       where
85         z' = adjustPrecision z
86
87 dwell :: Double -> Integer -> X.Complex Double -> Double
88 dwell er2 n z = fromIntegral n - logBase 2 (log (magnitude2' z) / log er2)
89
90 iterd :: NaturalNumber p => VComplex p -> Integer -> VComplex p -> VComplex p -> (VComplex p, VComplex p)
91 iterd !c !m !z !dz
92   | m == 0 = (z, dz)
93   | otherwise = iterd c (m - 1) (sqr z + c) (scaleComplex 1 (z * dz) + 1)
94
95 rayTraceReverseStep :: NaturalNumber p => RayTraceReverseContext p -> RayTraceReverseCallback a -> a
96 rayTraceReverseStep rtrc go
97   | escaped   = rayTraceReverseCallback go (Just (stop, Nothing `asTypeOf` Just (c, Nothing)))
98   | otherwise = rayTraceReverseCallback go (Just (continue, Just (newC, newB)))
99   where
100     stop     = RayTraceReverse{ rayTraceReverse = \go' -> rayTraceReverseCallback go' (Nothing `asTypeOf` Just (undefined, Just (c, Nothing))) }
101     continue = RayTraceReverse
102       { rayTraceReverse = rayTraceReverseStep rtrc
103         { rtrcEpsilon2 = scaleFloat (negate (2 * rtrcAccuracy rtrc)) (magnitude2 (newC - c))
104         , rtrcLastC = newC
105         , rtrcLastN = newN
106         , rtrcLastZ = newZ
107         , rtrcLastD = newD
108         }
109       }
110     er2 = er * er
111     er = rtrcEscapeR rtrc
112     eps2 = rtrcEpsilon2 rtrc
113     c = rtrcLastC rtrc
114     n = rtrcLastN rtrc
115     z = rtrcLastZ rtrc
116     d = rtrcLastD rtrc
117     d' = d - 1 / fromIntegral (rtrcSharpness rtrc)
118     m = ceiling d'
119     r = er ** (2 ** (fromIntegral m - d'))
120     a = X.phase z / (2 * pi)
121     t = a - fromIntegral (floor a :: Int)
122     k0 = adjustPrecision . fromComplexDouble $ X.mkPolar r (2 * pi *  t     )
123     k1 = adjustPrecision . fromComplexDouble $ X.mkPolar r (    pi *  t     )
124     k2 = adjustPrecision . fromComplexDouble $ X.mkPolar r (    pi * (t + 1))
125     step !k !cc =
126       let (f, df) = iterd cc m 0 0
127           dc = (f - k) / df
128           cc' = cc - dc
129       in  cc : if not (magnitude2 (cc' - cc) > eps2) then [] else step k cc'
130     steps k = step k c
131     c0 = last (steps k0)
132     (c1, c2) = last (steps k1 `zip` steps k2)
133     dc1 = magnitude2 (c1 - c)
134     dc2 = magnitude2 (c2 - c)
135     (n0, z0) = iter er2 c0 0 0
136     (n1, z1) = iter er2 c1 0 0
137     (n2, z2) = iter er2 c2 0 0
138     d0 = dwell er2 n0 z0
139     d1 = dwell er2 n1 z1
140     d2 = dwell er2 n2 z2
141     (newC, newN, newZ, newD, newB)
142       | m == n    = (c0, n0, z0, d0, Nothing)
143       | dc1 < dc2 = (c1, n1, z1, d1, if n1 == n then Nothing else Just False)
144       | dc2 < dc1 = (c2, n2, z2, d2, if n2 == n then Nothing else Just True )
145       | otherwise = error $ show (dc1, dc2, m, n, n1, n2)
146     escaped = m <= 0