283 lines
6.2 KiB
ObjectPascal
283 lines
6.2 KiB
ObjectPascal
program CircularPrimes;
|
|
//nearly the way it is done:
|
|
//http://www.worldofnumbers.com/circular.htm
|
|
|
|
{$IFDEF FPC}
|
|
{$Release}
|
|
{$MODE DELPHI}{$OPTIMIZATION ON,ALL}{$Coperators ON}
|
|
// {$O+,R+}
|
|
uses
|
|
Sysutils,gmp;
|
|
{$ENDIF}
|
|
{$IFDEF Delphi}
|
|
uses
|
|
System.Sysutils,?gmp?;
|
|
{$ENDIF}
|
|
|
|
{$IFDEF WINDOWS}
|
|
{$APPTYPE CONSOLE}
|
|
{$ENDIF}
|
|
const
|
|
Pot10 :array[0..19] of Uint64 =
|
|
(1,10,100,1000,10000,100000,1000000,
|
|
10000000,100000000,1000000000,10000000000,100000000000,1000000000000,10000000000000,
|
|
100000000000000,1000000000000000,10000000000000000,100000000000000000,1000000000000000000,10000000000000000000);
|
|
MAXCNTOFDIGITS = 19;
|
|
|
|
type
|
|
tDigits = 0..9;
|
|
tUsedDigits = set of tDigits;
|
|
tmyNum = record
|
|
Dgts: array[0..MAXCNTOFDIGITS] of tDigits;
|
|
num : Uint64;
|
|
UsedDgt : tUsedDigits;
|
|
end;
|
|
|
|
var
|
|
CheckNum : array[0..MAXCNTOFDIGITS] of Uint64;
|
|
Found : array[0..23] of Uint64;
|
|
mpz : mpz_t;
|
|
cntPrmTest,
|
|
cntRot : Uint64;
|
|
SolCount : Int32;
|
|
|
|
procedure ExtendRepUnit(var mpz:mpz_t;idx1,idx2:Int32);
|
|
begin
|
|
inc(idx1);
|
|
dec(idx2,9);
|
|
while idx1 < idx2 do
|
|
begin
|
|
mpz_mul_ui(mpz,mpz,1000000000);// 2020 Windows only uses Uint32
|
|
mpz_add_ui(mpz,mpz,0111111111);
|
|
inc(idx1,9);
|
|
end;
|
|
inc(idx2,9);
|
|
For idx1 := idx1 to idx2 do
|
|
begin
|
|
mpz_mul_ui(mpz,mpz,10);
|
|
mpz_add_ui(mpz,mpz,1);
|
|
end;
|
|
end;
|
|
|
|
procedure CheckOne(MaxIdx:integer);
|
|
begin
|
|
MaxIdx -= 1;
|
|
Found[SolCount] := CheckNum[MaxIdx];
|
|
repeat
|
|
mpz_set_ui(mpz,CheckNum[MaxIdx]);
|
|
If mpz_probab_prime_p(mpz,3)=0then
|
|
EXIT;
|
|
inc(cntPrmTest);
|
|
dec(MaxIdx);
|
|
until MaxIdx < 0;
|
|
inc(SolCount);
|
|
end;
|
|
|
|
procedure InitNum(var myNum:tmyNum;maxDgt:Int32);
|
|
var
|
|
i: Int32;
|
|
begin
|
|
fillchar(myNum,SizeOf(myNum),#0);
|
|
with myNum do
|
|
Begin
|
|
num := 0;
|
|
For i := 0 to MaxDgt-1 do
|
|
begin
|
|
Dgts[i] := 1;
|
|
num += Pot10[i];
|
|
end;
|
|
end;
|
|
end;
|
|
|
|
procedure CorrectNum(var myNum:tmyNum;maxIdx,idx:Int32);
|
|
var
|
|
n : Uint64;
|
|
dgt :int32;
|
|
begin
|
|
with myNum do
|
|
begin
|
|
dec(MaxIdx);
|
|
dgt := Dgts[maxIdx];
|
|
n := 0;
|
|
while maxIdx>idx do
|
|
Begin
|
|
n += Dgts[maxIdx]*Pot10[maxIdx];
|
|
dec(maxIdx);
|
|
end;
|
|
For Idx := Idx downto 0 do
|
|
Begin
|
|
Dgts[idx]:= dgt;
|
|
n += Dgt*Pot10[Idx];
|
|
end;
|
|
num := n;
|
|
end;
|
|
end;
|
|
|
|
function IncNum(var myNum:tmyNum;maxDgt:Int32):boolean;
|
|
//Next number with only digits of 1,3,7,9
|
|
var
|
|
n :Uint64;
|
|
i,dgt: Int64;
|
|
Udgt : tusedDigits;
|
|
begin
|
|
with myNum do
|
|
Begin
|
|
n := num;
|
|
i := 0;
|
|
repeat
|
|
dgt := Dgts[i];
|
|
dgt +=2;
|
|
n += 2*Pot10[i];
|
|
if dgt = 5 then
|
|
begin
|
|
dgt := 7;
|
|
n += 2*Pot10[i];
|
|
end;
|
|
if dgt < 10 then
|
|
begin
|
|
Dgts[i] := dgt;
|
|
Break;
|
|
end;
|
|
//correct values
|
|
Dgts[i] := 1;
|
|
dec(n,Pot10[i]*10);
|
|
inc(i);
|
|
until i >= maxDgt;
|
|
if i >= maxDgt then
|
|
EXIT(false);
|
|
|
|
dgt := Dgts[maxDgt-1];
|
|
if dgt > 1 then
|
|
Begin
|
|
i := maxDgt-2;
|
|
while i >= 0 do
|
|
begin
|
|
if Dgts[i]< dgt then
|
|
begin
|
|
//311->333 // 91111 -> 99999
|
|
CorrectNum(myNum,maxDgt,i);
|
|
EXIT(true);
|
|
end;
|
|
dec(i);
|
|
end;
|
|
end;
|
|
num := n;
|
|
end;
|
|
result := true;
|
|
end;
|
|
|
|
function rotateNum(n:Uint64;maxDgt:Int32):boolean;
|
|
//rotate number and check divisibility by 3 and 7
|
|
var
|
|
Pot,q ,dgt,n0: Uint64;
|
|
begin
|
|
n0 := n;
|
|
if n0 mod 3 = 0 then
|
|
exit(false);
|
|
if n0 mod 7 = 0 then
|
|
exit(false);
|
|
|
|
dec(maxDgt);
|
|
CheckNum[maxDgt] := n;
|
|
Pot := Pot10[maxDgt];
|
|
while maxDgt>0 do
|
|
Begin
|
|
inc(cntRot);
|
|
q := n div 10;
|
|
//last dgt = n mod 10
|
|
dgt := n - 10*q;
|
|
n := dgt*Pot+q;
|
|
//tested elsewhere before
|
|
if n < n0 then
|
|
exit(false);
|
|
|
|
if n mod 7 = 0 then
|
|
exit(false);
|
|
dec(maxDgt);
|
|
CheckNum[maxDgt] := n;
|
|
end;
|
|
result := true;
|
|
end;
|
|
|
|
var
|
|
myNum : TmyNum;
|
|
s :AnsiString;
|
|
T0: Int64;
|
|
idx,idx2 : NativeInt;
|
|
begin
|
|
T0 := GetTickCount64;
|
|
mpz_init(mpz);
|
|
|
|
SolCount := 0;
|
|
//one digit primes including 5
|
|
For idx := 2 to 10 do
|
|
if idx in[2,3,5,7] then
|
|
begin
|
|
Found[SolCount]:= Idx;
|
|
inc(SolCount);
|
|
end;
|
|
|
|
writeln(' search for circular Primes');
|
|
writeln(' digits found rotated numbers to test time in ms ');
|
|
cntPrmTest := 0;
|
|
cntRot := 0;
|
|
For idx := 2 to 16 do
|
|
begin
|
|
InitNum(myNum,idx);
|
|
repeat
|
|
if rotateNum(myNum.num,idx) then
|
|
CheckOne(idx);
|
|
until Not(IncNUm(myNum,idx));
|
|
writeln(idx:7,Solcount:7,cntRot:17,cntPrmTest:12,GetTickCount64-T0:8);
|
|
end;
|
|
|
|
writeln;
|
|
writeln('Found these ',solcount,' circular primes: ');
|
|
For idx := 0 to solCount-2 do
|
|
write(Found[idx],',');
|
|
writeln(Found[solCount-1]);
|
|
|
|
mpz_set_ui(mpz,1);
|
|
For idx := 2 to 1031 do
|
|
begin
|
|
mpz_mul_ui(mpz,mpz,10);
|
|
mpz_add_ui(mpz,mpz,1);
|
|
if mpz_probab_prime_p(mpz,1)=1 then
|
|
writeln('Found prime RepUnit(',idx,')');
|
|
end;
|
|
(*
|
|
idx := 1031;idx2 := 5003;
|
|
ExtendRepUnit(mpz,idx,idx2);
|
|
write(' Found prime RepUnit(',idx2,')');
|
|
writeln(boolean(Ord(mpz_probab_prime_p(mpz,1))));
|
|
idx := idx2;idx2 := 9887;
|
|
ExtendRepUnit(mpz,idx,idx2);
|
|
write(' Found prime RepUnit(',idx2,')');
|
|
writeln(boolean(Ord(mpz_probab_prime_p(mpz,1))));
|
|
idx := idx2;idx2 := 15073;
|
|
ExtendRepUnit(mpz,idx,idx2);
|
|
write(' Found prime RepUnit(',idx2,')');
|
|
writeln(boolean(Ord(mpz_probab_prime_p(mpz,1))));
|
|
idx := idx2;idx2 := 15073;
|
|
ExtendRepUnit(mpz,idx,idx2);
|
|
write(' Found prime RepUnit(',idx2,')');
|
|
writeln(boolean(Ord(mpz_probab_prime_p(mpz,1))));
|
|
idx := idx2;idx2 := 25031;
|
|
ExtendRepUnit(mpz,idx,idx2);
|
|
write(' Found prime RepUnit(',idx2,')');
|
|
writeln(boolean(Ord(mpz_probab_prime_p(mpz,1))));
|
|
idx := idx2;idx2 := 35317;
|
|
ExtendRepUnit(mpz,idx,idx2);
|
|
write(' Found prime RepUnit(',idx2,')');
|
|
writeln(boolean(Ord(mpz_probab_prime_p(mpz,1))));
|
|
idx := idx2;idx2 := 49081;
|
|
ExtendRepUnit(mpz,idx,idx2);
|
|
write(' Found prime RepUnit(',idx2,')');
|
|
writeln(boolean(Ord(mpz_probab_prime_p(mpz,1))));
|
|
// real 8m9,733s
|
|
*)
|
|
mpz_clear(mpz);
|
|
{$IFDEF WINDOWS}
|
|
readln;
|
|
{$ENDIF}
|
|
end.
|