数论算法ch31

31.1 基础数论概念

先简要回顾一下书中内容：

整除性与约数：d|a 表示为d整除a，存在整数k,使得a=kd

                         若d≥0，则称d是a的约数。

素数与合数：素数：如果能被平凡约数1和自身整除即为素数。

                      合数：如果整数a>1且不是素数，则称之为合数。

除法定理，余数和等模

                     除法定理： 对于任何整数a和任何正整数n,存在唯一整数q和r，满足0≤r<n且a=qn+r，称q=(向下取整)(a/n)为除法的商，值r=a mod n为除法的余数。

                     等模：对于整数模n的余数，可以划分为n个等价类。包含整数a的模n等价类为[a]n={a+kn:k∈Z}

公约数与最大公约数：

                      公约数：d|a,a|b,则d是a与b的公约数。

               最大公约数：两个不同时为0的整数a与b的公约数中最大的称为最大公约数。记做gcd(a,b).

               相关定理：1.若任意整数a和b不都为0，则gcd(a,b)是a与b的线性组合集{ax+by:x,y∈Z}中的最小正元素。

                                  2.对任意整数a与b，如果d|a且d|b,，则d|gcd(a,b).

                                  3.对所有整数a与b以及任意非负整数n.有gcd(an,bn)=ngcd(a,b)

                                  4.对于任意正整数n,a,b .如果n|ab且gcd(a,n)=1,则n|b。

互质数：如果两个数的只有公约数1，则a与b称为互质数。

                相关定理：5.对于任意整数a,b和p,如果gcd(a,p)=1且gcd(b,p)=1则 gcd(ab,p)=1.

唯一因子分解定理

             相关定理：1.对所有素数p和所有整数a,b，如果p|ab,则p|a或p|b(或者两者都成立)。

                               2.合数a技能以一种方法写成如下乘积形式：a=p(e1,1)p(e2,2)....p(er,r) 其中pi为素数，p1<p2<...<pr,且er为正整数。

 31.1-1  证明：若a>b>0,且c=a+b,则c mod a=b.

设c mod a=x,则存在整数k,则c=ak+x=a+b  若k=1,则x=b,若k≠1 ,由于a>b>0,c=a+b>0 ,k=向下取整(c/a)>0.所以k≥2,所以ak+x=a+b≥2a+x 所以b≥a+x≥a 这样与题目假设矛盾。

所以k=1,x=b,c mod a=b.

31.1-2证明有无穷多个素数。(提示：证明素数p1,p2....,pk都不能整除(p1p2....pk)+1)

首先证明提示部分：假设pi(i=1,2,...k)能整除(p1p2....pk)+1),则(p1p2....pk)+1)=k*pi=>pi(k-p1p2..p(i-1)...p(i+1)...pk)=1 pi=1/(k-p1p2..p(i-1)...p(i+1)...pk)>1所以k-p1p2..p(i-1)...p(i+1)...pk<1又因为pi(k-p1p2..p(i-1)...p(i+1)...pk)=1>0且pi>0 所以(k-p1p2..p(i-1)...p(i+1)...pk)>0 所以k-1<(p1p2..p(i-1)...p(i+1)...pk)<k ,又因为(p1p2..p(i-1)...p(i+1)...pk)为整数，但是它在(k-1,k)这个没有整数的区间里，所以与假设矛盾。

现在证明原题：假设只有有限个素数，p1,p2...pk.则(p1p2....pk)+1)比p1,p2...pk有限个素数都大的数是合数。但是由提示知：p1,p2....,pk任何一个素数都不能整除(p1p2....pk)+1

根据定理31.8。由于合数能分解成一组素数，那么(p1p2....pk)+1也不能被合数整除，所以(p1p2....pk)+1不能被任何除1和它本身的数整除，所以(p1p2....pk)+1是素数。得证！

31.1-3 证明：如果a|b且b|c，则a|c.

a|b =>存在整数k1,b=k1*a. 同理 存在整数k2,c=k2*b所以c=k2*b=k2*(k1*a)=(k1*k2)*a 所以a|c

31.1-4 证明：如果p是素数并且0<k<p,则gcd(k,p)=1.

若k为素数，因为0<k<p且p是素数,由于k的约数为1和k,p的约数为1和p,因为k≠p，所以两个素数只有1为公约数，所以gcd(k,p)=1.

若k为合数，因为0<k<p且p是素数,由于k的约数是1和数个小于k的数，而p的约数为1和p,因为k的所有约数都小于p，所以gcd(k,p)=1.

31.1-5 证明：对于任意正整数n,a,b .如果n|ab且gcd(a,n)=1,则n|b。

因为n|ab,所以存在整数k,使得ab=kn ,gcd(a,n)=1  =>gcd(a,ab/k)=1 =>kgcd(a,ab/k)=k  =>gcd(ak,ab)=k  => a*gcd(k,b)=k

所以ab=kn => ab= a*gcd(k,b)n  =>  b=gcd(k,b)n  => n|b

31.1-6 证明：如果p是素数且0<k<p,则.证明对所有整数a,b和素数p,有(a+b)^p≡a^p+b^p(mod p).

是组合数，所以是整数，所以k!(p-k)!|p!,p是素数且对于任意i∈(1,p-1),有gcd(p,i)=1 所以gcd(p,k!(p-k)!)=1(这里有个定理我不得不说，因为书上没有相关定理,gcd(m,a)=1,则有gcd(ab,m)=gcd(b,m))则对于任意k∈[1,p-1],k!(p-k)!)|(p-1)!(这里也有个定理，需要特别说明，gcd(m,a)=1,m|ab,则m|b)。所以存在整数z,使(p-1)!=z*k!(p-k)!.两边同乘以p,则有p!=zk!(p-k)!p => p!/(k!(p-k)!)=zp =>p| p!/(k!(p-k)!)=>

下面证明第二个问题 因为所以所以存在整数k,使(a+b)^p=a^p+b^p+kp. 因为b^p+kp≡ b^p(mod p)，所以可证结论。

31.1-7证明：如果a和b是任意正整数，且满足a|b,则对任意x,(x mod b) mod a= xmod a 在相同的假设下，证明对任意整数x和y,若果x≡y(mod b),则x≡y(mod a).

设x mod b=y则存在整数k1,使得x=bk1+y.设y mod a=z,则存在整数k2使得y=ak2+z 所以x=bk1+ak2+z 又因为a|b,则存在整数k. b=ak 所以x=(ak)k1+ak2+z=a(kk1+k2)+z所以x mod a=z=y mod a=(x mod b) mod a.

在相同假设下，x≡y(mod b) =>存在整数k1使得y≡bk1+x. 又因为a|b,则存在整数k. b=ak  使得y≡a(kk1)+x,所以x≡y(mod a).

31.1-8 对任意整数k>0,如果存在一个整数a,满足a^k=n,则称整数n是一个k次幂。如果对于某个整数k>1,n>1是一个k次幂，则称n是非平凡幂。说明如何在关于β的多项式时间内判定一个β位整数n是否是非平凡幂。

以下是代码：

[cpp] view plaincopyprint?

1. //最朴素的多项式时间内判断一个数是否为某个数的幂的形式：就是用枚举法，挨个找,但是这个是关于n的多项式，关于β的多项式暂时没有想出。  
2. #include <iostream>  
3. #include <math.h>  
4. using namespace std;  
5. void main()  
6. {  
7.     int n=64,flag=1;  
8.     //O(√nlgn)  
9.     for (int i=2;i<=sqrt(n);i++)//O(a=√n)  
10.     {  
11.         int m=n,k=0;  
12.         while (m%i==0)//O(k=lgn)  
13.         {  
14.             m=m/i;  
15.             k++;  
16.         }  
17.         if (m==1&&k>1)  
18.         {  
19.             cout<<n<<"是非平凡幂,存在一个整数"<<i<<"它的"<<k<<"次幂="<<n<<endl;  
20.             flag=0;  
21.             break;  
22.         }  
23.     }  
24.     if (flag)  
25.     {  
26.         cout<<"n="<<n<<"不存在非平凡幂"<<endl;  
27.     }  
28. }  

31.1-9证明等式(31.6)-31.10

需要证明的等式已用加粗

31.6  gcd(a,b)=gcd(b,a) 根据整数的交换率便知

31.7 设d是a与b的约数，则d|a,d|b.存在整数k有a=dk =>|a|=d(±k) =>d||a|,所以d|(±a),可见d|(-a)又因为d|b且-a与b的所有约数都一样，那么最大约数也一样故gcd(-a,b)=gcd(a,b)

31.8 根据31.7知：设d是a与b的约数,若d|a=>d||a| 同理d|b=>d||b| ,所以d是|a|与|b|的约数，既然所有约数都一样，那么最大的约数当然也一样了。gcd(a,b)=gcd(|a|,|b|)得证！

31.9设d是a的约数，则显然d也是0的约数（因为任何整数乘以0都得0,）所以d也是a与0的公约数，从中找出最大的,a的最大约数为它本身也就是a(当然a>0)所以a=gcd(a,0),若a<0,-a=gcd(-a,0)=gcd(a,0)(根据31.7) 所以对于任意整数a,有|a|=gcd(|a|,0)=gcd(a,0)(根据31.8)

31.10 若a>0,则 gcd(a,ka)=a*gcd(1,k)(根据推论31.4)，又因为任意整数与1的最大公约数肯定是1，gcd(a,ka)=a

          若a<0,则 gcd(a,ka)=gcd(-a,-ka)=(-a)*gcd(1,k)(根据推论31.4)，又因为任意整数与1的最大公约数肯定是1，gcd(a,ka)=-a

          若a=0,则gcd(0,0)=0，所以gcd(a,ka)=a 总结：gcd(a,ka)=|a|

31.1-10 证明：最大公约数运算满足结合律，即证明对所有整数a,b和c。gcd(a,gcd(b,c))=gcd(gcd(a,b),c)

在证明之前，先证明：gcd(a,b,c)=gcd(a,gcd(b,c))

设d=gcd(a,gcd(b,c)),则d|a,d|gcd(b,c).所以d|a,d|b,d|c.所以d是a,b,c的约数，所以gcd(a,b,c)≥gcd(a,gcd(b,c))。。。。（1）

设d'=gcd(a,b,c).则d'|a,d'|b,d'|c 由推论31.3知d'|gcd(b,c) 又因为d'|a,所以再次使用d'|gcd(a,gcd(b,c)) .因为最大公约数一定是正整数，所以存在正整数k,使得gcd(a,gcd(b,c))=kd'

d'≤gcd(a,gcd(b,c)) 即gcd(a,b,c)≤gcd(a,gcd(b,c)).。。.。。（2） 由（1）与（2）知gcd(a,b,c)=gcd(a,gcd(b,c)) 同理可证：gcd(a,b,c)=gcd(gcd(a,b),c)

所以gcd(a,b,c)=gcd(a,gcd(b,c))=gcd(gcd(a,b),c) 。

31.2-11 证明定理31.8

注：这个定理的证明可参考北师大《初等数论》第6讲，里面有详细解答。这里略过。

31.1-12 试写出计算β位整除除以短整数的高效算法，以及计算β位整数除以短整数的余数的高效算法。所给出的算法运行时间应为θ(β^2).(感觉31.1-12和31.1-13应该用FFT算法解决。)

既然用高效的算法，那就用位运算。

[cpp] view plaincopyprint?

