这个素数生成器的执行时间可以改进吗？

Question

我在写这篇文章时的最初目标是尽可能减少占用空间。我可以自信地说，这个目标已经实现。不幸的是，这给我带来了相当缓慢的实施。要生成所有低于 200 万的素数，在 3Ghz 英特尔芯片上大约需要 8 秒。

有没有什么办法可以在减少内存占用的情况下缩短这段代码的执行时间？或者，从功能的角度来看，我是不是走错了路？

代码

/// 6.5s for max = 2,000,000
let generatePrimeNumbers max =    
    let rec generate number numberSequence =
        if number * number > max then numberSequence else
        let filteredNumbers = numberSequence |> Seq.filter (fun v -> v = number || v % number <> 0L)
        let newNumberSequence = seq { for i in filteredNumbers -> i }
        let newNumber = newNumberSequence |> Seq.find (fun x -> x > number)
        generate newNumber newNumberSequence                
    generate 2L (seq { for i in 2L..max -> i })

更新

我调整了算法并设法缩短了 2 秒，但内存消耗增加了一倍。

/// 5.2s for max = 2,000,000
let generatePrimeNumbers max =    
    let rec generate number numberSequence =
        if number * number > max then numberSequence else
        let filteredNumbers = numberSequence |> Seq.filter (fun v -> v = number || v % number <> 0L) |> Seq.toArray |> Array.toSeq
        let newNumber = filteredNumbers |> Seq.find (fun v -> v > number)
        generate newNumber filteredNumbers                
    generate 2L (seq { for i in 2L..max -> i })

更新

显然，我使用的是旧编译器。使用最新版本，我的原始算法需要 6.5 秒而不是 8 秒。这是一个很大的进步。

Answer 1

Tomas Petricek's function 相当快，但我们可以让它更快一点。

比较以下：

let is_prime_tomas n =
    let ms = int64(sqrt(float(n)))
    let rec isPrimeUtil(m) =
        if m > ms then true
        elif n % m = 0L then false
        else isPrimeUtil(m + 1L)
    (n > 1L) && isPrimeUtil(2L)

let is_prime_juliet n =
    let maxFactor = int64(sqrt(float n))
    let rec loop testPrime tog =
        if testPrime > maxFactor then true
        elif n % testPrime = 0L then false
        else loop (testPrime + tog) (6L - tog)
    if n = 2L || n = 3L || n = 5L then true
    elif n <= 1L || n % 2L = 0L || n % 3L = 0L || n % 5L = 0L then false
    else loop 7L 4L

is_prime_juliet 的内部循环稍快一些。它是一种众所周知的素数生成策略，它使用“切换”以 2 或 4 的增量跳过非素数。比较：

> seq { 2L .. 2000000L } |> Seq.filter is_prime_tomas |> Seq.fold (fun acc _ -> acc + 1) 0;;
Real: 00:00:03.628, CPU: 00:00:03.588, GC gen0: 0, gen1: 0, gen2: 0
val it : int = 148933

> seq { 2L .. 2000000L } |> Seq.filter is_prime_juliet |> Seq.fold (fun acc _ -> acc + 1) 0;;
Real: 00:00:01.530, CPU: 00:00:01.513, GC gen0: 0, gen1: 0, gen2: 0
val it : int = 148933

我的版本快了大约 2.37 倍，而且它也非常接近于最快的命令式版本的速度。我们可以让它更快，因为我们不需要过滤2L .. 2000000L 的列表，我们可以在应用过滤器之前使用相同的策略来生成更优化的可能素数序列：

> let getPrimes upTo =
    seq {
        yield 2L;
        yield 3L;
        yield 5L;
        yield! (7L, 4L) |> Seq.unfold (fun (p, tog) -> if p <= upTo then Some(p, (p + tog, 6L - tog)) else None)
    }
    |> Seq.filter is_prime_juliet;;
Real: 00:00:00.000, CPU: 00:00:00.000, GC gen0: 0, gen1: 0, gen2: 0

val getPrimes : int64 -> seq<int64>

> getPrimes 2000000L |> Seq.fold (fun acc _ -> acc + 1) 0;;
Real: 00:00:01.364, CPU: 00:00:01.357, GC gen0: 36, gen1: 1, gen2: 0
val it : int = 148933

我们从 1.530 秒降至 01.364 秒，因此我们的速度提高了大约 1.21 倍。太棒了！

Answer 2

只是为了好玩，让我们看看this page。

pi(x) 是素数计数函数，它返回 x 以下素数的个数。您可以使用以下不等式近似 pi(x)：

(x/log x)(1 + 0.992/log x) < pi(x) < (x/log x)(1 + 1.2762/log x) 
// The upper bound holds for all x > 1

p(x) 是第 n 次的素数函数，可以近似为：

n ln n + n ln ln n - n < p(n) < n ln n + n ln ln n
// for n >= 6

考虑到这一点，here is a very fast generator 计算前 n 个素数，其中索引 i 处的每个元素等于 p(i)。因此，如果我们想将数组的质数限制在 2,000,000 以下，我们可以使用质数计数函数的上界不等式：

let rec is_prime (primes : int[]) i testPrime maxFactor =
    if primes.[i] > maxFactor then true
    else
        if testPrime % primes.[i] = 0 then false
        else is_prime primes (i + 1) testPrime maxFactor

let rec prime_n (primes : int[]) primeCount testPrime tog =
    if primeCount < primes.Length then
        let primeCount' =
            if is_prime primes 2 testPrime (float testPrime |> sqrt |> int) then
                primes.[primeCount] <- testPrime
                primeCount + 1
            else
                primeCount
        prime_n primes primeCount' (testPrime + tog) (6 - tog)

