1
#include<iostream>
2
#include<stdio.h>
3
#include<cmath>
4
#include<cstdlib>
5
#include<algorithm>
6
using namespace std;
7
long long mp(long long a,long long b,long long n)
8![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedBlockStart.gif)
![](http://www.cppblog.com/Images/OutliningIndicators/ContractedBlock.gif)
{
9
long long d=1,t=a;
10
while (b>0)
11![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{
12
if (b%2==1) d=d*t%n;
13
b/=2;
14
t=t*t%n;
15
}
16
return d;
17
}
18
long long tu, tv;
19
long long Gau(double x)
20![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedBlockStart.gif)
![](http://www.cppblog.com/Images/OutliningIndicators/ContractedBlock.gif)
{
21
if (x+1e-8>(long long)x) return (long long)x;
22
return (long long)x-1;
23
}
24
int main()
25![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedBlockStart.gif)
![](http://www.cppblog.com/Images/OutliningIndicators/ContractedBlock.gif)
{
26
long long p, M, a, B, r, A, x, y,ta,tb;
27
while (scanf("%lld", &p) == 1)
28![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{
29
if (p % 4 != 1)
30![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{
31
printf("Illegal\n");
32
continue;
33
}
34
B = 1;
35
a = 2;
36
A = mp(a , (p - 1)/ 4,p);
37
while ((A * A + 1) % p != 0)
38![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{
39
a++;
40
A = mp(a , (p - 1)/ 4,p);
41
}
42
M = ((A * A + 1) / p);
43
tu = A%M;
44
tv = B%M;
45
while (M>1)
46![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{
47
tu=A+M*Gau(0.5-double(A)/double(M));
48
tv=B+M*Gau(0.5-double(B)/double(M));
49
r=(tu*tu+tv*tv)/M;
50
ta=abs((tu*A+tv*B)/M);
51
tb=abs((tv*A-tu*B)/M);
52
A=ta;
53
B=tb;
54
M=r;
55
}
56
if(A>B)
57
swap(A,B);
58
printf("Legal %lld %lld\n", A, B);
59
}
60
}
61![](http://www.cppblog.com/Images/OutliningIndicators/None.gif)
62![](http://www.cppblog.com/Images/OutliningIndicators/None.gif)
posted on 2010-09-02 10:34
Lancelot 阅读(245)
评论(0) 编辑 收藏 引用