HDU 3542 （费马降阶法）

probelm

Given an odd prime p.
You are to tell whether p can be written into form a^2+b^2 (0

Input

There will be up to 2*10^5 test cases.
An odd prime p (3<=p<2*10^9) in each line.
End with EOF.

Output

For each number p,
if p can be written into form a^2+b^2 (0<=b),
output one line that contains “Legal a b”,
and if not just output one line that contains “Illegal”.
For a certain p that can be written into the form described above, any pair of a,b (0<=b) will do.

Sample Input

3
5
7
11
13
17
19
23
29

Sample Output

Illegal
Legal 1 2
Illegal
Illegal
Legal 2 3
Legal 1 4
Illegal
Illegal
Legal 2 5

思路

费马降阶法

//将素数表示成平方之和
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef pair<ll,ll> pll;
#define mp make_pair

ll mod_mul(ll a,ll b,ll c){//a*b %c 乘法改加法 防止超long long
    ll res=0;
    a=a%c;
    assert(b>=0);
    while(b)
    {
        if(b&1) res=(res+a)%c;
        b>>=1;
        a=(a+a)%c;
    }
    return res;
}

ll fast_exp(ll a,ll b,ll c){
    ll res=1;
    a=a%c;
    assert(b>=0);
    while(b>0)
    {
        if(b&1) res=a*res%c;
        b=b>>1;
        a=a*a%c;
    }
    return res;
}

//求解x^2 \equiv -1 (mod p)
ll ran(ll n){
    //要保证n%4=1
    assert(n%4==1);
    srand(time(0));
    for(;;){
        ll a=rand()%(n-1)+1;//1~ n-1
        ll b=fast_exp(a,(n-1)/4,n);
        if(b*b%n==n-1) return b<=(n-1)/2 ? b:n-b;
    }
}

pll solve(ll p)
{
    //算法保证M<p
    if(p==2) return mp(1,1);
    if(p%4!=1) return mp(-1,-1);
    ll A,B,u,v,M;
    A=ran(p); B=1;
    //__int128 tmp=1;
    ll tmp=1;
    M=ll((tmp*A*A+tmp*B*B)/p);
    while(M>1)
    {
        u=(A%M+M)%M;v=(B%M+M)%M;//注意可能是负数
        if(2*u>M) u-=M;
        if(2*v>M) v-=M;
        assert(u<=M/2 && u>=-M/2);
        assert(v<=M/2 && v>=-M/2);
        ll ta=A;
        A=ll((tmp*u*A+tmp*v*B)/M);
        B=ll((tmp*v*ta-tmp*u*B)/M);
        M=ll((tmp*u*u+tmp*v*v)/M);
    }
    return make_pair(min(abs(A),abs(B)),max(abs(A),abs(B)));
}

int main()
{
//    freopen("1.in","r",stdin);
//    freopen("1.out","w",stdout);
    ll p;
    while(cin>>p){
        pll ans=solve(p);
        if(ans.first==-1) cout<<"Illegal"<<endl;
        else{
            ll x=ans.first;ll y=ans.second;
            cout<<"Legal "<<x<<" "<<y<<endl;
            assert(x<y);
            assert((x*x+y*y)==p);
        }
    }
    return 0;
}