无题noname
Description
给定一个N，求出所有1到N之间的x，使得x^2=1(mod N)。
Input
一行一个正整数，没有多余字符，表示N。
Output
从小到大输出所有的x，一行一个。
Sample Input
1 | 7 |
Sample Output
1 | 1 |
Data Constraint
【数据范围】
30%的数据N<=20000；
100%的数据N<=2000000000。
Analysis
这道题有两种解题方法
***
第一种方法(扩展GCD)
$x^2\equiv1 (mod N)\longleftrightarrow x^2-1=N*T\longleftrightarrow (x+1)(x-1) = a*b*x0*(-y0)$
将N分解为ab，$T\in N$ 将T分解为x0(-y0),这样可以得到两个等式。
$$ \left\{
\begin{aligned}
x+1 = a*x0 \quad①\
x-1 = b*(-y0) \quad②\
\end{aligned}
\right.
$$
$a*x0+b*y0=2 \qquad ③=①-②$
***
下面是扩展GCD的一些知识
③可以用扩展GCD解得一组特解x0,y0，然后能通过特解的出一般解。
一般解为
$$ \left\{
\begin{aligned}
x0 + \frac {b}{p}*k \quad①\
y0 - \frac {a}{p}*k \quad②\
\end{aligned}
\right.
$$
***
扩展GCD$a*x+b*y=d$解得一组整数解x,y的充分必要条件是gcd（a,b）| d。
所以我们可以通过解出x，y，从而计算x的值
$2*x = a*x0-b*y0\qquad ④=①+②$
$x = \frac {a*x0-b*y0}{2}$
将$x$化为一般解，即：
$x = \frac {a*(x0 + \frac {b}{p}*k)-b*(y0- \frac {a}{p}*k)}{2}$
因为 $1\leq x\leq N$
所以可以得到k的取值范围
![取值范围]()
最后所求的x可化为：
$x = \frac {a*x0-b*y0}{2}+\frac {a*b*k}{p}$
详情代码1
***
第二种方法(中国剩余定理)
详情请参见References列表中的链接
详情代码2
Code
代码1
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
using namespace std;
typedef long long LL;
LL n,k,tmp,x3,y3,a,b,p;
int ans[1000];
LL gcd(LL a,LL b)
{
if (!b)
{
x3 = 2/a;
y3 = 0;
return a;
}
LL ret = gcd(b,a%b),t = x3;
x3 = y3;
y3 = t-a/b*y3;
return ret;
}
void work(LL a,LL b)
{
p = gcd(a,b);
if (2%p!=0) return;
LL x = a*x3-b*y3;
if (x%2==1) return;
LL t = 2*a*b/p;
double kmin = 1.0*(2-a*x3+b*y3)/t,
kmax = 1.0*(2*n-a*x3+b*y3)/t;
for (int k=ceil(kmin);k<=floor(kmax);k ++)
{
ans[++ ans[0]] = x/2+a*b*k/p;
}
}
int main()
{
//freopen("1460.in","r",stdin);
//freopen("1460.out","w",stdout);
scanf("%lld",&n);
if (n==1) return 0;
for (int i=1;i<=sqrt(n);i ++)
{
if (n%i==0)
{
a = i;b = n/i;
work(a,b);work(b,a);
}
}
sort(ans+1,ans+1+ans[0]);
printf("%d\n",ans[1]);
for (int i=2;i<=ans[0];i ++)
if (ans[i]!=ans[i-1]) printf("%d\n",ans[i]);
return 0;
}
***
代码2
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
var n,m,t,t1,t2,tt,tt2,i,j:longint;
mo:int64;
b,d,g,mm:array[0..10001] of int64;
c:array[0..10001,0..4] of int64;
p:array[0..100001] of boolean;
function exgcd(x,y:longint;var tx,ty:longint):longint;
var t1:int64;
begin
if y=0 then begin
tx:=1;ty:=0;
exit(x);
end;
exgcd:=exgcd(y,x mod y,tx,ty);
t1:=tx;tx:=ty;ty:=t1-(x div y)*ty
end;
procedure qsort(x,y:longint);
var i,j,k:longint;
begin
i:=x;j:=y;
k:=g[(i+j) div 2];
while i<=j do begin
while g[i]<k do inc(i);
while g[j]>k do dec(j);
if i<=j then begin
g[0]:=g[i];g[i]:=g[j];g[j]:=g[0];
inc(i);dec(j);
end;
end;
if i<y then qsort(i,y);
if x<j then qsort(x,j);
end;
procedure sy;
var i:longint;
ans:int64;
begin
ans:=0;
for i:=1 to tt do begin
exgcd(mm[i],b[i],t1,t2);
ans:=(ans+((d[i]*t1) mod mo*mm[i]) mod mo+mo*10) mod mo;
end;
inc(tt2);
g[tt2]:=ans;
end;
procedure dfs(x:longint);
var i:longint;
begin
if x=tt+1 then begin
sy;
exit;
end;
for i:=1 to c[x,0] do begin
d[x]:=c[x,i];
dfs(x+1);
end;
end;
begin
readln(n);
mo:=n;
t:=trunc(sqrt(n));
tt:=0;
m:=n;
for i:=2 to t do
if not p[i] then begin
if m mod i=0 then begin
inc(tt);
t2:=i;
b[tt]:=1;
while m mod i=0 do begin
m:=m div i;
b[tt]:=b[tt]*t2;
end;
end;
j:=i;
while j<=t do begin
p[j]:=true;
j:=j+i;
end;
end;
if m<>1 then begin
inc(tt);
b[tt]:=m;
end;
for i:=1 to tt do begin
if b[i] mod 2=0 then begin
if b[i]=2 then begin
c[i,0]:=1;
c[i,1]:=1;
end else if b[i]=4 then begin
c[i,0]:=2;
c[i,1]:=1;
c[i,2]:=3;
end else begin
c[i,0]:=4;
c[i,1]:=1;
c[i,2]:=b[i] div 2-1;
c[i,3]:=b[i] div 2+1;
c[i,4]:=b[i]-1;
end;
end else begin
c[i,0]:=2;
c[i,1]:=1;
c[i,2]:=b[i]-1;
end;
end;
for i:=1 to tt do mm[i]:=n div b[i];
dfs(1);
qsort(1,tt2);
for i:=1 to tt2 do writeln(g[i]);
end.
References
1 |
|
1 | var n,m,t,t1,t2,tt,tt2,i,j:longint; |