Sumdiv POJ

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.OutputThe 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).

题意：求A^B的所有约数的和  % MOD （9901）    题意中有点问题，我们知道0是没有约数的，我觉得A、B应该都是>0的

思路：我们可以把A分解质因数（p₁^c1  *  p₂^c2  *  .... * p_n^cn）^B

约数和：（1 + p₁ + p₁² + ... + p₁^B*c1）* （1 + p₂ + p₂² + ... + p₂^B*c2）* .... * （1 + p_n + p_n² + ... + p_n^B*cn）    （ 排列组合问题）

这样我们可以看出这是多个等比数列乘积，可以用等比数列求和公式  （a₁ *（1-qⁿ））/（1-q）,我们注意到这里有除法，但是同余模定理是对于加减乘的，那么我们可以利用费马小定理，

求出（1-q）的逆元，然后把除变成乘逆元

坑点：应为 9901 这个质数较小，很容易找到一个数x，（x-1）% MOD == 0 ,就说明这个数是没有逆元的（例217823），那么对于这种情况，我们不能用逆元算，你会发现这种情况下，

pⁿ % MOD == 1  （（p-1）% MOD == 0）），这样（1 + p_n + p_n² + ... + p_n^B*cn） == （1 % MOD + p_n %MOD + p_n² %MOD + ... + p_n^B*cn %MOD） == B*cn+1

#include<iostream>
#include<cstdio>
#include<math.h>
using namespace std;

const int maxn = 1e4;
const int mod = 9901;
int a,b;
int p[maxn];
int c[maxn];
int calc(int x)
{
    int m= 0;
    int up = sqrt(x);
    for(int i=2;i<=up;i++)
    {
        if(x % i == 0)p[++m] = i,c[m] = 0;
        while(x % i == 0)x/= i,c[m]++;
    }
    if(x > 1)p[++m] = x,c[m] = 1;
    return m;
}

typedef long long ll;

ll qpow(ll a,ll b)
{
    ll ans = 1;
    ll base = a;
    while(b)
    {
        if(b&1)ans = (ans * base)%mod;
        base = (base * base)%mod;
        b >>= 1;
    }
    return ans;
}
int main()
{
    scanf("%d%d",&a,&b);
    int n = calc(a);
    ll ans = 1;
    for(int i=1;i<=n;i++)
    {
        if((1-p[i])%mod == 0)
        {
            ans = (ans * (b * c[i]+ 1))%mod;
            continue;
        }
        ll Ni = qpow(1-p[i],mod-2);
        ll tmp = 1-qpow(p[i],c[i]*b+1);
        ans = (ans * (tmp*Ni%mod+mod)%mod)%mod;
    }
    printf("%lld
",ans);
}

还有一种写法，就是不用公式计算等比数列和，这样就避免了逆元的问题

sum（p，c） = （1 + p + p² + ... + p^k）

（1）c为奇数，sum（p，k）= sum（p，（k-1 ）/2）*（1+p^（k+1）/2）   

sum（p，c） = （1 + p + p² + ... + p^（k-1）/2）+ （p^（k+1）/2 + ... + p^k）           （c为奇数，加上0次幂，变成偶数，刚好可以分成两个等长的数列）

（2）c为偶数，sum（p，k）= sum（p，k/2-1）*（1+p^k/2）+ p^k 

sum（p，c） = （1 + p + p² + ... + p^k/2-1）+ （p^k/2 + ... + p^k-1）+ p^k                         （c+1是奇数）

#include<iostream>
#include<cstdio>
#include<math.h>
using namespace std;

const int maxn = 1e4;
const int mod = 9901;
int a,b;
int p[maxn];
int c[maxn];
int calc(int x)
{
    int m= 0;
    for(int i=2;i*i<=x;i++)
    {
        if(x % i == 0)p[++m] = i,c[m] = 0;
        while(x % i == 0)x/= i,c[m]++;
    }
    if(x > 1)p[++m] = x,c[m] = 1;
    return m;
}

typedef long long ll;
ll qpow(ll a,ll b)
{
    ll ans = 1;
    ll base = a;
    while(b)
    {
        if(b&1)ans = (ans * base)%mod;
        base = (base * base)%mod;
        b >>= 1;
    }
    return ans;
}
ll sum(ll p,ll c)
{
    if(c == 0)return 1;
    if(c&1)return ((1+qpow(p,(c+1)/2))%mod*(sum(p,(c-1)/2)%mod))%mod;
    else return ((1+qpow(p,c/2))%mod*(sum(p,c/2-1))%mod+qpow(p,c))%mod;
}

int main()
{
    scanf("%d%d",&a,&b);
    int n = calc(a);
    ll ans = 1;
    for(int i=1;i<=n;i++)
    {
        ans = (ans * sum(p[i],c[i]*b))%mod;
    }
    printf("%lld
",ans);
}