求解二次同余式

求解形如 x^2 = a (mod p) 这样的同余式

/*模P平方根：    求 X ^2 = a (mod p)    定理：当P为！！！奇素数 ！！！的时候    先判断（a / p )的勒让德符号， 若为-1则无解，若为1则有解    分解P-1,然后求B，然后求出X(t-1),和a的逆元    然后开始求 ans = （ a的逆元 * 上一个X的平方（t-k））的(t-k-1)次方对P取模    如果  ans == -1 则 J = 1；    如果  ans ==  1 则 J = 0；    然后开始求现在的 X = （上一个X * B的（J*2的k次方）次方    直到求出X0，也就是最后的解*/#include<bits/stdc++.h>using namespace std;void Divide(int p, int &t, int &s){    t = 0, s= 0;    while(p % 2 == 0)    {        t++;        p /= 2;    }    s = p;    return ;}int Pow(int a, int b, int mod){    int ans = 1, base = a;    while(b != 0)    {        if(b & 1)            ans = (ans * base) % mod;        base = (base * base) % mod;        b >>= 1;    }    return ans;}int Legendre(int a, int p){    if(a == 2)    {        int x = (p+1)*(p-1)/8;        if(x % 2 == 0)            return 1;        else            return -1;    }    else    {        int ans = Pow(a, (p-1)/2, p);        if(ans == p-1)            return -1;        else            return 1;    }}int FindN(int p){    for(int i = 1; i < p; i++)    {        if(Legendre(i, p) == -1)            return i;    }}int e_gcd(int a, int b, int &x, int &y){    if(b == 0)    {        x = 1; y = 0;        return a;    }    int q = e_gcd(b, a%b, y, x);    y -= a/b*x;    return q;}int Inverse(int a, int m){    int x, y;    int gcd = e_gcd(a, m, x, y);    if(1 % gcd != 0)        return -1;    x *= 1/gcd;    m = abs(m);    int ans = x % m;    if(ans <= 0)        ans += m;    return ans;}int JudgeJ(int A, int x, int t, int p){    cout << A << " " << x << " " << t << " " << p << endl;    x = ((x % p) * (x % p) % p);    int xx = (A * x) % p;    int ans = Pow(xx, pow(2, t), p);    if(ans == p-1)        return 1;    else        return 0;}int main(){    int a, p;    printf("请输入所求算式的 a 和 p：\n");    while( scanf("%d %d", &a, &p) != EOF)    {        if(Legendre(a, p) == -1)        {            cout << "无解" << endl;            continue;        }        int t, s;        Divide(p-1, t, s); //求t和s        int n = FindN(p);  //找到那个不符合条件的n        int b = Pow(n, s, p);        int *X;        X = (int *) malloc(sizeof(int) * (t+5));        X[t-1] = Pow(a, (s+1)/2, p);        t--;        int A = Inverse(a, p);          //求A的逆元        int k = 0;        while(t > 0)        {            int j = JudgeJ(A, X[t], t-1, p);            int B = Pow(b, j * pow(2, k), p);            X[t-1] = ((X[t] % p) * (B % p)) % p;            t--; k++;        }        printf("所求的解为：");        cout << X[0] << endl;    }    return 0;}