Haskell:Julia 集合函数

Haskell: Julia Set function

我正在尝试编写一个函数,returns Julia Set 的长度在它通过某个阈值 2 之前,最大迭代次数为 255。

我有函数f z = z * z - (a :+ b)我在一个复数上迭代这个函数直到a + bimagnitude越过2,然后计算它所进行的迭代次数。

我的第一次尝试是

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

但速度慢得令人痛苦。

我的第二次尝试是

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

考虑到我将不得不进行数百万次迭代,这仍然很慢。

有人可以帮助提高速度吗?

此外,Prelude 列表函数中是否内置了 list fusion

对于这样的循环,手工制作的递归通常是最好的方法。

下面是三种实现的比较:

  • iter1 - 原始代码
  • iter2 - 手工递归
  • iter3 - @Alec 的解决方案

时间是(-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)

我认为iter2更好的主要原因是因为它避免了 分配复数值。在任何情况下,它的分配都会少很多。

代码如下:

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
    ]
  ]