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

福大大架构师每日一题

发布于 2023-12-03 14:16:42

1440

发布于 2023-12-03 14:16:42

文章被收录于专栏：福大大架构师每日一题

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，3
func 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 Method
func 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
}