演算法導論第三十一(31)章數論演算法

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是否是非平凡冪。

以下是程式碼：

//最樸素的多項式時間內判斷一個數是否為某個數的冪的形式：就是用列舉法，挨個找,但是這個是關於n的多項式，關於β的多項式暫時沒有想出。
#include <iostream>
#include <math.h>
using namespace std;
void main()
{
    int n=64,flag=1;
	//O(√nlgn)
    for (int i=2;i<=sqrt(n);i++)//O(a=√n)
	{
		int m=n,k=0;
		while (m%i==0)//O(k=lgn)
		{
			m=m/i;
			k++;
		}
		if (m==1&&k>1)
		{
			cout<<n<<"是非平凡冪,存在一個整數"<<i<<"它的"<<k<<"次冪="<<n<<endl;
			flag=0;
			break;
		}
	}
	if (flag)
	{
		cout<<"n="<<n<<"不存在非平凡冪"<<endl;
	}
}

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演算法解決。)

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

//位運算的乘法與除法
#include <iostream>
using namespace std;
//位運算的乘法
int bit_Multiplication(int a,int b)
{
     int ans=0;
	 for (int i=1;i;i<<=1,a<<=1)
	 {
		 if (b&i)
		 {
			 ans+=a;
		 }
	 }
	 return ans;
}
//位運算的除法
int bit_Division1(int x,int y)
{
    int ans=0;
	for (int i=31;i>=0;i--)
	{
		if ((x>>i)>=y)
		{
			ans+=(1<<i);
			x-=(y<<i);
		}
	}
	return ans;
}
//計算整數的二進位制位數
int bit_num(int d)
{
   int i=0;
   while (d)
   {
	   d>>=1;
	   i++;
   }
   return i;
}
//位運算的除法 計算商
int bit_Division2_quotient(int x,int y)
{
	int c2=bit_num(x),c1=bit_num(y),quotient=0;
	for (int i=c2-c1;i>=0;i--)//i=c2-c1防止除數y移位後超過無符號整數最大值 時間複雜度O(c2-c1)
	{
		unsigned int a=(y<<i);//有了i=c2-c1保證了y<<i不會溢位 a有c1+c2-c1=c2位
		if (a<=x)
		{
			quotient+=(1<<i);
			x-=a;
		}
	}
	//總的時間複雜度為 O(c2)=O(x的二進位制位數)=O(b^2) b為除數的十進位制位數
	return quotient;
}
//位運算的除法 計算餘數 與計算商一樣，只是返回值不同
int bit_Division2_Remainder(int x,int y)
{
	int c2=bit_num(x),c1=bit_num(y),quotient=0;
	for (int i=c2-c1;i>=0;i--)//i=c2-c1防止除數y移位後超過無符號整數最大值 時間複雜度O(c2-c1)
	{
		unsigned int a=(y<<i);//有了i=c2-c1保證了y<<i不會溢位 a有c1+c2-c1=c2位
		if (a<=x)
		{
			quotient+=(1<<i);
			x-=a;
		}
	}
	//總的時間複雜度為 O(c2)=O(x的二進位制位數)=O(b^2) b為除數的十進位制位數
	return x;
}
void main()
{
	cout<<bit_Multiplication(350,43)<<endl;
	cout<<bit_Division1(350,43)<<endl;
	cout<<"商："<<bit_Division2_quotient(350,43)<<endl;
	cout<<"餘數："<<bit_Division2_Remainder(350,43)<<endl;
}

31.1-13 寫出一個高效演算法，用於將β位二進位制整數轉化為響應的十進位制表示。證明：如果長度至多為β的整數的乘法或除法運算所需時間為M(β)，則執行二進位制到十進位制轉換所需時間為θ(M(β)lgβ)。(提示：應用分治法，分別使用獨立的遞迴計算結果的前段和後段)(感覺31.1-12和31.1-13應該用FFT演算法解決。)
這次要用位運算+分治思想。我對我寫的這段程式碼持懷疑態度。如果牛人看到實際沒用到分治思想，那麼給出理由和比較好的建議。謝謝！

