立方根模 P——我该怎么做？

Question

我正在尝试用 Python 计算数百位数字模 P 的立方根，但惨败。

我找到了 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

来源：在 Python 中计算模平方根

Answer 1

稍后添加的注释：在 Tonelli-Shanks 算法中以及这里，假设 p 是素数。一般来说，如果我们能够快速计算模平方根以合成模，我们就可以快速分解数字。我很抱歉假设你知道 p 是素数。

请参阅此处或此处。请注意，以 p 为模的数字是具有 p 个元素的有限域。

编辑：另请参阅此（这是这些论文的祖父。 ）

最简单的部分是当 p = 2 mod 3 时，那么一切都是立方体，a 的立方根只是 a**((2*p-1)/3) %p

添加：这是执行所有操作的代码，但the primes 1 mod 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"

Answer 2

这是纯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 modulo 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]

Answer 3

我将上面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]

Answer 4

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