分治-约数之和

假设现在有两个自然数A和B，S是AB的所有约数之和。
请你求出S mod 9901的值是多少。

输入格式
在一行中输入用空格隔开的两个整数A和B。

输出格式
输出一个整数，代表S mod 9901的值。

数据范围
0≤A,B≤5×107

输入样例：
2 3

输出样例：
15
注意: A和B不会同时为0。

思路：
1、先将a分解质因数，即写成$a = (p1 ** k1) * (p2 ** k2) * (p3 ** k3) …(pn ** kn) $，那么a**b就变成了 $a ** b = (p1 ** (b * k1)) *(p2 ** (b * k2)) * (p3 ** (b * k3)) ……(pn ** (b * kn))$a∗∗b=(p1∗∗(b∗k1))∗(p2∗∗(b∗k2))∗(p3∗∗(b∗k3))……(pn∗∗(b∗kn))(这里因为markdown语法的原因把乘方写成了两个乘号)

2、分解完质因数之后，由乘法分配律我们可以知道$a**b$a∗∗b的约数之和就是
$（1 + p1 + p1 ** 2 + p1 ** 3 + …… p1 ** (b * k1) * (1 + p2 + p2 ** 2 + p2 ** 3 + p2 ** (b * k2)) ……(1 + pn + pn2 + pn3 + pn**(b*kn)) $

3、我们用sum(p,k)表示1 + p + p**2 + p ** 3 + …… + p * *(k)
那么如果k是奇数
$sum(p,k) = (1 + p + …… + p * * ((k-1)/2)) + (p ** ((k + 1) / 2) + …… + p**k)$sum(p,k)=(1+p+……+p∗∗((k−1)/2))+(p∗∗((k+1)/2)+……+p∗∗k)
进一步$sum(p,k) = (1 + p + …… + p**(k+1)/2) + p**((k+1)/2)*(1 + p + …… + p**(k+1)/2)$sum(p,k)=(1+p+……+p∗∗(k+1)/2)+p∗∗((k+1)/2)∗(1+p+……+p∗∗(k+1)/2)
所以$sum(p,k) = (1+ p**(k + 1)/2)*sum(p,(k-1)/2)$sum(p,k)=(1+p∗∗(k+1)/2)∗sum(p,(k−1)/2)

同理 当k为偶数的时候
$sum(p,k) = (1 + p ** (k / 2)) * sum(p , k/2 - 1) + p **k$sum(p,k)=(1+p∗∗(k/2))∗sum(p,k/2−1)+p∗∗k

因此我们可以想到用分治的思想来解决

代码：

#include<iostream>
using namespace std;

typedef long long ll;
const int mod = 9901;

ll a,b;
ll res;
ll fastpow(ll x , ll y)//快速幂函数
{
    ll ans = 1;
    while(y)
    {
        if(y & 1)
            ans = ans * x % mod;
        y >>= 1;
        x = (x * x) % mod;
    }
    return ans;
}

ll sum(ll p,ll k)
{
    if(k == 0)
        return 1;
    else if(k % 2 == 0)
        return ((1 + fastpow(p , k / 2)) * sum(p , k / 2 - 1) + fastpow(p,k)) % mod;
    else
        return((1 + fastpow(p , (k + 1) / 2)) * sum(p , (k - 1) / 2)) % mod;
}

int main(void)
{
    while(cin  >> a >> b)
    {
        res = 1;
        for(int i = 2 ; i <= a ; i ++)
        {
            ll cnt = 0;
            while(a % i ==0)//分解质因数
            {
                cnt += 1;
                a /= i;
            }
            if(cnt)
            {
                res = res * sum(i,cnt * b) % mod;
            }
        }
        
        if(a == 0)//特判
            cout << 0 << endl;
        else
            cout << res << endl;
    }
}