Cube root modulo P -- how do I do this?

Question

I am trying to calculate the cube root of a many-hundred digit number modulo P in Python, and failing miserably.

I found code for the Tonelli-Shanks algorithm which supposedly is simple to modify from square roots to cube roots, but this eludes me. I've scoured the web and math libraries and a few books to no avail. Code would be wonderful, so would an algorithm explained in plain English.

Here is the Python (2.6?) code for finding square roots:

def modular_sqrt(a, p):
    """ Find a quadratic residue (mod p) of 'a'. p
        must be an odd prime.

        Solve the congruence of the form:
            x^2 = a (mod p)
        And returns x. Note that p - x is also a root.

        0 is returned is no square root exists for
        these a and p.

        The Tonelli-Shanks algorithm is used (except
        for some simple cases in which the solution
        is known from an identity). This algorithm
        runs in polynomial time (unless the
        generalized Riemann hypothesis is false).
    """
    # Simple cases
    #
    if legendre_symbol(a, p) != 1:
        return 0
    elif a == 0:
        return 0
    elif p == 2:
        return n
    elif p % 4 == 3:
        return pow(a, (p + 1) / 4, p)

    # Partition p-1 to s * 2^e for an odd s (i.e.
    # reduce all the powers of 2 from p-1)
    #
    s = p - 1
    e = 0
    while s % 2 == 0:
        s /= 2
        e += 1

    # Find some 'n' with a legendre symbol n|p = -1.
    # Shouldn't take long.
    #
    n = 2
    while legendre_symbol(n, p) != -1:
        n += 1

    # Here be dragons!
    # Read the paper "Square roots from 1; 24, 51,
    # 10 to Dan Shanks" by Ezra Brown for more
    # information
    #

    # x is a guess of the square root that gets better
    # with each iteration.
    # b is the "fudge factor" - by how much we're off
    # with the guess. The invariant x^2 = ab (mod p)
    # is maintained throughout the loop.
    # g is used for successive powers of n to update
    # both a and b
    # r is the exponent - decreases with each update
    #
    x = pow(a, (s + 1) / 2, p)
    b = pow(a, s, p)
    g = pow(n, s, p)
    r = e

    while True:
        t = b
        m = 0
        for m in xrange(r):
            if t == 1:
                break
            t = pow(t, 2, p)

        if m == 0:
            return x

        gs = pow(g, 2 ** (r - m - 1), p)
        g = (gs * gs) % p
        x = (x * gs) % p
        b = (b * g) % p
        r = m

def legendre_symbol(a, p):
    """ Compute the Legendre symbol a|p using
        Euler's criterion. p is a prime, a is
        relatively prime to p (if p divides
        a, then a|p = 0)

        Returns 1 if a has a square root modulo
        p, -1 otherwise.
    """
    ls = pow(a, (p - 1) / 2, p)
    return -1 if ls == p - 1 else ls

Source: Computing modular square roots in Python

Answer 1

Here is a complete code in pure python. By considering special cases first, it is almost as fast as the Peralta algoritm.

