## Lenstra's algorithm for divisors in residue classes

For the discussion of math. Duh.

Moderators: gmalivuk, Moderators General, Prelates

LucasBrown
Posts: 299
Joined: Thu Apr 15, 2010 2:57 am UTC
Location: Poway, CA

### Lenstra's algorithm for divisors in residue classes

I'm attempting to implement Lenstra's algorithm for finding divisors in residue classes as described as Algorithm 9.1.29 in Cohen's A Course in Computational Algebraic Number Theory (also available here in PDF form).

My initial translation of the description into Python-with-goto is as follows:

Code: Select all

def isqrt(n):   # Shamelessly stolen from https://codegolf.stackexchange.com/a/9088.
if n < 0: return int(n)
c = n*4//3
d = c.bit_length()
a = d>>1
if d&1:
x = 1 << a
y = (x + (n >> a)) >> 1
else:
x = (3 << a) >> 2
y = (x + (c >> a)) >> 1
if x != y:
x, y = y, (y + n//y) >> 1
while y < x: x, y = y, (y + n//y) >> 1
return x

def xgcd_goto(a, b): # Algorithm 1.3.6
assert a >= 0 <= b
# 1: Initialize
u = 1
d = a
if b == 0: v = 0; return (u, v, d)
v1 = 0
v3 = b
# 2: Finished?
if v3 == 0: v = (d - a*u) // b; return (u, v, d)
# 3: Euclidean step
q, t3 = divmod(d, v3)
t1 = u - q * v1
u = v1
d = v3
v1 = t1
v3 = t3
goto 2

def moddiv_goto(r, s, n): # Algorithm 9.1.29
assert 0 <= r < s < n and gcd(r,s) == 1 and s**3 > n
# 1: Initialization
u, v, _ = xgcd(r, s)
rprime = (u * n) % s
assert 0 <= rprime < s
a0 = s
b0 = 0
c0 = 0
a1 = (u * rprime) % s
b1 = 1
c1s = u * (n - r * rprime)
assert c1s % s == 0
c1 = (c1s // s) % s
j = 1
if a1 == 0: a1 = s
assert 0 < a1 <= s
alist, blist, clist = [a0, a1], [b0, b1], [c0, c1]
# 2: Compute c
aj, bj, cj = alist[j], blist[j], clist[j]
if j % 2 == 0: c = cj
else: c = cj + s * ( (n + s**2 * (aj*bj-cj) ) // s**3 )
if c < 2 * aj * bj: goto 6  # TODO: This may belong under the "else" in the previous line.
A = c*s + aj*r + bj*rprime
B = aj * bj * n
D = A*A - 4*B
d = isqrt(D)
if d**2 != D: goto 5
t12, t22 = A + d, A - d
assert t12 % 2 == t22 % 2 == 0
t1, t2 = t12//2, t22//2
assert t1**2 - A * t1 + B == t2**2 - A * t2 + B == 0
# 4: Divisor found?
if t1 % aj == 0 and t2 % bj == 0 and (t1 // aj - r) % s == 0 and (t2 // bj - rprime) % s == 0: yield t1 // aj
# 5: Other value of c
if j % 2 == 0 and c > 0: c = c - s; goto 3
# 6: Next j
if aj == 0: return
j = j + 1
if j % 2 == 0: qj = alist[j-2] // alist[j-1]
else:          qj = (alist[j-2] - 1) // alist[j-1]
aj = alist[j-2] - qj * alist[j-1]
bj = blist[j-2] - qj * blist[j-1]
cj = clist[j-2] - qj * clist[j-1]
alist.append(aj)
blist.append(bj)
clist.append(cj)
goto 2

Translating this into proper Python yields

Code: Select all

def isqrt(n):   # Shamelessly stolen from https://codegolf.stackexchange.com/a/9088.
if n < 0: return int(n)
c = n*4//3
d = c.bit_length()
a = d>>1
if d&1:
x = 1 << a
y = (x + (n >> a)) >> 1
else:
x = (3 << a) >> 2
y = (x + (c >> a)) >> 1
if x != y:
x, y = y, (y + n//y) >> 1
while y < x: x, y = y, (y + n//y) >> 1
return x

def xgcd(a, b): # Algorithm 1.3.6
assert a >= 0 <= b
if b == 0: return (1, 0, d)
u, d, v, w = 1, a, 0, b
while w != 0:
q, r = divmod(d, w)
u, d, v, w = v, w, u - q * v, r
return (u, (d - a*u) // b, d)

def moddiv(r, s, n): # Algorithm 9.1.29
assert 0 <= r < s < n < s**3 and xgcd(r,s)[2] == 1
# 1: Initialization
u = xgcd(r,s)[0]
rprime = (u * n) % s
j, a1, c1 = 1, (u * rprime) % s, u * (n - r * rprime)
assert c1 % s == 0
if a1 == 0: a1 = s
alist, blist, clist = [s, a1], [0, 1], [0, (c1 // s) % s]
while True:
# 2: Compute c
aj, bj, cj, flag = alist[j], blist[j], clist[j], True
c = cj if j % 2 == 0 else cj + s * ( (n + s**2 * (aj*bj-cj) ) // s**3 )
if c < 2 * aj * bj and (j == 1): flag = False  # TODO: The second condition should perhaps be omitted.
while flag:
A, B = c*s + aj*r + bj*rprime, aj * bj * n
D = A*A - 4*B
d = isqrt(D)
if d**2 == D: #goto 5
t1, t2 = (A + d) // 2, (A - d) // 2
assert (A - t1) * t1 == (A - t2) * t2 == B
# 4: Divisor found?
if t1 % aj == t2 % bj == (t1//aj - r) % s == (t2//bj - rprime) % s == 0: yield t1//aj
# 5: Other value of c
if j % 2 == 0 < c: c -= s
else: flag = False
# 6: Next j
if aj == 0: return
j += 1
qj = (alist[j-2] - (j % 2)) // alist[j-1]
aj = alist[j-2] - qj * alist[j-1]
bj = blist[j-2] - qj * blist[j-1]
cj = clist[j-2] - qj * clist[j-1]
alist.append(aj)
blist.append(bj)
clist.append(cj)
#goto 2

I tested this by evaluating list(moddiv(19, 101, 40320)), which should return [120], but it either produces a ZeroDivisionError (if run as shown) or enters an infinite loop (if run with the suggested omission). Clearly I did something wrong, but I can't figure out what. Can anybody help?

Xanthir
My HERO!!!
Posts: 5321
Joined: Tue Feb 20, 2007 12:49 am UTC
Contact:

### Re: Lenstra's algorithm for divisors in residue classes

The two chunks of code are, as far as i can tell, equivalent. I haven't checked the paper to see if your transcription was correct, tho.
(defun fibs (n &optional (a 1) (b 1)) (take n (unfold '+ a b)))

eta oin shrdlu
Posts: 450
Joined: Sat Jan 19, 2008 4:25 am UTC

### Re: Lenstra's algorithm for divisors in residue classes

LucasBrown wrote:I'm attempting to implement Lenstra's algorithm for finding divisors in residue classes as described as Algorithm 9.1.29 in Cohen's A Course in Computational Algebraic Number Theory (also available here in PDF form).

[...]

I tested this by evaluating list(moddiv(19, 101, 40320)), which should return [120], but it either produces a ZeroDivisionError (if run as shown) or enters an infinite loop (if run with the suggested omission). Clearly I did something wrong, but I can't figure out what. Can anybody help?

Comparing your code with the Cohen text I see at least one difference; your conditional is

Code: Select all

if c < 2 * aj * bj and (j == 1): flag = False
while in step 2 I think the text is saying

Code: Select all

if c < 2 * aj * bj and (j % 2 == 1): flag = False

The ZeroDivisionError is probably because there's (I'm guessing) an implicit check in step 4 that the divisors aj and bj are nonzero, i.e. the condition "aj | t1" expands to Python "aj and t1 % aj == 0":

Code: Select all

if aj and bj and t1 % aj == t2 % bj == (t1//aj - r) % s == (t2//bj - rprime) % s == 0: yield t1//aj

But with those changes it still misses 120.

I glanced through the referenced Lenstra paper and saw one glaring difference with the Cohen text (I haven't checked it in detail, though the algorithm description is pretty short): the recurrence relation for the cj, at the top of p.333 in the linked version (don't worry, it's only 10 pages long, pp.331-340), two equations before equation (1.3), is not

Code: Select all

cj = clist[j-2] - qj * clist[j-1]
but

Code: Select all

cj = (clist[j-2] - qj * clist[j-1]) % s

With this change the code correctly finds 120 as a factor in your example, and more generally finds 120 for all s in range(37,120) with r=120%s and (r,s)=1. It typically finds all or almost all of the divisors in the specified residue class; the exceptions are that it sometimes doesn't find r if r is a divisor (e.g. moddiv(2,59,40320) finds 120 but not 2; however, moddiv(7,113,40320) finds both 7 and 120) and it sometimes doesn't find large factors (e.g. moddiv(37,83) finds 120 but not 10080). These may be problems with endpoints (e.g. maybe there should be a < instead of <= or vice versa) but I haven't looked carefully enough at the Lenstra paper to say.

(Also I've only looked at n=40320; looking at some other numbers would help to see if this pattern continues.)

[ETA missing "== 0" in a condition above.]