1. //位运算的乘法与除法  
2. #include <iostream>  
3. using namespace std;  
4. //位运算的乘法  
5. int bit_Multiplication(int a,int b)  
6. {  
7.      int ans=0;  
8.      for (int i=1;i;i<<=1,a<<=1)  
9.      {  
10.          if (b&i)  
11.          {  
12.              ans+=a;  
13.          }  
14.      }  
15.      return ans;  
16. }  
17. //位运算的除法  
18. int bit_Division1(int x,int y)  
19. {  
20.     int ans=0;  
21.     for (int i=31;i>=0;i--)  
22.     {  
23.         if ((x>>i)>=y)  
24.         {  
25.             ans+=(1<<i);  
26.             x-=(y<<i);  
27.         }  
28.     }  
29.     return ans;  
30. }  
31. //计算整数的二进制位数  
32. int bit_num(int d)  
33. {  
34.    int i=0;  
35.    while (d)  
36.    {  
37.        d>>=1;  
38.        i++;  
39.    }  
40.    return i;  
41. }  
42. //位运算的除法 计算商  
43. int bit_Division2_quotient(int x,int y)  
44. {  
45.     int c2=bit_num(x),c1=bit_num(y),quotient=0;  
46.     for (int i=c2-c1;i>=0;i--)//i=c2-c1防止除数y移位后超过无符号整数最大值 时间复杂度O(c2-c1)  
47.     {  
48.         unsigned int a=(y<<i);//有了i=c2-c1保证了y<<i不会溢出 a有c1+c2-c1=c2位  
49.         if (a<=x)  
50.         {  
51.             quotient+=(1<<i);  
52.             x-=a;  
53.         }  
54.     }  
55.     //总的时间复杂度为 O(c2)=O(x的二进制位数)=O(b^2) b为除数的十进制位数  
56.     return quotient;  
57. }  
58. //位运算的除法 计算余数 与计算商一样，只是返回值不同  
59. int bit_Division2_Remainder(int x,int y)  
60. {  
61.     int c2=bit_num(x),c1=bit_num(y),quotient=0;  
62.     for (int i=c2-c1;i>=0;i--)//i=c2-c1防止除数y移位后超过无符号整数最大值 时间复杂度O(c2-c1)  
63.     {  
64.         unsigned int a=(y<<i);//有了i=c2-c1保证了y<<i不会溢出 a有c1+c2-c1=c2位  
65.         if (a<=x)  
66.         {  
67.             quotient+=(1<<i);  
68.             x-=a;  
69.         }  
70.     }  
71.     //总的时间复杂度为 O(c2)=O(x的二进制位数)=O(b^2) b为除数的十进制位数  
72.     return x;  
73. }  
74. void main()  
75. {  
76.     cout<<bit_Multiplication(350,43)<<endl;  
77.     cout<<bit_Division1(350,43)<<endl;  
78.     cout<<"商："<<bit_Division2_quotient(350,43)<<endl;  
79.     cout<<"余数："<<bit_Division2_Remainder(350,43)<<endl;  
80. }  

31.1-13 写出一个高效算法，用于将β位二进制整数转化为响应的十进制表示。证明：如果长度至多为β的整数的乘法或除法运算所需时间为M(β)，则执行二进制到十进制转换所需时间为θ(M(β)lgβ)。(提示：应用分治法，分别使用独立的递归计算结果的前段和后段)(感觉31.1-12和31.1-13应该用FFT算法解决。)
这次要用位运算+分治思想。我对我写的这段代码持怀疑态度。如果牛人看到实际没用到分治思想，那么给出理由和比较好的建议。谢谢！

[cpp] view plaincopyprint?

1. //用分治思想进行进制转换(2->10),写得不好，凑合看吧  
2. /*#include <iostream>  
3. #include <string>  
4. using namespace std;  
5. #define BIT 6//二进制整数的位数，可根据所要输入的二进制位数设置BIT。  
6. int t=-1;  
7. int Bit_merge(int a[],int p,int r)  
8. {  
9.     static ans=0;  
10.     p=(p>t)?p:(t+1);  
11.     for (int i=p;i<=r;i++)  
12.     {  
13.         t++;  
14.         ans+=a[i]<<i;       
15.     }  
16.     return ans;  
17. }  
18. int bit_Multiplication(int a[],int p,int r,int flag)//x为二进制数  
19. {  
20.     int q;  
21.     if (p<r)  
22.     {  
23.        q=(p+r)/2;  
24.        bit_Multiplication(a,p,q,0);  
25.        bit_Multiplication(a,q+1,r,1);  
26.     }  
27.     if (p<=r&&t!=BIT)return flag==0?Bit_merge(a,p,q):Bit_merge(a,q+1,r);  
28. }  
29.   
30. void main()  
31. {  
32.     int a[BIT]={0};  
33.     string x;  
34.     cin>>x;  
35.     int j=0;  
36.     while (j!=BIT)  
37.     {  
38.        a[j]=x[BIT-1-j]-'0';  
39.        j++;  
40.     }  
41.     cout<<endl;  
42.     cout<<bit_Multiplication(a,0,BIT-1,0)<<endl;  
43. }  

如果β位整数的乘法或除法的运行时间为M(β)，那么用分治法的递归式为T(β)=2T(β/2)+θ(M(β))=>T(β)=θ(M(β)lgβ)

31.2最大公约数

欧几里得递归算法：

定理31.9（GCD递归定理）对任意非负整数a和任意正整数b. gcd(a,b)=gcd(b,a mod b)

引理31.10 如果a>b≥0,并且EUCLID(a,b)执行了k≥1次递归调用，则a≥F(k+2),b≥F(k+1)

定理31.11(Lame定理) 对任意整数k≥1,如果a>b≥1,且b<F(k+1),则EUCLID(a,b)的递归调用次数少于k次。

代码如下：

[cpp] view plaincopyprint?

1. //欧几里得算法递归形式  
2. #include <iostream>  
3. using namespace std;  
4. int Euclid(int a,int b)  
5. {  
6.     cout<<"gcd("<<a<<","<<b<<")=";  
7.     if (b==0)  
8.     {  
9.         return a;  
10.     }  
11.     else  
12.     {  
13.         return Euclid(b,a%b);  
14.     }  
15. }  
16. void main()  
17. {  
18.    cout<<Euclid(30,21)<<endl;  
19. }  

欧几里得算法扩展形式

代码如下：

[cpp] view plaincopyprint?

1. #include <iostream>  
2. using namespace std;  
3. struct a  
4. {  
5.     int d,x,y;  
6. }s;  
7. struct a extended_eucild(int a,int b)  
8. {  
9.     if(b==0)  
10.     {  
11.         s.d=a,s.x=1,s.y=0;  
12.         return s;   
13.     }  
14.     else  
15.     {  
16.         struct a ss=extended_eucild(b,a%b);  
17.         s.d=ss.d;  
18.         s.x=ss.y;  
19.         s.y=ss.x-(a/b)*ss.y;  
20.         return s;  
21.     }  
22. }  
23. int Fibonacci(int n)  
24. {  
25.     if (n==1)  
26.     {  
27.         return 1;  
28.     }  
29.     else if (n==0)  
30.     {  
31.         return 0;  
32.     }  
33.     else  
34.     {  
35.         return Fibonacci(n-2)+Fibonacci(n-1);  
36.     }  
37. }  
38. void main()  
39. {  
40.   struct a s=extended_eucild(99,78);  
41.   cout<<s.d<<" "<<s.x<<" "<<s.y<<endl;  
42. }  
43. <strong><span style="font-size:18px;color:#000000;">31.2-1 证明：由式(31.11)和式(31.12)可推得式(31.13）</span></strong>  

31.2-2计算调用过程EXTENDED-EUCLID(899,493)的返回值为(d,x,y).

ab向下取整a/bdxy899493129-6114934061295-640687429-1587581291-1582922901290一2910                  

31.2-3 证明：对所有整数a,k和n，gcd(a,n)=gcd(a+kn,n)

gcd(a,n)=gcd(n,a mod n)=gcd(a mod n,n)   对所有整数a,k',n和x,设a mod n=x,则a=nk'+x=>x=a-nk，所以gcd(a,n)=gcd(a-nk',n),因为k'为所有整数，则令k=-k'

则gcd(a,n)=gcd(a+nk,n).

31.2-4仅用常数大小的存储空间(即仅存储常数个整数值)把过程EUCLID改写成迭代形式。

代码如下：

[cpp] view plaincopyprint?

1. //欧几里得算法迭代形式  
2. #include <iostream>  
3. using namespace std;  
4. int Euclid(int a,int b)  
5. {  
6.     while (b)  
7.     {  
8.         int c=b;  
9.         b=a%b;  
10.         a=c;  
11.     }  
12.     return a;  
13. }  
14. void main()  
15. {  
16.    cout<<Euclid(89,55)<<endl;  
17. }  

31.2-5 如果a>b≥0，证明：EUCLID(a,b)至多执行1+logb次递归调用。把这个界改进为1+log(b/gcd(a,b)).

改进的界我暂时不会。

31.2-6 过程EXTENDED-EUCLID(F(k+1),Fk)返回什么值？证明答案的正确性。

[cpp] view plaincopyprint?

1. #include <iostream>  
2. using namespace std;  
3. struct a  
4. {  
5.     int d,x,y;  
6. }s;  
7. struct a extended_eucild(int a,int b)  
8. {  
9.     if(b==0)  
10.     {  
11.         s.d=a,s.x=1,s.y=0;  
12.         return s;   
13.     }  
14.     else  
15.     {  
16.         struct a ss=extended_eucild(b,a%b);  
17.         s.d=ss.d;  
18.         s.x=ss.y;  
19.         s.y=ss.x-(a/b)*ss.y;  
20.         return s;  
21.     }  
22. }  
23. int Fibonacci(int n)  
24. {  
25.     if (n==1)  
26.     {  
27.         return 1;  
28.     }  
29.     else if (n==0)  
30.     {  
31.         return 0;  
32.     }  
33.     else  
34.     {  
35.         return Fibonacci(n-2)+Fibonacci(n-1);  
36.     }  
37. }  
38. int power(int k)  
39. {  
40.     int i=1;  
41.     int t=1;  
42.     while (i!=k+1)  
43.     {  
44.        t*=-1;  
45.        i++;  
46.     }  
47.     return t;  
48. }  
49. void main()  
50. {  
51.   const int k=10;  
52.   struct a s=extended_eucild(Fibonacci(k+1),Fibonacci(k));  
53.   cout<<s.d<<" "<<s.x<<" "<<s.y<<endl;  
54.   cout<<"d="<<s.d<<" x="<<power(k-1)*Fibonacci(k-2)<<" y="<<power(k)*Fibonacci(k-1)<<endl;  
55. }  

31.2-7利用递归等式gcd(a0,a1,...,an)=gcd(a0,gcd(a1,...,an))定义多于两个变量的gcd函数。说明gcd函数的返回值与其参数次序无关。同时说明如何找出满足gcd(a0,a1,...,an)=a0x0+a1x1+...+anxn的整数x0,x1...,xn。证明所给出的算法执行除法运算次数为O(n+lg(max{a0,a1,....,an})).

在证明 说明gcd函数的返回值与其参数次序无关之前，有个定理必须说明：若D=gcd(a0,a1,...an) 则若d|ai(i=0,1,...n)，则d|D.还要说明的gcd(a0,a1,..an)=gcd(|a0|,|a1|,..|an|)也就是只需要求出ai的所有正整数的最大公约数即可。

下面需要证明如下等式：gcd(a0,a1,...,an)=gcd(ai,gcd(a1,.a(i-1),a(i+1)..,an)).

1.设d=gcd(a0,a1,...,an) 则d|ai(i=0,1,...n),d是ai的公约数，所以d|aj(j=0,1..n且(j≠i))且d|ai

所以d|gcd(a0,a1..a(i-1),a(i+1)..an),由于且d|ai ，

所以d|gcd(ai,gcd(a1,.a(i-1),a(i+1)..,an)),由此可知：gcd(ai,gcd(a1,.a(i-1),a(i+1)..,an))≥d=gcd(a0,a1,...,an) ..(1)

2.设d‘=gcd(ai,gcd(a1,.a(i-1),a(i+1)..,an))则d'|ai,d'|gcd(a1,.a(i-1),a(i+1)..,an)=>d'|ai(i=0,1,..n)=>d'是ai的公约数

所以的d'≤d，gcd(ai,gcd(a1,.a(i1),a(i+1)..,an))≤gcd(a0,a1,...,an)...(2)  

由(1)和(2)知gcd(ai,gcd(a1,.a(i-1),a(i+1)..,an))=d=gcd(a0,a1,...,an) 其中i=0,1....n 也就是说a0,a1...an的最大公约数就是等于其中任意一个ai与剩下的gcd(a1,.a(i-1),a(i+1)..,an)的最大公约数.

下面给出满足gcd(a0,a1,...,an)=a0x0+a1x1+...+anxn的算法代码如下：

[cpp] view plaincopyprint?

1. #include <iostream>  
2. #include <time.h>  
3. using namespace std;  
4. #define n 5//根据输入确定元素个数。  
5. struct t  
6. {  
7.     int d;  
8.     int x;  
9.     int y;  
10.     int z[n];//存放题目中的x0,x1...xn  
11. }s;  
12. struct t extended_eucild(int a,int b)  
13. {  
14.     if (b==0)  
15.     {  
16.         s.d=a;  
17.         s.x=1;  
18.         s.y=0;  
19.         return s;  
20.     }  
21.     else  
22.     {         
23.         struct t ss=extended_eucild(b,a%b);  
24.         s.d=ss.d;  
25.         s.x=ss.y;  
26.         s.y=ss.x-(a/b)*ss.y;  
27.         return s;  
28.     }  
29. }  
30. struct t extended_eucild1(int a[])  
31. {  
32.     struct t s3,s1,s2;  
33.     for (int i=0;i<n;i++)  
34.     {  
35.         s3.z[i]=0;  
36.     }  
37.     s1=extended_eucild(a[0],a[1]);//O(lgb)  
38.     s3.z[0]=s1.x;  
39.     s3.z[1]=s1.y;  
40.     int k=3;  
41.     while (k<n+1)  
42.     {     
43.         s2=extended_eucild(s1.d,a[k-1]);//O(lg(min{s1.d,a[k-1]}))........(3)  
44.         s3.d=s1.d=s2.d;  
45.         for (int j=0;j<k-1;j++)  
46.         {  
47.             s3.z[j]*=s2.x;  
48.         }  
49.         s3.z[k-1]=s2.y;  
50.         k++;  
51.     }  
52.     return s3;  
53. }  
54. void main()  
55. {  
56.     int a[n]={788,460,286,257,5689};  
57. //  int a[n]={340,510,1700,189,3388};  
58.     /*srand( (unsigned)time( NULL ) );// 
59.     int a[n]; 
60.     for (int i=0;i<n;i++) 
61.     { 
62.         a[i]=rand()%10000; 
63.     }*/  
64.     struct t s=extended_eucild1(a);  
65.     cout<<s.d<<" ";  
66.     for (int i=0;i<n;i++)  
67.     {  
68.         cout<<s.z[i]<<" ";  
69.     }  
70.     cout<<endl;  
71. }  