let getPrimes upTo =
    let x = float upTo
    let arraySize = int <| (x / log x) * (1.0 + 1.2762 / log x)
    let primes = Array.zeroCreate (max arraySize 3)
    primes.[0] <- 2
    primes.[1] <- 3
    primes.[2] <- 5

    prime_n primes 3 7 4
    primes

酷！那么它有多快呢？在我的 3.2ghz 四核上，我在 fsi 中得到以下信息：

> let primes = getPrimes 2000000;;
Real: 00:00:00.534, CPU: 00:00:00.546, GC gen0: 0, gen1: 0, gen2: 0

val primes : int [] =
  [|2; 3; 5; 7; 11; 13; 17; 19; 23; 29; 31; 37; 41; 43; 47; 53; 59; 61; 67; 71;
    73; 79; 83; 89; 97; 101; 103; 107; 109; 113; 127; 131; 137; 139; 149; 151;
    157; 163; 167; 173; 179; 181; 191; 193; 197; 199; 211; 223; 227; 229; 233;
    239; 241; 251; 257; 263; 269; 271; 277; 281; 283; 293; 307; 311; 313; 317;
    331; 337; 347; 349; 353; 359; 367; 373; 379; 383; 389; 397; 401; 409; 419;
    421; 431; 433; 439; 443; 449; 457; 461; 463; 467; 479; 487; 491; 499; 503;
    509; 521; 523; 541; ...|]

> printfn "total primes: %i. Last prime: %i" (primes.Length - 1) primes.[primes.Length - 1];;
total primes: 149973. Last prime: 2014853

所以这就是不到半秒内大约 200 万个素数 :)

Answer 3

编辑：以下更新版本，使用更少的内存并且速度更快

有时能够改变一些东西是件好事。这是一个以内存换速度的，无可否认的相当必要的版本。由于这个线程结果在 F# 中托管了一个很好的素数生成函数集合，我认为无论如何添加我的会很好。使用 BitArray 可以控制内存占用。

open System.Collections

let getPrimes nmax =
    let sieve = new BitArray(nmax+1, true)
    let result = new ResizeArray<int>(nmax/10)
    let upper = int (sqrt (float nmax))   
    result.Add(2)

    let mutable n = 3
    while n <= nmax do
       if sieve.[n] then
           if n<=upper then 
               let mutable i = n
               while i <= nmax do sieve.[i] <- false; i <- i + n
           result.Add n
       n <- n + 2
    result

更新：

上面的代码可以进一步优化：因为它只使用筛子中的奇数索引，通过将奇数 n 索引为 2m+1 可以将 BitArray 的大小减少到一半。新版本也更快：

let getPrimes2 nmax =
    let sieve = new BitArray((nmax/2)+1, true)
    let result = new ResizeArray<int>(nmax/10)
    let upper = int (sqrt (float nmax))   
    if nmax>1 then result.Add(2) //fixes a subtle bug for nmax < 2

    let mutable m = 1
    while 2*m+1 <= nmax do
       if sieve.[m] then
           let n = 2*m+1
           if n <= upper then 
               let mutable i = m
               while 2*i < nmax do sieve.[i] <- false; i <- i + n
           result.Add n
       m <- m + 1
    result

时序（intel core duo 2.33GHz）：

> getPrimes 2000000 |> Seq.length;;
Real: 00:00:00.037, CPU: 00:00:00.046, GC gen0: 0, gen1: 0, gen2: 0
val it : int = 148933
> getPrimes2 2000000 |> Seq.length;;
Real: 00:00:00.026, CPU: 00:00:00.031, GC gen0: 0, gen1: 0, gen2: 0
val it : int = 148933

Answer 4

尹发布的命令式版本非常快。在我的机器上也是大约 0.5 秒。 但是，如果你想写一个简单的函数式解决方案，你可以这样写：

let isPrime(n) =
  let ms = int64(sqrt(float(n)))
  let rec isPrimeUtil(m) =
    if m > ms then true
    elif n % m = 0L then false
    else isPrimeUtil(m + 1L)
  (n > 1L) && isPrimeUtil(2L)

[ 1L .. 2000000L ] |> List.filter isPrime

这只是测试一个数字是否是 100 万以内所有数字的质数。它不使用任何复杂的算法（有趣的是，最简单的解决方案通常就足够了！）。在我的机器上，您的更新版本运行大约 11 秒，而这大约运行 2 秒。

更有趣的是，这很容易并行化。如果您使用 PLINQ，您可以编写下面的版本，它在双核上的运行速度几乎快 2 倍。这意味着在四核上，它可能与此处所有答案中最快的解决方案一样快，但只需最少的编程工作:-)（当然，使用四核不是生态的，但是......好吧）

[ 1L .. 2000000L ] |> PSeq.ofSeq |> PSeq.filter isPrime |> List.ofSeq

PSeq 函数是我为我的书创建的 PLINQ 包装器（它使使用 F# 中的 PLINQ 更加自然）。它们在source code for Chapter 14 中可用。

Answer 5

我写了一个命令式版本，速度更快。以纯函数方式编写 Sieve of Eratosthenes 以达到相同的速度可能是不可能的，因为每个数字都必须有一个二进制状态。

let generatePrimes max=
    let p = Array.create (max+1) true
    let rec filter i step = 
        if i <= max then 
            p.[i] <- false
            filter (i+step) step
    {2..int (sqrt (float max))} |> Seq.map (fun i->filter (i+i) i) |> Seq.length |> ignore
    {2..max} |> Seq.filter (fun i->p.[i]) |> Seq.toArray