//用分治思想進行進位制轉換(2->10),寫得不好，湊合看吧
/*#include <iostream>
#include <string>
using namespace std;
#define BIT 6//二進位制整數的位數，可根據所要輸入的二進位制位數設定BIT。
int t=-1;
int Bit_merge(int a[],int p,int r)
{
	static ans=0;
	p=(p>t)?p:(t+1);
	for (int i=p;i<=r;i++)
	{
		t++;
		ans+=a[i]<<i;		
	}
	return ans;
}
int bit_Multiplication(int a[],int p,int r,int flag)//x為二進位制數
{
	int q;
	if (p<r)
	{
	   q=(p+r)/2;
	   bit_Multiplication(a,p,q,0);
	   bit_Multiplication(a,q+1,r,1);
	}
	if (p<=r&&t!=BIT)return flag==0?Bit_merge(a,p,q):Bit_merge(a,q+1,r);
}

void main()
{
	int a[BIT]={0};
	string x;
	cin>>x;
    int j=0;
    while (j!=BIT)
	{
	   a[j]=x[BIT-1-j]-'0';
	   j++;
    }
	cout<<endl;
	cout<<bit_Multiplication(a,0,BIT-1,0)<<endl;
}

如果β位整數的乘法或除法的執行時間為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次。

程式碼如下：

//歐幾里得演算法遞迴形式
#include <iostream>
using namespace std;
int Euclid(int a,int b)
{
	cout<<"gcd("<<a<<","<<b<<")=";
	if (b==0)
	{
		return a;
	}
	else
	{
		return Euclid(b,a%b);
	}
}
void main()
{
   cout<<Euclid(30,21)<<endl;
}

歐幾里得演算法擴充套件形式

程式碼如下：

#include <iostream>
using namespace std;
struct a
{
	int d,x,y;
}s;
struct a extended_eucild(int a,int b)
{
	if(b==0)
	{
        s.d=a,s.x=1,s.y=0;
		return s; 
	}
	else
	{
		struct a ss=extended_eucild(b,a%b);
		s.d=ss.d;
		s.x=ss.y;
		s.y=ss.x-(a/b)*ss.y;
		return s;
	}
}
int Fibonacci(int n)
{
	if (n==1)
	{
		return 1;
	}
	else if (n==0)
	{
		return 0;
	}
	else
	{
		return Fibonacci(n-2)+Fibonacci(n-1);
	}
}
void main()
{
  struct a s=extended_eucild(99,78);
  cout<<s.d<<" "<<s.x<<" "<<s.y<<endl;
}
<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).

a	b	向下取整a/b	d	x	y
899	493	1	29	-6	11
493	406	1	29	5	-6
406	87	4	29	-1	5
87	58	1	29	1	-1
58	29	2	29	0	1
29	0	一	29	1	0

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改寫成迭代形式。

程式碼如下：

//歐幾里得演算法迭代形式
#include <iostream>
using namespace std;
int Euclid(int a,int b)
{
	while (b)
	{
		int c=b;
		b=a%b;
		a=c;
	}
	return a;
}
void main()
{
   cout<<Euclid(89,55)<<endl;
}

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)返回什麼值？證明答案的正確性。

#include <iostream>
using namespace std;
struct a
{
	int d,x,y;
}s;
struct a extended_eucild(int a,int b)
{
	if(b==0)
	{
        s.d=a,s.x=1,s.y=0;
		return s; 
	}
	else
	{
		struct a ss=extended_eucild(b,a%b);
		s.d=ss.d;
		s.x=ss.y;
		s.y=ss.x-(a/b)*ss.y;
		return s;
	}
}
int Fibonacci(int n)
{
	if (n==1)
	{
		return 1;
	}
	else if (n==0)
	{
		return 0;
	}
	else
	{
		return Fibonacci(n-2)+Fibonacci(n-1);
	}
}
int power(int k)
{
	int i=1;
	int t=1;
    while (i!=k+1)
    {
       t*=-1;
	   i++;
    }
	return t;
}
void main()
{
  const int k=10;
  struct a s=extended_eucild(Fibonacci(k+1),Fibonacci(k));
  cout<<s.d<<" "<<s.x<<" "<<s.y<<endl;
  cout<<"d="<<s.d<<" x="<<power(k-1)*Fibonacci(k-2)<<" y="<<power(k)*Fibonacci(k-1)<<endl;
}

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的演算法程式碼如下：

#include <iostream>
#include <time.h>
using namespace std;
#define n 5//根據輸入確定元素個數。
struct t
{
	int d;
	int x;
	int y;
	int z[n];//存放題目中的x0,x1...xn
}s;
struct t extended_eucild(int a,int b)
{
	if (b==0)
	{
		s.d=a;
		s.x=1;
		s.y=0;
		return s;
	}
	else
	{		
		struct t ss=extended_eucild(b,a%b);
		s.d=ss.d;
		s.x=ss.y;
		s.y=ss.x-(a/b)*ss.y;
		return s;
	}
}
struct t extended_eucild1(int a[])
{
	struct t s3,s1,s2;
	for (int i=0;i<n;i++)
	{
		s3.z[i]=0;
	}
	s1=extended_eucild(a[0],a[1]);//O(lgb)
	s3.z[0]=s1.x;
	s3.z[1]=s1.y;
	int k=3;
	while (k<n+1)
	{	
		s2=extended_eucild(s1.d,a[k-1]);//O(lg(min{s1.d,a[k-1]}))........(3)
		s3.d=s1.d=s2.d;
		for (int j=0;j<k-1;j++)
		{
			s3.z[j]*=s2.x;
		}
		s3.z[k-1]=s2.y;
		k++;
	}
	return s3;
}
void main()
{
	int a[n]={788,460,286,257,5689};
//	int a[n]={340,510,1700,189,3388};
	/*srand( (unsigned)time( NULL ) );//
	int a[n];
	for (int i=0;i<n;i++)
	{
		a[i]=rand()%10000;
	}*/
	struct t s=extended_eucild1(a);
	cout<<s.d<<" ";
	for (int i=0;i<n;i++)
	{
		cout<<s.z[i]<<" ";
	}
	cout<<endl;
}

根據(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)(最小公倍數)

程式碼如下：