根据(3)知： 每次while循环都要调用一次extended_eucild函数，而调用一次函数需要经过O(1+lg(min{gcd(ai,aj),ak})) ，再根据31.2-5的更紧确的界,实际只需经过O(1+lg(min{gcd(ai,aj),ak}/gcd(gcd(ai,aj),ak))). 例如第一次调用extended_eucild函数时，需要O(1+lg(a0/gcd(a0,a1)))时间,(这里设ai是非负的).所以经过n-1次while循环.可以求出这n个数的最大公约数。

31.2-8说明如何使用(具有两个自变量的)gcd函数作为子程序才能高效计算出lcm(a1,a2,...an)(最小公倍数)

代码如下：

[cpp] view plaincopyprint?

1. //求N个数的最小公倍数  
2. #include <iostream>  
3. #include <time.h>  
4. using namespace std;  
5. #define n 4//根据输入确定元素个数。  
6. struct t  
7. {  
8.     int d;  
9.     int x;  
10.     int y;  
11.     int z[n];//存放题目中的x0,x1...xn  
12. }s;  
13. struct t extended_eucild(int a,int b)  
14. {  
15.     if (b==0)  
16.     {  
17.         s.d=a;  
18.         s.x=1;  
19.         s.y=0;  
20.         return s;  
21.     }  
22.     else  
23.     {         
24.         struct t ss=extended_eucild(b,a%b);  
25.         s.d=ss.d;  
26.         s.x=ss.y;  
27.         s.y=ss.x-(a/b)*ss.y;  
28.         return s;  
29.     }  
30. }  
31. struct t extended_eucild1(int a[])  
32. {  
33.     struct t s3,s1,s2;  
34.     for (int i=0;i<n;i++)  
35.     {  
36.         s3.z[i]=0;  
37.     }  
38.     s1=extended_eucild(a[0],a[1]);//O(lgb)  
39.     s3.z[0]=s1.x;  
40.     s3.z[1]=s1.y;  
41.     int k=3;  
42.     s1.d=a[0]*(a[1]/s1.d);  
43.     while (k<n+1)  
44.     {     
45.         s2=extended_eucild(s1.d,a[k-1]);  
46.         s3.d=s2.d;  
47.         s1.d=s1.d*(a[k-1]/s2.d);  
48.         for (int j=0;j<k-1;j++)  
49.         {  
50.             s3.z[j]*=s2.x;  
51.         }  
52.         s3.z[k-1]=s2.y;  
53.         k++;  
54.     }  
55.     s3.d=s1.d;  
56.     return s3;  
57. }  
58. void main()  
59. {  
60.     int a[n]={756,4400,19845,9000};  
61.     //int a[n]={12,30,50};  
62.     struct t s=extended_eucild1(a);  
63.     cout<<s.d<<" ";  
64.     for (int i=0;i<n;i++)  
65.     {  
66.         cout<<s.z[i]<<" ";  
67.     }  
68.     cout<<endl;  
69. }  

31.2-9证明：n1,n2,n3,n4是两两互质的当且仅当gcd(n1n2,n3n4)=gcd(n1n3,n2n4)=1,更一般地，证明：n1,n2,...nk,当仅当从ni中导出的lgk对整数互为质数。

证明：

(1)n1,n2,n3,n4是两两互质=>gcd(n1n2,n3n4)=gcd(n1n3,n2n4)=1

因为n1,n2,n3,n4是两两互质,所以gcd(n1,n3)=gcd(n2,n3)=1,根据定理31.6 gcd(n1n2,n3)=1.同理gcd(n1n2,n4)=1,再次利用31.6得：gcd(n1n2,n3n4)=1,同理可证：gcd(n1n3,n2n4)=1

(2)因为gcd(n1n2,n3n4)=gcd(n1n3,n2n4)=1=>n1,n2,n3,n4是两两互质,存在n1(n2x)+n2(n3y)=1.所以可以看出这是一个n1,n2的线性组合，gcd(n1,n2)是这个线性组合的最小正整数，那么1刚好是最小正整数，所以gcd(n1,n2)=1,类似的由定理31.2可导出gcd(n1,n3)=gcd(n2,n3)=gcd(n3,n4)=1，所以n1,n2,n3,n4是两两互质。得证！

第二个证明不太会

31.3模运算

群的定义：群(S,(+))是一个集合S和定义在S上的二进制运算(+)

群的性质：封闭性，单位元，结合律，逆元

子群：如果(S,(+))是一个群，S'‘⊆S,并且(S',(+))称为(S,(+))的子群。

31.3-1画出群(Z4,+4)和群(Z5*,X5)的运算表。通过找这两个群的元素间的一一对应关系α。满足a+b ≡c(mod 4)当仅当α(a)Xα(b)≡α(c)(mod 5),来证明这两个群是同构的。

+401230012311230223013301

2

x5123411234224133314244321

可以证明满足a+b ≡c(mod 4)当仅当α(a)Xα(b)≡α(c)(mod 5).考察两个群中的每一个元素以及对应的二元运算结果。a(行)与b(列)代表群(Z4,+4)中的元素，c代表a与b的二元运算结果。并且经过关系α ，有α(a)(行)与α(b)(列)是相同位置下的群(Z5*,X5)的元素，α(c)是α(a)与α(b)经过群(Z5*,X5)二元运算结果。举3个例子来说明：当a=0,b=0时.c=0 则α(a)=1,α(b)=1, α(c)=1 => 0+0≡0mod4 当仅当1*1≡1mod5

           当a=3,b=2时,c=1.则α(a)=4,α(b)=3, α(c)=2.=>3+2≡1mod4, 当仅当4*3≡2mod5

           当a=2,b=1,c=3时 。则α(a)=3,α(b)=2, α(c)=1.=>2+1≡3mod4, 当仅当3*2≡1mod5.

类似地，运算表中所有数据均可以得出证明结论。但是关系α是一一对应的。这点我还不太懂。当然如果证明关系α是一一对应的,那么自然就证明了他们是同构的。

31.3-2列举出Z9和Z13*的所有子群。

Z9各元素={0,1,2,3,4,5,6,7,8} 各元素生成的子群<0>={0}. <1>={0,1,2,3,4,5,6,7,8}. <1>=<2>=<4>=<5>=<7>=<8>. <3>=<6>={0,3,6}

Z13*各元素={1,2,3,4,5,6,7,8,9,10,11,12} 各元素生成的子群 <1>={1}. <2>=<6>=<7>=<11>={1,2,3,4,5,6,7,8,9,10,11,12}.<3>=<9>={1,3,9}.<4>={1,3,4,9,10,12},<5>=<8>={1,5,8,12},<10>={1,3,4,9,10,12} <12>={1,12}

31.3-3证明定理31.14

