codepad

Language:

# Square Form Factorization
# Jason E Gower and Samuel S Wagstaff Jr
# AMS Mathematics of Computation
# Volume 77, Number 261, January 2008, Pages 551-588
# S 0025-5718(07)02010-8
# http://homes.cerias.purdue.edu/~ssw/squfof.pdf
# Algorithm 3.3 SQUFOF (continued fraction method)

# see also
# SQUFOF notes by Daniel Shanks
# collected by Steve McMath
# http://www.usna.edu/Users/math/wdj/_files/
#     documents/mcmath/shanks_squfof.pdf

# note the typo in Gower and Wagstaff section 5.2
# in the description of step 2.b which should read
#     put (g, P mod g) onto the QUEUE
# instead of
#     put (g, B mod g) onto the QUEUE

verbose = True

def squfof(n, m=1): # with optional multiplier
    # useful multipliers are 3, 5, 7, 11, 15, 21,
    #   33, 35, 55, 77, 105, 165, 231, 385, 1155
    
    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. Initialise
    
    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

print squfof(42854447)

Private [?]

Run code