算法导论31（数论算法）

31.1 基础数论概念

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

$$\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 \Rightarrow d|\gcd (a,b)$
证：$\gcd (a,b) = ax + by$

推论31.5：$n|ab,\gcd (a,n) = 1 \Rightarrow n|b$
证：$ab = kn,ax + by = 1 \Rightarrow n(kx + by) = b$

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

$$\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)$
证：

$$\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 \ge 1$并且${\rm{EUCLID}}\left( {a,b} \right)$执行了$k \ge 1$次递归调用，则$a \ge {F_{k + 2}},b \ge {F_{k + 1}}$。
证：数学归纳法

$$\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 \ge 1$，如果$a > b \ge 1$，且$b < {F_{k + 1}}$，则${\rm{EUCLID}}\left( {a,b} \right)$的递归调用次数少于$k$次。
证：反证法
如果${\rm{EUCLID}}\left( {a,b} \right)$的递归调用次数不少于$k$次，则根据引理31.10，$b \ge {F_{k + 1}}$，与$b < {F_{k + 1}}$矛盾。

欧几里得算法的扩展形式

$$\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 \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 元素的幂

${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 素数的测试

伪素数测试过程
${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;
        }
    }
}