求解形如 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;
}