//求N個數的最小公倍數
#include <iostream>
#include <time.h>
using namespace std;
#define n 4//根據輸入確定元素個數。
struct t
{
	int d;
	int x;
	int y;
	int z[n];//存放題目中的x0,x1...xn
}s;
struct t extended_eucild(int a,int b)
{
	if (b==0)
	{
		s.d=a;
		s.x=1;
		s.y=0;
		return s;
	}
	else
	{		
		struct t ss=extended_eucild(b,a%b);
		s.d=ss.d;
		s.x=ss.y;
		s.y=ss.x-(a/b)*ss.y;
		return s;
	}
}
struct t extended_eucild1(int a[])
{
	struct t s3,s1,s2;
	for (int i=0;i<n;i++)
	{
		s3.z[i]=0;
	}
	s1=extended_eucild(a[0],a[1]);//O(lgb)
	s3.z[0]=s1.x;
	s3.z[1]=s1.y;
	int k=3;
	s1.d=a[0]*(a[1]/s1.d);
	while (k<n+1)
	{	
		s2=extended_eucild(s1.d,a[k-1]);
		s3.d=s2.d;
		s1.d=s1.d*(a[k-1]/s2.d);
		for (int j=0;j<k-1;j++)
		{
			s3.z[j]*=s2.x;
		}
		s3.z[k-1]=s2.y;
		k++;
	}
	s3.d=s1.d;
	return s3;
}
void main()
{
	int a[n]={756,4400,19845,9000};
	//int a[n]={12,30,50};
	struct t s=extended_eucild1(a);
	cout<<s.d<<" ";
	for (int i=0;i<n;i++)
	{
		cout<<s.z[i]<<" ";
	}
	cout<<endl;
}

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),來證明這兩個群是同構的。

+4	0	1	2	3
0	0	1	2	3
1	1	2	3	0
2	2	3	0	1
3	3	0	1	2

x5	1	2	3	4
1	1	2	3	4
2	2	4	1	3
3	3	1	4	2
4	4	3	2	1

可以證明滿足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)

程式碼如下：

<span style="color:#000000;">#include <iostream>
using namespace std;
struct t
{
	int d,x,y,z;
}s;
struct t extended_eucild(int a,int b)
{
	if (b==0)
	{
		s.d=a;
		s.x=1;
		s.y=0;
		return s;
	}
	else
	{
		
		struct t ss=extended_eucild(b,a%b);
		s.d=ss.d;
		s.x=ss.y;
		s.y=ss.x-(a/b)*ss.y;
		return s;
	}
}
void MODULAR_LINEAR_EQUATION_SOLVER(int a,int b,int n)
{
   struct t ss=extended_eucild(a,n);
   if (b%ss.d==0)
   {
	   int x=(ss.x*(b/ss.d))%(n);//100k-105
	   if (x<0)
	   {
		   x+=n;
	   }
	   for (int i=0;i<ss.d;i++)
	   {
		   cout<<(x+i*(n/ss.d))%n<<" ";
	   }
   }
   else
   {
	   cout<<"no solutions"<<endl;
   }
}
void main()
{
     MODULAR_LINEAR_EQUATION_SOLVER(25,100,50);
}</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)

程式碼如下：

<span style="color:#000000;">//孫子定理(中國剩餘定理)
#include <iostream>
using namespace std;
#define num 3//方程個數
#define LEN 10
struct t
{
	int d,x,y,z;
}s;
struct t extended_eucild(int a,int b)
{
	if (b==0)
	{
		s.d=a;
		s.x=1;
		s.y=0;
		return s;
	}
	else
	{
		
		struct t ss=extended_eucild(b,a%b);
		s.d=ss.d;
		s.x=ss.y;
		s.y=ss.x-(a/b)*ss.y;
		return s;
	}
}
int  MODULAR_LINEAR_EQUATION_SOLVER(int a,int b,int n)
{
	struct t ss=extended_eucild(a,n);
	int *c=new int[LEN];
	if (b%ss.d==0)
	{
		int x=(ss.x*(b/ss.d))%(n);//100k-105
		if (x<0)
		{
			x+=n;
		}
		return x;
	}
	else
	{
		cout<<"no solutions"<<endl;
		return -1;
	}
}
void Chinese_remainder_theorem(int m[],int a[])
{
	int s=1;
    for (int i=0;i<num;i++)
    {
      s*=m[i];
    }
	int M[num]={0},MM[num]={0};
	for (int j=0;j<num;j++)
	{
       M[j]=s/m[j];
	}
	for (int k=0;k<num;k++)
	{
       MM[k]=MODULAR_LINEAR_EQUATION_SOLVER(M[k],1,m[k]);
	}
	cout<<"x≡";
	int sum=0;
	i=0;
	while (i!=num)
	{
		sum+=MM[i]*M[i]*a[i];
		i++;
	}
	cout<<sum%s<<"(mod"<<s<<")"<<endl;
}
void main()
{
	//num值根據方程個數調整
	//int a[num]={4,5};
	//int m[num]={5,11};
	int a[num]={1,2,3};//輸入餘數
	int m[num]={9,8,7};//輸入除數
	Chinese_remainder_theorem(m,a);
}</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的各位的順序為從右向左，而非從左向右

#include <iostream>
using namespace std;
#define len 10
//int c[len]={0};
int* BIT(int bb,int b[])
{
   int i=0;   
   while (i!=len&&b!=0)
   {
    b[i]=bb%2;
    bb=bb/2;
    i++;
   }
   return b;
}
//從右向左檢查b的各位順序
int MODULAR_EXPONENTIATION(int a,int b[],int n)
{
 int c=0;
 int d=a;
 int t=1;
 for (int i=0;i<len;i++)
 {
  if (b[i]==1)
  {
   c++;
   t=(t*d)%n;
  }
  c=2*c;
  d=(d*d)%n;
 }
 return t;
}
void main()
{
 int *b=new int[len];
 cout<<MODULAR_EXPONENTIATION(100,BIT(560,b),560)<<endl;
 
}

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)即可得到乘法逆元。

