模P的立方根 -- 如何实现？
我正在尝试在Python中计算一个几百位数字的立方根，并且遇到了很多困难。
我找到了一种叫做Tonelli-Shanks算法的代码，听说这个算法从计算平方根改成计算立方根很简单，但我还是搞不明白。我在网上、数学库和几本书里找了很多资料，但都没有找到解决办法。如果能有代码示例就太好了，或者用简单的语言解释一下这个算法也行。
以下是用于计算平方根的Python（2.6？）代码：
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
4 个回答
2
我把上面Rolandb的代码转换成了python3。如果你把这个放到一个文件里，你可以在python3中导入并运行它。如果你单独运行这个文件，它会验证代码是否正常工作。
#! /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]
4
这里有一段用纯Python写的完整代码。通过先考虑一些特殊情况，它的速度几乎和Peralta算法一样快。
#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)
这是对Tonelli-Shank算法的扩展。
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 []
你可以用以下方式来测试它：
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)
结果应该是（很快）：
y^3=17 模 1459
p%3=1 [483, 329, 647] ---> [17, 17, 17]
y^3=17 模 1000003
p%3=1 [785686, 765339, 448981] ---> [17, 17, 17]
y^3=17 模 10000019
p%3=2 [5188997] ---> [17]
y^3=17 模 1839598566765178548164758165715596714561757494507845814465617175875455789047
p%3=1 [753801617033579226225229608063663938352746555486783903392457865386777137044, 655108821219252496141403783945148550782812009720868259303598196387356108990, 430688128512346825798124773706784225426198929300193651769561114101322543013] ---> [17, 17, 17]
13
后面补充说明：在Tonelli-Shanks算法中，以及这里的讨论中，假设p是一个质数。如果我们能快速计算出复合模数的模平方根，那我们就能快速分解数字。我为假设你知道p是质数而感到抱歉。
可以查看这里或这里。注意，模p的数字构成了一个有p个元素的有限域。
编辑：也可以看看这个（这是那些论文的祖先）。
简单的部分是当p等于2模3时，所有东西都是立方体，a的立方根就是a**((2*p-1)/3) %p。
补充：这里有代码可以处理除了1模9的质数以外的所有情况。我会尽量在这个周末完成。如果没有其他人先做的话。
#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"