## SQUFOF, Continued Fraction Version

### July 8, 2014

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