//找乘法逆元
#include <iostream>
#include <math.h>
using namespace std;
#define len 10
#define nn 120
struct t
{
 int d,x,y,z;
}s;
struct t extended_eucild(int a,int b)
{
 if (b==0)
 {
  s.d=a;
  s.x=1;
  s.y=0;
  return s;
 }
 else
 {
  
  struct t ss=extended_eucild(b,a%b);
  s.d=ss.d;
  s.x=ss.y;
  s.y=ss.x-(a/b)*ss.y;
  return s;
 }
}
int* BIT(int bb,int b[])
{
 int i=0;   
 while (i!=len&&b!=0)
 {
  b[i]=bb%2;
  bb=bb/2;
  i++;
 }
 return b;
}
//從左向右檢查b的各位順序
int MODULAR_EXPONENTIATION(int a,int b[],int n)
{
 int c=0;
 int d=1;
 for (int i=len-1;i>=0;i--)
 {
  c=2*c;
  d=(d*d)%n;
  if (b[i]==1)
  {
   c++;
   d=(d*a)%n;
  }
 }
 return d;
}
//計算F(n) 尤拉常數
int f(int n)
{
 double temp=n; 
 for (int i=2;i<=n;i++)
 {
  int flag=0;
  if (n%i==0)
  {
   for (int j=2;j<=sqrt(i);j++)
   {
    if (i%j==0)
    {
     flag=1;
    }
   }
   
   if (flag==0)
   {
    double k=i;
    temp*=((k-1.0)/k);
   }
  }
 }
 return temp;
}
//求Zn*中所有元素a的集合 也是對於模n的乘法群
int ZnX(int a[],int n)
{
 int t=0;
 for (int i=1;i<n;i++)
 {
  if(extended_eucild(i,n).d==1)
  {
   a[t++]=i;
   cout<<i<<" ";
  }
 }
   return t;
}
void main()
{
 int *b=new int[len];
 int a[nn]={0};
 cout<<"a的所有元素如下："<<endl;
 int t=ZnX(a,nn);
 cout<<endl;
 cout<<"對應的乘法逆元如下："<<endl;
 for (int i=0;i<t;i++)
 {
  cout<<MODULAR_EXPONENTIATION(a[i],BIT(f(nn)-1,b),nn)<<" ";
 }
}

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數過程以及簡單利用費馬定理求素數的程式碼：

//查詢carmichael數:這個程式碼利用樸素的試除法與簡單的費馬定理作對比，兩者結果相反，則是carmichael數
#include <iostream>
#include <math.h>
using namespace std;
#define len 15
#define COMPOSITE 0
#define PRIME 1
int* BIT(int bb,int b[])
{
	int i=0;   
	while (i!=len&&b!=0)
	{
		b[i]=bb%2;
		bb=bb/2;
		i++;
	}
	return b;
}
//從左向右檢查b的各位順序
int MODULAR_EXPONENTIATION(int a,int b[],int n)
{
	int c=0;
	int d=1;
	for (int i=len-1;i>=0;i--)
	{
		c=2*c;
		d=(d*d)%n;
		if (b[i]==1)
		{
			c++;
			d=(d*a)%n;
		}
	}
	return d;
}
bool PSEUDOPRIME(int n,int b[])
{
	if (MODULAR_EXPONENTIATION(2,BIT(n-1,b),n)!=1)
	{
		return COMPOSITE;
	}
	else return PRIME;
}
int Primer(int n)
{
	int Flag=0;
	for (int i=2;i<=sqrt(n);i++)
	{
		if (n%i==0)
		{
			Flag=1;
			cout<<i<<" ";
			break;
		}
	}
	if (Flag==0)
	{
		return 1;
	}
	else return 0;
}
void main()
{
	int *b=new int[len];
    for (int i=2;i<10000;i++)
    {
		if (PSEUDOPRIME(i,b)==1&&Primer(i)==0)
		{
			cout<<i<<" "<<endl;
		}
    }
	cout<<endl;
}

MILLE-RABIN素數測試：選取多個基值a,測試素數，雖然也可能出錯，但是出錯機率大大降低，出錯原因也不是依賴待測數本身以及壞的輸入，而是選取基值a抽籤運氣。但是這裡我們不用隨機函式進行選取a,而是利用從網上看到的在一定數範圍內選取固定的幾個a，雖然我不懂為何選取這幾個值作為基值a,但是經過我大量資料測試，這種選取特定a值對素數測試結果還是比較可信的。這是我經過3-4天的努力，終於解決了大資料的素數測試程式碼，唯一有瑕疵的地方就是使用了除法，但是是較小整數的除法，不會產生溢位的，所以還是可行的。