定理31.14(一个有限群的非空封闭子集是一个子群) 如果(S,(+))是一个有限群，S'是S的任意一个非空子集并满足对所有a,b∈S'，有a(+)b∈S‘，则(S',(+))是(S,(+))的一个子群。

证明：因为S’是S的一个非空子集，所以有S‘⊆S。由所有a,b∈S'，有a(+)b∈S‘，由群的定义知：(S’,(+))也是一个群。又因为：(S,(+))是一个群,则根据子群定义知：(S’,(+))是(S,(+))一个子群。

31.3-4证明：如果p是素数且e是正整数，则φ(p^e)=(p^(e-1))(p-1)

对于整数p^e的素因数只有一个p,所以φ(p^e)=(p^e)π(1-1/p)=(p^e)(1-1/p)=(p^(e-1))(p-1)

31.3-5 证明：对任意n>1和任意a∈Zn*,由式f(x)=ax mod n 所定义的函数f: Zn*->Zn*是Zn*的一个置换。

对于每一个a和n，都有唯一的一个x对应唯一的一个y，也就是f(x)是一个双射函数。所以函数f: Zn*->Zn*是Zn*的一个置换。

31.4 求解模线性方程

定理31.20 对任意正整数a和n，如果d=gcd(a,n),则在Zn中，<a>=<d>={0,d,2d,....((n/d)-1)d},因此，|<a>|=n/d;

推论31.21 当且仅当d|b,方程ax≡b(mod n)对于未知量x有解，这里d=gcd(a,n).

推论31.22 方程ax≡b(mod n)或者对模n有d个不同的解，或者无解，这里d=gcd(a,n).

定理31.23 令d=gcd(a,n).假设对某些整数x'和y'，有 d=ax'+ny',如果d|b，则方程ax≡b(mod n)有一个解的值为x0,这里 x0=x'(b/d)mod n

定理31.24 假设方程ax≡b(mod n)有解(d|b,d=gcd(a,n)),且x0是该方程任意一个解。因此，该方程对模n恰有n个不同解，分别为xi=x0+i(n/d)

代码如下：

[cpp] view plaincopyprint?

1. <span style="color:#000000;">#include <iostream>  
2. using namespace std;  
3. struct t  
4. {  
5.     int d,x,y,z;  
6. }s;  
7. struct t extended_eucild(int a,int b)  
8. {  
9.     if (b==0)  
10.     {  
11.         s.d=a;  
12.         s.x=1;  
13.         s.y=0;  
14.         return s;  
15.     }  
16.     else  
17.     {  
18.           
19.         struct t ss=extended_eucild(b,a%b);  
20.         s.d=ss.d;  
21.         s.x=ss.y;  
22.         s.y=ss.x-(a/b)*ss.y;  
23.         return s;  
24.     }  
25. }  
26. void MODULAR_LINEAR_EQUATION_SOLVER(int a,int b,int n)  
27. {  
28.    struct t ss=extended_eucild(a,n);  
29.    if (b%ss.d==0)  
30.    {  
31.        int x=(ss.x*(b/ss.d))%(n);//100k-105  
32.        if (x<0)  
33.        {  
34.            x+=n;  
35.        }  
36.        for (int i=0;i<ss.d;i++)  
37.        {  
38.            cout<<(x+i*(n/ss.d))%n<<" ";  
39.        }  
40.    }  
41.    else  
42.    {  
43.        cout<<"no solutions"<<endl;  
44.    }  
45. }  
46. void main()  
47. {  
48.      MODULAR_LINEAR_EQUATION_SOLVER(25,100,50);  
49. }</span>  

31.4-1 找出方程35x≡10(mod50)

带入上面代码后得：

31.4-2 证明：只要gcd(a,n)=1,方程ax≡ay(mod n) 就意味着x≡y(mod n).通过一个过gcd(a,n)>1情况下的反例，证明条件gcd(a,n)=1是必要的。

由ax≡ay(mod n) 得ax=ay+nk => a(x-y)=nk  => n|a(x-y)且gcd(a,n)=1,根据推论31.5 得n|(x-y)  =>x-y=nk => x=y+nk => x≡y(mod n),若gcd(a,n)>1,举个例子，设a=12,n=3 gcd(a,n)=3  由a(x-y)=nk  设k=4 ,x-y=1,x=6,y=5 但是 x-y≠nk 所以 x≡y(mod n).不成立。

31.4-3 考察下列对过程MODULAR-LINEAR-EQUATION-SOLVER的第3行修改：3 x0=x'(b/d)mod (n/d) 能正确运行吗？ 解释能或不能的原因。

能，因为x'(b/d)mod (n/d) 与x'(b/d)mod n 同余。

具体来说一下为什么以上两者同余。因为x0≡x'(b/d)mod (n/d)，设x0'≡x'(b/d)mod n 则n|(x0'-x'(b/d)) =>x0'-x'(b/d)=nk=>d|n(n=dk') => x0'-x'(b/d)=dk'k=>k'|x0'-x'(b/d)即(n/d)|x0'-x'(b/d) 根据同余定义知：x0'≡x'(b/d)(mod n/d),所以x0=x0'。也就是x0≡x'(b/d)mod n => x0≡x'(b/d)mod (n/d),x0如果是x'(b/d)mod n 的余数，那么也就是x'(b/d)mod (n/d)的余数.x0=x'(b/d)mod (n/d) 与x0'=x'(b/d)mod n 所得余数x0与x0'结果一样。

由此引出一条普遍又常用的定理：如果a≡b(mod m)， 若m1|m，则 a≡b(mod m1)

31.4-4 略。

31.5中国余数定理

定理31.27(中国余数定理) 令n=n1n2...nk，其中因子ni两两互质.考虑一下对应关系：a<->(a1,a2,....,ak) ..(1)这里a∈Zn，ai∈Zn,而且对i=1，2,....k,ai= a mod ni 因此，映射(1)是一个在Zn与笛卡尔积Zn1XZn2X...XZnk之间的一一对应没(双射)。通过在合适的系统中对每个坐标位置独立地执行操作。对Zn中元素所执行的运算可以等价地作用于对应的可k元祖，也就是说，如果a<->(a1,a2,....,ak) ,b<->(b1,b2,....,bk)  那么 (a+b)mod n<->((a1+b1)mod n1,.....(ak+bk)mod nk) 类似地a-b,ab也有类似对应关系。

推论31.28如果n1,n2...nk两两互质，且n=n1n2...nk,则对任意整数a1,a2...ak,关于未知量x的联立方程组 x≡ai(mod ni),i=1,2,...k 对模n有唯一解。

推论31.29 如果n1,n2...nk两两互质，且n=n1n2...nk,则对所有整数x和a， x≡a(mod ni)(其中i=1,2...k)当仅当  x≡a(mod n)

代码如下：

[cpp] view plaincopyprint?

1. <span style="color:#000000;">//孙子定理(中国剩余定理)  
2. #include <iostream>  
3. using namespace std;  
4. #define num 3//方程个数  
5. #define LEN 10  
6. struct t  
7. {  
8.     int d,x,y,z;  
9. }s;  
10. struct t extended_eucild(int a,int b)  
11. {  
12.     if (b==0)  
13.     {  
14.         s.d=a;  
15.         s.x=1;  
16.         s.y=0;  
17.         return s;  
18.     }  
19.     else  
20.     {  
21.           
22.         struct t ss=extended_eucild(b,a%b);  
23.         s.d=ss.d;  
24.         s.x=ss.y;  
25.         s.y=ss.x-(a/b)*ss.y;  
26.         return s;  
27.     }  
28. }  
29. int  MODULAR_LINEAR_EQUATION_SOLVER(int a,int b,int n)  
30. {  
31.     struct t ss=extended_eucild(a,n);  
32.     int *c=new int[LEN];  
33.     if (b%ss.d==0)  
34.     {  
35.         int x=(ss.x*(b/ss.d))%(n);//100k-105  
36.         if (x<0)  
37.         {  
38.             x+=n;  
39.         }  
40.         return x;  
41.     }  
42.     else  
43.     {  
44.         cout<<"no solutions"<<endl;  
45.         return -1;  
46.     }  
47. }  
48. void Chinese_remainder_theorem(int m[],int a[])  
49. {  
50.     int s=1;  
51.     for (int i=0;i<num;i++)  
52.     {  
53.       s*=m[i];  
54.     }  
55.     int M[num]={0},MM[num]={0};  
56.     for (int j=0;j<num;j++)  
57.     {  
58.        M[j]=s/m[j];  
59.     }  
60.     for (int k=0;k<num;k++)  
61.     {  
62.        MM[k]=MODULAR_LINEAR_EQUATION_SOLVER(M[k],1,m[k]);  
63.     }  
64.     cout<<"x≡";  
65.     int sum=0;  
66.     i=0;  
67.     while (i!=num)  
68.     {  
69.         sum+=MM[i]*M[i]*a[i];  
70.         i++;  
71.     }  
72.     cout<<sum%s<<"(mod"<<s<<")"<<endl;  
73. }  
74. void main()  
75. {  
76.     //num值根据方程个数调整  
77.     //int a[num]={4,5};  
78.     //int m[num]={5,11};  
79.     int a[num]={1,2,3};//输入余数  
80.     int m[num]={9,8,7};//输入除数  
81.     Chinese_remainder_theorem(m,a);  
82. }</span>  

31.5-1 找出所有解，使方程x≡4(mod 5) x≡5(mod 11)同时成立。

根据以上面代码，num调整为2，这样可得结果为x≡49(mod 55)

31.5-2 找出被9,8,7除时，余数分别为1,2,3的所有整数x.

根据以上面代码，num调整为3，这样可得结果为x≡10(mod 504)

31.5-3 略。不太懂a<->(a1,a2,....,ak) 这种式子的含义。

31.5-4 在定理31.27的定义下，证明：对于任意的多项式f,方程f(x)≡0(mod n)....(1)的根的个数等于 f(x)≡0(mod n1),f(x)≡0(mod n2)...f(x)≡0(mod nk)...(2)每个方程根的个数的积。

证明：(1)的解一定是(2)的解：由于f(x)≡0(mod n) 所以n|f(x),又ni|n=>ni|f(x)=>f(x)≡0(mod ni)

           (2)的解一定是(1)的解： f(x)≡0(mod ni)  => ni|f(x) => 又因为ni两两互质，n1n2..nk|f(x) n|f(x)=>f(x)≡0(mod n)

所以(1)与(2)等价。

现在设 ： f(x)≡0(mod ni)的解为 x=biti (ti=1,2...Ti) 一共有Ti个解

从这k个方程中，每个方程任意取一个解组成了一个方程组 x≡b1t1(mod m1),x≡b2t2(mod m2),.....x≡bktk(mod mk),根据中国剩余定理知：有唯一解，

x≡∑M1M1'biti(mod m)...(3) 方程组每个方程分别有T1,T2...Tk个解，所以有T1*T2...*Tk个排列选取方法，也就是这么多个根，又根据上面的(1)与(2)等价.所以(1)的根个数就是(2)的根个数乘积。

31.6元素的幂

欧拉定理：对于任意整数n>1,a^φ(n)≡1(mod n)对所有a∈Zn*都成立。

费马定理：如果p是素数，则a^(p-1))≡1(mod p)对所有a∈Zp*都成立.

定理31.32 对所有的素数p>2和所有正整数e,使得Zn*是循环群的n>1的值为2,4，p^e 2p^e.

如果ordn(g)=|Zn*|,则对模n,Zn*中的每个元素都是g的一个幂，且g是Zn*的一个原根或生成元。

如果g是原根且a是Zn*中的任意元素，则存在一个z,使得g^z≡a(mod n),这个z称为对模n到基g上的a的一个离散对数或指数，这个值用ind(a)

定理31.34 如果p是一个奇素数且e≥1，则方程x^2≡1(mod p^e)仅有两个解，即x=1和x=-1，这两个解就是平凡平方根，如果x≠这两个根(1,-1)，那么则x是一个以n为模的1的非平凡平方根。

推论31.35 如果对模n存在1的非平凡平方根，则n是合数。

31.6-1 画出一张表，展示Z11*中每个元素的阶。找出最小的原根g,并画出一张表，对所有x∈Z11*,给出相应的ind(x)的值。

元素=
1
2
3
4
5
6
7
8
9
10
i次=1
1
2
3
4
5
6
7
8
9
10
2
1
4
9
5
3
3
5
9
4
1
3
1
8
5
9
4
7
2
6
3
10
4
1
5
4
3
9
9
3
4
5
1
5
1
10
1
1
1
10
10
10
1
10
6
1
9
3
4
5
5
4
3
9
1
7
1
7
9
5
3
8
6
2
4
10
8
1
3
5
9
4
4
9
5
3
1
9
1
6
4
3
9
2
8
7
5
10
10
1
1
1
1
1
1
1
1
1
1
阶ord
1
10
5
5
5
10
10
10
5
2

最小原根g=2

x=
1
2
3
4
5
6
7
8
9
10
ind
10
1
8
2
4
9
7
3
6
5

31.6-2写出一个模取幂算法，要求该算法检查b的各位的顺序为从右向左，而非从左向右

[cpp] view plaincopyprint?

1. #include <iostream>  
2. using namespace std;  
3. #define len 10  
4. //int c[len]={0};  
5. int* BIT(int bb,int b[])  
6. {  
7.    int i=0;     
8.    while (i!=len&&b!=0)  
9.    {  
10.     b[i]=bb%2;  
11.     bb=bb/2;  
12.     i++;  
13.    }  
14.    return b;  
15. }  
16. //从右向左检查b的各位顺序  
17. int MODULAR_EXPONENTIATION(int a,int b[],int n)  
18. {  
19.  int c=0;  
20.  int d=a;  
21.  int t=1;  
22.  for (int i=0;i<len;i++)  
23.  {  
24.   if (b[i]==1)  
25.   {  
26.    c++;  
27.    t=(t*d)%n;  
28.   }  
29.   c=2*c;  
30.   d=(d*d)%n;  
31.  }  
32.  return t;  
33. }  
34. void main()  
35. {  
36.  int *b=new int[len];  
37.  cout<<MODULAR_EXPONENTIATION(100,BIT(560,b),560)<<endl;  
38.    
39. }  

31.6-3 假设已知φ(n),说明如何运用过程MODULAR-EXPONENTIATION,对任意a∈Zn*,计算出a^(-1) mod n的值。

