[Delphi]快速计算质数序列

       本程序(uPrimes)的算法是从Julia语言的Primes软件包算法转换过来，能实现1秒内找出0到1亿范围内的所有质数，这是编译为64bit的测试结果， 如果编译为32bit的话稍慢些。

       常用的找质数的方法如uPrimesSimple单元所示，方法很简单，但0到1亿范围内找出所有的质数要超过7秒，性能上与前者相比，有7倍差距。

unit uPrimes;

interface

uses SysUtils, Math;

type
  TIntArray  = array of Int64;

function Primes(Lo, Hi:Int64): TIntArray; overload;
function Primes(N: Int64): TIntArray; overload;
function IsPrime(N: Int64): Boolean;

implementation

const
  wheel        : array[1..8 ] of Int64 = (4,  2,  4,  2,  4,  6,  2,  6);
  wheel_primes : array[1..8 ] of Int64 = (7, 11, 13, 17, 19, 23, 29, 31);
  wheel_indices: array[1..31] of Int64 = (0,  0,  0,  0,  0,  0,  0,  1,
                                          1,  1,  1,  2,  2,  3,  3,  3,
                                          3,  4,  4,  5,  5,  5,  5,  6,
                                          6,  6,  6,  6,  6,  7,  7);

  bases: array[1..256] of UInt16 =(
    15591,  2018,   166,  7429,  8064, 16045, 10503,  4399,  1949,  1295,  2776,  3620,
      560,  3128,  5212,  2657,  2300,  2021,  4652,  1471,  9336,  4018,  2398, 20462,
    10277,  8028,  2213,  6219,   620,  3763,  4852,  5012,  3185,  1333,  6227,  5298,
     1074,  2391,  5113,  7061,   803,  1269,  3875,   422,   751,   580,  4729, 10239,
      746,  2951,   556,  2206,  3778,   481,  1522,  3476,   481,  2487,  3266,  5633,
      488,  3373,  6441,  3344,    17, 15105,  1490,  4154,  2036,  1882,  1813,   467,
     3307, 14042,  6371,   658,  1005,   903,   737,  1887,  7447,  1888,  2848,  1784,
     7559,  3400,   951, 13969,  4304,   177,    41, 19875,  3110, 13221,  8726,   571,
     7043,  6943,  1199,   352,  6435,   165,  1169,  3315,   978,   233,  3003,  2562,
     2994, 10587, 10030,  2377,  1902,  5354,  4447,  1555,   263, 27027,  2283,   305,
      669,  1912,   601,  6186,   429,  1930, 14873,  1784,  1661,   524,  3577,   236,
     2360,  6146,  2850, 55637,  1753,  4178,  8466,   222,  2579,  2743,  2031,  2226,
     2276,   374,  2132,   813, 23788,  1610,  4422,  5159,  1725,  3597,  3366, 14336,
      579,   165,  1375, 10018, 12616,  9816,  1371,   536,  1867, 10864,   857,  2206,
     5788,   434,  8085, 17618,   727,  3639,  1595,  4944,  2129,  2029,  8195,  8344,
     6232,  9183,  8126,  1870,  3296,  7455,  8947, 25017,   541, 19115,   368,   566,
     5674,   411,   522,  1027,  8215,  2050,  6544, 10049,   614,   774,  2333,  3007,
    35201,  4706,  1152,  1785,  1028,  1540,  3743,   493,  4474,  2521, 26845,  8354,
      864, 18915,  5465,  2447,    42,  4511,  1660,   166,  1249,  6259,  2553,   304,
      272,  7286,    73,  6554,   899,  2816,  5197, 13330,  7054,  2818,  3199,   811,
      922,   350,  7514,  4452,  3449,  2663,  4708,   418,  1621,  1171,  3471,    88,
    11345,   412,  1559,   194);


type
  TBoolArray = array of Boolean; //为与Julia Primes保持一致,下标范围实际只使用1:n

procedure push(var aArray: TIntArray; N: Int64);
var
  Len: Int64;
begin
  Len := Length(aArray);
  SetLength(aArray, Len + 1);
  aArray[Len] := N;
end;

function max(a, b: Int64): Int64;
begin
  if a > b then
    Result := a
  else
    Result := b;
end;

function isqrt(x: Int64): Int64;
var
  s: Int64;
begin
  if x = 0 then
  begin
    Result := x;
    Exit;
  end;
  s := Int64(trunc(sqrt(x)));
  s := (s + (x div s)) shr 1;
  if s*s > x then
    result := s-1
  else
    Result := s;
end;

function wheel_index(n: Int64): Int64; inline;
var
  d, r: Int64;
begin
  d := (n-1) div 30;
  r := n-1 - d*30;
  Result :=  8*d + wheel_indices[r + 2];
end;

function wheel_prime(n: Int64): Int64; inline;
var
  d, r: Int64;
