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

╰半夏微凉° 2021-10-29 23:10 235阅读 0赞

[![知识共享许可协议][80x15.png]][80x15.png 1]  
本作品采用[知识共享署名-相同方式共享 4.0 国际许可协议][80x15.png 1]进行许可。  
题目链接：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

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 ^ (c1*b)) \* …… \* (1 + pm + pm ^ 2 + …… + pm ^ (cm*b))

所以我们先把a的所有质因子分解出来。  
每一个括号都是一个等比数列，自然我们可以用等比数列来求值，不过没有这么简单。  
等比数列求和公式：  
(p1^(b*c1+1) - 1) / (p1- 1)  
对于p1^(b*c1+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);
        }
    }

[80x15.png]: https://i.creativecommons.org/l/by-sa/4.0/80x15.png
[80x15.png 1]: http://creativecommons.org/licenses/by-sa/4.0/