#include <iostream>
using namespace std;
#define len 80
#define COMPOSITE 0
#define PRIME 1
unsigned __int64 bit_Multiplication(__int64 a, __int64 b,unsigned __int64 m)//防止溢位__int64
{//思想是將a與b儘量轉化成較小的數再求餘數.轉化過程中遇＞m的情況，則對m求餘
// 比如求(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 
	__int64 rem=0;//初始化餘數
	while (b)
	{
	__int64	t=a;//將乘數a用臨時變數t儲存起來
		if ((a<<=1)>=m)//如果在移位時＞m了，
		{
			a-=(a/m)*m;//那麼只保留餘數
		}
		if (b&1)//如果乘數b的二進位制位上是1，則執行if語句
		{
			rem+=t;//對分解a與b過程中產生的剩餘部分求和
			if (rem>=m)//如果這個和＞m
			{
				rem-=(rem/m)*m;//那麼對剩餘部分求餘
			}
		}
    	b>>=1;
	}
	 return rem;
}
int* BIT(unsigned __int64 bb, int b[])
{
	int i = 0;
	while (i != len&&bb != 0)
	{
		b[i] = bb % 2;
		bb>>=1;
		i++;
	}
	return b;
}
//從左向右檢查b的各位順序
unsigned __int64 MODULAR_EXPONENTIATION(unsigned __int64 a, int b[], unsigned __int64 n)
{
	int c = 0;
	unsigned __int64 d = 1;
	for (int i = len - 1; i >= 0; i--)
	{
		d=bit_Multiplication(d,d,n);
		if (d>n/2)
		{
			d = n-d;//利用同餘原理，取較小的餘數
		}
		if (b[i] == 1)
		{
			d=bit_Multiplication(d,a,n);
			if (d>n/2)
			{
				d = n-d;//利用同餘原理，取較小的餘數
			}
		}
	}
	return d;
}
bool WITNESS(unsigned __int64 a, int b[], unsigned __int64 n)//由於n為偶數時比為合數，所以我們忽略偶數情況，只對奇數進行判斷！
{
	unsigned __int64 temp = n - 1, s = 0;
	int t=0;
	while (!bit_Multiplication(temp,1,2))
	{
		temp /= 2;
		t++;
	}
	unsigned __int64 u = temp;
	__int64 x = MODULAR_EXPONENTIATION(a, BIT(u, b), n);
    __int64 y=1;
	for (int i = 0; i<t; i++)
	{
		if (x==-1)
		{
			y=1;
		}
		else if (x<0&&x!=-1)
		{
			y=bit_Multiplication(-x,-x,n);
		}
		else
		{
            y=bit_Multiplication(x,x,n);
		}
		if (y>n / 2)
		{
			y -= n;//利用同餘原理，取較小的餘數
		}
		if (y == 1 && x != 1 &&x!=-1&& x != n - 1)
		{
			return true;
		}
		x=y;
	}
	if (y != 1)
	{
		return true;
	}
	return false;
}
int MILLER_RABIN(unsigned __int64 n, unsigned __int64 s, int b[])
{
	unsigned __int64 t = bit_Multiplication(n,1,2);
	if (n==2||n==3)
	{
		return PRIME;
	}
	if ((!t))
	{
		return COMPOSITE;
	}
		if( n < 1373653 ) //下面這一組if-else是網上找到選取基值a方法 
		{  
			if( WITNESS(2, b, n)   
				|| WITNESS(3, b, n)  )  
				return COMPOSITE;  
		}  
		else if( n < 9080191 )  
		{  
			if( WITNESS(31, b, n)    
				|| WITNESS(73, b, n)  )  
				return COMPOSITE; 
		}    
		else if( n < 4759123141 ) 
    { 
        if( WITNESS(2, b, n) 
                ||WITNESS(3, b, n) 
                ||WITNESS(5, b, n) 
                ||WITNESS(11, b, n) ) 
            return COMPOSITE; 
    } 
		else if( n < 2152302898747 ) 
		{ 
			if( WITNESS(7, b, n)
                ||WITNESS(3, b, n)
                || WITNESS(5, b, n)
                || WITNESS(2, b, n) 
                || WITNESS(11, b, n) ) 
				return COMPOSITE; 
		} 
		else 
		{ 
			if( WITNESS(7, b, n) 
                || WITNESS(3, b, n)
                || WITNESS(5, b, n)
                || WITNESS(2, b, n)
                || WITNESS(11, b, n)
                || WITNESS(31, b, n)
                || WITNESS(61, b, n)
                || WITNESS(73, b, n) ) 
				return COMPOSITE; 
    } 
     return PRIME;
}
int main()
{
	int j = 0;

	for (double i = 11111111900; i<11111112000; i++)
	//for (double i=2;i<100;i++)
	{
		int *b = new int[len];
		if (MILLER_RABIN(i, 1, b))
		{
			j++;
			cout << i << " ";
		}
	}
	cout << "共" << j << endl;
	return 0;
}

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啟發式方法程式碼如下：

