16

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

5
  • 1
    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
    – g.d.d.c
    Commented Jul 19, 2011 at 18:55
  • 1
    Have you looked into this: scipy.special.cbrt (docs.scipy.org/doc/scipy/reference/generated/…)?
    – yasouser
    Commented Jul 19, 2011 at 19:40
  • 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.
    – deinst
    Commented Jul 20, 2011 at 14:00
  • Ummmmmmmmmmmmmmmmmmmmmmmmmmmmm... it may be a composite of 3 primes. Is that bad? I need to check my algorithms.
    – Seth
    Commented Jul 21, 2011 at 5:03
  • I'm an idiot, you're a genius, and your code works perfectly (once I used the right number)!!! Thank you again!
    – Seth
    Commented Jul 21, 2011 at 5:28

4 Answers 4

13

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"
14
  • 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!
    – Seth
    Commented Jul 20, 2011 at 2:18
  • The python pow() function can take a 3rd argument for a modulus, making it equivalent to your powerMod().
    – phkahler
    Commented Jul 20, 2011 at 2:50
  • 2
    This code provides one of the three cube roots. How to get the other two? Commented Jan 17, 2018 at 22:31
  • 2
    @SheblaTsama multiply by a primitive cube root of 1.
    – deinst
    Commented Jan 18, 2018 at 20:46
  • 1
    @Denis Yes. In most cases when p is 1 mod 3 there is no cube root.
    – deinst
    Commented Jun 21, 2018 at 12:31
4

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]

5
  • This answer has some deficiencies. For example it cannot calculate cube roots mod 19 correctly.
    – PouJa
    Commented Mar 8, 2022 at 15:41
  • @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?
    – Rolandb
    Commented Mar 26, 2022 at 15:02
  • 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.
    – PouJa
    Commented Apr 3, 2022 at 17:41
  • 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? Commented May 7, 2023 at 23:17
  • @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.
    – Rolandb
    Commented May 9, 2023 at 6:28
2

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]
1
  • 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? Commented May 7, 2023 at 23:19
2

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.