首先证明求乘法逆元的数学原理：根据欧拉定理 :a^φ(n)≡1(mod n)  根据乘法逆元定义：gcd(a,n)=1,n>1 ax≡1(mod n)(x为乘法逆元=a^(-1)(mod n)) 由于a^φ(n)与1同余，ax与1同余，所以根据同余的传递性，a^φ(n)与ax同余，∴a^φ(n)≡ax(mod n) ∴ a^(φ(n)-1)≡x(mod n) 即为∴ x≡a^(φ(n)-1)(mod n),利用MODULAR-EXPONENTIATION函数求a^(φ(n)-1)(mod n)即可得到乘法逆元。

[cpp] view plaincopyprint?

1. //找乘法逆元  
2. #include <iostream>  
3. #include <math.h>  
4. using namespace std;  
5. #define len 10  
6. #define nn 120  
7. struct t  
8. {  
9.  int d,x,y,z;  
10. }s;  
11. struct t extended_eucild(int a,int b)  
12. {  
13.  if (b==0)  
14.  {  
15.   s.d=a;  
16.   s.x=1;  
17.   s.y=0;  
18.   return s;  
19.  }  
20.  else  
21.  {  
22.     
23.   struct t ss=extended_eucild(b,a%b);  
24.   s.d=ss.d;  
25.   s.x=ss.y;  
26.   s.y=ss.x-(a/b)*ss.y;  
27.   return s;  
28.  }  
29. }  
30. int* BIT(int bb,int b[])  
31. {  
32.  int i=0;     
33.  while (i!=len&&b!=0)  
34.  {  
35.   b[i]=bb%2;  
36.   bb=bb/2;  
37.   i++;  
38.  }  
39.  return b;  
40. }  
41. //从左向右检查b的各位顺序  
42. int MODULAR_EXPONENTIATION(int a,int b[],int n)  
43. {  
44.  int c=0;  
45.  int d=1;  
46.  for (int i=len-1;i>=0;i--)  
47.  {  
48.   c=2*c;  
49.   d=(d*d)%n;  
50.   if (b[i]==1)  
51.   {  
52.    c++;  
53.    d=(d*a)%n;  
54.   }  
55.  }  
56.  return d;  
57. }  
58. //计算F(n) 欧拉常数  
59. int f(int n)  
60. {  
61.  double temp=n;   
62.  for (int i=2;i<=n;i++)  
63.  {  
64.   int flag=0;  
65.   if (n%i==0)  
66.   {  
67.    for (int j=2;j<=sqrt(i);j++)  
68.    {  
69.     if (i%j==0)  
70.     {  
71.      flag=1;  
72.     }  
73.    }  
74.      
75.    if (flag==0)  
76.    {  
77.     double k=i;  
78.     temp*=((k-1.0)/k);  
79.    }  
80.   }  
81.  }  
82.  return temp;  
83. }  
84. //求Zn*中所有元素a的集合 也是对于模n的乘法群  
85. int ZnX(int a[],int n)  
86. {  
87.  int t=0;  
88.  for (int i=1;i<n;i++)  
89.  {  
90.   if(extended_eucild(i,n).d==1)  
91.   {  
92.    a[t++]=i;  
93.    cout<<i<<" ";  
94.   }  
95.  }  
96.    return t;  
97. }  
98. void main()  
99. {  
100.  int *b=new int[len];  
101.  int a[nn]={0};  
102.  cout<<"a的所有元素如下："<<endl;  
103.  int t=ZnX(a,nn);  
104.  cout<<endl;  
105.  cout<<"对应的乘法逆元如下："<<endl;  
106.  for (int i=0;i<t;i++)  
107.  {  
108.   cout<<MODULAR_EXPONENTIATION(a[i],BIT(f(nn)-1,b),nn)<<" ";  
109.  }  
110. }  

31.7 RSA公钥加密系统

在RSA公钥加密系统中，一个参与者按下列过程创建他的公钥和密钥。

1.随机选取两个大素数p和q，使得p≠q，例如，素数p和q可能各有1024位。

2.计算n=pq

3.选取一个与φ(n)互质的小奇数e,其中由等式(31.20),φ(n)=(p-1)(q-1)

4.对模φ(n)，计算出e的乘法逆元d的值。

5.将对P=(e,n)公开，并作为参与者的RSA公钥。

6.使对S=(d,n)保密，并作为参与者的RSA密钥。

31.7-1 考虑一个RSA密钥集合，其中p=11,q=29,n=319,e=3.在密钥中用到的d值应当是多少？对消息M=100加密后得到什么消息？

φ(n)=(p-1)(q-1)=280 3d≡1(mod 280) =>d=187

M=100时，密文C=P(M)=M^e mod n=100^3 mod 319=254
31.7-2和31.7-3以我现有知识无法解决，略去.

31.8素数的测试

素数的密度：定理31.37（素数定理） lim(π(n)/(n/lnn))=1

我可以用试除法测试素数，但是经过我实际测试发现，效率没有MILLER-RABIN高。

如果n是一个合数，且a^(n-1)≡1(mod n) 则称n是一个基为a的伪素数。

Carmichael数：就是用简单的利用费马定理测试素数过程中会产生合数，这个合数是符合费马定理的，但是多数情况，用费马定理测试素数还是有效的。因为这种数极少，1亿里面只有255个。下面就是检测Carmichael数过程以及简单利用费马定理求素数的代码：

[cpp] view plaincopyprint?

1. //查找carmichael数:这个代码利用朴素的试除法与简单的费马定理作对比，两者结果相反，则是carmichael数  
2. #include <iostream>  
3. #include <math.h>  
4. using namespace std;  
5. #define len 15  
6. #define COMPOSITE 0  
7. #define PRIME 1  
8. int* BIT(int bb,int b[])  
9. {  
10.     int i=0;     
11.     while (i!=len&&b!=0)  
12.     {  
13.         b[i]=bb%2;  
14.         bb=bb/2;  
15.         i++;  
16.     }  
17.     return b;  
18. }  
19. //从左向右检查b的各位顺序  
20. int MODULAR_EXPONENTIATION(int a,int b[],int n)  
21. {  
22.     int c=0;  
23.     int d=1;  
24.     for (int i=len-1;i>=0;i--)  
25.     {  
26.         c=2*c;  
27.         d=(d*d)%n;  
28.         if (b[i]==1)  
29.         {  
30.             c++;  
31.             d=(d*a)%n;  
32.         }  
33.     }  
34.     return d;  
35. }  
36. bool PSEUDOPRIME(int n,int b[])  
37. {  
38.     if (MODULAR_EXPONENTIATION(2,BIT(n-1,b),n)!=1)  
39.     {  
40.         return COMPOSITE;  
41.     }  
42.     else return PRIME;  
43. }  
44. int Primer(int n)  
45. {  
46.     int Flag=0;  
47.     for (int i=2;i<=sqrt(n);i++)  
48.     {  
49.         if (n%i==0)  
50.         {  
51.             Flag=1;  
52.             cout<<i<<" ";  
53.             break;  
54.         }  
55.     }  
56.     if (Flag==0)  
57.     {  
58.         return 1;  
59.     }  
60.     else return 0;  
61. }  
62. void main()  
63. {  
64.     int *b=new int[len];  
65.     for (int i=2;i<10000;i++)  
66.     {  
67.         if (PSEUDOPRIME(i,b)==1&&Primer(i)==0)  
68.         {  
69.             cout<<i<<" "<<endl;  
70.         }  
71.     }  
72.     cout<<endl;  
73. }  

MILLE-RABIN素数测试：选取多个基值a,测试素数，虽然也可能出错，但是出错几率大大降低，出错原因也不是依赖待测数本身以及坏的输入，而是选取基值a抽签运气。但是这里我们不用随机函数进行选取a,而是利用从网上看到的在一定数范围内选取固定的几个a，虽然我不懂为何选取这几个值作为基值a,但是经过我大量数据测试，这种选取特定a值对素数测试结果还是比较可信的。这是我经过3-4天的努力，终于解决了大数据的素数测试代码，唯一有瑕疵的地方就是使用了除法，但是是较小整数的除法，不会产生溢出的，所以还是可行的。

[cpp] view plaincopyprint?

1. #include <iostream>  
2. using namespace std;  
3. #define len 80  
4. #define COMPOSITE 0  
5. #define PRIME 1  
6. unsigned __int64 bit_Multiplication(__int64 a, __int64 b,unsigned __int64 m)//防止溢出__int64  
7. {//思想是将a与b尽量转化成较小的数再求余数.转化过程中遇＞m的情况，则对m求余  
8. // 比如求(56*121)mod 200=(56*(4*30+1))mod 200=((56*4-200)*30+56*1)mod 200=(24*(9*3+3)+56*1)mod 200=((24*9-200)*3+3*24+56)mod 200=(48+73+56) mod200   
9.     __int64 rem=0;//初始化余数  
10.     while (b)  
11.     {  
12.     __int64 t=a;//将乘数a用临时变量t保存起来  
13.         if ((a<<=1)>=m)//如果在移位时＞m了，  
14.         {  
15.             a-=(a/m)*m;//那么只保留余数  
16.         }  
17.         if (b&1)//如果乘数b的二进制位上是1，则执行if语句  
18.         {  
19.             rem+=t;//对分解a与b过程中产生的剩余部分求和  
20.             if (rem>=m)//如果这个和＞m  
21.             {  
22.                 rem-=(rem/m)*m;//那么对剩余部分求余  
23.             }  
24.         }  
25.         b>>=1;  
26.     }  
27.      return rem;  
28. }  
29. int* BIT(unsigned __int64 bb, int b[])  
30. {  
31.     int i = 0;  
32.     while (i != len&&bb != 0)  
33.     {  
34.         b[i] = bb % 2;  
35.         bb>>=1;  
36.         i++;  
37.     }  
38.     return b;  
39. }  
40. //从左向右检查b的各位顺序  
41. unsigned __int64 MODULAR_EXPONENTIATION(unsigned __int64 a, int b[], unsigned __int64 n)  
42. {  
43.     int c = 0;  
44.     unsigned __int64 d = 1;  
45.     for (int i = len - 1; i >= 0; i--)  
46.     {  
47.         d=bit_Multiplication(d,d,n);  
48.         if (d>n/2)  
49.         {  
50.             d = n-d;//利用同余原理，取较小的余数  
51.         }  
52.         if (b[i] == 1)  
53.         {  
54.             d=bit_Multiplication(d,a,n);  
55.             if (d>n/2)  
56.             {  
57.                 d = n-d;//利用同余原理，取较小的余数  
58.             }  
59.         }  
60.     }  
61.     return d;  
62. }  
63. bool WITNESS(unsigned __int64 a, int b[], unsigned __int64 n)//由于n为偶数时比为合数，所以我们忽略偶数情况，只对奇数进行判断！  
64. {  
65.     unsigned __int64 temp = n - 1, s = 0;  
66.     int t=0;  
67.     while (!bit_Multiplication(temp,1,2))  
68.     {  
69.         temp /= 2;  
70.         t++;  
71.     }  
72.     unsigned __int64 u = temp;  
73.     __int64 x = MODULAR_EXPONENTIATION(a, BIT(u, b), n);  
74.     __int64 y=1;  
75.     for (int i = 0; i<t; i++)  
76.     {  
77.         if (x==-1)  
78.         {  
79.             y=1;  
80.         }  
81.         else if (x<0&&x!=-1)  
82.         {  
83.             y=bit_Multiplication(-x,-x,n);  
84.         }  
85.         else  
86.         {  
87.             y=bit_Multiplication(x,x,n);  
88.         }  
89.         if (y>n / 2)  
90.         {  
91.             y -= n;//利用同余原理，取较小的余数  
92.         }  
93.         if (y == 1 && x != 1 &&x!=-1&& x != n - 1)  
94.         {  
95.             return true;  
96.         }  
97.         x=y;  
98.     }  
99.     if (y != 1)  
100.     {  
101.         return true;  
102.     }  
103.     return false;  
104. }  
105. int MILLER_RABIN(unsigned __int64 n, unsigned __int64 s, int b[])  
106. {  
107.     unsigned __int64 t = bit_Multiplication(n,1,2);  
108.     if (n==2||n==3)  
109.     {  
110.         return PRIME;  
111.     }  
112.     if ((!t))  
113.     {  
114.         return COMPOSITE;  
115.     }  
116.         if( n < 1373653 ) //下面这一组if-else是网上找到选取基值a方法   
117.         {    
118.             if( WITNESS(2, b, n)     
119.                 || WITNESS(3, b, n)  )    
120.                 return COMPOSITE;    
121.         }    
122.         else if( n < 9080191 )    
123.         {    
124.             if( WITNESS(31, b, n)      
125.                 || WITNESS(73, b, n)  )    
126.                 return COMPOSITE;   
127.         }      
128.         else if( n < 4759123141 )   
129.     {   
130.         if( WITNESS(2, b, n)   
131.                 ||WITNESS(3, b, n)   
132.                 ||WITNESS(5, b, n)   
133.                 ||WITNESS(11, b, n) )   
134.             return COMPOSITE;   
135.     }   
136.         else if( n < 2152302898747 )   
137.         {   
138.             if( WITNESS(7, b, n)  
139.                 ||WITNESS(3, b, n)  
140.                 || WITNESS(5, b, n)  
141.                 || WITNESS(2, b, n)   
142.                 || WITNESS(11, b, n) )   
143.                 return COMPOSITE;   
144.         }   
145.         else   
146.         {   
147.             if( WITNESS(7, b, n)   
148.                 || WITNESS(3, b, n)  
149.                 || WITNESS(5, b, n)  
150.                 || WITNESS(2, b, n)  
151.                 || WITNESS(11, b, n)  
152.                 || WITNESS(31, b, n)  
153.                 || WITNESS(61, b, n)  
154.                 || WITNESS(73, b, n) )   
155.                 return COMPOSITE;   
156.     }   
157.      return PRIME;  
158. }  
159. int main()  
160. {  
161.     int j = 0;  
162.   
163.     for (double i = 11111111900; i<11111112000; i++)  
164.     //for (double i=2;i<100;i++)  
165.     {  
166.         int *b = new int[len];  
167.         if (MILLER_RABIN(i, 1, b))  
168.         {  
169.             j++;  
170.             cout << i << " ";  
171.         }  
172.     }  
173.     cout << "共" << j << endl;  
174.     return 0;  
175. }  

