Sumdiv poj1845（逆元+快速幂+快乘）

本作品采用知识共享署名-相同方式共享 4.0 国际许可协议进行许可。
题目链接：http://poj.org/problem?id=1845
Description
Consider two natural numbers A and B. Let S be the sum of all natural divisors of A^B. Determine S modulo 9901 (the rest of the division of S by 9901).

Input
The only line contains the two natural numbers A and B, (0 <= A,B <= 50000000)separated by blanks.

Output
The only line of the output will contain S modulo 9901.

Sample Input

2 3

Sample Output

15

Hint
2^3 = 8.
The natural divisors of 8 are: 1,2,4,8. Their sum is 15.
15 modulo 9901 is 15 (that should be output).

Source
Romania OI 2002

分析：
我好菜呀。又wa了好久。结果还是要用快乘。 搞得我现在看到啥都想打个快乘了。

约数和定理：
a所有的正约数的和为:
(1 + p1 + p1 ^ 2 + …… + p1 ^ c1) * …… * (1 + pm + pm ^ 2 + …… + pm ^ cm)
观察可以得知，每一个括号里面的值是一个等比数列。
可得a^b 的正约数之和为：
(1 + p1 + p1 ^ 2 + …… + p1 ^ (c1b)) * …… * (1 + pm + pm ^ 2 + …… + pm ^ (cmb))

所以我们先把a的所有质因子分解出来。
每一个括号都是一个等比数列，自然我们可以用等比数列来求值，不过没有这么简单。
等比数列求和公式：
(p1^(bc1+1) - 1) / (p1- 1)
对于p1^(bc1+1)，我们要用到快速幂 （一定要快乘，不然爆longlong 绝望）
对于除上(p1- 1)，我们应该要用逆元来进行计算。
此时引入费马小定理：
若mod为质数，则p^(m -2) 为p的逆元。
但是需要注意一点，当p1-1为mod的倍数时，逆元不存在。
但是我们可以的值（p1 - 1） % mod == 1.故：
(1 + p1 + p1 ^ 2 + …… + p1 ^ (c1*b)) = 1 + 1 + …… + 1 = (b * c1 + 1)% mod;

#include"stdio.h"
#include"string.h"
#include"algorithm"
using namespace std;
typedef long long ll;
const ll mod = 9901;

ll a,b,top;
ll p[10100],c[10010];

void divide()
{
    for(int i = 2; i * (ll)i <= a; i ++)
    {
        if(a % i)
            continue;
        p[++ top] = i;
        while(a % i == 0)
        {
            a = a / i;
            c[top] ++;
        }
    }
    if(a != 1)
    {
        p[++ top] = a;
        c[top] ++;
    }
}
ll multi(ll a, ll b)
{
    ll ans = 0;
    while(b)
    {
        if(b & 1)
            ans = (ans + a) % mod;
        a = (a << 1) % mod;
        b = b >> 1;
    }
    return ans % mod;
}
ll Quick_Pow(ll a,ll b)
{
    ll res = 1;
    while(b)
    {
        if(b & 1)
            res = multi(res,a);
        a = multi(a,a);
        b >>= 1;
    }
    return res % mod;
}

int main()
{
    while(~scanf("%lld%lld",&a,&b))
    {
        top = 0;
        divide();
        ll sum = 1;
        for(int i = 1; i <= top; i ++)
        {
            if((p[i] - 1)% mod == 0)
            {
                sum = (sum % mod * ((b * c[i] + 1) % mod)) % mod;
                continue;
            }
            ll x = Quick_Pow(p[i],c[i] * b + 1);
            x = (x - 1 + mod) % mod;
            ll y = p[i] - 1;
            y = Quick_Pow(y,mod - 2);
            sum = ((sum % mod * (x % mod)) % mod * (y % mod)) % mod;
        }
        printf("%lld\n",sum);
    }
}

标签：p1,9901,res,ll,Sumdiv,逆元,poj1845,pm,mod
来源： https://blog.csdn.net/qq_43506138/article/details/99693320