//POLLARD_RHO大整數因子分解，對於19位以內的整數有效。
#include <iostream>
#include <stdlib.h>
using namespace std;
#define len 80
#define COMPOSITE 0
#define PRIME 1
unsigned __int64 bit_Multiplication(__int64 a, __int64 b,unsigned __int64 m)//防止溢位__int64
{//思想是將a與b儘量轉化成較小的數再求餘數.轉化過程中遇＞m的情況，則對m求餘
	// 比如求(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 
	__int64 rem=0;//初始化餘數
	while (b)
	{
		__int64	t=a;//將乘數a用臨時變數t儲存起來
		if ((a<<=1)>=m)//如果在移位時＞m了，
		{
			a-=(a/m)*m;//那麼只保留餘數
		}
		if (b&1)//如果乘數b的二進位制位上是1，則執行if語句
		{
			rem+=t;//對分解a與b過程中產生的剩餘部分求和
			if (rem>=m)//如果這個和＞m
			{
				rem-=(rem/m)*m;//那麼對剩餘部分求餘
			}
		}
		b>>=1;
	}
	return rem;
}
__int64 Euclid(__int64 a,__int64 b)
{
	if (b==0)
	{
		return a;
	}
	else
	{
		return Euclid(b,a%b);
	}
}
int* BIT(unsigned __int64 bb, int b[])
{
	int i = 0;
	while (i != len&&bb != 0)
	{
		b[i] = bb % 2;
		bb>>=1;
		i++;
	}
	return b;
}
//從左向右檢查b的各位順序
unsigned __int64 MODULAR_EXPONENTIATION(unsigned __int64 a, int b[], unsigned __int64 n)
{
	int c = 0;
	unsigned __int64 d = 1;
	for (int i = len - 1; i >= 0; i--)
	{
		d=bit_Multiplication(d,d,n);
		if (d>n/2)
		{
			d = n-d;//利用同餘原理，取較小的餘數
		}
		if (b[i] == 1)
		{
			d=bit_Multiplication(d,a,n);
			if (d>n/2)
			{
				d = n-d;//利用同餘原理，取較小的餘數
			}
		}
	}
	return d;
}
bool WITNESS(unsigned __int64 a, int b[], unsigned __int64 n)//由於n為偶數時比為合數，所以我們忽略偶數情況，只對奇數進行判斷！
{
	unsigned __int64 temp = n - 1, s = 0;
	int t=0;
	while (!bit_Multiplication(temp,1,2))
	{
		temp /= 2;
		t++;
	}
	unsigned __int64 u = temp;
	__int64 x = MODULAR_EXPONENTIATION(a, BIT(u, b), n);
    __int64 y=1;
	for (int i = 0; i<t; i++)
	{
		if (x==-1)
		{
			y=1;
		}
		else if (x<0&&x!=-1)
		{
			y=bit_Multiplication(-x,-x,n);
		}
		else
		{
            y=bit_Multiplication(x,x,n);
		}
		if (y>n / 2)
		{
			y -= n;//利用同餘原理，取較小的餘數
		}
		if (y == 1 && x != 1 &&x!=-1&& x != n - 1)
		{
			return true;
		}
		x=y;
	}
	if (y != 1)
	{
		return true;
	}
	return false;
}
int MILLER_RABIN(unsigned __int64 n, unsigned __int64 s, int b[])
{
	unsigned __int64 t = bit_Multiplication(n,1,2);
	if (n==2||n==3)
	{
		return PRIME;
	}
	if ((!t))
	{
		return COMPOSITE;
	}
		if( n < 1373653 )  
		{  
			if( WITNESS(2, b, n)   
				|| WITNESS(3, b, n)  )  
				return COMPOSITE;  
		}  
		else if( n < 9080191 )  
		{  
			if( WITNESS(31, b, n)    
				|| WITNESS(73, b, n)  )  
				return COMPOSITE; 
		}    
		else if( n < 4759123141 ) 
    { 
        if( WITNESS(2, b, n) 
                ||WITNESS(3, b, n) 
                ||WITNESS(5, b, n) 
                ||WITNESS(11, b, n) ) 
            return COMPOSITE; 
    } 
		else if( n < 2152302898747 ) 
		{ 
			if( WITNESS(7, b, n)
                ||WITNESS(3, b, n)
                || WITNESS(5, b, n)
                || WITNESS(2, b, n) 
                || WITNESS(11, b, n) ) 
				return COMPOSITE; 
		} 
		else 
		{ 
			if( WITNESS(7, b, n) 
                || WITNESS(3, b, n)
                || WITNESS(5, b, n)
                || WITNESS(2, b, n)
                || WITNESS(11, b, n)
                || WITNESS(31, b, n)
                || WITNESS(61, b, n)
                || WITNESS(73, b, n) ) 
				return COMPOSITE; 
    } 
     return PRIME;
}
unsigned __int64 POLLARD_RHO(unsigned __int64 n,int b[])//分解成某個n的因子，該因子可能不是最小的。
{
   int i=1;
   if (MILLER_RABIN(n,1,b))return NULL;
   __int64 x;
   if (n<32767) x=rand()%n;
   else if(n<1073676289) x=rand()*rand();
   else if (n<35181150961663)x=rand()*rand()*rand();
   else if (n<1152780773560811521)x=rand()*rand()*rand()*rand();
   else return NULL;
   x=bit_Multiplication(x,1,n);//對於隨機超過待測試數n的x，我們對其求餘，保證它隨機數在(1,n)之間的數
   __int64 y=x;
   __int64 k=2;
   while (1)
   {
	   i++;
	   x=bit_Multiplication(x-1,x+1,n);
	   __int64 d=Euclid(y-x,n);
	   d=d>0?d:-d;
	   if (d!=1&&d!=n)
	   {
		   return d;
	   }
	   if (i==k)
	   {
		   y=x;
		   k<<=2;
	   }
   }
}
void main()
{
	int *b = new int[len];
	__int64 t=POLLARD_RHO(1152780773560811517,b);
}

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的方法：