31.8-1 证明：如果一个奇整数n>1不是素数或素数的幂，则存在一个以n为模的1的非平凡平方根。

设n=(p1^e1)*(p2^e2)*...(pi^ei) 其中p1,p2...pn≥3的两两不同的素数。由于奇整数n不是素数或者素数的幂，则n必然是奇合数，至少存在两个不同的奇素因子，i≥2，且p1≠p2。

由于p1,p2...pn≥3的两两不同的素数,所以p1^e1,p2^e2,.....pi^ei是两两互质。 故根据孙子定理：求x² ≡1(mod n)的解<=>即为求方程组x²≡1(mod pi^ei)的解 (i=1,2...k)，根据定理31.34 对于奇素数pi,ei≥1方程x²≡1(mod pi^ei)的解仅有2个解，又根据书上题目31.5-4的结论知：由于存在至少两个素因子p1,p2,所以至少存在两个方程f(x)=x²-1≡0(mod p1^e1)和f(x)=x²-1≡0(mod p2^e2)他们的解的乘积=2X2=4即为原方程f(x)=x²-1≡0(mod n)的解，解数≥4，除去平凡根±1外，应该还有至少2个非平凡平方根。

31.8-2 可以把欧拉定理稍微加强为如下形式：对所有a∈Zn*,a^λ(n)≡1(mod n) 其中n=(p1^e1)*(p2^e2).....(pr^er),且λ(n)定义为λ(n)=lcm(φ(p1^e1),....,φ(pr^er)).证明：λ(n)|φ(n).如果 λ(n)|n-1,则合数n为Carmichael数。最小的Carmichael数为561=3*11*17；这里，λ(n)=lcm(2,10,16)=80,它可以整除560.证明Carmichael数必须既是“无平方数”(不能被任何素数的平方整除)，又是至少三个素数的积。(因此，Carmichael数不是很常见)。

(1)由φ(pr^er))=pr^(er-1)(pr-1) 得φ(n)=n*(1-1/p1)(1-1/p2)...(1-1/pr)=(p1^e1)*(p2^e2).....(pr^er)(1-1/p1)(1-1/p2)...(1-1/pr)=φ(p1^e1)*φ(p2^e2)*.....φ(pr^er)，又因为λ(n)=lcm(φ(p1^e1),....,φ(pr^er)),根据公倍数一定是最小公倍数的倍数知：λ(n)|φ(n).

(2)若λ(n)|n-1,则可以写成n-1=λ(n)k,k∈Z,对所有a∈Zn*,a^λ(n)≡1(mod n) =>根据同余是可以相乘的：a^(λ(n)k)≡1^k(mod n) =>a^(n-1)≡1(mod n)...(1) 由于n是合数且满足(1)式，则n为Carmichael数

(3) Carmichael数对于每个a∈(1,n)都满足 a^(n-1)≡1(mod n) 则设p是n的一个素因子，并且它的最高次幂是p^(e+1)且e>0  令a=p^e∈(1,n) 则满足p^(ne)≡p^e(mod n)=>n|p^(ne)-p^e又因为p^(e+1)|n,则根据整除的传递性：p^(e+1)|p^(ne)-p^e =>p^(ne)-p^e=p^(e+1)k => p^(ne-e)-1=pk=>p|p^(ne-e)-1,由p|p^(ne-e)得 p|1,这样推出矛盾，所以只有e=0才可能成立。那么Carmichael数的每个素因子的幂肯定只有≤1次。

(4)λ(n)是最小公倍数，所以φ(pi^ei)|λ(n) 由(2)知：是Carmichael数的条件是λ(n)|n-1,那么φ(pi^ei)|n-1,所以pi-1|n-1.n为合数，假设它只有两个素因子n=p1p2(p1≠p2)=> p1-1|p1p2-1 => p1-1|p2(p1-1)+p2-1 =>p1-1|p2-1=>同理 p2-1|p1-1 那么p1=p2 和假设矛盾，所以至少有3个作为Carmichael数的素因子。

31.8-3 证明：如果x是以n为模的1的非平凡平方根，则gcd(x-1,n) 和gcd(x+1,n)都是n的非平凡约数。

x是以n为模的1的非平凡平方根=> n|x²-1.=> n|(x-1)(x+1)  假设gcd(x-1,n)=1,则n|(x+1) => x≡(-1)(mod n) =>则x是-1的平凡平方根，与已知矛盾，同理可证gcd(x+1,n)>1。所以原题得证！

31.9 整数的因子分解

Pollard的rho启发式方法代码如下：

[cpp] view plaincopyprint?

1. //POLLARD_RHO大整数因子分解，对于19位以内的整数有效。  
2. #include <iostream>  
3. #include <stdlib.h>  
4. using namespace std;  
5. #define len 80  
6. #define COMPOSITE 0  
7. #define PRIME 1  
8. unsigned __int64 bit_Multiplication(__int64 a, __int64 b,unsigned __int64 m)//防止溢出__int64  
9. {//思想是将a与b尽量转化成较小的数再求余数.转化过程中遇＞m的情况，则对m求余  
10.     // 比如求(56*121)mod 200=(56*(4*30+1))mod 200=((56*4-200)*30+56*1)mod 200=(24*(9*3+3)+56*1)mod 200=((24*9-200)*3+3*24+56)mod 200=(48+73+56) mod200   
11.     __int64 rem=0;//初始化余数  
12.     while (b)  
13.     {  
14.         __int64 t=a;//将乘数a用临时变量t保存起来  
15.         if ((a<<=1)>=m)//如果在移位时＞m了，  
16.         {  
17.             a-=(a/m)*m;//那么只保留余数  
18.         }  
19.         if (b&1)//如果乘数b的二进制位上是1，则执行if语句  
20.         {  
21.             rem+=t;//对分解a与b过程中产生的剩余部分求和  
22.             if (rem>=m)//如果这个和＞m  
23.             {  
24.                 rem-=(rem/m)*m;//那么对剩余部分求余  
25.             }  
26.         }  
27.         b>>=1;  
28.     }  
29.     return rem;  
30. }  
31. __int64 Euclid(__int64 a,__int64 b)  
32. {  
33.     if (b==0)  
34.     {  
35.         return a;  
36.     }  
37.     else  
38.     {  
39.         return Euclid(b,a%b);  
40.     }  
41. }  
42. int* BIT(unsigned __int64 bb, int b[])  
43. {  
44.     int i = 0;  
45.     while (i != len&&bb != 0)  
46.     {  
47.         b[i] = bb % 2;  
48.         bb>>=1;  
49.         i++;  
50.     }  
51.     return b;  
52. }  
53. //从左向右检查b的各位顺序  
54. unsigned __int64 MODULAR_EXPONENTIATION(unsigned __int64 a, int b[], unsigned __int64 n)  
55. {  
56.     int c = 0;  
57.     unsigned __int64 d = 1;  
58.     for (int i = len - 1; i >= 0; i--)  
59.     {  
60.         d=bit_Multiplication(d,d,n);  
61.         if (d>n/2)  
62.         {  
63.             d = n-d;//利用同余原理，取较小的余数  
64.         }  
65.         if (b[i] == 1)  
66.         {  
67.             d=bit_Multiplication(d,a,n);  
68.             if (d>n/2)  
69.             {  
70.                 d = n-d;//利用同余原理，取较小的余数  
71.             }  
72.         }  
73.     }  
74.     return d;  
75. }  
76. bool WITNESS(unsigned __int64 a, int b[], unsigned __int64 n)//由于n为偶数时比为合数，所以我们忽略偶数情况，只对奇数进行判断！  
77. {  
78.     unsigned __int64 temp = n - 1, s = 0;  
79.     int t=0;  
80.     while (!bit_Multiplication(temp,1,2))  
81.     {  
82.         temp /= 2;  
83.         t++;  
84.     }  
85.     unsigned __int64 u = temp;  
86.     __int64 x = MODULAR_EXPONENTIATION(a, BIT(u, b), n);  
87.     __int64 y=1;  
88.     for (int i = 0; i<t; i++)  
89.     {  
90.         if (x==-1)  
91.         {  
92.             y=1;  
93.         }  
94.         else if (x<0&&x!=-1)  
95.         {  
96.             y=bit_Multiplication(-x,-x,n);  
97.         }  
98.         else  
99.         {  
100.             y=bit_Multiplication(x,x,n);  
101.         }  
102.         if (y>n / 2)  
103.         {  
104.             y -= n;//利用同余原理，取较小的余数  
105.         }  
106.         if (y == 1 && x != 1 &&x!=-1&& x != n - 1)  
107.         {  
108.             return true;  
109.         }  
110.         x=y;  
111.     }  
112.     if (y != 1)  
113.     {  
114.         return true;  
115.     }  
116.     return false;  
117. }  
118. int MILLER_RABIN(unsigned __int64 n, unsigned __int64 s, int b[])  
119. {  
120.     unsigned __int64 t = bit_Multiplication(n,1,2);  
121.     if (n==2||n==3)  
122.     {  
123.         return PRIME;  
124.     }  
125.     if ((!t))  
126.     {  
127.         return COMPOSITE;  
128.     }  
129.         if( n < 1373653 )    
130.         {    
131.             if( WITNESS(2, b, n)     
132.                 || WITNESS(3, b, n)  )    
133.                 return COMPOSITE;    
134.         }    
135.         else if( n < 9080191 )    
136.         {    
137.             if( WITNESS(31, b, n)      
138.                 || WITNESS(73, b, n)  )    
139.                 return COMPOSITE;   
140.         }      
141.         else if( n < 4759123141 )   
142.     {   
143.         if( WITNESS(2, b, n)   
144.                 ||WITNESS(3, b, n)   
145.                 ||WITNESS(5, b, n)   
146.                 ||WITNESS(11, b, n) )   
147.             return COMPOSITE;   
148.     }   
149.         else if( n < 2152302898747 )   
150.         {   
151.             if( WITNESS(7, b, n)  
152.                 ||WITNESS(3, b, n)  
153.                 || WITNESS(5, b, n)  
154.                 || WITNESS(2, b, n)   
155.                 || WITNESS(11, b, n) )   
156.                 return COMPOSITE;   
157.         }   
158.         else   
159.         {   
160.             if( WITNESS(7, b, n)   
161.                 || WITNESS(3, b, n)  
162.                 || WITNESS(5, b, n)  
163.                 || WITNESS(2, b, n)  
164.                 || WITNESS(11, b, n)  
165.                 || WITNESS(31, b, n)  
166.                 || WITNESS(61, b, n)  
167.                 || WITNESS(73, b, n) )   
168.                 return COMPOSITE;   
169.     }   
170.      return PRIME;  
171. }  
172. unsigned __int64 POLLARD_RHO(unsigned __int64 n,int b[])//分解成某个n的因子，该因子可能不是最小的。  
173. {  
174.    int i=1;  
175.    if (MILLER_RABIN(n,1,b))return NULL;  
176.    __int64 x;  
177.    if (n<32767) x=rand()%n;  
178.    else if(n<1073676289) x=rand()*rand();  
179.    else if (n<35181150961663)x=rand()*rand()*rand();  
180.    else if (n<1152780773560811521)x=rand()*rand()*rand()*rand();  
181.    else return NULL;  
182.    x=bit_Multiplication(x,1,n);//对于随机超过待测试数n的x，我们对其求余，保证它随机数在(1,n)之间的数  
183.    __int64 y=x;  
184.    __int64 k=2;  
185.    while (1)  
186.    {  
187.        i++;  
188.        x=bit_Multiplication(x-1,x+1,n);  
189.        __int64 d=Euclid(y-x,n);  
190.        d=d>0?d:-d;  
191.        if (d!=1&&d!=n)  
192.        {  
193.            return d;  
194.        }  
195.        if (i==k)  
196.        {  
197.            y=x;  
198.            k<<=2;  
199.        }  
200.    }  
201. }  
202. void main()  
203. {  
204.     int *b = new int[len];  
205.     __int64 t=POLLARD_RHO(1152780773560811517,b);  
206. }  

