求大素数原根算法（python代码）

定义：

代码可直接看下面求素数p原根的最终方法

对于素数 p，如果存在一个正整数 1<a<p，
使得 a1,a2,…,ap−1a1,a2,…,ap−1 模 p 的值
取遍 1,2,…,p−11,2,…,p−1 的所有整数，
称 a 是 p 的一个原根（primitive root），就是循环群的生成元。

下面求原根是基于素数来说的

原理

欧拉定理

欧拉定理表明，若n,a为正整数，且互素，则
$a^{\phi(n)} = 1\ \ (mod\ n)\\ \phi(n)是欧拉函数，素数的\phi(n)=n-1$aϕ(n)=1  (mod n)ϕ(n)是欧拉函数，素数的ϕ(n)=n−1
所以对于素数来说

a的n-1次方模n 一定是1

对于 a^j ≡ a^i(mod p)，则 i≡j(mod p−1)。这里有两个例子：

• 3是7的原根，因为3–>2–>6–>4–>5–>1，然后开始循环
• 2不是7的原根，因为2–>4–>1–>2–>4–>1，过早的循环了

$对于素数p\\a^{p-1} ≡ 1 （mod\ \ p）\\$对于素数pap−1≡1（mod  p）

原根在p-1之前不会有循环

也就是原根的循环节为 p−1，非原根有较小的循环节，且是 p−1 的约数（因为元素的阶整除群的阶）

算法

（参考别人的博客）

所以我们可以求p-1非自身的约数m1,m2,m3…，如果出现
$a^{m_n} ≡ 1（mod\ \ p）$amn​≡1（mod  p）
则a不是原根，否则找到一个原根

考虑优化这个算法。

注意到若
$$
g^m≡1 (mod p)\

当m’=2m,3m,…时，均有g^{m’}≡1 (mod p)\
因此我们并不需要枚举p-1的所有约数\
只要质因数即可
$$

（我写完代码验证后出错）

不过我搜的几篇博客都这样说

验证出错

经验证，上面的优化算法是有问题的。
比如p=32683079412992931673的质因子为
m=[2, 3, 149, 9139563594237397]
s^m %p都不为1，
得到结果原根为2，
但其实p的原根为5，10等

而且
2^4085384926624116459% p = 1
4085384926624116459=9139563594237397 * 149 * 3

$也就是当m'=2m,3m,…时，均有g^{m'}≡1 (mod p)\\ 有可能g^{3m}≡1 (mod p) 我们只能省去小于p-1的 6m,9m\\而不能省去m,2m .\ \ \ \ \ \ \ \ \ 因为循环节不一定是质数\\ 所以稳妥一点，还是对所有的约数进行检测$也就是当m′=2m,3m,…时，均有gm′≡1(modp)有可能g3m≡1(modp)我们只能省去小于p−1的6m,9m而不能省去m,2m.         因为循环节不一定是质数所以稳妥一点，还是对所有的约数进行检测

我的思考

之后，我想出了一种更方便的检测方法
$如果a是原根,一定有 \ \ a^n\ !=\ 1\ \ (mod\ p) \\ 这里n是小于p-1的任何数$如果a是原根,一定有  an != 1  (mod p)这里n是小于p−1的任何数
p的循环节一定是p-1的约数
那么也就是说我们可以验证最大的p-1的约数（即p-1除以最小的p-1的质因数）
$如果\\ a^x\ !=\ 1\ \ (mod\ p) \ \ \ x是p-1最大约数\\ 那么a就一定是原根$如果ax != 1  (mod p)   x是p−1最大约数那么a就一定是原根

这么一来我们在分解p-1因数时只需要分解出来除1以外的最小因子就可以
而由于我们探讨是p是素数
p-1即为偶数
p-1的最小质因子是2
我们只需要验证 x=（p-1）/2即可

则求大素数p原根的最终方法

$for\ \ a\ in\ range(2,p-1):\\ \ \ \ 如果\ \ a^{\frac{(p-1)}2}\ =\ 1\ \ (mod\ p)\\ 则a不是原根，否则是素数p的一个原根$for  a in range(2,p−1):   如果  a2(p−1)​ = 1  (mod p)则a不是原根，否则是素数p的一个原根

不管怎么说这个我没有搜到，
是自己想出来，挺开心的，不晓得是不是原创算法呢。（想想还有点小激动）

python代码

def yg(n):		# 这样默认求最小原根
    k=(n-1)//2
    for i in range(2,n-1):
        if multimod(i,k,n)!=1:
            return i

print(yg(n))

multimod是快速幂求模的算法

python代码

def multimod(a,k,n):    #快速幂取模
    ans=1
    while(k!=0):
        if k%2:         #奇数
            ans=(ans%n)*(a%n)%n
        a=(a%n)*(a%n)%n
        k=k//2          #整除2
    return ans

具体原理可以看我另一篇博客

快速幂求模