加快特殊的取模运算

空木秋人 · 2020-10-19 20:53:45 · 个人记录

在使用 字符串哈希 和 哈希表 的时候，我们常常会遇到取模运算，当然这里并不是讨论快速乘法取模，或者是对于所有的a%b都有用，目前我知道的是如果选择合适的模可以利用一些位运算技巧加快计算。

insight 1.

首先我们要知道的一个很显然的取模技巧：if b=p+1 ,then (a b)%p=(ap+a)%p=a. 我把它叫做 insight 1.

有了上面的基础就可以解决这个问题: 如何快速计算 a%p，if p=2^m-1 .

我们把a用成 2^m 进制的数表示出来：

a= a0 + a1*(p+1)+a2*(p+1) +... +an*(p+1).

那么 a%p= a0 + a1 +a2+..+an ，即把a用p+1进制表示后base前面系数的和。例子如下：

127=001 111 001
// 8进制，3 位对1位
121%7=(1* 8^2 + 7 * 8 +1*8^0) %7 =(1 + 7+ 1)%7

insight 2

基于此，如果你知道16进制是二进制4位对1位，那么其它进制同理可以写出一段代码, 类似（p=2^{31}-1 ,因为x最多64位，因此计算两次就够了）：

x = (x >> 31) + (x & 2147483647);

x = (x >> 31) + (x & 2147483647);

return x;

但这段代码有问题，当x=mod的时候，无法得到正确的输出0，实际上输出范围是[1,2^{31}-1]. 因此要对它修正。我们可以在计算之前先让x+1,最后输出的时候再输出x-1. 这样可以对x=2^{31}-1修正，而又不影响其它数的答案。这是一个很有趣的技巧（我把它叫做insight 2）. 但是这样会引入另一个bug，当x=2^{64}-1时先加一会溢出得到0,最后减1的时候又会溢出而得到错误的结果。因此实际上我的代码改成如下：

//  x%(2^m-1)
 ull modulo(ull x, int m)
{
    ull mod = (1LL << m) - 1;
    while (x > mod)
    {
        x = (x >> m) + (x & mod);
    }
    if(x==mod)
        return 0;
    else {
        return x;
    }
}

insight 3

最后我们再来看如何计算(a* b)%(2^61-1). a和b都是unsigned long long类型，且小于2^61-1。如果我们不用__int128这样的类型存储（实际上很多oj不允许用这个类型），直接相乘会有溢出的问题。下面一代码段充分利用了位操作，可以准确得到结果。

const uint64_t mod = (1ull<<61) - 1;
uint64_t modmul(uint64_t a, uint64_t b){
    uint64_t l1 = (uint32_t)a, h1 = a>>32, l2 = (uint32_t)b, h2 = b>>32;
    uint64_t l = l1*l2, m = l1*h2 + l2*h1, h = h1*h2;
    uint64_t ret = (l&mod) + (l>>61) + (h << 3) + (m >> 29) + (m << 35 >> 3) + 1;
    ret = (ret & mod) + (ret>>61);
    ret = (ret & mod) + (ret>>61);
    return ret-1;
}

如前面所说这段代码运用了大量技巧，我花了好几个小时才得到结论（但是我相信有经验的选手可以在10分钟内得到相同的结论），解释如下：

把a，b拆成低32位和高29位, 取低位可以通过截断来实现，就像代码那样，也可以用&实现。我不清楚这两者哪个更好^-^
高位相乘之后应该左移64位，但是按照insight 1 ，我们只需要左移三位就够了
低位相乘之后不需要移位，但是可能结果太大，l1* l2<= (2^{32}-1)* (2^{32}-1)=2^{64}-2^{33}+1.所以按照insight 2,我们可以先对它做一次模运算。
中间位相乘之后需要左移32位，但是这可能导致溢出，我们实际上可以根据insight 1右移29位.

但是代码中的m<<35>>3 是我觉得最难以理解的。根据我几个小时的研究（事实上我写了一份去掉它的代码，然后看这两份代码有什么区别，在过了一个小时后我的机器还没有找出可以产生出差别的a，b，因此我放弃了暴力搜索，开始坐下来分析），最终得到如下解释：

m可能是一个比较小的数，在m>>29后得0导致错误。我们可以按照前面的想法改进它（原谅我不能找到更好的方法吧）。我们把m分成低29位和高位。这样根据insight 1,我们只需要计算高位加上低位左移32位的和就行了。为了得到低位，我们需要先移35位（总共64位，左移后先挤掉高位），再右移3位。

你可能想要一份完整的推导，就像这样：

$b= high2 * 2^{32} + low2$ \ $ a*b =high1 * high2 * 2^{64} +high1*low2* 2^{32}+ high2* low1* 2^{32} +low1*low2

$ \qquad m=(high1*low2+high2*low1)*2^{32}

$ \qquad mod=2^{61}-1

令 $ m/2^{32}= x* 2^{29} +y$ , 那么 $ m%mod=(x* 2^{61} + y * 2^32)%mod =(x+ y*2^{32})%mod

最后再讲一个很简单的吧，如果mod=2^m, 求a%mod. a%mod= a&(mod-1). 这同样可以用基于进制的思想来解释。

reference:

https://codeforces.com/blog/entry/60442

https://stackoverflow.com/questions/763137/computing-ab-mod-c-quickly-for-c-2n-1/763226#763226