Data update
This commit is contained in:
Ingy döt Net
parent
35bcdeebf8
commit
74c69a0df6
2427 changed files with 31826 additions and 3468 deletions
|
|
@ -0,0 +1,293 @@
|
|||
program SquFoF_console;
|
||||
|
||||
{$mode objfpc}{$H+}
|
||||
|
||||
uses SquFoF_utils;
|
||||
|
||||
type TResultKind =
|
||||
(rkPrelim, // a factor was found by the preliminary routine
|
||||
rkMain, // a factor was found by the main algorithm
|
||||
rkFail); // no factor was found
|
||||
type TAlgoResult = record
|
||||
Kind : TResultKind;
|
||||
Mult : word;
|
||||
Factor : UInt64;
|
||||
end;
|
||||
|
||||
// Preliminary to G-W algorithm. Returns D and S of the algorithm.
|
||||
// Also returns a non-trivial factor if found, else returns factor = 1.
|
||||
procedure GWPrelim( N : UInt64; // number to be factorized
|
||||
m : word; // multiplier
|
||||
out D, S, factor : UInt64);
|
||||
var
|
||||
sqflag : boolean;
|
||||
begin
|
||||
D := m*N;
|
||||
sqflag := SquFoF_utils.IsSquare( D, S);
|
||||
if m = 1 then
|
||||
if sqflag then factor := S
|
||||
else factor := 1
|
||||
else
|
||||
factor := GCD( N,m);
|
||||
end;
|
||||
|
||||
// Tries to factorize N by applying Gower-Wagstaff alsorithm to m*N.
|
||||
function GW_with_multiplier( N : UInt64;
|
||||
m : word) : TAlgoResult;
|
||||
var
|
||||
D, S, P, P_prev, Q, L, B: Uint64;
|
||||
r : UInt64;
|
||||
i, j, k : integer;
|
||||
f, g : UInt64;
|
||||
sqrtD : double;
|
||||
endCycle : boolean;
|
||||
// Queue is not much used, so make it a simple array.
|
||||
type TQueueItem = record
|
||||
Left, Right : UInt64;
|
||||
end;
|
||||
const QUEUE_CAPACITY = 50; // as suggested by Gower & Wagstaff
|
||||
var
|
||||
queue : array [0..QUEUE_CAPACITY - 1] of TQueueItem;
|
||||
queueCount : integer;
|
||||
begin
|
||||
result.Mult := m;
|
||||
|
||||
// Filter out special cases (differs from Gower & Wagstaff). Note:
|
||||
// (1) multiplier m is assumed to be squarefree;
|
||||
// (2) if we proceed to the main algorithm, mN must not be square
|
||||
// (otherwise Q = 0 and division by Q causes an error).
|
||||
GWPrelim( N, m, {out} D, S, f);
|
||||
if f > 1 then begin
|
||||
result.Kind := rkPrelim;
|
||||
result.Factor := f;
|
||||
exit;
|
||||
end;
|
||||
// Not special, proceed to main algorithm
|
||||
result.Kind := rkMain;
|
||||
result.Mult := m;
|
||||
result.Factor := 1;
|
||||
queueCount := 0; // Clear queue
|
||||
P := S;
|
||||
Q := 1;
|
||||
i := -1; // keep i same as in G & W; algo fails if i > B
|
||||
sqrtD := SquFoF_utils.FSqrt( D);
|
||||
// L := Trunc( 2.0*Sqrt( 2.0*sqrtD));
|
||||
L := Trunc( Sqrt( 2.0*sqrtD)); // as in Section 5.2 of G&W paper
|
||||
B := 2*L;
|
||||
|
||||
// Start forward cycle
|
||||
endCycle := false;
|
||||
while not endCycle do begin
|
||||
// We update Q here, P at the end of the loop
|
||||
Q := (D - P*P) div Q;
|
||||
if (not Odd(i)) and SquFoF_utils.IsSquare( Q, r) then begin
|
||||
// Q is square for even i.
|
||||
// Possibly (?probably?) ends the forward cycle,
|
||||
// but we need to inspect the queue first.
|
||||
endCycle := true;
|
||||
j := queueCount; // working backwards down the queue
|
||||
if r = 1 then begin // the method may fail
|
||||
while (j > 0) and (result.Kind <> rkFail) do begin
|
||||
dec( j);
|
||||
if queue[j].Left = 1 then result.Kind := rkFail;
|
||||
end;
|
||||
if result.Kind = rkFail then exit;
|
||||
end
|
||||
else begin // if r > 1
|
||||
while (j > 0) and (endCycle) do begin
|
||||
dec( j);
|
||||
if (queue[j].Left = r)
|
||||
and ((P - queue[j].Right) mod r = 0) then begin
|
||||
// Deleting up to the *last* item in the list that
|
||||
// satisfies the condition, Should it be the *first* ?
|
||||
// Delete queue items 0..j inclusive
|
||||
inc(j); k := 0;
|
||||
while j < queueCount do begin
|
||||
queue[k] := queue[j];
|
||||
inc(j); inc(k);
|
||||
end;
|
||||
queueCount := k;
|
||||
endCycle := false;
|
||||
end; // if
|
||||
end;
|
||||
end;
|
||||
end; // if i even and Q square
|
||||
if not endCycle then begin
|
||||
g := Q div SquFoF_utils.GCD( Q, 2*m);
|
||||
if g <= L then begin
|
||||
if queueCount < QUEUE_CAPACITY then begin
|
||||
with queue[queueCount] do begin
|
||||
Left := g; Right := P mod g;
|
||||
end;
|
||||
inc( queueCount);
|
||||
end
|
||||
else begin // queue overflow, fail
|
||||
result.Kind := rkFail;
|
||||
exit;
|
||||
end;
|
||||
end;
|
||||
inc(i);
|
||||
if i > B then begin
|
||||
result.Kind := rkFail;
|
||||
exit;
|
||||
end;
|
||||
end;
|
||||
P := S - ((S + P) mod Q);
|
||||
end; // while not endCycle
|
||||
Assert( (D - P*P) mod r = 0); // optional check
|
||||
P := S - ((S + P) mod r);
|
||||
Q := r;
|
||||
// Start backward cycle
|
||||
endCycle := false;
|
||||
while not endCycle do begin
|
||||
P_prev := P;
|
||||
Q := (D - P*P) div Q;
|
||||
P := S - ((S + P) mod q);
|
||||
endCycle := (P = P_prev);
|
||||
end; // while not endCycle
|
||||
// Finished
|
||||
result.Factor := Q div SquFoF_utils.GCD( Q, 2*m);
|
||||
end;
|
||||
|
||||
const NR_RC_VALUES = 28;
|
||||
RC_VALUES : array [0..NR_RC_VALUES - 1] of UInt64 =
|
||||
( 2501, 12851, 13289, 75301, 120787, 967009, 997417, 7091569, 13290059,
|
||||
42854447, 223553581, 2027651281, 11111111111, 100895598169, 1002742628021,
|
||||
60012462237239, 287129523414791, 9007199254740931, 11111111111111111,
|
||||
314159265358979323, 384307168202281507, 419244183493398773,
|
||||
658812288346769681, 922337203685477563, 1000000000000000127,
|
||||
1152921505680588799, 1537228672809128917, 4611686018427387877);
|
||||
|
||||
type TMultAndMaxN = record
|
||||
Mult : word; // small multiplier
|
||||
MaxN : UInt64; // maximum N for that multiplier (N*multiplier < 2^64)
|
||||
end;
|
||||
const NR_MULTIPLIERS = 16;
|
||||
const MULTIPLIERS : array [0..NR_MULTIPLIERS - 1] of TMultAndMaxN =
|
||||
((Mult: 1; MaxN: 18446744073709551615),
|
||||
(Mult: 3; MaxN: 6148914691236517205),
|
||||
(Mult: 5; MaxN: 3689348814741910323),
|
||||
(Mult: 7; MaxN: 2635249153387078802),
|
||||
(Mult: 11; MaxN: 1676976733973595601),
|
||||
(Mult: 15; MaxN: 1229782938247303441),
|
||||
(Mult: 21; MaxN: 878416384462359600),
|
||||
(Mult: 33; MaxN: 558992244657865200),
|
||||
(Mult: 35; MaxN: 527049830677415760),
|
||||
(Mult: 55; MaxN: 335395346794719120),
|
||||
(Mult: 77; MaxN: 239568104853370800),
|
||||
(Mult: 105; MaxN: 175683276892471920),
|
||||
(Mult: 165; MaxN: 111798448931573040),
|
||||
(Mult: 231; MaxN: 79856034951123600),
|
||||
(Mult: 385; MaxN: 47913620970674160),
|
||||
(Mult: 1155; MaxN: 15971206990224720));
|
||||
|
||||
function GowerWagstaff( N : UInt64) : TAlgoResult;
|
||||
var
|
||||
j : integer;
|
||||
begin
|
||||
j := 0;
|
||||
result.Kind := rkFail;
|
||||
while (result.Kind = rkFail)
|
||||
and (j < NR_MULTIPLIERS)
|
||||
and (N <= MULTIPLIERS[j].MaxN) do
|
||||
begin
|
||||
result := GW_with_multiplier( N, MULTIPLIERS[j].Mult);
|
||||
if result.Kind = rkFail then inc(j);
|
||||
end;
|
||||
end;
|
||||
|
||||
// Main program
|
||||
var
|
||||
j : integer;
|
||||
ar : TAlgoResult;
|
||||
kindStr : string;
|
||||
N, cofactor : UInt64;
|
||||
begin
|
||||
WriteLn( ' Number Mult M/P Factorization');
|
||||
for j := 0 to NR_RC_VALUES - 1 do begin
|
||||
N := RC_VALUES[j];
|
||||
ar := GowerWagstaff( N);
|
||||
if ar.Kind = rkFail then
|
||||
WriteLn( N:20, ' No factor found')
|
||||
else begin
|
||||
case ar.Kind of
|
||||
rkPrelim: kindStr := 'Prelim';
|
||||
rkMain : kindStr := 'Main ';
|
||||
end;
|
||||
cofactor := N div ar.Factor;
|
||||
Assert( cofactor * ar.Factor = N); // check that all has gone well
|
||||
WriteLn( N:20, ar.Mult:5, ' ',
|
||||
kindStr:6, ' ', ar.Factor, ' * ', cofactor);
|
||||
end;
|
||||
end;
|
||||
end.
|
||||
|
||||
unit SquFoF_utils;
|
||||
|
||||
{$mode objfpc}{$H+}
|
||||
|
||||
interface
|
||||
|
||||
// Returns floating-point square root of 64-bit unsigned integer.
|
||||
function FSqrt( x : UInt64) : double;
|
||||
|
||||
// Returns whether a 64-bit unsigned integer is a perfect square.
|
||||
// In either case, returns floor(sqrt(x)) in the out parameter.
|
||||
function IsSquare( x : UInt64; out iroot : UInt64) : boolean;
|
||||
|
||||
// Returns g.c.d. of 64-bit and 16-bit unsigned integer.
|
||||
function GCD( u : UInt64; x : word) : word;
|
||||
|
||||
implementation
|
||||
|
||||
function FSqrt( x : UInt64) : double;
|
||||
// Both Free Pascal and Delphi 7 seem unreliable when casting
|
||||
// a 64-bit integer to floating point. We use a workaround.
|
||||
type TSplitUint64 = packed record case boolean of
|
||||
true: (All : UInt64);
|
||||
false: (Lo, Hi : longword); // longword is 32-bit unsigned
|
||||
end;
|
||||
var
|
||||
temp : TSplitUInt64;
|
||||
begin
|
||||
temp.All := x;
|
||||
result := Sqrt( 1.0*temp.Lo + 4294967296.0*temp.Hi);
|
||||
end;
|
||||
|
||||
// Based on Rosetta Code ISqrt, solution for Modula-2.
|
||||
// Trunc of the f.p. square root won't do, bacause of rounding errors..
|
||||
function IsSquare( x : UInt64; out iroot : UInt64) : boolean;
|
||||
var
|
||||
Xdiv4, q, r, s, z : UInt64;
|
||||
begin
|
||||
Xdiv4 := X shr 2;
|
||||
q := 1;
|
||||
while q <= Xdiv4 do q := q shl 2;
|
||||
z := x;
|
||||
r := 0;
|
||||
repeat
|
||||
s := q + r;
|
||||
r := r shr 1;
|
||||
if z >= s then begin
|
||||
z := z - s;
|
||||
r := r + q;
|
||||
end;
|
||||
q := q shr 2;
|
||||
until q = 0;
|
||||
iroot := r;
|
||||
result := (z = 0);
|
||||
end;
|
||||
|
||||
function GCD( u : UInt64; x : word) : word;
|
||||
var
|
||||
y, t : word;
|
||||
begin
|
||||
y := u mod x;
|
||||
while y <> 0 do begin
|
||||
t := x mod y;
|
||||
x := y;
|
||||
y := t;
|
||||
end;
|
||||
result := x;
|
||||
end;
|
||||
end.
|
||||
|
|
@ -0,0 +1,56 @@
|
|||
const multipliers = divisors(3*5*7*11).grep { _ > 1 }
|
||||
|
||||
func sff(N) {
|
||||
|
||||
N.is_prime && return 1 # n is prime
|
||||
N.is_square && return N.isqrt # n is square
|
||||
|
||||
multipliers.each {|k|
|
||||
|
||||
var P0 = isqrt(k*N) # P[0]=floor(sqrt(N)
|
||||
var Q0 = 1 # Q[0]=1
|
||||
var Q = (k*N - P0*P0) # Q[1]=N-P[0]^2 & Q[i]
|
||||
var P1 = P0 # P[i-1] = P[0]
|
||||
var Q1 = Q0 # Q[i-1] = Q[0]
|
||||
var P = 0 # P[i]
|
||||
var Qn = 0 # P[i+1]
|
||||
var b = 0 # b[i]
|
||||
|
||||
while (!Q.is_square) { # until Q[i] is a perfect square
|
||||
b = idiv(isqrt(k*N) + P1, Q) # floor(floor(sqrt(N+P[i-1])/Q[i])
|
||||
P = (b*Q - P1) # P[i]=b*Q[i]-P[i-1]
|
||||
Qn = (Q1 + b*(P1 - P)) # Q[i+1]=Q[i-1]+b(P[i-1]-P[i])
|
||||
(Q1, Q, P1) = (Q, Qn, P)
|
||||
}
|
||||
|
||||
b = idiv(isqrt(k*N) + P, Q) # b=floor((floor(sqrt(N)+P[i])/Q[0])
|
||||
P1 = (b*Q0 - P) # P[i-1]=b*Q[0]-P[i]
|
||||
Q = (k*N - P1*P1)/Q0 # Q[1]=(N-P[0]^2)/Q[0] & Q[i]
|
||||
Q1 = Q0 # Q[i-1] = Q[0]
|
||||
|
||||
loop {
|
||||
b = idiv(isqrt(k*N) + P1, Q) # b=floor(floor(sqrt(N)+P[i-1])/Q[i])
|
||||
P = (b*Q - P1) # P[i]=b*Q[i]-P[i-1]
|
||||
Qn = (Q1 + b*(P1 - P)) # Q[i+1]=Q[i-1]+b(P[i-1]-P[i])
|
||||
break if (P == P1) # until P[i+1]=P[i]
|
||||
(Q1, Q, P1) = (Q, Qn, P)
|
||||
}
|
||||
with (gcd(N,P)) {|g|
|
||||
return g if g.is_ntf(N)
|
||||
}
|
||||
}
|
||||
|
||||
return 0
|
||||
}
|
||||
|
||||
[ 11111, 2501, 12851, 13289, 75301, 120787, 967009, 997417, 4558849, 7091569, 13290059,
|
||||
42854447, 223553581, 2027651281, 11111111111, 100895598169, 1002742628021, 60012462237239,
|
||||
287129523414791, 11111111111111111, 384307168202281507, 1000000000000000127, 9007199254740931,
|
||||
922337203685477563, 314159265358979323, 1152921505680588799, 658812288346769681,
|
||||
419244183493398773, 1537228672809128917
|
||||
].each {|n|
|
||||
var v = sff(n)
|
||||
if (v == 0) { say "The number #{n} is not factored." }
|
||||
elsif (v == 1) { say "The number #{n} is a prime." }
|
||||
else { say "#{n} = #{[n/v, v].sort.join(' * ')}" }
|
||||
}
|
||||
Loading…
Add table
Add a link
Reference in a new issue