写了一个超快的求余算法，求验证

前几天，在做“有道难题”的第一次练习赛中的第二题时，最终写了个自认为正确且速度超快的算法，但是我不敢肯定其正确性，也无法进行有效的验证，故现发帖广征大家的意见，并帮我验证或修改此算法，若已有类似的算法也望大家能贴出来以共勉。

此题大概要求如下：

1、a、b为正整数，其中，a > 0,b < 2^31

2、求a^b mod 9907

此例可用下面2种算法进行求解：

1、a^b = 9907*n+m，其中，n >= 0，0 <= m < 9907，m，n均为自然数，m即为所求余数。

2.、此处涉及两个公式：

1）a^b = (9907+m)^n，其中，m，n为自然数，a >= 9907，m^n mode 9907即为所求余数。

2）mn mode 9907 = [(m mod 9907)(n mod 9907)]mod 9907，其中，m，n为正整数。

用第一种原理求解，运算时间很长，而且当b很大时，a^b是难以直接求解的。故本算法采用算法2。

算法2似乎存在局限性，因为仅讨论了a >= 9907的情况，若a < 9907该怎么办？是的，a < 9907时用算法2不行，但是我可以将已知条件处理一下，取一正整数x，使a^x >= 9907，然后，将b赋值为b与x的商，最后得到等式a^b = ((a^x)^(b/x))*(a^(b%x))（b/x,b%x为C语言中的求商和求余）。这样(a^x)^(b/x)不就可以用算法2中的公式1求解了吗？最后再用公式2就可以求解a^b mod 9907了。

当b很大时，就需要用公式1不断降阶直到阶小于x时算法终止。用此算法可以实现快速降阶，因为每次x都不一样，而每次b = b/x，从而使b越来越小，最后使b < x，假设b = 2^30=1024*1024*1024，而a=2，则首次降阶x=14，b = (2^30)/14 < 1024*1024*74，a =(2^14)%9907 = 6477，看不出来有什么变化吗？那算一下1024*1024*1024 - 1024*1024*74等于多少吧。这可直接减掉了大约1024*1024*950的“赘肉”呀！接下来就只用求6477^(1024*1024*74) mod 9907了，如法炮制便可很快实现降阶。

所以，最后的算法概貌如下：

tmp=1; do{ if(a == 1) break; if(a % 9907 == 0){ tmp = 0; break; } x = log(a,9907); //向上取整，即取大于或等于log(a,9907)的最小整数 mod = a^x % 9907; //mod为中间结果的余数 tmp *= a^(b%x); tmp %= 9907; //随时计算余数，防止累乘结果过大 b /= x; a = mod; }while(b != 0); //b<x时结束循环 return tmp; //tmp即为所求余数

下面是我用C语言写的具体实现：

#include <stdio.h> #include <malloc.h> #include <math.h> unsigned int Power(int a,int b); int main() { int num = 0; int *a = NULL,*b = NULL; int m,n,i; printf("input num:"); scanf("%d",&num); if(num <= 0 || num > 10000) return -1; a = (int *)malloc(num*sizeof(int)); b = (int *)malloc(num*sizeof(int)); if(a == NULL || b == NULL) return -1; for(i = 0; i < num; i++){ printf("input %dth:",i); scanf("%d %d",&m,&n); if(m < 0 || n > (unsigned long(-1) - 1)) return -1; a[i] = m; b[i] = n; } for(i = 0; i < num; i++){ printf("%d/n",Power(a[i],b[i])); } return 0; } unsigned int Power(int a,int b) { unsigned int tmp = 1; unsigned int mod; unsigned int x; double y; do{ //中间结果的余为1，则tmp即为所求余数 if(a == 1) break; //中间结果的余为0，则a^b与9907的余也为0 if(a % 9907 == 0){ tmp = 0; break; } y = log(9907) / log(a); x = (int)floor(y); //取大于或等于log(9907)/log(a)的最大整数 if(y - x != 0.0) x += 1; mod = ((unsigned int)pow(a,x)) % 9907; //求tmp*a^(b%x)与9907的余数 tmp *= (unsigned int)pow(a,b%x); tmp %= 9907; b /= x; //降阶 a = mod; //求中间余数与9907的余 }while(b != 0); mod = tmp; return mod; }

恩，就是这些了。请大家帮我检查一下，该算法是否有问题，或是提供一些验证实例或方法。当然，若有更好的算法，也希望您能与大家分享，谢谢！