abmod与powmod

原文链接：http://blog.sina.com.cn/s/blog_66ad7bba0100hm8n.html

对于搞过竞赛算法的人来说，powmod可能不会陌生，它是一个计算a^b mod m的函数，

但abmod你可能不知道，它其实意思更简单，是计算a*b mod m的函数

powmod的出现看起来很自然，因为a^b可能非常巨大，但a*b的结果很小，有存在的必要吗？

 

比如a,b,m都是int，那么a*b的结果有可能越int，你会说，那用__int64或者long long保存就行了

但，如果a,b,m都是long long，你又不想使用麻烦的大整数呢？

 

其实，如果你懂了计算a^b mod m的原理，那解决这个也相当容易，不过还要先说一下：

a^b * a^c mod m == a^(b+c) mod m

这恒等式貌似是初中的内容吧，这个要是明白了，那来再下一步：

举例子：2^9 mod 5

== 2^8 * 2^1 mod 5

== ( (2^8 mod 5) * (2^1 mod 5) ) mod 5

== ( (2^4 mod 5)^2 * (2^1 mod 5) ) mod 5

== ( (2^2 mod 5)^4 * (2^1 mod 5) ) mod 5

== ( (2^1 mod 5)^8 * (2^1 mod 5) ) mod 5

从这一系列等式中你看出些什么？

现在，我们来假设一下，我们只懂算乘法，去通过这些等式算结果。

首先2^1 mod 5 == 2，于是有

(2^8 mod 5) * (2 mod 5)

2^8不会算，但我们会把2^8变为(2^2)^4

于是有4^4 mod 5 * 2

然后有(4^2 mod 5)^2 mod 5 * 2

得到1^2 mod 5 * 2

最后结果为2，而实际上，2^9 = 512， 512 mod 5 == 2

 

可能你会问，这样好像程序不好实现啊？不会的，我们换个角度看

 

a^b，把b拆成二进制，像刚刚的a^9 mod m，可以拆为a^(2^3 + 2^0) mod m

于是，可以变形为(a^(2^3) mod m) * (a^(2^0) mod m) mod m 间接算得

并且，a^(2^3) mod m的结果，可以由 (a^(2^2) mod m) ^ 2 mod m 获得

于是，我们只要保证平方运算不会越界，就能通过递推得到a^b mod m的结果

 

示例代码如下：

int powmod( int a, int n, int k )
{
    int d = 1;
    for (a %= k; n > 0; n >>= 1)
    {
        if(n & 1)
            d = (d*a)%k;
        a = (a*a)%k;
    }
    return d;
}

 

好了，可能有人未必看的懂前面的，现在这里就讲一下简单一些的abmod

a*b mod m，可能会越界，这是刚刚解释过的，我们先个个假设这三个数都是byte

现在要算99 * 66 mod 100

我们把式子改写为 (99 * (64 + 2)) mod 100

== ( 99 * 64 + 99 * 2 ) mod 100

其中99*64的可以由99*32递推出，然后又可以从99*16的递推出。。。。

 

于是，我们可以知道，只要m不大于char最大值的一半，就可以保证结果是正确的

参考代码：

INT abmod(INT a, INT b, INT m)
{
    a %= m; b %= m;
    INT s = 0;
    for (INT i=b; i>0; a = (a<<1)%m,i>>=1)
        if (i&1) s = (s+a) % m;
    return s;
} 

 

讲解的最后，再回答一个问题：为什么采用二进制分解？

答，为了提高可适用范围。如果是三进制，前面的powmod要保证m^3不越界，后面的abmod要保证3*m不越界，

于是m的范围就小了，精确性低了。

 

好了，admod的问题解决了，而这时候你发现，前面的powmod是因为有乘法，所以令m的平方不能越界，

如果，我们把这两个函数同时用上呢？

同时用上的话，我们可以把powmod扩展成superpowmod！可以在__int64下不借助大整数能正常工作！！

因为最大的瓶颈是m*2不能越界，并且计算过程我们不需要负数，于是，参数类型我们用unsigned

就能保证__int64下能正常工作！比传统的powmod可计算范围大大增加哦！

 

终极代码：

typedef unsigned __int64 INT;
INT abmod(INT a, INT b, INT m)
{
    a %= m; b %= m;
    INT s = 0;
    for (INT i=b; i>0; a = (a<<1)%m,i>>=1)
        if (i&1) s = (s+a) % m;
    return s;
}

INT superPowmod( INT a, INT n, INT k )
{
    INT d = 1;
    for (a %= k; n > 0; n >>= 1)
    {
        if(n & 1)
            d = abmod(d,a,k); //(d*a)%k;
        a = abmod(a,a,k); //(a*a)%k;
    }
    return d;
}