POJ 1845（Sumdiv）

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

题意：给出A和B，求A^B这个数所有约数的和。

思路：设A^B=(p1^k1)*...(pn^kn)，那么约数和就是ans=(1+p1+p1^2+...p1^k1)*....(1+pn+pn^2+...pn^kn)。然后对于每一个素因子用一个等比数列求和公式算出来就行了，但是因为涉及到除法取模，直接求逆元wa了，毕竟费马小定理要求 a^(p-1)≡1(mod p)中a和p互素，且p是素数。所以用之前见过的另一个知识点。

(A/B)%mod=(A%(mod*B))/B%mod。对B*mod取余，剩余A的值必定是B的倍数。

等比求和公式应该就不用证明A是B的倍数了吧。。剩下就直接代代进去就行了。

#include <set>
#include <map>
#include <deque>
#include <stack>
#include <queue>
#include <time.h>
#include <vector>
#include <string>
#include <math.h>
#include <cstring>
#include <cstdlib>
#include <stdio.h>
#include <iomanip>
#include <iostream>
#include <algorithm>
#define PI acos(-1)
#define ll long long
#define LL long long
#define inf 0x3f3f3f3f
#define ull unsigned long long
using namespace std;

const int MAXN = 50005;
bool flag[MAXN];
int primes[MAXN/3],pi;
void Prime()
{
    int i,j;
    pi=0;
    memset(flag,false,sizeof(flag));
    for(i=2;i<MAXN;i++)
    {
        if(!flag[i]) primes[++pi] = i;
        for(j=1;(j<=pi)&&(i*primes[j]<MAXN);j++)
        {
            flag[i*primes[j]]=true;
            if (i%primes[j]==0)//这句保证每个非素数只被筛去一次、、
                break;
        }
    }
}
long long mul_mod(long long a,long long b,long long mod)
{
    a=a%mod;
    long long res=0;
    while(b)
    {
        if(b&1)
        {
            res+=a;
            res%=mod;
        }
        a+=a;
        a%=mod;
        b>>=1;
    }//
    return res;
}
long long qpow(long long a,long long b,long long mod)
{
    a=a%mod;
    long long res=1;
    while(b)
    {
        if(b&1)
        {
            res=mul_mod(res,a,mod);
        }
        a=mul_mod(a,a,mod);
        b>>= 1;
    }
    return res;
}
void solve(ll a,ll b)
{
    ll ans=1;
    ll cnt;
    for(ll i=1;i<=pi&&primes[i]*primes[i]<=a;i++)
    {
        if(a%primes[i]==0)
        {
            cnt=0;
            while(a%primes[i]==0)
            {
                cnt++;
                a=a/primes[i];
            }
            long long m = (primes[i] - 1) * 9901;
			ans*=(qpow(primes[i],cnt*b+1,m)+m-1)/(primes[i]-1);
			ans%= 9901;
        }
    }
    if(a>1)
    {
        long long m = 9901 * (a - 1);
		ans*=(qpow(a,b+1,m)+m-1)/(a-1);
		ans%= 9901;
    }
    printf("%lld\n",ans);
}
int main()
{
    Prime();
    ll a,b;
    while(scanf("%lld%lld",&a,&b)!=EOF)
    {
        if(a<=1||b==0)
        {
            cout<<1<<endl;
            continue;
        }
        solve(a,b);
    }
}