prime.pas 1.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081
  1. unit Prime;
  2. interface
  3. function isPrime(X : Longint) : Boolean;
  4. function isPrimeQ(X : Longint) : Boolean;
  5. {
  6. History :
  7. Created : Dec. 2004, Weidong Hu
  8. }
  9. implementation
  10. function isPrime(X : Longint) : Boolean;
  11. var
  12. i : longint;
  13. begin
  14. isPrime := false;
  15. if x < 2 then exit;
  16. for i := 2 to trunc(sqrt(x)) do
  17. if x mod i = 0 then
  18. exit;
  19. isPrime := true;
  20. end;
  21. function getExpMod(base, exp, _mod : longint) : longint;
  22. var
  23. tmp : int64;
  24. begin
  25. if exp = 0 then getExpMod := 1 else
  26. begin
  27. tmp := getExpMod(base, exp shr 1, _mod);
  28. tmp := tmp * tmp mod _mod;
  29. if odd(exp) then
  30. tmp := tmp * base mod _mod;
  31. getExpMod := tmp;
  32. end;
  33. end;
  34. function isPrimeQ(X : Longint) : Boolean;
  35. { Rabin-Miller Strong Pseudoprime Test }
  36. const
  37. S : array[1..4] of longint = (2, 3, 5, 7{, 11});
  38. var
  39. d, m : longint;
  40. tmp : int64;
  41. a : longint;
  42. i : integer;
  43. OK : boolean;
  44. begin
  45. isPrimeQ := false;
  46. if x < 2 then exit;
  47. d := X - 1;
  48. while d and 1 = 0 do
  49. d := d shr 1;
  50. for i := Low(S) to High(S) do
  51. begin
  52. OK := false;
  53. a := S[i];
  54. if a >= X then break;
  55. tmp := getExpMod(a, d, X);
  56. if (tmp = 1) or (tmp = x - 1) then OK := true else
  57. begin
  58. m := d shl 1;
  59. while m < X do
  60. begin
  61. tmp := tmp * tmp mod x;
  62. if tmp = x - 1 then
  63. begin
  64. OK := true;
  65. break;
  66. end;
  67. m := m shl 1;
  68. end;
  69. end;
  70. if not OK then exit;
  71. end;
  72. isPrimeQ := true;
  73. end;
  74. end.