Haskell: Julia Set fonction

Question

Je suis en train d'écrire une fonction qui retourne la longueur d'un ensemble de Julia où avant qu'il ne passe un certain seuil, 2, avec un nombre maximum d'itérations de 255.

j'ai fonction f z = z * z - (a :+ b) J'itère cette fonction sur un nombre complexe jusqu'à ce que le magnitude de a + bi croise 2 puis compte le nombre d'itérations qu'il a fallu.

ma première tentative a été

iter1 = genericLength . take 255 . takeWhile ((<=2) . magnitude) . iterate f

mais il est très lent.

ma deuxième tentative est

iters2 x y = iters2' f ((<=2) . magnitude) (a :+ b) 
iters2' f' p c = len c 0 
    where 
    len c acc = if p c && acc < 255 then len (f' c) (acc + 1) else acc 

Ceci est encore lent, étant donné que je vais devoir faire des millions d'itérations.

Est-ce que quelqu'un peut vous aider à accélérer le processus?

En outre, y a-t-il list fusion intégré dans les fonctions de liste ?

Answer 1

Pour une boucle comme celle-ci, la récursion artisanale est généralement la meilleure approche.

Voici une comparaison de trois implémentations:

• iter1 - code original
• iter2 - récursion fabriqués à la main
• iter3 - @ solution d'Alec

Les timings sont (avec -O2) :

benchmarking julia/1 
time     17.30 μs (17.06 μs .. 17.56 μs) 
        0.995 R² (0.991 R² .. 0.998 R²) 
mean     17.79 μs (17.37 μs .. 18.51 μs) 
std dev    1.839 μs (1.008 μs .. 3.285 μs) 
variance introduced by outliers: 86% (severely inflated) 

benchmarking julia/2 
time     1.456 μs (1.448 μs .. 1.466 μs) 
        0.999 R² (0.999 R² .. 1.000 R²) 
mean     1.457 μs (1.444 μs .. 1.470 μs) 
std dev    42.03 ns (34.14 ns .. 54.37 ns) 
variance introduced by outliers: 38% (moderately inflated) 

benchmarking julia/3 
time     13.53 μs (13.35 μs .. 13.69 μs) 
        0.999 R² (0.998 R² .. 0.999 R²) 
mean     13.42 μs (13.26 μs .. 13.58 μs) 
std dev    550.8 ns (445.3 ns .. 768.8 ns) 
variance introduced by outliers: 50% (moderately inflated) 

Je pense que la principale raison pour laquelle iter2 est meilleur est parce qu'il évite l'allocation de valeurs complexes. En tout cas, il fait beaucoup moins d'allocations.

Code suit:

import Data.List 
import Data.Ratio 
import Data.Complex 
import Criterion 
import Criterion.Main 

iter1 :: Double -> Double -> Complex Double -> Int 
iter1 a b = genericLength . take 255 . takeWhile ((<=2) . magnitude) . iterate f 
    where 
    f z = z * z - (a :+ b) 

iter2 :: Double -> Double -> Complex Double -> Int 
iter2 a b (x :+ y) = loop 0 x y 
    where 
    loop n x y 
     | n > 255 || x*x + y*y > 4 = n 
    loop n x y = loop (n+1) (x*x-y*y-a) (2*x*y-b) 

iter3 :: Double -> Double -> Complex Double -> Int 
iter3 a b z = fst $ until (\(n,z) -> magnitude z > 2 || n == 255) 
         (\(n,z) -> (n+1,f z)) 
         (0,z) 
    where 
    f z = z * z - (a :+ b) 

-- should give 101 iterations 
z0 = 7.396147400000001e-3 :+ 0.6418972592000001 

main = defaultMain [ 
    bgroup "julia" 
    [ 
     bench "3" $ whnf (iter3 (negate 0.40) 0.65) z0 
    , bench "2" $ whnf (iter2 (negate 0.40) 0.65) z0 
    , bench "1" $ whnf (iter1 (negate 0.40) 0.65) z0 
    ] 
    ]