## 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 xdef 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 2def 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 xdef 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) == 1    # 1: Initialization    u = xgcd(r,s)    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 , 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? https://pypi.python.org/pypi/labmath
3239×21412529+1 is prime!

Xanthir
My HERO!!!
Posts: 5426
Joined: Tue Feb 20, 2007 12:49 am UTC
Location: The Googleplex
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: 451
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 , 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.]

Return to “Mathematics”

### Who is online

Users browsing this forum: No registered users and 18 guests