Корень кубика по модулю P - как мне это сделать?
Я пытаюсь вычислить кубический корень из числа многих сотен по модулю P в Python, и с треском провалился.
Я нашел код для алгоритма Тонелли-Шенкса, который, предположительно, легко изменить с квадратных корней на кубические, но это ускользает от меня. Я просмотрел веб-и математические библиотеки и несколько книг безрезультатно. Код был бы замечательным, так же как и алгоритм, объясненный простым языком.
Вот код 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 ответа
Примечание добавлено позже: в алгоритме Тонелли-Шенкса и здесь предполагается, что p прост. Если бы мы могли быстро вычислить модульные квадратные корни для составных модулей, мы могли бы быстро вычислять числа. Я прошу прощения за предположение, что вы знали, что p было простым
Смотрите здесь или здесь. Обратите внимание, что числа по модулю p являются конечным полем с p элементами.
Редактировать: Смотрите также (это дедушка из этих бумаг.)
Легкая часть - когда p = 2 mod 3, тогда все является кубом, и корень куба 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"
Вот полный код на чистом питоне.
#assumes p prime, it returns all cube roots of a mod p
def cuberoots(a, p):
#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 xrange(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 xrange(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):
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 [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 xrange(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 []
Вы можете проверить это, используя:
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]
у ^3=17 по модулю 10000019
р%3=2 [5188997] ---> [17]
у ^3=17 по модулю 1839598566765178548164758165715596714561757494507845814465617175875455789047
р%3=1 [753801617033579226225229608063663938352746555486783903392457865386777137044, 655108821219252496141403783945148550782812009720868259303598196387356108990, 430688128512346825798124773706784225426198929300193651769561114101322543013] ---> [17, 17, 17]
Я преобразовал код 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]
Sympy имеет хорошую реализацию для произвольного целочисленного модуля и произвольной мощности: 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