Data.Complex.magnitude is very slow sometimes
The magnitude
function tries very hard to support accurate computations for large exponents:
magnitude (x:+y) = scaleFloat k
(sqrt (sqr (scaleFloat mk x) + sqr (scaleFloat mk y)))
where k = max (exponent x) (exponent y)
mk = - k
sqr z = z * z
However this can lead to more than 10x slower programs. Here is an example program:
import qualified Data.ByteString.Lazy as B
import Data.Complex
import System.Environment
iter = 50
limit = 2
main = do
[arg] <- getArgs
putStr $ "P4\n" ++ arg ++ " " ++ arg ++ "\n"
B.putStr . B.pack . makePBM . map fractal $ points (read arg) (read arg)
points w h =
[ (2 * x / w - 1.5) :+ (2 * y / h - 1)
| y <- [0 .. h - 1]
, x <- [0 .. w - 1]
]
fractal c = length . takeIter $ iterate (\z -> z * z + c) c
where
takeIter = take iter . takeWhile (\z -> magnitude z < limit)
makePBM =
map
( foldl (\n x -> n * 2 + if x == iter then 1 else 0) 0
. rightPad 8 0
)
. chunksOf 8
chunksOf _ [] = []
chunksOf n xs = let (l, r) = splitAt n xs in l : chunksOf n r
rightPad n x xs = foldr (\x go n -> x : go (n - 1)) (`replicate` x) xs n
When compiled with ghc -O2 -fllvm
and running with ./SimpleMandelbrot 1024 >/dev/null +RTS -s
takes 2.136s. In comparison replacing magnitude z
with sqrt (x * x + y * y)
where x :+ y = z
, makes the program run in 0.184s. This difference gets even larger for larger inputs such as 16000
for which I've measured about 700 seconds without vs. about 40 seconds with this simple change.
Note that using LLVM is quite critical, I think LLVM is able to compile that sqrt x < c
to x < c * c
, but even if I apply that manually LLVM still outperforms the native back end.
Is this enough reason to change the default implementation of magnitude
? Should we add a second magnitudeFast
function?
Environment
- GHC version used: 8.10.7