fork(2) download
copy 
  1. {
  2.   Copyright (C) 2020 I. Kakoulidis, <ioulianos.kakoulidis@hotmail.com>
  3.   https://g...content-available-to-author-only...b.com/JulStrat/primesieve-pas
  4.  
  5.   This file is distributed under the BSD 2-Clause License.
  6. }
  7.  
  8. program hprimes;
  9. {$IF Defined(FPC)}{$MODE Delphi}{$ASMMODE Intel}{$ENDIF}
  10.  
  11. function MulModASM(a, b: UInt64; m: UInt64): UInt64;
  12. label
  13.   MUL_MOD;
  14.  
  15. asm
  16.   MOV RCX, m
  17.   MOV RAX, a
  18.   CMP RAX, RCX
  19.   JB MUL_MOD
  20.   XOR RDX, RDX
  21.   DIV RCX
  22.   MOV RAX, RDX
  23.  
  24. MUL_MOD:
  25.   MUL b
  26.   DIV RCX
  27.   MOV @Result, RDX
  28. end;
  29.  
  30. function PowMod(b, e: UInt64; m: UInt64): UInt64;
  31. begin
  32.   if b >= m then b := b mod m;
  33.   Result := 1;
  34.  
  35.   while e <> 0 do
  36.   begin
  37.     if (e and 1) = 1 then Result := MulModASM(Result, b, m);
  38.     e := e shr 1;
  39.     b := MulModASM(b, b, m);
  40.   end;
  41. end;
  42.  
  43. function MillerRabin(n: UInt64): boolean;
  44. var
  45.   a, d, x: UInt64;
  46.   i, r: integer;
  47.   witness: boolean;
  48.  
  49. begin
  50.   case n of
  51.     2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37: Exit(True);
  52.   end;
  53.   if ((n and 1) = 0) or (n = 1) then Exit(False);
  54.   d := n - 1;
  55.   r := 0;
  56.  
  57.   while (d and 1) = 0 do
  58.   begin
  59.     Inc(r);
  60.     d := d shr 1;
  61.   end;
  62.  
  63.   for a in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37] do
  64.   begin
  65.     x := PowMod(a, d, n);
  66.     if (x = 1) or (x = n - 1) then continue;
  67.     witness := false;
  68.  
  69.     for i := 1 to r - 1 do
  70.     begin
  71.       x := MulModASM(x, x, n);
  72.       if x = n - 1 then
  73.       begin
  74.         witness := true;
  75.         break;
  76.       end;
  77.     end;
  78.  
  79.     if witness = false then Exit(false);
  80.   end;
  81.   Result := true;
  82. end;
  83.  
  84. var
  85.   m: UInt64;
  86.   c: integer;
  87.  
  88. begin
  89.   m := $FFFFFFFFFFFFFFFF;
  90.   c := 0;
  91.  
  92.   while c < 10 do
  93.   begin
  94.     if MillerRabin(m) then
  95.     begin
  96.       WriteLn(($FFFFFFFFFFFFFFFF - m) + 1);
  97.       Inc(c);
  98.     end;
  99.     m := m - 1;
  100.   end;
  101.  
  102. end.
Success #stdin #stdout 0s 4360KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
59
83
95
179
189
257
279
323
353
363
https://ideone.com/oO83MV
language:
Pascal (fpc 3.0.4)
created:
3 years ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!