2023-12-02：用 go 语言，如何求模立方根？ x^3=a mod p， p 是大于等于 3 的大质数， a 是 1 到 p-1 范围的整数常数， x 也是 1 到 p-1 范围的整数，求 x。 p 过大，x 不能从 1 到 p-1 遍

作者：福大大架构师每日一题

2023-12-02
北京
本文字数：3126 字
阅读完需：约 10 分钟

2023-12-02：用 go 语言，如何求模立方根？

x^3=a mod p，

p 是大于等于 3 的大质数，

a 是 1 到 p-1 范围的整数常数，

x 也是 1 到 p-1 范围的整数，求 x。

p 过大，x 不能从 1 到 p-1 遍历。

答案 2023-12-02：

灵捷3.5

大体步骤如下：

1.判断是否存在模立方根。有 0，1，3 个根这三种情况。

1.1.求 p-1 和 3 的最大公约数 gcd(p-1,3)。最后结果要么是 1，要么是 3。如果是 1，那肯定模立方根，但只有 1 个根。如果是 3，进行下一步。

1.2.欧拉判别法。a**[(p-1)/3]==1 mod p。如果等于 1，那就有 3 个根。如果不等于 1，那就是 0 个根。

2.Peralta 算法。求 y。

2.1.当只有 0 个根时，直接返回。

2.2.当只有 1 个根时，a ^ ((p-1)/3) mod p 就是答案。

2.3.当有 3 个根时，这个很难描述，具体见代码。

2.3.1.定义复数乘法和复数的快速幂。这虽然叫复数，但跟传统意义上的复数是不一样的。

2.3.2.确定一个常数 r（r>=1 并且 r<p），使得 x ^ 3=r ^ 3 - a mod p 无根。

2.3.3.确定一个复数根，对这个复数根作复数的快速幂运算，指数是(p^2+p+1)/3，最终结果就是需要的根。

时间复杂度为 O((log p)^3)。

额外空间复杂度为 O(1)。

go 完整代码如下：

package main
import (  "fmt"  "math/big")
func main() {  if true {
    if false {      p := big.NewInt(0)      p.SetString("FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFC2F", 16)      for c := big.NewInt(20000); c.Cmp(big.NewInt(30000)) <= 0; c.Add(c, big.NewInt(1)) {        fmt.Println("c = ", c, "-------------")        r := ModCbrt(c, p)        fmt.Println("答案：", r)        for i := 0; i < len(r); i++ {          if big.NewInt(0).Exp(r[i], big.NewInt(3), p).Cmp(c) == 0 {
          } else {            fmt.Println("答案错误", r[i], "，c = ", big.NewInt(0).Exp(r[i], big.NewInt(3), p))            return          }        }      }      return    }    if true {      p := big.NewInt(0)      p.SetString("FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEBAAEDCE6AF48A03BBFD25E8CD0364141", 16)      for c := big.NewInt(20000); c.Cmp(big.NewInt(30000)) <= 0; c.Add(c, big.NewInt(1)) {        fmt.Println("c = ", c, "-------------")        r := ModCbrt(c, p)        fmt.Println("答案：", r)        for i := 0; i < len(r); i++ {          if big.NewInt(0).Exp(r[i], big.NewInt(3), p).Cmp(c) == 0 {
          } else {            fmt.Println("答案错误", r[i], "，c = ", big.NewInt(0).Exp(r[i], big.NewInt(3), p))            return          }        }      }      return    }
    if true {      p := big.NewInt(997)      for c := big.NewInt(1); c.Cmp(big.NewInt(0).Add(p, big.NewInt(-1))) <= 0; c.Add(c, big.NewInt(1)) {        fmt.Println("c = ", c, "-------------")        r := ModCbrt(c, p)        fmt.Println("答案：", r)        for i := 0; i < len(r); i++ {          if big.NewInt(0).Exp(r[i], big.NewInt(3), p).Cmp(c) == 0 {
          } else {            fmt.Println("答案错误", r[i], "，c = ", big.NewInt(0).Exp(r[i], big.NewInt(3), p))            return          }        }      }    }    return  }  fmt.Println("")}
