轉載自https://blog.csdn.net/new_ke_2014/article/details/21829921

#include <iostream>
using namespace std;
#define LL __int64
LL w;
struct Point//x + y*sqrt(w)
{
    LL x;
    LL y;
};

LL mod(LL a, LL p)
{
    a %= p;
    if (a < 0)
    {
        a += p;
    }
    return a;
}

Point point_mul(Point a, Point b, LL p)
{
    Point res;
    res.x 
 
= mod(a.x * b.x, p);
    res.x += mod(w * mod(a.y * b.y, p), p);
    res.x = mod(res.x, p);
    res.y = mod(a.x * b.y, p);
    res.y += mod(a.y * b.x, p);
    res.y = mod(res.y , p);
    return res;
}

Point power(Point a, LL b, LL p)
{
    Point res;
    res.x = 1;
    res.y = 0;
    while(b)
    {
        
 
if (b & 1)
        {
            res = point_mul(res, a, p);
        }
        a = point_mul(a, a, p);
        b = b >> 1;
    }

    return res;
}

LL quick_power(LL a, LL b, LL p)//(a^b)%p
{
    LL res = 1;
    while(b)
    {
        if (b & 1)
        {
            res = (res * a) % p;
        }

        a 
 
= (a * a) % p;
        b = b >> 1;
    }

    return res;
}


LL Legendre(LL a, LL p) // a^((p-1)/2)
{
    return quick_power(a, (p - 1) >> 1, p);
}

LL equation_solve(LL b, LL p)//求解x^2=b(%p)方程解
{
    if ((Legendre(b, p) + 1) % p == 0)
    {
        return -1;//表示沒有解
    }
    LL a, t;
    while(true)
    {
        a = rand() % p;
        t = a * a - b;
        t = mod(t, p);
        if ((Legendre(t, p) + 1) % p == 0)
        {
            break;
        }
    }

    w = t;
    Point temp, res;
    temp.x = a;
    temp.y = 1;
    res = power(temp, (p + 1) >> 1, p);
    return res.x;
}

int main()
{
    LL b, p;
    scanf("%I64d %I64d", &b, &p);
    printf("%I64d\n", equation_solve(b, p));//輸出為-1表示，不存在解
    return 0;

}