begin
  d := (n - 1) shr 3;
  r := (n - 1) and 7;
  Result := 30*d + wheel_primes[r + 1];
end;

function _primesmask(limit:Int64): TBoolArray; overload;
var
  m, n, p, q, i, j: Int64;
begin
  if (limit < 7) then
      raise Exception.Create('The condition limit ≥ 7 must be met.');
  n := wheel_index(limit);
  m := wheel_prime(n);
  SetLength(Result, n+1);
  for i := 1 to n do Result[i] := True;
  for i := 1 to wheel_index(isqrt(limit)) do
  begin
    if Result[i] = True then
    begin
      p := wheel_prime(i);
      q := p*p;
      j := ((i - 1) and 7) + 1;
      while q <= m do
      begin
        Result[wheel_index(q)] := False;
        q := q + wheel[j] * p;
        j := (j and 7) + 1;
      end;
    end;
  end;
end;

function _primesmask(lo: Int64; hi: Int64): TBoolArray; overload;
var
  r, q: Int64;
  i, j, wlo, whi, m, p: Int64;
  small_sieve: TBoolArray;
begin
  if not ((7 <= lo) and (lo <= hi)) then
    Raise Exception.Create('The condition 7 ≤ lo ≤ hi must be met.');
  if lo = 7 then
  begin
    Result := _primesmask(hi);
    Exit;
  end;
  wlo := wheel_index(lo - 1);
  whi := wheel_index(hi);
  m := wheel_prime(whi);
  SetLength(Result, whi - wlo + 1);
  for i := 1 to whi - wlo do Result[i] := True;
  if hi < 49 then Exit;
  small_sieve := _primesmask(isqrt(hi));
  for i := 1 to length(small_sieve) - 1 do // don't use eachindex here
  begin
    if small_sieve[i] = True then
    begin
      p := wheel_prime(i);
      j := wheel_index(2 * ((lo - p - 1) div (2*p)) + 1);
      r := p * wheel_prime(j + 1);
      if r > m then Continue; // use widemul to avoid r <= m caused by overflow
      j := (j and 7) + 1;
      q := r;
      // q < 0 indicates overflow when incrementing q inside loop
      while (0 <= q) and (q <= m) do
      begin
        Result[wheel_index(q) - wlo] := False;
        q := q + wheel[j] * p ;
        j := j and 7 + 1;
      end;
    end;
  end;
end;

function primesmask(lo: Int64; hi: Int64): TBoolArray; overload;
var
  i, lsi, lwi: Int64;
  wheel_sieve: TBoolArray;
begin
  if not ((0 < lo) and (lo <= hi)) then
     raise Exception.Create('The condition 0 < lo ≤ hi must be met.');
  SetLength(Result, hi - lo + 1 + 1);
  for i := 1 to hi - lo + 1 do Result[i] := False;

  if (lo <= 2) and (2 <= hi) then Result[3 - lo] := True;
  if (lo <= 3) and (3 <= hi) then Result[4 - lo] := True;
  if (lo <= 5) and (5 <= hi) then Result[6 - lo] := True;
  if (hi < 7) then Exit;
  wheel_sieve := _primesmask(max(7, lo), hi);
  lsi := lo - 1;
  lwi := wheel_index(lsi);
  for i := 1 to length(wheel_sieve)-1 do   // don't use eachindex here
  begin
    Result[wheel_prime(i + lwi) - lsi] := wheel_sieve[i];
  end;
end;

function primesmask(limit: Int64): TBoolArray; overload;
begin
  Result := primesmask(1, limit);
end;

function Primes(Lo, Hi: Int64): TIntArray;
var
  i, lwi: Int64;
  Sieve: TBoolArray;
begin
  SetLength(Result, 0);
  if not (Lo <= Hi) then
      Raise Exception.Create('The condition lo ≤ hi must be met.');
  if (Lo <= 2) and (Hi >= 2) then push(Result, 2);
  if (Lo <= 3) and (Hi >= 3) then push(Result, 3);
  if (Lo <= 5) and (Hi >= 5) then push(Result, 5);
  if Hi < 7 then Exit;
  Lo := max(2, Lo);
  Sieve := _primesmask(max(7, Lo), Hi);
  lwi := wheel_index(Lo - 1);
  for i := 1 to Length(Sieve)-1 do   // don't use eachindex here
  begin
    if Sieve[i] = True then
      push(Result, wheel_prime(i + lwi));
  end;
end;

function Primes(N: Int64): TIntArray;
begin
  Result := Primes(1, N);
end;

function trailing_zeros(n: Int64): Int64;
begin
  if n = 0 then
  begin
    Result := 64;
    Exit;
  end;
  Result := 0;
  while (n and 1) = 0 do
  begin
    Inc(Result);
    n := (n shr 1);
  end;
end;

