Cube root modulo P -- how do I do this?
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]
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"
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]