Haskell:Julia 集合函数
Haskell: Julia Set function
我正在尝试编写一个函数,returns Julia Set 的长度在它通过某个阈值 2 之前,最大迭代次数为 255。
我有函数f z = z * z - (a :+ b)
我在一个复数上迭代这个函数直到a + bi
的magnitude
越过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
]
]
我正在尝试编写一个函数,returns Julia Set 的长度在它通过某个阈值 2 之前,最大迭代次数为 255。
我有函数f z = z * z - (a :+ b)
我在一个复数上迭代这个函数直到a + bi
的magnitude
越过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
]
]