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

Precedence

I'm curious to look over the full algorithm, but haven't had a chance to yet. As a shorthand for roots other than the square root, you can always use 6 ** (1.0/3). Raising to the 1/3 power -> cube root. It's possible you may lose a bit of precision though: (5 ** (1.0/3)) ** 3 -> 4.9999999999999982 – Anastigmatic

Have you looked into this: scipy.special.cbrt (docs.scipy.org/doc/scipy/reference/generated/…)? – Yorkist

are you sure p is prime? pow(2,p-1,p) != 1 so either the pow function is broken (doubtful) or p is not prime. Pari also thinks that p is composite. – Exhilarate

Ummmmmmmmmmmmmmmmmmmmmmmmmmmmm... it may be a composite of 3 primes. Is that bad? I need to check my algorithms. – Precedence

I'm an idiot, you're a genius, and your code works perfectly (once I used the right number)!!! Thank you again! – Precedence

Answer 1

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"

Exhilarate

You are amazing! I had printed those papers out and was planning on a long, frustrating night trying to understand them, and probably failing. I used your code on some smaller numbers and it worked perfectly! When I used larger numbers I exceeded the recursion limit in Python. I figured out how to increase the limit, and the result is Python 3.2 crashes, and Python 2.6 returns a result which seemed incorrect and when cubed doesn't become my original number. Do I need to try a new version of Python or a different computer? Should I give you my number? Thanks again! – Precedence

The python pow() function can take a 3rd argument for a modulus, making it equivalent to your powerMod(). – Helsell

@seth The second paper has a fairly straghtforward algorithm written out. If you replace the powerMod with pow and in phkahler's suggestion the recursion problems should go away. I'm not sure of the errors, I only tested it on small primes. I've got to go to work. – Exhilarate

@seth could you give a large example that fails? – Exhilarate

I added a set of numbers to my question. Unless I really screwed something up, which is possible, that number is a cube followed by the prime/finite field. Your latest version has no issues with recursion. I haven't really been able to look at it today but I will go over it again tomorrow and Friday and re-check all of my assumptions. Thank you for all of your help! – Precedence

Never mind, once I used the correct number it worked perfectly, thank you! – Precedence

I'll get to the complicated part this weekend if I have time. – Exhilarate

This code provides one of the three cube roots. How to get the other two? – Gader

@SheblaTsama multiply by a primitive cube root of 1. – Exhilarate

@Exhilarate please, can you explain what to do if this function returns None. Is that mean that number has no cube roots? – Plasticizer

@Denis Yes. In most cases when p is 1 mod 3 there is no cube root. – Exhilarate

Does any one know how to get the other two roots?! As I brute force for p = 31, I see there are 3 roots in some cases – Hachure

@PouJa, if r is 1 of the 3 cubic roots and m is the multiplicative order of r, the other 2 cubic roots would be r^(1 + m/3) and r^(1 + 2*m/3) – Regenerator

Still waiting for the last part of the code. Would you please complete that part? – Hachure

Answer 2

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

Bestead

This answer has some deficiencies. For example it cannot calculate cube roots mod 19 correctly. – Hachure

@pouJa I spent an afternoon checking if the algorithm works, but I can't find any fault. Do you have an example showing that the algorithm is incorrect? – Bestead

I checked it again and sorry my mistake. I had to correct indentation of your code after pasting the code to my editor and maybe I have made an incorrect change. Please correct the indentations in the code block. For example in the first code block every thing after #---MAIN--- should be indented by one tab. – Hachure

Hi @Rolandb, I have re-implemented this in Julia, and compared it to a brute force approach. I have a discrepancy for p=307 and p=739 (so the case where p mod 27 = 10). For example, for p=307 and a=2 I get no solutions with your algorithm, but 52, 270, 292 are solutions. If I remove the special-case-speed-up it works (using Tonelli-Shanks). Did I make a mistake? – Bloodstream

@ReinerMartin I have checked your examples (2.307) and (2.739), and the result is correct. user3435487's Python3 version also works correctly. It is possible that the indentations are not proper when copying the code. – Bestead

Answer 3

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://mcmap.net/q/730956/-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]

Footstep

Hi @user3435487, I have re-implemented this in Julia, and compared it to a brute force approach. I have a discrepancy for p=307 and p=739 (so the case where p mod 27 = 10). For example, for p=307 and a=2 I get no solutions with your algorithm, but 52, 270, 292 are solutions. If I remove the special-case-speed-up it works (using Tonelli-Shanks). Did I make a mistake? – Bloodstream

Answer 4

Sympy has a nice implementation for arbitrary integer modulo and arbitrary power: https://docs.sympy.org/latest/modules/ntheory.html#sympy.ntheory.residue_ntheory.nthroot_mod

from sympy.ntheory.residue_ntheory import nthroot_mod

a = 17
n = 3
modulo = 10000019
roots = nthroot_mod(a, n, modulo)
print(roots)

# 5188997

Demantoid