初探Euclid算法

学过编程的人应该对Euclid算法（也叫辗转相除法）不陌生，许多书籍都把它作为递归的范例。
由于最近写到了一些涉及数论的题目，因此本人对其进行了更加深入的探究。

用Euclid算法求两个自然数的最大公约数
这是Euclid算法最常见的应用。它基于这么一个事实：
gcd(a,b)=gcd(b,a mod b)
《算法导论》中把这称作GCD递归定理，具体证明可以参考《算法导论》31.2小节。
31.2小节也讨论了Euclid算法的运行时间，结论Euclid执行中的递归调用次数是O(lgb)。

用Euclid求解二元一次不定方程的整数解
定理31.2：如果a和b不是都为0的任意整数，则gcd(a,b)是a与b的线性组合集合{ax+by}(x,y均为整数)中的最小正元素。
我们可以使用一个改进的Euclid算法（也叫做扩展Euclid算法）求解ax+by=gcd(a,b) 中x和y的值。
算法框架如下（其中d是a和b的最大公约数，x、y分别是a、b的系数）

ext-euclid(a,b)
  if(b = 0)
    then return (a,1,0)
(d',x',y')<--ext-euclid(b,a mod b)
(d,x,y)<--(d',y',x' - (a div b) * y')
return (d,x,y)

倒二行是Euclid算法的精髓。可以这样推导：
由GCD递归定理知d = gcd(a,b) = gcd(a, a % b)
那么就有d = ax+by = bx' + (a mod b)y'
又a mod b = a - (a div b) * b
因此d = ax + by = bx' + (a - (a div b) * b)y' = ay' + b(x' - a div b * y')
可知x = y'   y = x' - a div b * y'
当b为0时，显然x=1,y=0。

如果要求解不定方程ax+by=c（a,b,c均为整数）的整数解呢？

求解之前，要判断这个方程是否有整数解，这个方程有整数解的充要条件是gcd(a,b) | c。
然后，用扩展Euclid求解出ax+by=gcd(a,b)的一组解(x1,y1)。
那么(x1 * c / gcd(a,b),y1 * c / gcd(a,b))就显然是一组解了。
根据裴属定理，(x+b/gcd(a,b) * k,y + a/gcd(a,b) * k)，k为整数，也是这个方程的一组解。

一个计数问题——求散列函数碰撞次数

题目来源：NOI导刊
题目大意：某人设计了一个hash函数hash(x,y) = x * y + x + y,把平面上某个点(x,y)(x,y均为非负整数)映射到一个非负整数。现在输入H，要求你给出符合hash(x,y)=h的点的个数。
输入数据：第一行为测试数据个数T，后有T行H。
输出数据：T行，每一行为H对应的点的个数。
数据范围：0
解法一：
把函数变形为y=(h - x)/(x + 1)，然后对[1,h]这个范围内的x逐个测试(h-x) mod (x+1)是否为0。若为0，则ans++。
一个优化：容易知道hash(x,y) = hash(y,x)，据此可以把时间优化到一半。
这个算法整体时间上限是O(MaxT * MaxH)，只能通过3/10的测试点。

解法二：
把函数变形为H+1 = x * y + x + y + 1；因式分解得H+1 = (x + 1) * (y + 1)。
于是，原问题就转化成求H+1能表示成几对正整数相乘了。
故，把H+1质因数分解，然后用乘法原理计算即可。
由于H最大只有100,000,000，因此只要用筛法预处理[1,10000]内的质数就足够了。
需要注意的是，质因数分解后还可能剩下一个大于10000的质数，它的指数一定是1，故把答案×2即可。
解法二的时间复杂度近似T * sqrt(MaxH)，可以通过所有数据。

C语言: 高亮代码由发芽网提供

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define MAXP 20000

long totp = 0,p[10000];
int idx[10000];  //指数

void init_prime()   //普通筛法求素数
{
   static _Bool pb[MAXP];
   pb[1] = pb[0] = 1;
   long i,j;
   for(i = 2;i < MAXP;i ++)
       if(pb[i] == 0)
       {
           p[totp ++] = i;
           for(j = i + i;j < MAXP;j += i)
               pb[j] = 1;
       }
}

inline void solve(long x)
{
   x ++;
   long ans = 1;
   memset(idx,0,sizeof(idx));
   long i;
   for(i = 0;i < totp;i ++)
   {
       if(x == 1)    break;
       while(x % p[i] == 0)
       {
           idx[i] ++;
           x /= p[i];
       }
       ans *= idx[i] + 1;
   }
   if(x != 1)    ans *= 2;   //剩余一个大的质因数
   printf("%ld\n",ans);
}

int main()
{
   init_prime();
   freopen("hash.in","r",stdin);
   freopen("hash.out","w",stdout);
   int T;
   long x;
   scanf("%d\n",&T);
   while(T --)
   {
       scanf("%ld",&x);
       solve(x);
   }
   fclose(stdout);
   return 0;
}