算法导论31（数论算法）

31.1 基础数论概念

定理31.2：gcd(a,b)=min{ax+by|ax+by>0}\gcd (a,b) = \min \{ ax + by|ax + by > 0\}

证：

s=ax+bygcd(a,b)|s⇒gcd(a,b)≤sa%s=a−s⌊a/s⌋=a−(ax+by)⌊a/s⌋=a(1−⌊a/s⌋x)+b(−⌊a/s⌋y)0≤a%s<s⇒a%s=0⇒s|a,s|b⇒gcd(a,b)≥s\begin{array}{l}
s = ax + by\\
\\
\gcd (a,b)|s \Rightarrow \gcd (a,b) \le s\\
\\
a\% s = a - s\left\lfloor {a/s} \right\rfloor = a - (ax + by)\left\lfloor {a/s} \right\rfloor = a(1 - \left\lfloor {a/s} \right\rfloor x) + b( - \left\lfloor {a/s} \right\rfloor y)\\
0 \le a\% s < s \Rightarrow a\% s = 0 \Rightarrow s|a,s|b \Rightarrow \gcd (a,b) \ge s\\
\end{array}

推论31.3：d|a,d|b⇒d|gcd(a,b)d|a,d|b \Rightarrow d|\gcd (a,b)

证：gcd(a,b)=ax+by\gcd (a,b) = ax + by

推论31.5：n|ab,gcd(a,n)=1⇒n|bn|ab,\gcd (a,n) = 1 \Rightarrow n|b

证：ab=kn,ax+by=1⇒n(kx+by)=bab = kn,ax + by = 1 \Rightarrow n(kx + by) = b

定理31.6：gcd(a,p)=gcd(b,p)=1⇒gcd(ab,p)=1\gcd (a,p) = \gcd (b,p) = 1 \Rightarrow \gcd (ab,p) = 1

证：

ax+py=1,bx′+py′=11=(ax+py)(bx′+py′)=ab(xx′)+p(axy′+bx′y+pyy′)\begin{array}{l}
ax + py = 1,bx' + py' = 1\\
1 = (ax + py)(bx' + py') = ab(xx') + p(axy' + bx'y + pyy')
\end{array}

31.2 最大公约数

定理31.9：gcd(a,b)=gcd(b,a%b)\gcd (a,b) = \gcd (b,a\% b)

证：

d=gcd(a,b)⇒d|a,d|ba%b=a−b⌊a/b⌋⇒d|(a%b)⇒d|gcd(b,a%b)⇒gcd(a,b)|gcd(b,a%b)d=gcd(b,a%b)⇒d|b,d|(a%b)a=b⌊a/b⌋+a%b⇒d|a⇒d|gcd(a,b)⇒gcd(b,a%b)|gcd(a,b)\begin{array}{l}
d = \gcd (a,b) \Rightarrow d|a,d|b\\
a\% b = a - b\left\lfloor {a/b} \right\rfloor \Rightarrow d|(a\% b)\\
\Rightarrow d|\gcd (b,a\% b) \Rightarrow \gcd (a,b)|\gcd (b,a\% b)\\
\\
d = \gcd (b,a\% b) \Rightarrow d|b,d|(a\% b)\\
a = b\left\lfloor {a/b} \right\rfloor + a\% b \Rightarrow d|a\\
\Rightarrow d|\gcd (a,b) \Rightarrow \gcd (b,a\% b)|\gcd (a,b)
\end{array}

引理31.10：如果a>b≥1a > b \ge 1并且EUCLID(a,b){\rm{EUCLID}}\left( {a,b} \right)执行了k≥1k \ge 1次递归调用，则a≥Fk+2,b≥Fk+1a \ge {F_{k + 2}},b \ge {F_{k + 1}}。

证：数学归纳法

a≥F3=2,b≥F2=1b≥Fk+1,a%b≥Fka>b⇒a=b⌊a/b⌋+a%b≥b+a%b=Fk+1+Fk=Fk+2\begin{array}{l}
a \ge {F_3} = 2,b \ge {F_2} = 1\\
\\
b \ge {F_{k + 1}},a\% b \ge {F_k}\\
a > b \Rightarrow a = b\left\lfloor {a/b} \right\rfloor + a\% b \ge b + a\% b = {F_{k + 1}} + {F_k} = {F_{k + 2}}
\end{array}