31.9-1 在图31-7(a)所示的执行过程中，过程POLLARD-RHO在何时输出1387的因子73

在x=84,y=814 时，输出73.

31.9-2 假设给定函数f:Zn->Zn和一个初值x0∈Zn。定义xi=f(x(i-1)),i=1,2....令t和u>0是满足x(t+1)=x(t+u+i)(i=0,1...)的最小值。在Pollard的rho算法的术语中，t为rho的尾的长度，u是rho的回路长度。试写出一个计算t和u的值的有效算法，并分析其运行时间。

由于上面已给出完整的Pollard启发式方法的代码，下面仅给出计算t和u的方法：

[cpp] view plaincopyprint?

1. <span style="color:#000000;">unsigned __int64 POLLARD_RHO(unsigned __int64 n,int b[])//分解成某个n的因子，该因子可能不是最小的。  
2. {  
3.    int i=1,flag=0,u=0,t=0;  
4.    List *x=NULL,*head=NULL;//由于需要用数组来保存每一个x以便查找重复的数据来确定rho回路大小，但是我们不知道该数组大小，所以我才用动态数组。  
5.    if (MILLER_RABIN(n,1,b))return NULL;  
6.    x=new List[LEN];  
7.    head=x;  
8.    x->num=i;  
9.    if (n<32767) x->key=rand()%n;//这个ifelse结构用来随机(1,n)范围内的数据初始化x  
10.    else if(n<1073676289) x->key=rand()*rand();  
11.    else if (n<35181150961663)x->key=rand()*rand()*rand();  
12.    else if (n<1152780773560811521)x->key=rand()*rand()*rand()*rand();  
13.    else return NULL;  
14.    x->key=bit_Multiplication(x->key,1,n);//对于随机超过待测试数n的x，我们对其求余，保证它随机数在(1,n)之间的数  
15.    __int64 y=x->key;  
16.    __int64 k=2;  
17.    while (1)  
18.    {  
19.        i++;  
20.        int z=x->key;  
21.        x->next=new List[LEN];  
22.        x=x->next;  
23.        x->next=NULL;  
24.        x->key=bit_Multiplication(z-1,z+1,n);  
25.        x->num=i;  
26.        __int64 d=Euclid(y-x->key,n);  
27.        d=d>0?d:-d;  
28.        if (d!=1&&d!=n&&flag<1)  
29.        {  
30.            printf("n的因子=%d\n",d);  
31.            flag++;  
32.        }  
33.        struct List*yy=head;  
34.        while (yy->num<x->num)  
35.        {  
36.            if (yy->key==x->key)  
37.            {  
38.                u=(x->num)-(yy->num);  
39.                t=yy->num;  
40.                cout<<"rho回路长度u="<<u<<endl;  
41.                cout<<"rho尾长度t="<<t<<endl;  
42.                return d;  
43.            }  
44.            yy=yy->next;  
45.        }  
46.        if (i==k)  
47.        {  
48.            y=x->key;  
49.            k=2*k;  
50.        }  
51.    }  
52. }</span>  

31.9-3为了发现形如p^e的数(其中p是素数，e>1)的一个因子，POLLARD-RHO要执行多少步？ θ(√p)步

31.9-4 略。

33-1(二进制的gcd算法) 与计算余数的执行速度相比，大多数计算机执行减法运算，测试一个二进制整数的奇偶性运算以及折半运算的执行速度都要更快些。本题所讨论的二进制gcd算法中避免了欧几里得算法中对余数的计算过程。

a.证明：如果a和b都是偶数，则gcd(a,b)=2gcd(a/2,b/2).

设gcd(a,b)=d1,gcd(a/2,b/2)=d2.由a与b为偶数，2|a,2|b,2|gcd(a,b)=2|d1，由d1|a,d1|b,所以(d1/2)|(a/2),(d2/2)|(b/2) 所以d1|gcd(a/2,b/2) => (d1/2)|d2 => d1/2≤d2....(1) d2|(a/2)且d2|(b/2)=>2d2|a且2d2|b=>2d2|gcd(a,b) => 2d2|d1 => 2d2≤d1...(2) 由(1)与(2)得：d1=2d2.得证！

b.证明：如果a是奇数，b是偶数，则gcd(a,b)=gcd(a,b/2).

设gcd(a,b)=d1,gcd(a,b/2)=d2. 由d1|a且d1|b => gcd(d1,2)=1(原因：若2|d1=>2|a=>因为a是奇数,所以推出矛盾,则gcd(2,d1)=1),d1|(b/2)=>d1|gcd(a,b/2)=>d1|d2=> d1≤d2...(1) 由d2|a,d2|(b/2)=> b/2=d2*k=>b=d2*(2k)=>d2|b =>d2|gcd(a,b) =>d2|d1 =>d2≤d1...(2) 由(1)与(2)知：d1=d2,得证！

c.证明：如果a和b都是奇数，则gcd(a,b)=gcd((a-b)/2,b).

设gcd(a,b)=d2 ,gcd((a-b)/2,b)=d1 则d1|(a-b)/2 d1|b =>d1|((a-b)/2)*2+b =>d1|a =>d1|gcd(a,b) =>d1|d2 =>d1≤d2...(1)  由d2|a且d2|b =>d2|a-b 因为gcd(d2,2)=1(原因：若2|d2 => 2|b =>因为b是奇数=>推出矛盾=>gcd(d2,2)=1) 所以d2|(a-b)/2 =>d2|gcd((a-b)/2,b) =>d2|d1 => d2≤d1...(2) 由(1)与(2)知： d1=d2.得证！

d.设计一个有效的二进制算法，输入整数为a和b(a≥b)，并且算法的运行时间为O(lg a).假定每个减法运算，测试奇偶性运算以及折半运算都能在单位时间行。

代码如下：

[cpp] view plaincopyprint?

1. //二进制的gcd算法  
2. #include <iostream>  
3. using namespace std;  
4. int Abs(int x)  
5. {  
6.     return x<0?-x:x;  
7. }  
8. int Binary_gcd(int a,int b)  
9. {  
10.     if (a==0)  
11.     {  
12.         return b;  
13.     }  
14.     else if(!(a&1)&&!(b&1))  
15.     {  
16.         return 2*Binary_gcd(a>>1,b>>1);  
17.     }  
18.     else if ((a&1)&&(b&1))  
19.     {  
20.         if (a<b)  
21.         {  
22.             a = a^b;    
23.             b = b^a;    
24.             a = a^b;   
25.         }  
26.         return Binary_gcd(Abs((a-b)>>1),b);  
27.     }  
28.     else if (!(a&1)&&(b&1))  
29.     {  
30.         return Binary_gcd(a>>1,b);  
31.     }  
32.     else   
33.     {  
34.         return Binary_gcd(a,b>>1);  
35.     }  
36. }  
37. void main()  
38. {  
39.     cout<<Binary_gcd(2*3*7*13*5,17*3*14*2*5)<<endl;  
40.     cout<<Binary_gcd(17,15)<<endl;  
41. }  

31-2(对欧几里得算法中位操作的分析)

a.考虑用普通的“纸和笔”算法来实现长除法的运算，用a除以b，得到商q和余数r。证明，这种算法需要执行O((1+lgq)lgb)次位操作。

以下代码完全按照用“纸和笔”手写计算除法的方式计算二进制除法：

[cpp] view plaincopyprint?

1. <span style="color:#000000;">//我不清楚如何达到O((1+lgq)lgb，但是除去log函数可能用到的位操作，这里只使用了O(lgqlgb)次位操作  
2. #include <iostream>  
3. #include<math.h>  
4. using namespace std;  
5. __int64 Dev(__int64 a, __int64 b)  
6. {  
7.     __int64 s = 1, i = 0, ans = 0,j=1;  
8.     __int64 alength = log(a) / log(2);  
9.     while (alength>i )//O(lgqlgb)  
10.     {//内存循环执行次数与外层循环执行次数乘积就是移位操作符执行的总次数  
11.         while (s < b&&alength>i++)//O(lgb)  
12.         {  
13.             s <<= 1;//位操作  
14.             if (a&j << alength - i )s++;  
15.             ans <<= 1;//位操作  
16.         }  
17.         if (s >= b)//O(lgq)  
18.         {//if语句执行的次数完全取决于商的二进制1的个数，假设商每位都是1，那么就达到了if语句执行次数最大值，也就是lgq  
19.             ans++;  
20.             s -= b;                                                              
21.         }  
22.     }  
23.     printf("商=%I64d\n",ans);  
24.     return s;  
25. }  
26. void main()  
27. {  
28.       
29.     printf("余数=%I64d", Dev(2561, 147));  
30. }</span>  

c.证明：EUCILD(a,b)通常需要执行O(μ(a,b))次位操作；当其输入为两个β位数时，需要执行的位操作次数为O(β²).

由31.2-5结论知：EUCILD函数执行了(1+lgb)次递归调用，由于a%b取余操作执行了lga+1次位运算，所以总的位运算次数为O((1+lga)(1+lgb))=O(μ(a,b)) 当输入两个β位时，具有β位的数a与b，它们的位数就是lga=β与lgb=β，所以运算次数为O((1+β)²) 去掉低次数项就得到结论。

其中的求余数的位操作代码如下：

[cpp] view plaincopyprint?

