Shanks's square forms factorization: Difference between revisions

From testwiki
Jump to navigation Jump to search
No edit summary
 
(No difference)

Latest revision as of 12:13, 16 December 2023

Template:More footnotes

Shanks' square forms factorization is a method for integer factorization devised by Daniel Shanks as an improvement on Fermat's factorization method.

The success of Fermat's method depends on finding integers x and y such that x2y2=N, where N is the integer to be factored. An improvement (noticed by Kraitchik) is to look for integers x and y such that x2y2(modN). Finding a suitable pair (x,y) does not guarantee a factorization of N, but it implies that N is a factor of x2y2=(xy)(x+y), and there is a good chance that the prime divisors of N are distributed between these two factors, so that calculation of the greatest common divisor of N and xy will give a non-trivial factor of N.

A practical algorithm for finding pairs (x,y) which satisfy x2y2(modN) was developed by Shanks, who named it Square Forms Factorization or SQUFOF. The algorithm can be expressed in terms of continued fractions or in terms of quadratic forms. Although there are now much more efficient factorization methods available, SQUFOF has the advantage that it is small enough to be implemented on a programmable calculator. Shanks programmed it on an HP-65, made in 1974, which has storage for only nine digit numbers and allows only 100 steps/keystrokes of programming. There are versions of the algorithm that use little memory and versions that store a list of values that run more quickly.

In 1858, the Czech mathematician Václav Šimerka used a method similar to SQUFOF to factor (10171)/9 = 11111111111111111 = 20717235363222357.[1]

Algorithm

Note This version of the algorithm works on some examples but often gets stuck in a loop.

This version does not use a list.

Input: N, the integer to be factored, which must be neither a prime number nor a perfect square, and a small positive integer, k.

Output: a non-trivial factor of N.

The algorithm:

Initialize i=0,P0=kN,Q1=1,Q0=kNP02.

Repeat

i=i+1,bi=P0+Pi1Qi1,Pi=biQi1Pi1,Qi=Qi2+bi(Pi1Pi)

until Qi is a perfect square at some odd value of i.

Start the second phase (reverse cycle).

Initialize b0=P0PiQi, Q1=Qi, and P0=b0Qi+Pi, where P0,Pi, and Qi are from the previous phase. The b0 used in the calculation of P0 is the recently calculated value of b0.

Set i=0 and Q0=kNP02Q1, where P0 is the recently calculated value of P0.

Repeat

i=i+1,bi=P0+Pi1Qi1,Pi=biQi1Pi1,Qi=Qi2+bi(Pi1Pi)

until Pi=Pi1.Template:Cn

Then if f=gcd(N,Pi) is not equal to 1 and not equal to N, then f is a non-trivial factor of N. Otherwise try another value of k.Template:Cn

Shanks' method has time complexity O(N4).[2]

Stephen S. McMath wrote a more detailed discussion of the mathematics of Shanks' method, together with a proof of its correctness.[3]

Example

Let N=11111

Q1=1

Cycle forward
i bi Pi Qi
0 105 86
1 2 67 77
2 2 87 46
3 4 97 37
4 5 88 91
5 2 94 25

Here Q5=25 is a perfect square, so the first phase ends.

For the second phase, set Q1=25=5. Then:

Reverse cycle
i bi Pi Qi
0 2 104 59
1 3 73 98
2 1 25 107
3 1 82 41
4 4 82

Here P3=P4=82, so the second phase ends. Now calculate gcd(11111,82)=41, which is a factor of 11111.

Thus, N=11111=41271.

Example implementation

Below is an example of C function for performing SQUFOF factorization on unsigned integer not larger than 64 bits, without overflow of the transient operations. Template:Citation needed

#include <inttypes.h>
#define nelems(x) (sizeof(x) / sizeof((x)[0]))

const int multiplier[] = {1, 3, 5, 7, 11, 3*5, 3*7, 3*11, 5*7, 5*11, 7*11, 3*5*7, 3*5*11, 3*7*11, 5*7*11, 3*5*7*11};

uint64_t SQUFOF( uint64_t N )
{
    uint64_t D, Po, P, Pprev, Q, Qprev, q, b, r, s;
    uint32_t L, B, i;
    s = (uint64_t)(sqrtl(N)+0.5);
    if (s*s == N) return s;
    for (int k = 0; k < nelems(multiplier) && N <= UINT64_MAX/multiplier[k]; k++) {
        D = multiplier[k]*N;
        Po = Pprev = P = sqrtl(D);
        Qprev = 1;
        Q = D - Po*Po;
        L = 2 * sqrtl( 2*s );
        B = 3 * L;
        for (i = 2 ; i < B ; i++) {
            b = (uint64_t)((Po + P)/Q);
            P = b*Q - P;
            q = Q;
            Q = Qprev + b*(Pprev - P);
            r = (uint64_t)(sqrtl(Q)+0.5);
            if (!(i & 1) && r*r == Q) break;
            Qprev = q;
            Pprev = P;
        };
        if (i >= B) continue;
        b = (uint64_t)((Po - P)/r);
        Pprev = P = b*r + P;
        Qprev = r;
        Q = (D - Pprev*Pprev)/Qprev;
        i = 0;
        do {
            b = (uint64_t)((Po + P)/Q);
            Pprev = P;
            P = b*Q - P;
            q = Q;
            Q = Qprev + b*(Pprev - P);
            Qprev = q;
            i++;
        } while (P != Pprev);
        r = gcd(N, Qprev);
        if (r != 1 && r != N) return r;
    }
    return 0;
}

References

Template:Reflist

Template:Number theoretic algorithms