定理31.11（Lame定理）：对任意整数k≥1k \ge 1，如果a>b≥1a > b \ge 1，且b<Fk+1b < {F_{k + 1}}，则EUCLID(a,b){\rm{EUCLID}}\left( {a,b} \right)的递归调用次数少于kk次。

证：反证法

如果EUCLID(a,b){\rm{EUCLID}}\left( {a,b} \right)的递归调用次数不少于kk次，则根据引理31.10，b≥Fk+1b \ge {F_{k + 1}}，与b<Fk+1b < {F_{k + 1}}矛盾。

欧几里得算法的扩展形式

d=gcd(a,b)=ax+byd′=gcd(b,a%b)=bx′+(a%b)y′d=d′=bx′+(a%b)y′=bx′+(a−b⌊a/b⌋)y′=ay′+b(x′−⌊a/b⌋y′)x=y′,y=x′−⌊a/b⌋y′\begin{array}{l}
d = \gcd (a,b) = ax + by\\
d' = \gcd (b,a\% b) = bx' + (a\% b)y'\\
d = d' = bx' + (a\% b)y' = bx' + (a - b\left\lfloor {a/b} \right\rfloor )y' = ay' + b(x' - \left\lfloor {a/b} \right\rfloor y')\\
x = y',y = x' - \left\lfloor {a/b} \right\rfloor y'
\end{array}

//欧几里得算法
unsigned Euclid(unsigned a,unsigned b)
{
if(b==0)return a;
else return Euclid(b,a%b);
}

//欧几里得算法的扩展形式
struct result
{
unsigned d;
int x,y;
result(unsigned i,int j,int k):d(i),x(j),y(k){}
};

result extendedEuclid(unsigned a,unsigned b)
{
if(b==0)return result(a,1,0);
else
{
result res=extendedEuclid(b,a%b);
return result(res.d,res.y,res.x-a/b*res.y);
}
}

31.4 求解模线性方程

ax≡b(modn)ax \equiv b(\bmod n)

void modularLinearEquationSolver(unsigned a,int b,unsigned n)
{
result res=extendedEuclid(a,n);
if(res.d%b==0)
{
unsigned x0=res.x*(b/res.d)%n;
for(int i=0;i<res.d;++i)cout<<(x0+i*(n/res.d))%n<<' ';
cout<<endl;
}
else cout<<"no solutions"<<endl;
}

31.6 元素的幂

abmodn{a^b}\bmod n

//反复平方法
unsigned modularExponentiation(unsigned a,unsigned b,unsigned n)
{
unsigned d=1;
for(int k=31;k>=0;--k)
{
d=d*d%n;
if(b&(1<<k))d=d*a%n;
}
return d;
}

31.8 素数的测试

伪素数测试过程

an−1≡1(modn){a^{n - 1}} \equiv 1(\bmod n)

bool pseudoprime(unsigned n)
{
return modularExponentiation(2,n-1,n)==1;
}

Miller-Rabin随机性素数测试方法

bool witness(unsigned a,unsigned n)
{
unsigned t=0;
for(;t<32;++t)
{
if((n-1)&(1<<t))break;
}
unsigned u=(n-1)>>t;
unsigned x0=modularExponentiation(a,u,n);
unsigned x1=0;
for(unsigned i=1;i<=t;++i)
{
x1=x0*x0%n;
if(x1==1&&x0!=1&&x0!=n-1)return true;
x0=x1;
}
if(x1!=1)return true;
return false;
}

bool MillerRabin(unsigned n,unsigned s)
{
for(unsigned j=1;j<=s;++j)
{
unsigned a=rand()%(n-1)+1;
if(witness(a,n))return false;
}
return true;
}

31.9 整数的因子分解

int gcd(int a,int b)
{
return b?gcd(b,a%b):a;
}

void PollardRho(int n)
{
int i=1,x=rand()%n,y=x,k=2;
while(true)
{
++i;
x=(x*x-1)%n;
int d=gcd(y-x,n);
if(d!=1&&d!=n)cout<<d<<' ';
if(i==k)
{
y=x;
k*=2;
}
}
}