1. <span style="color:#000000;">#include <iostream>  
2. using namespace std;  
3. int Dev(int a,int b)  
4. {  
5.    int ans=0,j=1,temp=a;  
6.    while ((temp>>=1))j++;//计算a的位数为lga  
7.    for (int i=j;i>=0;i--)//O(lga+1)  
8.    {  
9.        if ((a>>i)>=b)  
10.        {  
11.            ans+=(1<<i);  
12.            a-=(b<<i);  
13.        }  
14.    }  
15.    cout<<"余数="<<a<<endl;  
16.    return ans;  
17. }  
18. void main()  
19. {  
20.    cout<<"商="<<Dev(1157923,4713)<<endl;  
21. }</span>  

b.定义μ(a,b)=(1+lga)(1+lgb).证明：过程EUCLID在把计算gcd(a,b)的问题转化为计算gcd(b,a mod b)的问题时，所执行的位操作次数至多为c(μ(a,b)-μ(b,a mod b),其中c>0,为某一个足够大的常数。

根据c的结论很容易便知：EUCILD(a,b)通常需要执行O(μ(a,b))次位操作，同理EUCILD(b,a mod b)通常需要执行O(μ(b,a mod b))次位操作，则从EUCILD(a,b)->EUCILD(b,a mod b)过程 就是两者之差：O(μ(a,b))-O(μ(b,a mod b))=O(μ(a,b)-μ(b,a mod b))=c(μ(a,b)-μ(b,a mod b)
31-3(关于斐波那契数的三个算法) 在已知n的情况下，本题对计算第n个斐波那契数Fn的三种算法的效率进行了比较。

假定两个数的加法，减法和乘法的代价都是O(1),与数的大小无关。

a.证明：基于递归式(3.22)计算Fn的直接递归方法的运行时间为n的幂。

由下图可知，求Fn就是在求一颗递归树，这颗树的所有叶子的总和就是待求Fn的值.叶子数量一共有2^n,显然她是n的幂

                f(10)
               /        \
            f(9)         f(8)
          /     \       /    \
       f(8)     f(7)  f(7)   f(6)
      /   \     /   \ 
   f(7)  f(6)  f(6) f(5)

b.试说明如何运用记忆法在O(n)时间内计算Fn.

说白了，就是倒着求，预求Fn,先求F0与F1，然后保存记忆下F0与F1的值，再求F2=F1+F0的值，以此类推直到求得Fn，所以可以做一个循环

代码如下：

[cpp] view plaincopyprint?

1. #include <iostream>  
2. using namespace std;  
3. #define n 10  
4. double Fibonacci()  
5. {  
6.     double  F[n+1]={0,1};  
7.     for (int i=0;i<n-1;i++)  
8.     {  
9.         F[i+2]=F[i]+F[i+1];  
10.     }  
11.     return F[n];  
12. }  
13. void main()  
14. {  
15.     cout<<Fibonacci()<<endl;  
16. }  

c.试说明如何仅用整数加法和乘法运算，就可以在O(lgn)的时间内计算Fn.考虑2X2矩阵[0,1,1,1]。

Fn为[0,1,1,1]^n的第一行第一列，这个可用数学归纳法证明，下面是代码：

[cpp] view plaincopyprint?

1. #include<iostream>  
2. #include <assert.h>   
3. using namespace std;  
4. struct Matrix2X2  
5. {  
6.     Matrix2X2//2X2矩阵  
7.         (int F00 , int F01 , int F10 , int F11 ):F_00(F00), F_01(F01), F_10(F10), F_11(F11){}     
8.     int F_00;  
9.     int F_01;  
10.     int F_10;  
11.     int F_11;  
12. };  
13. Matrix2X2 MatrixMultiply(const Matrix2X2& matrix1, const Matrix2X2& matrix2)  
14. {//两个矩阵乘积  
15.     return Matrix2X2(  
16.         matrix1.F_00 * matrix2.F_00 + matrix1.F_01 * matrix2.F_10,  
17.         matrix1.F_00 * matrix2.F_01 + matrix1.F_01 * matrix2.F_11,  
18.         matrix1.F_10 * matrix2.F_00 + matrix1.F_11 * matrix2.F_10,  
19.         matrix1.F_10 * matrix2.F_01 + matrix1.F_11 * matrix2.F_11);  
20. }  
21. // n阶矩阵  
22. // 0  1  
23. // 1  1  
24. Matrix2X2 Fibonacci(unsigned int n)  
25. {  
26.     assert(n > 0);     
27.     Matrix2X2 matrix(0,1,1,1);  
28.     if(n==1)  
29.     {  
30.         matrix = Matrix2X2(0, 1, 1, 1);  
31.     }  
32.     else if(!(n&1))  
33.     {  
34.         matrix = Fibonacci(n>>1);  
35.         matrix = MatrixMultiply(matrix, matrix);  
36.     }  
37.     else if(n&1)  
38.     {  
39.         matrix = Fibonacci((n-1)>>1);  
40.         matrix = MatrixMultiply(matrix, matrix);  
41.         matrix = MatrixMultiply(matrix, Matrix2X2(0, 1, 1, 1));  
42.     }  
43.     return matrix;  
44. }  
45. void main()  
46. {  
47.     for (int i=0;i<=30;i++)  
48.     {  
49.         cout<<Fibonacci(i+1).F_00<<endl;  
50.     }  
51. }  

d.现在假设对两个β位数相加需要θ(β)时间，对两个β位数相乘需要θ(β²)时间。如果这样更合理地估计基本算数运算代价，这三种方法运行时间又是多少？

第一种所有叶子由0与1组成，那么相加需要θ(1*2^n),而lgn=β,所以θ(1*2^(2^β))=θ(2^(2^β))

第二种由于lgn=β，所以θ(2^β)
第三种T(n)=T(n/2)+θ(β²) 其中lg(n/2)=β 利用主方法.T(n)=(β²)

31-4(二次余数)设p是一个奇素数。如果关于未知量x的方程x²=a(mod p)有解，则数a∈Zp*就是一个二次余数。

a.证明：对模p,恰有(p-1)/2个二次余数。

考虑模p的绝对最小剩余系-(p-1)/2,-(p-1)/2+1,...-1,1,...(p-1)/2-1,(p-1)/2 d是模p的二次剩余当仅当d≡(-(p-1)/2)²,(-(p-1)/2+1)²,...(-1)²,1²...

((p-1)/2-1)²或((p-1)/2)²(mod p). 由于(-j)²≡j²(mod p),所以d是模p的二次剩余当且仅当d≡1²,...,((p-1)/2-1)²或((p-1)/2)²(mod p).1≤i<j≤p-1)/2

时，i²≠j²(mod p),所以d的二次剩余共有(p-1)/2个

b.如果p是素数，对a∈Zp*,定义勒让德符号(a/p),若a是对模p的二次余数，则它等于1，；否则它等于-1.证明：如果a∈Zp*,则(a/p)=a^((p-1)/2)(mod p)给出一个有效的算法，使其能确定一个给定的数a是否是对模p的二次余数。分析所给算法的效率。

(1)若a是对模p的二次余数,则x²=a(mod p)有解,且gcd(a,p)=1 =>p|x²-a 考虑x^p-x=x^p-(a^((p-1)/2))x+(a^((p-1)/2))x-x=x(x^(p-1)-a^((p-1)/2))+(a^(p-1)/2-1)x=x(x²-a)f(x)+(a^((p-1)/2)-1)x 又因为对数素数p 有x^p≡x(mod p) 即p|x^p-x  所以p|(a^((p-1)/2)-1)x 用反证法可以证明：x²≠0(mod p)=>同样继续用反证法可以证x≠0(mod p)=>gcd(x,p)=1 所以p|(a^((p-1)/2)-1) =>(a^((p-1)/2)≡1(mod p) 也可以写成1≡(a^((p-1)/2)(mod p)勒让德符号是二次剩余则(a/p)=1。所以(a/p)≡(a^((p-1)/2)(mod p)

(2)若a是对模p的二次非余数，则考虑对于素数p, a^(p-1)≡1(mod p)=> p|a^(p-1)-1 =>p-1是偶数，所以可以进行因式分解=〉p|(a^((p-1)/2)-1)(a^((p-1)/2)+1)，由于是二次非剩余，所以gcd(p,a^((p-1)/2)-1)=1 所以p|(a^((p-1)/2)+1),即-1≡(a^((p-1)/2)+1)(mod p) 勒让德符号无解表示为(a/p)=-1 所以(a/p)≡(a^((p-1)/2)(mod p)

由(1)与(2)知(a/p)≡(a^((p-1)/2)(mod p)

书中介绍的快速模取幂的方法是高效的计算余数的算法并且能计算大整数的模运算，其运行时间也在书中提及，具体在31.6 元素的幂中有介绍。

关于算法可以用模取幂的那个程序求模取幂的余数，判断余数与±1是否同余，若与1同余代表有二次剩余，否则就是非二次剩余。

以下是代码：

[cpp] view plaincopyprint?

1. //用反复平方法求数的幂 也称为模取幂  
2. #include <iostream>  
3. using namespace std;  
4. #define len 10  
5. #define p 7  
6. #define aa 34  
7. int* BIT(int bb,int b[])  
8. {  
9.    int i=0;     
10.    while (i!=len&&b!=0)  
11.    {  
12.        b[i]=bb%2;  
13.        bb=bb/2;  
14.        i++;  
15.    }  
16.    return b;  
17. }  
18. //从左向右检查b的各位顺序  
19. int MODULAR_EXPONENTIATION(int a,int b[],int n)  
20. {  
21.   int c=0;  
22.   int d=1;  
23.   for (int i=len-1;i>=0;i--)  
24.   {  
25.       c=2*c;  
26.       d=(d*d)%n;  
27.       if (b[i]==1)  
28.       {  
29.           c++;  
30.           d=(d*a)%n;  
31.       }  
32.   }  
33.   return d;  
34. }  
35. void main()  
36. {  
37.     int *b=new int[len];  
38.     //判断a^((p-1)/2)(mod p)是否为1，若为1，则有二次余数。  
39.     cout<<MODULAR_EXPONENTIATION(aa,BIT((p-1)/2,b),p)<<endl;//结果是-1，说明是二次非剩余  
40.     cout<<MODULAR_EXPONENTIATION(254,BIT((p-1)/2,b),p)<<endl;//结果是1，说明是二次剩余   
41. }  

c.证明：如果p是形如4k+3的素数，且a是Zp*中一个二次余数，则a^(k+1)mod p 是对模p的a的平方根。找出一个以p为模的二次余数a的平方根需要多长时间？

因为a是Zp*中一个二次余数,所以有(a^((p-1)/2)≡1(mod p) 将p=4k+3带入 a^(2k+1)≡1(mod p) => a^(2k+2)≡a(mod p) =>(a^(k+1))^2≡a(mod p)有解=>所以a^(k+1)mod p 是模p对a的平方根。用快速模取幂的方法算出a^(k+1)mod p,需要进行算术运算的总次数是O(β)，并且需要的位操作的总次数是O(β^3).(若a,(p-1)/2与p都是β位数)
d.试描述一个有效的随机算法，找出一个以任意素数p为模的非二次余数，也就是指Zp*中不是二次余数的成员。所给出的算法平均需要执行多少次算术运算？

a=random(1,p-1) 在(1,p-1)内循环随机选取一个数作为a的值来用快速模取幂的方法测试是否是二次非剩余直到找到一个非剩余为止，利用(a^((p-1)/2)≡(-1)(mod p)这个同余式，只要余数为-1(mod p),那么就是非剩余。这里做算术运算的主要过程就是快速模取幂函数的时间以及随机选取a值的时间乘积，书中31.6节已有介绍模取幂的运行时间。现在来看下选中非剩余a的运行时间，因为由题意知，非剩余有(p-1)/2个，剩余有(p-1)/2个,而选中的概率有1/2，平均选2个a值就能找到一个非剩余a,所以循环随机选a的次数是一个常数，所以只有快速模取幂才是程序主要运行时间，平均要执行的算术运算和模取幂的一样。

感觉31.1-12和31.1-13应该用FFT算法解决。