poj1845 数论 快速幂

Sumdiv

Time Limit: 1000MS		Memory Limit: 30000K
Total Submissions: 16466		Accepted: 4101

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

要求的是A^B的所有因子的和之后再mod 9901的值。

(1+a1+a1^2+...a1^n1)*(1+a2+a2^2+...a2^n2)*(1+a3+a3^2+...a3^n2)*...(1+am+am^2+...am^nm) mod 9901。

对于每一个(1+a1+a1^2+...a1^n1) mod 9901 

等于 (a1^(n1+1)-1)/(a1-1) mod 9901，这里用到逆元的知识：a/b mod c = (a mod (b*c))/ b 

所以就等于(a1^(n1+1)-1)mod (9901*(a1-1)) / (a1-1)。

至于前面的a1^(n1+1)，快速幂。

#include<cstdio>
#include<cmath>
#include<iostream>
#include<algorithm>
#include<cstring>
#define mod 9901
#define N 10007
#define ll long long
using namespace std;

int prime[10007];

void getPrime()
{
    for(int i=2;i<=10000;i++)
    {
        if(!prime[i])prime[++prime[0]]=i;
        for(int j=1;j<=prime[0]&&prime[j]*i<=10000;j++)
        {
            prime[prime[j]*i]=1;
            if(i%prime[j]==0) break;
        }
    }
}

long long factor[100][2];
int fatCnt;
int getFactors(long long x)
{
    fatCnt=0;
    long long tmp=x;
    for(int i=1;prime[i]<=tmp/prime[i];i++)
    {
        factor[fatCnt][1]=0;
        if(tmp%prime[i]==0)
        {
            factor[fatCnt][0]=prime[i];
            while(tmp%prime[i]==0)
            {
                factor[fatCnt][1]++;
                tmp/=prime[i];
            }
            fatCnt++;
        }
    }
    if(tmp!=1)
    {
        factor[fatCnt][0]=tmp;
        factor[fatCnt++][1]=1;
    }
    return fatCnt;
}
long long pow_m(long long a,long long n)//快速模幂运算
{
    ll ans=1;a%=mod;
    while(n)
    {
        if (n&1) ans=(ans*a)%mod;
        n>>=1;
        a=(a*a)%mod;
    }
    return ans;
}
long long sum(long long p,long long n)//计算1+p+p^2+````+p^n
{
    if(p==0)return 0;
    if(n==0)return 1;
    if(n&1) return ((1+pow_m(p,n/2+1))%mod*sum(p,n/2)%mod)%mod;
    else return ((1+pow_m(p,n/2+1))%mod*sum(p,n/2-1)+pow_m(p,n/2)%mod)%mod;

}
int main()
{
    int A,B;
    getPrime();
    while(~scanf("%d%d",&A,&B))
    {
        getFactors(A);
        long long ans=1;
        for(int i=0;i<fatCnt;i++)
            ans=(ans*sum(factor[i][0],B*factor[i][1])%mod)%mod;
        printf("%lld
",ans);
    }
}