POJ 1845 Sumdiv 求幂级数的因子和+二分

点击打开链接

Sumdiv

Time Limit: 1000MS		Memory Limit: 30000K
Total Submissions: 13055		Accepted: 3183

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

让你求A^B的所有约数之和。
可以对A进行素因子分解：A=p1^k1*p2^k2*...*pn^kn.
那么A^B=p1^(k1*B)*p2^(k2*B)*...*pn^(kn*B)
利用乘性函数的性质，所有的约数和就是：S=(1+p1+p1^2+...+p1^(k1*B))*(1+p2+p2^2+...+p2^(k2*B))*...*(1+pn+pn^2+...+pn^(kn*B))
然后分别求出每个等比数列的和，求等比数列的和可以用递归二分的形式。
如果n为奇数，那么就有偶数项：

1 + p + p^2 + p^3 +...+ p^n

 = (1+p^(n/2+1)) + p * (1+p^(n/2+1)) +...+ p^(n/2) * (1+p^(n/2+1))

 = (1 + p + p^2 +...+ p^(n/2)) * (1 + p^(n/2+1))

如：1 + p + p^2 + p^3 + p^4 + p^5 =
(1 + p + p^2) * (1 + p^3)

如果n为偶数，那么就有奇数项：

1 + p + p^2 + p^3 +...+ p^n

 = (1+p^(n/2+1)) + p * (1+p^(n/2+1)) +...+ p^(n/2-1) * (1+p^(n/2+1)) + p^(n/2)
 = (1 + p + p^2 +...+ p^(n/2-1)) * (1+p^(n/2+1)) + p^(n/2);

如：1 + p + p^2 + p^3 + p^4 = (1
+ p) * (1 + p^3) + p^2

//396K	0MS
#include<stdio.h>
#include<string.h>
#include<math.h>
#define LL long long
#define N 50010
#define M 9901
LL s
[2],len;
long long multi(long long a,long long b,long long m)//a*b%m
{
long long ret=0;
while(b>0)
{
if(b&1)ret=(ret+a)%m;
b>>=1;
a=(a<<1)%m;
}
return ret;
}
long long quick_mod(long long a,long long b,long long m)//a^b%m
{
long long ans=1;
a%=m;
while(b)
{
if(b&1)
{
ans=multi(ans,a,m);
b--;
}
b/=2;
a=multi(a,a,m);
}
return ans;
}

long long get(long long n)  //S[i][0]代表第i个素数,S[i][1]代表第i个素数的个数
{
len=0;
for(long long i=2;i*i<=n;i++)
{
if(n%i==0)
{
s[len][0]=i;s[len][1]=0;
do{n/=i;s[len][1]++;}while(n%i==0);
len++;
}
}
if(n>1){s[len][0]=n;s[len++][1]=1;}
}
LL solve(LL p,LL a)
{
LL count=0;
if(!a)return 1;
if(a&1)
{
LL ans=solve(p,a/2)%M;
count=(ans+quick_mod(p,a/2+1,M)*ans%M)%M;
return count;
}
else
{
LL ans=solve(p,a/2-1)%M;
count=(ans+quick_mod(p,a/2,M)*(1+(p*ans)%M))%M;
return count;
}
return count;
}
int main()
{
LL a,b,ans;
while(scanf("%I64d%I64d",&a,&b)!=EOF)
{
ans=1;
get(a);
for(LL i=0;i<len;i++)
ans=(ans*solve(s[i][0],s[i][1]*b)%M)%M;
printf("%I64d\n",ans);
}
return 0;
}