// 求模立方根的个数0，1，3func ModCbrtCount(c, p *big.Int) int {  t := big.NewInt(0)  t.Add(p, big.NewInt(-2))  t.Mod(t, big.NewInt(3))  if t.Cmp(big.NewInt(0)) == 0 {    return 1  }  t = big.NewInt(0).Add(p, big.NewInt(-1))  t.Div(t, big.NewInt(3))  if big.NewInt(0).Exp(c, t, p).Cmp(big.NewInt(1)) == 0 {    return 3  } else {    return 0  }}
// Peralta Methodfunc ModCbrt(a, p *big.Int) (ans []*big.Int) {  ans = make([]*big.Int, 0)  count := ModCbrtCount(a, p)  if count == 1 { //有1个解    t := big.NewInt(0).Lsh(p, 1)    t.Mod(t, p)    t = t.Add(t, big.NewInt(-1))    t.Mod(t, p)    t.Mul(t, big.NewInt(0).ModInverse(big.NewInt(3), p))    t.Mod(t, p)    ans = append(ans, big.NewInt(0).Exp(a, t, p))  } else if count == 3 { //有3个解，Peralta Method算法
    w := big.NewInt(0)    p3 := big.NewInt(0).Add(p, big.NewInt(-1)) //(p-1)/3    p3.Mul(p3, big.NewInt(0).ModInverse(big.NewInt(3), p))    p3.Mod(p3, p)    for i := big.NewInt(1); i.Cmp(p) < 0; i.Add(i, big.NewInt(1)) {      w.Exp(i, p3, p)      if w.Cmp(big.NewInt(1)) != 0 {        break      }    }    var x *big.Int    key := big.NewInt(0)    for x = big.NewInt(1); x.Cmp(p) < 0; x.Add(x, big.NewInt(1)) {      key.Exp(x, big.NewInt(3), p) //key=x^3-a      key.Add(key, big.NewInt(0).Neg(a))      key.Mod(key, p)      if key.Cmp(big.NewInt(0)) != 0 && ModCbrtCount(key, p) == 0 {        break      }    }    r := Ring{x, big.NewInt(0).Add(p, big.NewInt(-1)), big.NewInt(0), key}    pp := big.NewInt(0).Mul(p, p) // pp = (p*p+p+1)/3，注意pp是不能 mod p的，有点反直觉    pp.Add(pp, p)    pp.Add(pp, big.NewInt(1))    pp.Div(pp, big.NewInt(3))    ansr := powerModI(r, pp, p)    ans0 := ansr.a    ans1 := big.NewInt(0)    ans1.Mul(ans0, w)    ans1.Mod(ans1, p)    ans2 := big.NewInt(0)    ans2.Mul(ans1, w)    ans2.Mod(ans2, p)    ans = append(ans, ans0, ans1, ans2)  }  return}
type Ring struct {  a *big.Int  b *big.Int  c *big.Int  w *big.Int}
// 复数乘法func mulI(x Ring, y Ring, p *big.Int) Ring {  var res Ring  res.a = big.NewInt(0)  res.b = big.NewInt(0)  res.c = big.NewInt(0)  res.w = x.w  w := x.w
  a1 := big.NewInt(0)  a2 := big.NewInt(0)  a3 := big.NewInt(0)  a1.Mul(x.a, y.a) //x.a*y.a  a1.Mod(a1, p)  a2.Mul(x.b, y.c) //x.b*y.c*key  a2.Mod(a2, p)  a2.Mul(a2, w)  a2.Mod(a2, p)  a3.Mul(x.c, y.b) //x.c*y.b*key  a3.Mod(a3, p)  a3.Mul(a3, w)  a3.Mod(a3, p)  res.a.Add(a1, a2)  res.a.Mod(res.a, p)  res.a.Add(res.a, a3)  res.a.Mod(res.a, p)
  b1 := big.NewInt(0)  b2 := big.NewInt(0)  b3 := big.NewInt(0)  b1.Mul(x.a, y.b) //x.a*y.b  b1.Mod(b1, p)  b2.Mul(x.b, y.a) //x.b*y.a  b2.Mod(b2, p)  b3.Mul(x.c, y.c) //x.c*y.c*key  b3.Mod(b3, p)  b3.Mul(b3, w)  b3.Mod(b3, p)  res.b.Add(b1, b2)  res.b.Mod(res.b, p)  res.b.Add(res.b, b3)  res.b.Mod(res.b, p)
  c1 := big.NewInt(0)  c2 := big.NewInt(0)  c3 := big.NewInt(0)  c1.Mul(x.a, y.c) //x.a*y.c  c1.Mod(c1, p)  c2.Mul(x.b, y.b) //x.b*y.b  c2.Mod(c2, p)  c3.Mul(x.c, y.a) //x.c*y.a  c3.Mod(c3, p)  res.c.Add(c1, c2)  res.c.Mod(res.c, p)  res.c.Add(res.c, c3)  res.c.Mod(res.c, p)
  return res}
// 复数快速幂，注意b不能取模func powerModI(a Ring, b, p *big.Int) Ring {  res := Ring{big.NewInt(1), big.NewInt(0), big.NewInt(0), a.w}  for b.Cmp(big.NewInt(0)) != 0 {    if big.NewInt(0).Mod(b, big.NewInt(2)).Cmp(big.NewInt(1)) == 0 {      res = mulI(res, a, p)    }    a = mulI(a, a, p)    b.Rsh(b, 1)  }  return res}