第31章：数论算法

最大公约数：

欧几里得算法：

对于任意非负整数a和任意正整数b来说，当求它们的最大公约数时，有如下定理：gcd(a,b)=gcd(b,a mod b)；根据这个定理，可以写出求最大公约数的代码：

// a和b是任意非负整数；

//递归求解；
unsigned long Euclid_Rec(unsigned long a,unsigned long b)
{
        if(b==0)
                return a;
        else
                return Euclid_Rec(b,a%b);
}

//迭代求解：
unsigned long Euclid_Iter(unsigned long a,unsigned long b)
{
        while(b!=0){
                unsigned long temp=a%b;
                a=b;
                b=temp;
        }

        return a;
}

欧几里得算法的扩展形式：

我们知道gcd(a,b)是a与b的线性组合，可以写成d=gcd(a,b)=a*x+b*y的形式，我们可以推广欧几里得算法算出满足上述形式的d和系数x，y。

//计算满足d=gcd(a,b)=ax+by的d和系数x,y;
//tuple中存储的第一个元素为d=gcd(a,b)，第二个元素为x，第三个元素为y;
tuple<long,long,long> extend_Euclid(unsigned long a,unsigned long b)
{
        if(b==0)
                return make_tuple(a,1,0);

        auto temp=extend_Euclid(b,a%b);
        return make_tuple(get<0>(temp),get<2>(temp),get<1>(temp)-get<2>(temp)*(a/b));
}

在extend_Euclid的过程中，首先测试b是否为0，如果为0，则gcd(a,b)=d=a，x=1，y=0；如果b不为0，则计算出d’=b*x’+y’(a%b)。我们知道d=d’，则d=d’=b*x’+y’(a-b*(a/b))=a*y’+b*(x’-(a/b)*y’)=a*x+b*y，则x=y’，y=x’-(a/b)*y’。注：在这里a/b的结果是整数，如果商包含小数部分，会被截除掉。

求解模线性方程：

当我们要求解(ax)(mod n)=b(mod n)这个方程时，有如下定理可以可以应用：

假设方程 (ax)(mod n)=b(mod n) 有解（即b能够整除d，d=gcd(a,n)），并且x0（d=a*x’+n*y’，x0=(x’(b/d))(mod n)）是该方程的任意一个解，则该方程有d个不同的解，分别为xi=x0+i*(n/d)，这里i=0,1,2,…,d-1。

依据上述的定理，可以写出如下的程序来求模线性方程：

//计算满足(ax)(mod n)=b(mod n)的所有大于0的x;
//a和n为任意正整数，b为任意整数
void modular_linear_equation_solver(unsigned long a,long b,unsigned long n)
{
        auto temp=extend_Euclid(a,n);//d=a*x+n*y;

        auto d=get<0>(temp);
        auto x=get<1>(temp);
        auto y=get<2>(temp);

        if(b%d==0){
                long temp=x*(b/d);
                //确保x*(b/d)大于0，因为如果x*(b/d)小于0，则下面x*(b/d)%n的结果的正负号将是机器依赖的;
                //比如x*(b/d)=-30, n=100，(-30)%100有的机器上是-30，有的是70，但我们应该要使结果为70;   
                if(temp<0){
                        int k=0;
                        temp=-temp;
                        while(++k){
                                if(temp<k*n){
                                        temp=k*n-temp;
                                        break;
                                }
                        }
                }

                unsigned long x0=temp%n;
                for(int i=0;i!=d;++i)
                        cout<<(x0+i*(n/d))%n<<" ";
                cout<<endl;
        }
        else
               cout<<"no solutions to the equation!"<<endl;

        return;
}

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

// eval the value of (a^p)%n;
unsigned long evalRemainder(unsigned long a,unsigned long p,unsigned long n)
{
        unsigned long r=a%n;
        unsigned long tmp=1;

        while(p>1){
                if(p%2!=0)
                        tmp=(tmp*r)%n;
                        if((double)(r)<=sqrt(ULONG_MAX))
                                r=(r*r)%n;
                        else
                                throw runtime_error("the number is so large that it surpass the limit of unsiged long");
                p=p/2;
        }
        return (r*tmp)%n;
}

bool rabbinMiller(unsigned long a,unsigned long k,unsigned long q,unsigned long n)
{
        if(evalRemainder(a,q,n)!=1)
        {
                int e=1;
                for(int i=0;i!=k;++i)
                {
                        if(evalRemainder(a,q*e,n)==n-1)
                                return true;
                        e*=2;
                }
                return false;
        }

        return true;
}

bool isPrime(unsigned long n)
{
        if(n<2)
                return false;

        if(n==2||n==3)
                return true;

        if(n%2==0)
                return false;

        // n=2^k*q+1;
        unsigned long k=0,q=0;
        unsigned long tmp=n-1;
        while(tmp%2==0){
                tmp=tmp/2;
                k++;
        }
        q=tmp;
/*

        if(n<1373653){
                if(rabbinMiller(2,k,q,n)==false||rabbinMiller(3,k,q,n)==false)
                        return false;
                else
                        return true;
        }
        else if(n<9080191){
                if(rabbinMiller(31,k,q,n)==false||rabbinMiller(73,k,q,n)==false)
                        return false;
                else
                        return true;
        }
        else if(n<4759123141){
                if(rabbinMiller(2,k,q,n)==false||rabbinMiller(3,k,q,n)==false||rabbinMiller(5,k,q,n)==false||rabbinMiller(11,k,q,n)==false)
                        return false;
                else
                        return true;
        }
               else if(n<2152302898747){
                if(rabbinMiller(2,k,q,n)==false||rabbinMiller(3,k,q,n)==false||rabbinMiller(5,k,q,n)==false||rabbinMiller(7,k,q,n)==false||rabbinMiller(11,k,q,n)==false)
                        return false;
                else
                        return true;
        }
*/
// make use of random number in the range[2,n-2] to test if n is prime number;

        const int loopNum=10;

        srand((unsigned)time(NULL));
        for(int i=0;i!=loopNum;++i)
                if(!rabbinMiller(rand()%(n-3)+2,k,q,n))
                        return false;
        return true;
}