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

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

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

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

go完整代码如下:

packagemain
import(
"fmt"
"math/big"
)
funcmain(){
iftrue{
iffalse{
p:=big.NewInt(0)
p.SetString("FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFC2F",16)
forc:=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)
fori:=0;iifbig.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
}
iftrue{
p:=big.NewInt(0)
p.SetString("FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEBAAEDCE6AF48A03BBFD25E8CD0364141",16)
forc:=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)
fori:=0;iifbig.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
}
iftrue{
p:=big.NewInt(997)
forc:=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)
fori:=0;iifbig.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
funcModCbrtCount(c,p*big.Int)int{
t:=big.NewInt(0)
t.Add(p,big.NewInt(-2))
t.Mod(t,big.NewInt(3))
ift.Cmp(big.NewInt(0))==0{
return1
}
t=big.NewInt(0).Add(p,big.NewInt(-1))
t.Div(t,big.NewInt(3))
ifbig.NewInt(0).Exp(c,t,p).Cmp(big.NewInt(1))==0{
return3
}else{
return0
}
}
//PeraltaMethod
funcModCbrt(a,p*big.Int)(ans[]*big.Int){
ans=make([]*big.Int,0)
count:=ModCbrtCount(a,p)
ifcount==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))
}elseifcount==3{//有3个解,PeraltaMethod算法
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)
fori:=big.NewInt(1);i.Cmp(p)<0;i.Add(i,big.NewInt(1)){
w.Exp(i,p3,p)
ifw.Cmp(big.NewInt(1))!=0{
break
}
}
varx*big.Int
key:=big.NewInt(0)
forx=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)
ifkey.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是不能modp的,有点反直觉
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
}
typeRingstruct{
a*big.Int
b*big.Int
c*big.Int
w*big.Int
}
//复数乘法
funcmulI(xRing,yRing,p*big.Int)Ring{
varresRing
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)
returnres
}
//复数快速幂,注意b不能取模
funcpowerModI(aRing,b,p*big.Int)Ring{
res:=Ring{big.NewInt(1),big.NewInt(0),big.NewInt(0),a.w}
forb.Cmp(big.NewInt(0))!=0{
ifbig.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)
}
returnres
}

打开网易新闻 查看精彩图片