寻找原根的更快算法

Posted 2023-02-22

【中文标题】寻找原根的更快算法

【英文标题】：Faster algorithm for finding primitive roots

【发布时间】：2014-10-29 16:51:17

【问题描述】：

我正在尝试用这个算法找到素根：

std::vector<unsigned long long> Keyexchange::primroot(unsigned long long val) 

    std::vector<unsigned long long> res;

    for (unsigned long long i = 2; i<val - 1; i++) 

        unsigned long long start = 1;
        bool flag = 1;

        for (unsigned long long j = 0; j<val / 2; j++) 
            start = (start * i) % val;
            if (start % val == 1) 
                flag = 0;
                break;
            
        
        if (flag) 
            res.push_back(i);
        
    
    return res;

效果很好，但速度非常慢。 我想计算像 1073741789 这样的大数的原始根。如果有可能设置一个范围，那将是最好的，因为我现在正在计算整个集合。

所以基本上我正在寻找一种方法 [代码片段会很棒] 从给定的大数字中生成大约 100.000 个最大的原始根。

我知道使用 Eulersche φ 函数会快得多，但我不知道如何实现它。

非常感谢。

【问题讨论】：

如果输入数字是半素数并且你知道它的（两个）素数因子，那么使用 φ 函数是可行的。你输入的数字是这样的吗？

【参考方案1】：

首先，如果你选择一个从 2 到 p-1 的随机整数，那么它很有可能成为原始根。所以你选择一个随机整数（或者你从 2 开始），检查它，如果失败，你选择下一个等等。

检查 x 是否为原根：这意味着 x^(p-1) = 1 (模 p)，但 p 的幂没有更小。以 p = 31, p-1 = 30 = 2 x 3 x 5 为例。如果 p 不是原根，则 x^(30/2)、x^(30/3) 和 x^(30 /5) 必须为 1（模 p）。

在其素因数中分解p-1，对每一个素因数f计算x^((p-1)/f) (模p)，如果结果都不为1，则x为原根。

当然x^y（模p）需要通过重复平方/乘法来计算。例如，要计算 x^10，您将按该顺序计算 x^2、x^4、x^5、x^10。

一旦你找到了一个原根 g，如果 gcd (k, p-1) = 1，g^k 就是一个原根。但是你关心多个原根的情况很少见。

【参考方案2】：

如果输入数字是semi-prime，并且你手头有它的（两个）质因数，那么你可以使用这个：

vector<uint64> Roots(uint64 p,uint64 q)

    vector<uint64> roots;

    uint64 zstar = p*q;
    for (uint64 y=1; y<zstar; y++)
    
        if (GCD(zstar,y) == 1 && InQR(y,p,q))
        
            uint64 yp = PowMod(y,(p+1)/4,p);
            uint64 yq = PowMod(y,(q+1)/4,q);
            uint64 r1 = Map(0+yp,0+yq,p,q);
            uint64 r2 = Map(0+yp,q-yq,p,q);
            uint64 r3 = Map(p-yp,0+yq,p,q);
            uint64 r4 = Map(p-yp,q-yq,p,q);
            roots.push_back(r1);
            roots.push_back(r2);
            roots.push_back(r3);
            roots.push_back(r4);
        
    

    return roots;

以下是辅助功能：

uint64 GCD(uint64 a,uint64 b)

    uint64 c = a%b;
    if (c == 0)
        return b;
    return GCD(b,c);


uint64 PowMod(uint64 x,uint64 e,uint64 n)

    uint64 y = 1;
    while (e > 0)
    
        if (e & 1)
            y = (y*x)%n;
        x = (x*x)%n;
        e >>= 1;
    
    return y;


bool InQR(uint64 y,uint64 p)

    return PowMod(y,(p-1)/2,p) == 1;


bool InQR(uint64 y,uint64 p,uint64 q)

    return InQR(y,p) && InQR(y,q);


uint64 Map(uint64 u,uint64 v,uint64 p,uint64 q)

    uint64 a = q*Inverse(p,q);
    uint64 b = p*Inverse(q,p);
    return (u*a+v*b)%(p*q);


uint64 Inverse(uint64 n,uint64 a)

    int64  x1 = 1;
    int64  x2 = 0;
    int64  y1 = 0;
    int64  y2 = 1;
    uint64 r1 = n;
    uint64 r2 = a;

    while (r2 != 0)
    
        uint64 r3 = r1%r2;
        uint64 q3 = r1/r2;
        int64  x3 = x1-q3*x2;
        int64  y3 = y1-q3*y2;

        x1 = x2;
        x2 = x3;
        y1 = y2;
        y2 = y3;
        r1 = r2;
        r2 = r3;
    

    return (uint64)(y1>0? y1:y1+n);