## 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.

Step 2b (section 5.2) with multiplier m=1 does not reduce to the routine where multipliers are not used.

This says (m=1):

If Q is even: g=Q/2<L or Q<2L

If Q is odd: g=Q<L

In contrast to (no multipliers, main routine):

If Q is even: Q<L

If Q is odd: Q<L/2

So I wonder whether section 5.2 should read:

Change Step 2b to: g<-Q/gcd(Q,2m) g<=L/2 etc. instead of L.