function _witnesses(n:UInt64): TIntArray;
var
  i: UInt64;
begin
  i := ((n shr 16) xor n) * $45d9f3b;
  i := ((i shr 16) xor i) * $45d9f3b;
  i := ((i shr 16) xor i) and 255 + 1;
  SetLength(Result, 1);
  Result[0] := Int64(bases[i]);
end;

function witnesses(n: Int64): TIntArray;
begin
  if n < 4294967296 then
      Result := _witnesses(UInt64(n))
  else if n < 2152302898747 then
      Result := [2, 3, 5, 7, 11]
  else if n < 3474749660383 then
      Result := [2, 3, 5, 7, 11, 13]
  else
      Result := [2, 325, 9375, 28178, 450775, 9780504, 1795265022];
end;

function QuickPower(Base, N: Int64): Int64;
begin
  if N = 0 then
    Result := 1
  else if N = 1 then
    Result := Base
  else
  begin
    Result := 1;
    while N > 0 do
    begin
      if (N and 1) > 0 then
      begin
        Result := Result * Base;
      end;
      Base := Base * Base;
      N := (N shr 1);
    end;
  end;
end;

function _prevpow2(x: Int64): Int64;
var
  n: Integer;
begin
  if x < 0 then
    raise Exception.Create('`x` must be ≥ 0');
  if x <= 1 then
    Result := x
  else
  begin
    n := 0;
    while x > 0 do
    begin
      Result := x;
      x := x shr 1;
      if x = 0 then
      begin
         Result := Result shl n;
         Break;
      end;
      n := n + 1;
    end;
  end;
end;

function prevpow(a:Int64; x:Int64): Int64;
var
  n, p: Int64;
begin
  if x < 1 then
    raise Exception.Create('`x` must be ≥ 1');
  // See comment in nextpos() for a == special case.
  if a = 2 then
  begin
    Result :=_prevpow2(x);
    Exit
  end;
  if a <= 1 then raise Exception.Create('`a` must be greater than 1.');
  n := Floor(ln(x)/ln(a));
  p := QuickPower(a, n+1);
  if p <= x then
    Result := p
  else
    Result := QuickPower(a,n);
end;

function powermod(x: Int64; p: Int64; m: Int64): Int64;
var
  b, t, r: Int64;
begin
  //if p < 0 then Exit(powermod(invmod(x, m), -p, m));
  if p < 0 then raise Exception.Create('From function "powermod": p < 0 not allowed');
  if p = 0 then
  begin
    Result := (1 mod m);
    Exit;
  end;
  if (m = 1) or (m = -1) then
  begin
    Result := 0;
    Exit;
  end;
  b := (x mod m);   //this also checks for divide by zero
  t := prevpow(2, p);
  r := 1;
  while true do
  begin
    if p >= t then
    begin
      r := (r*b) mod m;
      p := p - t;
    end;
    t  := (t shr 1);
    if (t <= 0) then break;
    r := (r * r) mod m;
  end;
  Result := r;
end;

function IsPrime(N: Int64): Boolean;
var
  m, s, d, a, x, t: Int64;
begin
  // Small precomputed primes + Miller-Rabin for primality testing:
  //  https://en.wikipedia.org/wiki/Miller–Rabin_primality_test
  //  https://github.com/JuliaLang/julia/issues/11594
  for m in [2, 3, 5, 7, 11, 13, 17, 19, 23] do
  begin
    if (N mod m) = 0 then
    begin
      result := (N = m);
      Exit;
    end;
  end;
  if N < 841 then
  begin
    result := (N > 1);
    Exit;
  end;
  s := trailing_zeros(N - 1);
  d := (N - 1) shr s;
  for a in witnesses(N) do
  begin
    x := powermod(a, d, N);
    if x = 1 then continue;
    t := s;
    while x <> (N - 1) do
    begin
      t := t - 1;
      if t <= 0 then
      begin
        Result := False;
        Exit;
      end;
      x := (x * x) mod N;
      if x = 1 then
      begin
        Result := false;
        Exit;
      end;
    end;
  end;
  Result := True;
end;

end.

unit uPrimesSimple;

interface

function PrimesSimple(n: Int64): TArray<Int64>;

implementation

function PrimesSimple(n: Int64): TArray<Int64>;
var
  i, j, len: Int64;
  prime: TArray<Boolean>;
begin
  SetLength(prime, n+1);
  prime[1] := False;
  for i := 2 to n do
    prime[i] := True;

  for i := 2 to n do
  begin
    j := i*i;
    while j <= n do
    begin
      prime[J] := False;
      j := j + i;
    end;
  end;

  SetLength(Result, 0);
  len := 0;
  for i := 1 to n do
  begin
    if prime[i] = True then
    begin
      len := Len + 1;
      SetLength(result, Len);
      Result[Len-1] := i;
    end;
  end;
end;

end.