摘自数论概论的内容:

素数的两平方数之和定理：设p是素数，则p是两平方数之和的充要条件是p=1(mod 4) (或 p = 2).

两平方数之和定理实际上由两个陈述组成：

陈述1：如果p是两平方数之和，则p = 1(mod 4).

证明：设p = a^2 + b^2,p是奇数，所以a，b为一奇一偶，设a为奇数， b为偶数.比如 a= 2*n+1 b = 2*m.p = a^2 + b^2 = 4n^2+4n+1+4m^2 = 1 (mod 4).

陈述2：如果p=1(mod 4)，则p是两平方数之和. 这个的证明很麻烦，主要依据费马降阶法，可以参考数论概论第26章。

简单的说，如果p=1(mod 4)，不直接获得p是两平方数之和，而是将p的某个倍数表示成两个平方数之和。由二次互反律知x^2=-1(mod p)有一解，令x = a,b = 1,

a*a + b*b = Mp.利用费马降阶不断减小p的倍数使其可以表示两平方数之和，最终使p变成两平方数之和。如何利用已知的a,b,M来产生新的a,M.有恒等式：

(v^2+v^2)(a^2+b^2) = (ua+vb)^2 + (va-ub)^2.降阶程序有5个断言，只列出内容：1）a^2 + b^2 = Mp; 应用恒等式，我们选取的u，v满足u=a(mid M),v= b(mod M)

-M/2<= u,v,<= M/2. 于是有,u^2 + v^2 = a^2 + b^2= 0 (mod M),u^2 + v^2能被M整除，设u^2 + v^2 = Mr.其余四个断言陈述：2)r>=1; 3)r < M; 4)ua + vb能被M整除，

5）va-ub能被M整除。

#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <cstdio>
#include <cstdlib>
#include <cmath>
using namespace std;
 
typedef __int64 lint;
 
lint pow_mod(lint r,lint x,lint p) {
    lint pm = 1;
    while (x) {
        if (x&1)
            pm = (pm*r)%p;
        r = r*r%p;
        x >>= 1;
    }
    return pm;
}
 
int main()
{
    lint p,a,r,x,s,M,u,k,y,z,pm;
    while (scanf("%I64d",&p) != EOF) {
        if ((p-1)%4)
            printf("Illegal\n");
        else {
            b = 1;
            srand(NULL);
            r = rand()%(p-2)+2;
            x = (p-1)>>2;
            pm = pow_mod(r,p);
            while ((pm*pm)%p != p - 1) {
                r = rand()%(p-1)+1;
                pm = pow_mod(r,p);
            }
            a = pm;
            s = a*a + b*b;
            while (s != p) {
                M = s/p;
                k = M>>1;
                u = (a%M + M)%M;
                v = (b%M + M)%M;
                if (u > k)
                    u = M - u;
                if (v > k)
                    v = M - v;
                if ((u*a + v*b)%M)
                    swap(a,b);
                y = (u*a + v*b)/M;
                z = (v*a - u*b)/M;
                s = y*y + z*z;
                a = y;
                b = z;
            }
            if (a < 0)
                a = -a;
            if (b < 0)
                b = -b;
            if (a > b)
                swap(a,b);
            printf("Legal %I64d %I64d\n",b);
        }
    }
    return 0;
}