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

user_12386130 · Answer

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

user_3435487 · Answer

I converted the code by Rolandb above into python3. If you put this into a file, you can import it and run it in python3, and if you run it standalone it will validate that it works.

#! /usr/bin/python3

def ts_cubic_modular_roots (a, p):
  """ python3 version of cubic modular root code posted
      by Rolandb on stackoverflow.  With new formatting.
      https://stackoverflow.com/questions/6752374/cube-root-modulo-p-how-do-i-do-this

  """
  
  #Non-trivial solution of x**r = 1
  def onemod (p, r):
    t = p - 2
    while pow (t, (p - 1) // r, p) == 1: 
      t -= 1
    return pow (t, (p - 1) // r, p)

  def solution(p, root): 
    g = onemod (p, 3)
    return [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)] #Eén oplossing

  #There are 3 or no solutions 

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

  if p % 9 == 4:                                #[13, 31, 67]
    root = pow (a, ((2 * p) + 1) // 9, p)  
    if pow (root, 3, p) == a: 
      return solution (p, root)        
    else: 
      return []

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

  if p % 27 == 10:                              #[37, 199]
    root = pow (a, ((2 * p) + 7) // 27, p)         
    h = onemod (p, 9)
    for i in range (0,9):
      if pow (root, 3, p) == a: 
        return solution (p, root)                
      root *= h
    return []        

  if p % 27 == 19:                              #[19, 73, 127, 181]
    root = pow (a, (p + 8)//27, p)
    h = onemod (p, 9)
    for i in range (0, 9):
      if pow (root, 3, p) == a: 
        return solution (p, root)
      root *= h
    return []        

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


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

  #Non-trivial solution of x**r = 1
  def onemod (p, r):
    t = p - 2
    while pow (t, (p - 1) // r, p) == 1: 
      t -= 1
    return pow (t, (p - 1) // r, p)

  def solution (p, root):
    g = onemod (p, 3)
    return [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)] #Eén oplossing

  #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_mod(Integer(c), 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_mod(Integer(r), p)

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

if '__name__' == '__main__':
  import ts_cubic_modular_roots
  tscr = ts_cubic_modular_roots.ts_cubic_modular_roots
  test=[(17,1459),(17,1000003),(17,10000019),(17,1839598566765178548164758165715596714561757494507845814465617175875455789047)]
  for a,p in test:
    print ("y**3=%s modulo %s"%(a,p))
    sol = tscr (a,p)
    print ("p%s3=%s"%("%",p % 3), sol, [pow (t,3,p) for t in sol])

# results of the above
#y**3=17 modulo 1459
#p%3=1 [] []
#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]

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

3 Answers