Lenstra's algorithm for divisors in residue classes

[LucasBrown]

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.
    # 3: Solve quadratic equation
    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:
            # 3: Solve quadratic equation
            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]

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.

[eta oin shrdlu]

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