c语言寻找素数的一个生成元（本原元、本原根）

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

long long s, t;

long long llabs(long long x)
{
    return x < 0 ? -x : x;
}

long long multiply(long long a, long long b, long long mod)
{
    if (llabs(a) > mod / 2)
    {
        a = a < 0 ? a + mod : a - mod;
    }
    long long x = llabs(a), y = llabs(b);
    int a_0 = x % 2;
    long long result = a_0 * y % mod;
    x = x / 2;
    while (x > 0)
    {
        y = (y % mod * 2) % mod;
        if (x % 2 == 1)
        {
            result = (result + y) % mod;
        }
        x = x / 2;
    }
    if (a < 0)
    {
        result = -result;
    }
    if (b < 0)
    {
        result = -result;
    }

    if (llabs(result) > mod / 2)
    {
        result = result < 0 ? result + mod : result - mod;
    }
    return result;
}

long long power(int *flag, long long b, long long n, long long mod) // 这里n>=1,是奇数
{
    long long A = 1;
    int i;
    while (n > 1)
    {
        if (n % 2 == 1)
        {
            A = multiply(A, b, mod);
        }
        if (A == -1 || A == 1)
        {
            *flag = 0;
            return 0;
        }
        b = multiply(b, b, mod);
        if (b == -1 || b == 1)
        {
            *flag = 0;
            return 0;
        }
        n = n / 2;
    }
    A = multiply(A, b, mod);
    if (A == 1)
    {
        *flag = 0;
        return 0;
    }
    return A;
}

long long gcd(long long a, long long b)
{
    if (b == 0)
        return a;
    return gcd(b, a % b);
}

long long pollard_rho(long long n)
{
    if (n % 2 == 0)
        return 2;
    long long x = rand() % (n - 3) + 2;
    long long y = x;
    long long c = rand() % (n - 1) + 1;
    long long d = 1;
    while (d == 1)
    {
        x = (x * x + c) % n;
        y = (y * y + c) % n;
        y = (y * y + c) % n;
        d = gcd(abs(x - y), n);
    }
    return d;
}

int factorize(long long t, long long A2, long long p)
{
    int flag = 1;
    long long A3, A4, m, n, r1, r2;
    if (t > 7)
    {
        long long factor = pollard_rho(t);
        if (factor != t)
        {
            m = t / factor;
            n = factor;
            if (m < n) // 保证m>=n
            {
                long long temp = m;
                m = n;
                n = temp;
            }
            A3 = power(&flag, A2, n, p);
            if (A3 == -1)
            {
                return 0;
            }
            if (m == n)
            {
                A4 = A3;
            }
            else
            {
                r1 = m / n;
                r2 = m % n;
                A4 = power(&flag, A3, r1, p);
                if (r1 % 2 == 0)
                {
                    A4 = multiply(A4, power(&flag, A2, r2, p), p);
                }
                else if (r2 % 2 == 0 && r2 > 0)
                {
                    A4 = multiply(multiply(A4, power(&flag, A2, r2 - 1, p), p), A2, p);
                }
            }
            if (A4 == -1)
            {
                return 0;
            }
            factorize(gcd(m, n), A2, p);
        }
    }
    return 1;
}

long long findGenerator(long long p) // p是素数
{
    if (p == 2)
    {
        return 1;
    }
    if (p == 3 || p == 5 || p == 11)
    {
        return 2;
    }
    if (p == 7)
    {
        return 3;
    }
    s = 0;
    t = p - 1;
    while (t % 2 == 0)
    {
        t = t / 2;
        s++;
    } // 易知s>=1,t>=1
    long long bound;
    // s>1且t>1时phi(2^s*t)/2>=(2^(s-1)-2^(s-2))*(t-1)/2；
    // s=1时，t一定大于等于3，phi(2^s*t)/2>=(t-1)/2
    // t=1时，s一定大于2，phi(2^s*t)/2>=2^(s-1)-2^(s-2),
    long long r1 = (p - 1) / (4 * t), r2 = (t - 1) / 2;
    if (s > 1 && t > 1)
    {
        bound = r1 * r2;
    }
    else if (s == 1)
    {
        bound = r2;
    }
    else if (t == 1)
    {
        bound = r1;
    }
    bound = (p - 1) / 2 + 1 - bound; // 可以证明bound>=6
    long long A, i, f, g;
    int flag;
    long long sqrtbound = floor(sqrt(bound)); // 下取整
    long long begin, end = 0;
    long long Max;
    long long sqrtT = floor(sqrt(t));
    if (sqrtbound * sqrtbound == bound)
    {
        Max = sqrtbound;
    }
    else
    {
        Max = sqrtbound + 1;
    }
    for (f = 2; f <= Max; f++)
    {
        begin = end + 2;
        if (f == sqrtbound + 1)
        {
            end = bound;
        }
        else
        {
            end = begin + 2 * f - 3;
        }
        for (g = begin; g <= end; g++)
        {
            // printf("g=%lld\n", g);
            if (g == 2 && llabs(4 - (p % 8)) == 3)
            {
                continue;
            }
            if (g == 3 && llabs(6 - (p % 12)) == 5)
            {
                continue;
            }
            flag = 1;
            A = power(&flag, g, t, p);
            if (flag == 0)
            {

                continue;
            }
            for (i = 0; i < s; i++)
            {
                if (A == -1)
                {
                    if (i == s - 1)
                    {
                        if (t > 1)
                        {
                            long long A2 = g, A3, A4, j, k;
                            for (j = 0; j < s - 1; j++)
                            {
                                A2 = multiply(A2, A2, p);
                            }
                            if (A2 == -1)
                            {
                                goto sign;
                            }
                            if (!factorize(t, A2, p)) // 最占用运行时间的函数
                            {
                                goto sign;
                            }
                        }
                        return g;
                    }
                    else
                    {
                    sign:
                        break;
                    }
                }
                if (i < s - 1)
                {
                    A = multiply(A, A, p);
                }
            }
        }
    }
    return -1;
}

int main()
{
    long long p = 66666666666666601; // p不能超过2^64严格的说是10进制不要超过9223372036854775807，否则会溢出
    double dur;
    clock_t start, end;
    start = clock(); // 统计运行时间
    long long g = findGenerator(p);
    end = clock();
    dur = (double)(end - start) / CLOCKS_PER_SEC * 1000;
    printf("Use Time:%f ms\n", dur);
    printf("The generator of %lld(2^%lld*%lld+1) is %lld\n", p, s, t, g);
}

时间复杂度O(p(logp)^2)，空间复杂度O(1)，如有疑问，请私信我！