等比数列求和的二分法

常见问题:给定整数M,等比数列${{a_n}}$
求其前n项和对M的模
朴素的想法,我们可以用等比数列求和公式$S_n=\frac {a_1(1-q^n)}{1-q}$求得和再模除M,但是计算过程中取模导致分母下面的1-q不能直接除,要求1-q的乘法逆元;
这里介绍另一种思路,就是使用二分的思想,很像矩阵快速幂 首先,我们先讨论一个简单的问题:
求$a^1+a^2+a^3+a^4+a^5+a^6+a^7+a^8$的和
提取公因式,我们能把原式变为
$(a^1+a^2+a^3+a^4)(1+a^4)$
进一步变换为
$(a^1+a^2)(1+a^2)(1+a^4)$
$=a(a+a^1)(1+a^2)(1+a^4)$
我们把这种等比数列前$2^k$项和记为$S(k)$有递推式$S(k)=S(k-1)(1+a^{2^{k-1}})$
这是二分法的关键,也就是说我们能在log(n)的时间内求出等比数列前$2^k$项的和
我们要求的不是前2的整数次幂的时候我们将n分解为2的整次幂的和;
比如求前等比数列11项的和,我们将11(二进制表示1011)分解为1+2+8(二进制1+10+1000) 前11项的和可以表示为:
$(a^1)+(a^2+a^3)+(a^4+a^5+a^6+a^7+a^8+a^9+a^{10}+a^{11})$
可进一步表示为
$a^0(a^1)+a^1(a^1+a^2)+a^3(a^1+a^2+a^3+a^4+a^5+a^6+a^7+a^8)$
即:$a^0S(0)+a^1S(1)+a^3S(3)$
Well done!我们又把原式分解成了多个$a^xS(y)$的和 在这里可能可能会有一个小小的疑问,就是$a^xS(y)$中的x和y是怎么确定的,这里我是参考的矩阵快速幂的思想,x就是已经处理的指数的和,y根据当前处理的比特位置决定
一个快速幂大概这么写

int pow(int a,int b)
{
    int ans=1,power=a;
    while (b)
    {
        if (b&1)
            ans*=power;
        power*=power;
        b>>=1;
    }
    return ans;
}

注意,这里面的power一直充当$a^{2^k}$的角色,在等比数列求和中,他担任的功能还包括提供$S(k)=S(k-1)(1+a^{2^{k-1}})$中的$2^{k-1}$ 下面是完整的样例程序:

//等比数列求和
//author:comzyh
//等比数列求和
#include <cstdio>
const int MOD=9901;
int A,B;
int sum(int n,long long k)
{
    n%=MOD;
    int ans=0;
    int power=n;
    int powersum=n;//S(y)
    int mul=1;//a^x
    while (k)
    {
        if (k&1)
        {
            ans+=mul*powersum;  ans%=MOD;
            mul*=power;         mul%=MOD;
        }
        powersum*=(power+1);    powersum%=MOD;
        power*=power;           power%=MOD;
        k>>=1;
    }
    return ans;
}
int main()
{
    int n,k;
    while (scanf("%d%d",&n,&k)!=EOF)
    {
        printf("%d\n",sum(n,k));
    }
    return 0;
}