<span style="color:#000000;">unsigned __int64 POLLARD_RHO(unsigned __int64 n,int b[])//分解成某個n的因子，該因子可能不是最小的。
{
   int i=1,flag=0,u=0,t=0;
   List *x=NULL,*head=NULL;//由於需要用陣列來儲存每一個x以便查詢重複的資料來確定rho迴路大小，但是我們不知道該陣列大小，所以我才用動態陣列。
   if (MILLER_RABIN(n,1,b))return NULL;
   x=new List[LEN];
   head=x;
   x->num=i;
   if (n<32767) x->key=rand()%n;//這個ifelse結構用來隨機(1,n)範圍內的資料初始化x
   else if(n<1073676289) x->key=rand()*rand();
   else if (n<35181150961663)x->key=rand()*rand()*rand();
   else if (n<1152780773560811521)x->key=rand()*rand()*rand()*rand();
   else return NULL;
   x->key=bit_Multiplication(x->key,1,n);//對於隨機超過待測試數n的x，我們對其求餘，保證它隨機數在(1,n)之間的數
   __int64 y=x->key;
   __int64 k=2;
   while (1)
   {
	   i++;
	   int z=x->key;
	   x->next=new List[LEN];
	   x=x->next;
	   x->next=NULL;
	   x->key=bit_Multiplication(z-1,z+1,n);
	   x->num=i;
	   __int64 d=Euclid(y-x->key,n);
	   d=d>0?d:-d;
	   if (d!=1&&d!=n&&flag<1)
	   {
		   printf("n的因子=%d\n",d);
		   flag++;
	   }
	   struct List*yy=head;
	   while (yy->num<x->num)
	   {
		   if (yy->key==x->key)
		   {
			   u=(x->num)-(yy->num);
			   t=yy->num;
			   cout<<"rho迴路長度u="<<u<<endl;
			   cout<<"rho尾長度t="<<t<<endl;
			   return d;
		   }
		   yy=yy->next;
	   }
	   if (i==k)
	   {
		   y=x->key;
		   k=2*k;
	   }
   }
}</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).假定每個減法運算，測試奇偶性運算以及折半運算都能在單位時間行。

程式碼如下：

//二進位制的gcd演算法
#include <iostream>
using namespace std;
int Abs(int x)
{
	return x<0?-x:x;
}
int Binary_gcd(int a,int b)
{
	if (a==0)
	{
		return b;
	}
	else if(!(a&1)&&!(b&1))
	{
		return 2*Binary_gcd(a>>1,b>>1);
	}
	else if ((a&1)&&(b&1))
	{
        if (a<b)
        {
			a = a^b;  
	        b = b^a;  
			a = a^b; 
        }
        return Binary_gcd(Abs((a-b)>>1),b);
	}
	else if (!(a&1)&&(b&1))
	{
		return Binary_gcd(a>>1,b);
	}
	else 
	{
		return Binary_gcd(a,b>>1);
	}
}
void main()
{
	cout<<Binary_gcd(2*3*7*13*5,17*3*14*2*5)<<endl;
	cout<<Binary_gcd(17,15)<<endl;
}

31-2(對歐幾里得演算法中位操作的分析)

a.考慮用普通的“紙和筆”演算法來實現長除法的運算，用a除以b，得到商q和餘數r。證明，這種演算法需要執行O((1+lgq)lgb)次位操作。

以下程式碼完全按照用“紙和筆”手寫計算除法的方式計算二進位制除法：

<span style="color:#000000;">//我不清楚如何達到O((1+lgq)lgb，但是除去log函式可能用到的位操作，這裡只使用了O(lgqlgb)次位操作
#include <iostream>
#include<math.h>
using namespace std;
__int64 Dev(__int64 a, __int64 b)
{
	__int64 s = 1, i = 0, ans = 0,j=1;
	__int64 alength = log(a) / log(2);
	while (alength>i )//O(lgqlgb)
	{//記憶體迴圈執行次數與外層迴圈執行次數乘積就是移位操作符執行的總次數
		while (s < b&&alength>i++)//O(lgb)
		{
			s <<= 1;//位操作
			if (a&j << alength - i )s++;
			ans <<= 1;//位操作
		}
		if (s >= b)//O(lgq)
		{//if語句執行的次數完全取決於商的二進位制1的個數，假設商每位都是1，那麼就達到了if語句執行次數最大值，也就是lgq
			ans++;
			s -= b;                                                            
		}
	}
	printf("商=%I64d\n",ans);
	return s;
}
void main()
{
	
	printf("餘數=%I64d", Dev(2561, 147));
}</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+β)²) 去掉低次數項就得到結論。

其中的求餘數的位操作程式碼如下：

