Programming Praxis

Here is our version of SQUFOF, written in Python because it so clearly shows the structure of the algorithm described by Gower and Wagstaff:

def squfof(n, m=1):
    
    def gcd(a, b): # euclid's algorithm
        if b == 0: return a
        return gcd(b, a % b)
    
    def isqrt(n): # newton's method
        x = n; y = (x + 1) // 2
        while y < x:
            x = y; y = (x + n // x) // 2
        return x
    
    def a(i): return 1 if i % 2 == 0 else -1
    
    if n % 2 == 0: return 2
    g = gcd(n, m)
    if g > 1: return g
    
  # 1. Initialize:
    
    s = isqrt(n)
    if s * s == n: return s
    
    d = 2 * m * n if (m * n) % 4 == 1 else m * n
    s = isqrt(d)
    qHat = 1
    p = s
    q = d - p * p
    ell = 2 * isqrt(2 * isqrt(d))
    b = 2 * ell
    i = 0
    queue = []
    
    if verbose:
        print "n =", n
        print "d =", d
        print "s =", s
        print "ell =", ell
        print "b =", b
        print "-------------------------"
        print 0, -1*q, p, 1, queue
    
  # 2. Cycle forward to find a proper square form:
    
    while True:
        
        if verbose: print i+1, a(i)*qHat, p, -1*a(i)*q, queue
        
        # 2a
        littleQ = (s + p) // q
        pPrime = littleQ * q - p
        
        # 2b
        g = q / gcd(q, 2 * m)
        if g <= ell:
            queue.append(g)
            queue.append(p % g)
        
        # 2c
        t = qHat + littleQ * (p - pPrime)
        qHat = q
        q = t
        p = pPrime
        
        # 2d
        if i % 2 == 0:
            
            r = isqrt(q)
            if r * r == q:
                
                z = 0
                while z < len(queue):
                    if r == queue[z] and (p - queue[z+1]) % r == 0:
                        break
                    z = z + 2
                    
                if z < len(queue):
                    if r > 1:
                        queue = queue[z+2:]
                    elif r == 1:
                        return 0 # square form doesn't exist
                else: break # go to step 3
                
        # 2e
        i = i + 1
        if i > b:
            return 0 # timeout
        
  # 3. Compute an inverse square root of the square form:
    
    if verbose: print i+2, -1*a(i)*qHat, p, a(i)*q, queue
    if verbose: print "-------------------------"
    
    qHat = r
    p = p + r * ((s - p) // r)
    q = (d - p * p) / qHat
    
    if verbose: print 0, -1*a(i)*qHat, p, a(i)*q
    
  # 4. Cycle in the reverse direction to find a factor of N:
    
    for i in range(1, b):
        
        # 4a
        littleQ = (s + p) // q
        pPrime = littleQ * q - p
        
        # 4b
        if p == pPrime: break
        
        # 4c
        t = qHat + littleQ * (p - pPrime)
        qHat = q
        q = t
        p = pPrime
    
        if verbose: print i, -1*a(i)*qHat, p, a(i)*q
        
  # 5. Print the factor of N:
    
    if verbose: print i, -1*a(i)*q, p, a(i)*qHat
    if verbose: print "-------------------------"
    
    q = q / gcd(q, 2 * m)
    return q

Gower and Wagstaff show the calculations to find a factor of 22117019; here is the output of our program with verbose = True:

>>> squfof(22117019)
n = 22117019
d = 22117019
s = 4702
ell = 192
b = 384
-------------------------
0 -8215 4702 1 []
1 1 4702 -8215 []
2 -8215 3513 1190 []
3 1190 3627 -7531 []
4 -7531 3904 913 []
5 913 4313 -3850 []
6 -3850 3387 2765 []
7 2765 2143 -6338 []
8 -6338 4195 713 []
9 713 4361 -4346 []
10 -4346 4331 773 []
11 773 4172 -6095 []
12 -6095 1923 3022 []
13 3022 4121 -1699 []
14 -1699 4374 1757 []
15 1757 4411 -1514 []
16 -1514 4673 185 []
17 185 4577 -6314 [185, 48]
18 -6314 1737 3025 [185, 48]
-------------------------
0 -55 4652 8653
1 8653 4001 -706
2 -706 4471 3013
3 3013 4568 -415
4 -415 4562 3145
5 3145 1728 -6083
6 -6083 4355 518
7 518 4451 -4451
8 -4451 4451 518
-------------------------
4451

At http://programmingpraxis.codepad.org/izOZliiR you can see the factorization of 42854447, which demonstrates the use of the queue.

Pages: 1 2