#assumes p prime, it returns all cube roots of a mod p
def cuberoots(a, p):

    #Non-trivial solutions of x**r=1
    def onemod(p,r):
        sols=set()
        t=p-2
        while len(sols)<r:        
            g=pow(t,(p-1)//r,p)
            while g==1: t-=1; g=pow(t,(p-1)//r,p)
            sols.update({g%p,pow(g,2,p),pow(g,3,p)})
            t-=1
        return sols

    def solutions(p,r,root,a): 
        todo=onemod(p,r)
        return sorted({(h*root)%p for h in todo if pow(h*root,3,p)==a})

#---MAIN---
a=a%p

if p in [2,3] or a==0: return [a]
if p%3 == 2: return [pow(a,(2*p - 1)//3, p)] #One solution

#There are three or no solutions 

#No solution
if pow(a,(p-1)//3,p)>1: return []

if p%9 == 7:                                #[7, 43, 61, 79, 97, 151]   
    root = pow(a,(p + 2)//9, p)
    if pow(root,3,p) == a: return solutions(p,3,root,a)
    else: return []

if p%9 == 4:                                #[13, 31, 67, 103, 139]
    root = pow(a,(2*p + 1)//9, p) 
    print(root)
    if pow(root,3,p) == a: return solutions(p,3,root,a)        
    else: return []        
            
if p%27 == 19:                              #[19, 73, 127, 181]
    root = pow(a,(p + 8)//27, p)
    return solutions(p,9,root,a)

if p%27 == 10:                              #[37, 199, 307]
    root = pow(a,(2*p +7)//27, p)  
    return solutions(p,9,root,a) 

#We need a solution for the remaining cases
return tonelli3(a,p,True)

An extension of Tonelli-Shank algorithm.

def tonelli3(a,p,many=False):

    def solution(p,root):
        g=p-2
        while pow(g,(p-1)//3,p)==1: g-=1  #Non-trivial solution of x**3=1
        g=pow(g,(p-1)//3,p)
        return sorted([root%p,(root*g)%p,(root*g**2)%p])

#---MAIN---
a=a%p
if p in [2,3] or a==0: return [a]
if p%3 == 2: return [pow(a,(2*p - 1)//3, p)] #One solution

#No solution
if pow(a,(p-1)//3,p)>1: return []

#p-1=3**s*t
s=0
t=p-1
while t%3==0: s+=1; t//=3
 
#Cubic nonresidu b
b=p-2
while pow(b,(p-1)//3,p)==1: b-=1

c,r=pow(b,t,p),pow(a,t,p)    
c1,h=pow(c,3**(s-1),p),1    
c=pow(c,p-2,p) #c=inverse modulo p

for i in range(1,s):
    d=pow(r,3**(s-i-1),p)
    if d==c1: h,r=h*c,r*pow(c,3,p)
    elif d!=1: h,r=h*pow(c,2,p),r*pow(c,6,p)           
    c=pow(c,3,p)
    
if (t-1)%3==0: k=(t-1)//3
else: k=(t+1)//3

r=pow(a,k,p)*h
if (t-1)%3==0: r=pow(r,p-2,p) #r=inverse modulo p

if pow(r,3,p)==a: 
    if many: 
        return solution(p,r)
    else: return [r]
else: return [] 

You can test it using:

test=[(17,1459),(17,1000003),(17,10000019),(17,1839598566765178548164758165715596714561757494507845814465617175875455789047)]

for a,p in test:
    print "y^3=%s modulo %s"%(a,p)
    sol=cuberoots(a,p)
    print "p%s3=%s"%("%",p%3),sol,"--->",map(lambda t: t^3%p,sol)

which should yield (fast):

y^3=17 modulo 1459
p%3=1 [483, 329, 647] ---> [17, 17, 17]
y^3=17 modulo 1000003
p%3=1 [785686, 765339, 448981] ---> [17, 17, 17]
y^3=17 modulo 10000019
p%3=2 [5188997] ---> [17]
y^3=17 modulo 1839598566765178548164758165715596714561757494507845814465617175875455789047
p%3=1 [753801617033579226225229608063663938352746555486783903392457865386777137044, 655108821219252496141403783945148550782812009720868259303598196387356108990, 430688128512346825798124773706784225426198929300193651769561114101322543013] ---> [17, 17, 17]

Answer 2

Note added later: In the Tonelli-Shanks algorithm and here it is assumed that p is prime. If we could compute modular square roots to composite moduli quickly in general we could factor numbers quickly. I apologize for assuming that you knew that p was prime.

See here or here. Note that the numbers modulo p are the finite field with p elements.

Edit: See this also (this is the grandfather of those papers.)

The easy part is when p = 2 mod 3, then everything is a cube and athe cube root of a is just a**((2*p-1)/3) %p

Added: Here is code to do all but the primes 1 mod 9. I'll try to get to it this weekend. If no one else gets to it first

#assumes p prime returns cube root of a mod p
def cuberoot(a, p):
    if p == 2:
        return a
    if p == 3:
        return a
    if (p%3) == 2:
        return pow(a,(2*p - 1)/3, p)
    if (p%9) == 4:
        root = pow(a,(2*p + 1)/9, p)
        if pow(root,3,p) == a%p:
            return root
        else:
            return None
    if (p%9) == 7:
        root = pow(a,(p + 2)/9, p)
        if pow(root,3,p) == a%p:
            return root
        else:
            return None
    else:
        print "Not implemented yet. See the second paper"