<span style="color:#000000;">#include <iostream>
using namespace std;
int Dev(int a,int b)
{
   int ans=0,j=1,temp=a;
   while ((temp>>=1))j++;//計算a的位數為lga
   for (int i=j;i>=0;i--)//O(lga+1)
   {
	   if ((a>>i)>=b)
	   {
		   ans+=(1<<i);
		   a-=(b<<i);
	   }
   }
   cout<<"餘數="<<a<<endl;
   return ans;
}
void main()
{
   cout<<"商="<<Dev(1157923,4713)<<endl;
}</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的冪

                <wbr>f(10)
              <wbr><wbr><wbr><wbr><wbr><wbr><wbr><wbr><wbr><wbr><wbr><wbr><wbr><wbr> <wbr>/   <wbr><wbr> <wbr> <wbr> <wbr> <wbr> <wbr>\
  <wbr><wbr>          <wbr><wbr><wbr><wbr><wbr><wbr><wbr><wbr><wbr><wbr>f(9) <wbr> <wbr> <wbr> <wbr> <wbr> <wbr>  <wbr><wbr> <wbr>f(8)
       <wbr><wbr><wbr><wbr><wbr><wbr><wbr>   <wbr><wbr><wbr>/ <wbr>  <wbr><wbr> <wbr> <wbr>\ <wbr> <wbr> <wbr>  <wbr><wbr> <wbr> <wbr>/  <wbr> <wbr> <wbr>\
 <wbr>   <wbr><wbr><wbr>   <wbr><wbr>f(8)   <wbr><wbr>  <wbr><wbr>f(7)  <wbr>f(7) <wbr> <wbr> <wbr>f(6)
      <wbr><wbr><wbr><wbr><wbr><wbr>/   <wbr><wbr>\     <wbr><wbr><wbr><wbr><wbr>/  <wbr><wbr> <wbr>\<wbr><wbr><wbr><wbr><wbr><wbr><wbr><wbr><wbr><wbr><wbr><wbr><wbr><wbr><wbr><wbr><wbr><wbr> 
<wbr>   <wbr><wbr><wbr>f(7)  <wbr>f(6)  f(6) f(5)<wbr><wbr><wbr>

b.試說明如何運用記憶法在O(n)時間內計算Fn.

說白了，就是倒著求，預求Fn,先求F0與F1，然後儲存記憶下F0與F1的值，再求F2=F1+F0的值，以此類推直到求得Fn，所以可以做一個迴圈

程式碼如下：

#include <iostream>
using namespace std;
#define n 10
double Fibonacci()
{
	double  F[n+1]={0,1};
    for (int i=0;i<n-1;i++)
    {
		F[i+2]=F[i]+F[i+1];
    }
	return F[n];
}
void main()
{
	cout<<Fibonacci()<<endl;
}

c.試說明如何僅用整數加法和乘法運算，就可以在O(lgn)的時間內計算Fn.考慮2X2矩陣[0,1,1,1]。

Fn為[0,1,1,1]^n的第一行第一列，這個可用數學歸納法證明，下面是程式碼：

#include<iostream>
#include <assert.h> 
using namespace std;
struct Matrix2X2
{
	Matrix2X2//2X2矩陣
		(int F00 , int F01 , int F10 , int F11 ):F_00(F00), F_01(F01), F_10(F10), F_11(F11){}	
	int F_00;
	int F_01;
	int F_10;
	int F_11;
};
Matrix2X2 MatrixMultiply(const Matrix2X2& matrix1, const Matrix2X2& matrix2)
{//兩個矩陣乘積
	return Matrix2X2(
		matrix1.F_00 * matrix2.F_00 + matrix1.F_01 * matrix2.F_10,
		matrix1.F_00 * matrix2.F_01 + matrix1.F_01 * matrix2.F_11,
		matrix1.F_10 * matrix2.F_00 + matrix1.F_11 * matrix2.F_10,
		matrix1.F_10 * matrix2.F_01 + matrix1.F_11 * matrix2.F_11);
}
// n階矩陣
// 0  1
// 1  1
Matrix2X2 Fibonacci(unsigned int n)
{
	assert(n > 0);	
	Matrix2X2 matrix(0,1,1,1);
	if(n==1)
	{
        matrix = Matrix2X2(0, 1, 1, 1);
	}
	else if(!(n&1))
	{
		matrix = Fibonacci(n>>1);
		matrix = MatrixMultiply(matrix, matrix);
	}
	else if(n&1)
	{
		matrix = Fibonacci((n-1)>>1);
		matrix = MatrixMultiply(matrix, matrix);
		matrix = MatrixMultiply(matrix, Matrix2X2(0, 1, 1, 1));
	}
	return matrix;
}
void main()
{
    for (int i=0;i<=30;i++)
    {
		cout<<Fibonacci(i+1).F_00<<endl;
    }
}

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同餘代表有二次剩餘，否則就是非二次剩餘。

以下是程式碼：

//用反覆平方法求數的冪 也稱為模取冪
#include <iostream>
using namespace std;
#define len 10
#define p 7
#define aa 34
int* BIT(int bb,int b[])
{
   int i=0;   
   while (i!=len&&b!=0)
   {
	   b[i]=bb%2;
	   bb=bb/2;
	   i++;
   }
   return b;
}
//從左向右檢查b的各位順序
int MODULAR_EXPONENTIATION(int a,int b[],int n)
{
  int c=0;
  int d=1;
  for (int i=len-1;i>=0;i--)
  {
	  c=2*c;
	  d=(d*d)%n;
	  if (b[i]==1)
	  {
		  c++;
		  d=(d*a)%n;
	  }
  }
  return d;
}
void main()
{
	int *b=new int[len];
	//判斷a^((p-1)/2)(mod p)是否為1，若為1，則有二次餘數。
	cout<<MODULAR_EXPONENTIATION(aa,BIT((p-1)/2,b),p)<<endl;//結果是-1，說明是二次非剩餘
	cout<<MODULAR_EXPONENTIATION(254,BIT((p-1)/2,b),p)<<endl;//結果是1，說明是二次剩餘	
}

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演算法解決。