Square form factorization: Difference between revisions
Added solution for Pascal
(Added Easylang) |
(Added solution for Pascal) |
||
Line 1,965:
1537228672809128917 26675843 57626245319
4611686018427387877 343242169 13435662733</pre>
=={{header|Pascal}}==
{{works with|Free Pascal}}
A console program written in Free Pascal. The code is based on:
Jason E. Gower and Samuel S. Wagstaff, jr., "Square form factorization",
Mathematics of Computation Volume 77, Number 261, January 2008, Pages 551–588
S 0025-5718(07)02010-8
Article electronically published on May 14, 2007
I'm not sure whether this is the same as reference [1] in the task description; the words "a detailed analysis of SquFoF" appear in the abstract.
The Pascal program includes some small changes to the Gower-Wagstaff algorithm, as noted in the comments. The output shows the successful multiplier (if any) and whether the factor was found in the main or preliminary part of the algorithm.
Further to the advice (on the Discussion page) not to use the Wikipedia version of the algorithm, I tested the Gower-Wagstaff version for all odd composite numbers less than 10^9, and it found a factor in every case. The Wikipedia version failed in 790 cases.
<syntaxhighlight lang="pascal">
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.
</syntaxhighlight>
{{out}}
<pre>
Number Mult M/P Factorization
2501 3 Main 61 * 41
12851 1 Main 71 * 181
13289 1 Main 97 * 137
75301 3 Main 293 * 257
120787 1 Main 43 * 2809
967009 7 Main 1609 * 601
997417 1 Main 257 * 3881
7091569 1 Prelim 2663 * 2663
13290059 1 Main 3119 * 4261
42854447 3 Main 4423 * 9689
223553581 1 Main 11213 * 19937
2027651281 1 Main 44021 * 46061
11111111111 3 Main 21649 * 513239
100895598169 11 Main 898423 * 112303
1002742628021 No factor found
60012462237239 1 Main 6862753 * 8744663
287129523414791 5 Main 6059887 * 47381993
9007199254740931 1 Main 10624181 * 847801751
11111111111111111 1 Main 2071723 * 5363222357
314159265358979323 1 Main 317213509 * 990371647
384307168202281507 1 Main 415718707 * 924440401
419244183493398773 1 Main 48009977 * 8732438749
658812288346769681 1 Main 62222119 * 10588072199
922337203685477563 1 Main 110075821 * 8379108103
1000000000000000127 1 Main 111756107 * 8948056861
1152921505680588799 3 Main 139001459 * 8294312261
1537228672809128917 3 Main 26675843 * 57626245319
4611686018427387877 1 Main 343242169 * 13435662733
</pre>
=={{header|Perl}}==
| |||