POJ 1845 Sumdiv (因子和)

题目：http://poj.org/problem?id=1845

题意：给定A,B，S为A^B的因子和，求S%9901。

分析：这里给出三种方法。

方法一：首先将A进行素因子分解，S=p1^a1*p2^a2*p3^a3......pn^an,A^B=P1^(B*a1)*p2^(B*a2)+.....+pn^(B*an)由求因子和的公式http://blog.csdn.net/w20810/article/details/43702087因子和σ(n)=(1+p1+p1^2+p1^3+...+p1^（B*a1）)(1+p2+p2^3+p2^3+...p2^（B*a2）)...(1+pn+pn^2+pn^3+...+pn^（B*an）)。括号里面的为等比数列的和，http://blog.csdn.net/w20810/article/details/43458217二分就行了。

代码：

#include <iostream>
#include <cstdio>
#include <cstring>
#include <vector>

using namespace std;
typedef long long LL;
#define MOD 9901

LL nprime,prime[10000];
bool isprime[10000];

vector <LL > vt;

void doprime()
{
	nprime=0;
	LL i,j,k;
	for(i=1;i<10000;i+=2)
		isprime[i]=true;
	isprime[1]=false;
	isprime[2]=true;
	for(i=2;i<10000;i++)
	{
		if(isprime[i])
		{
			prime[nprime++]=i;
			for(j=i*i;j<7200;j+=i)
				isprime[j]=false;
		}
	}
}

void divd(LL n)
{
	LL t=n;
	vt.clear();
	for(LL i=0;i<nprime && prime[i]*prime[i]<=n;i++)
	{
		while(t%prime[i]==0)
		{
			vt.push_back(prime[i]);
			t/=prime[i];
		}
		if(t==1)
			break;
	}
	if(t>1)
		vt.push_back(t);
}

LL my_pow(LL a,LL n)
{
	LL ans=1;
	while(n)
	{
		if(n&1)
			ans=((ans%MOD)*(a%MOD))%MOD;
		n>>=1;
		a=((a%MOD)*(a%MOD))%MOD;
	 }
	 return ans;
}

LL bsum(LL q,LL n)
{
	if(n==1)
		return q%MOD;
	LL ans=bsum(q,n/2);
	if(n&1)
	{
		LL temp=my_pow(q,n/2+1);
		return ((1+temp)%MOD*(ans%MOD)%MOD+temp%MOD)%MOD;
	}
	else
	{
		LL temp=my_pow(q,n/2);
		return (1+temp)%MOD*(ans%MOD)%MOD;
	}
}

LL cnt(LL p)
{
	LL ans=0,i=0;
	while(i<vt.size() && vt[i]!=p) i++;
	for(;i<vt.size();i++)
		if(p==vt[i])
			ans++;
		else	
			break;
	return ans;
}

void solve(LL A,LL B)
{
	divd(A);
	LL i,j;
	LL p,count;
	LL ans=1;
	for(i=0;i<vt.size();i++)
	{
		p=vt[i];
		count=cnt(p);
		ans=ans*(bsum(p,count*B)+1)%MOD;
		i+=count-1;
	}
	printf("%lld\n",ans);
}

int main()
{
	doprime();
	LL A,B;
	while(scanf("%lld%lld",&A,&B)!=EOF)
	{
		if(B==0)
		{
			printf("1\n");
			continue;
		}
		solve(A,B);
	}
	return 0; 
} 

方法二：σ(n)=(1+p1+p1^2+p1^3+...+p1^（B*a1）)(1+p2+p2^3+p2^3+...p2^（B*a2）)...(1+pn+pn^2+pn^3+...+pn^（B*an）)。括号里面的可以用等比数列前n项和的公式求出。那么就是((p^(eB+1)-1)/(p-1))%MOD,这个式子可以转换成(p^(eB+1)-1)%(MOD*(p-1))/(p-1)%MOD。因为转换后的等式等价为[（p^(eB+1)-1）+K*MOD*(p-1)]/(p-1)%MOD={[(p^(eB+1)-1)]/(p-1)+K*MOD}%MOD=｛[(p^(eB+1)-1)]/(p-1)｝%MOD。

代码：

#include <iostream>
#include <cstdio>
#include <vector>
using namespace std;
typedef long long LL;
#define MOD 9901

int nprime,prime[10000];
bool isprime[10000];

vector <LL > vt; 

void doprime()
{
	nprime=0;
	LL i,j;
	for(i=1;i<10000;i+=2)
		isprime[i]=true;
	isprime[1]=false;
	isprime[2]=true;
	for(i=2;i<10000;i++)
	{
		if(isprime[i])
		{
			prime[nprime++]=i;
			for(j=i*i;j<10000;j+=i)
				isprime[j]=false;
		}
	}
}

void divd(LL n)
{
	LL i,j,t=n;
	vt.clear();
	for(i=0;i<nprime && prime[i]*prime[i]<=n;i++)
	{
		while(t%prime[i]==0)
		{
			vt.push_back(prime[i]);
			t/=prime[i];
		}
		if(t==1)
			break;
	}
	if(t>1)
		vt.push_back(t);
}

LL multi(LL a,LL b,LL m)
{
	LL ans=0;
	while(b)
	{
		if(b&1)
			ans=(ans+a)%m;
		b>>=1;
		a=(a<<1)%m;
	}
	return ans;
}

LL my_pow(LL a,LL n,LL m)
{
	LL ans=1;
	while(n)
	{
		if(n&1)
			ans=multi(ans,a,m);
		n>>=1;
		a=multi(a,a,m);
	}
	return ans;
}

LL cnt(LL p)
{
	LL ans=0,i=0;
	while(i<vt.size() && vt[i]!=p) i++;
	for(;i<vt.size();i++)
		if(p==vt[i])
			ans++;
		else	
			break;
	return ans;
}

void solve(LL A,LL B)
{
	divd(A);
	LL i,j,k,p,ans=1,count;
	for(i=0;i<vt.size();i++)
	{
		p=vt[i];
		count=cnt(p);
		ans=ans*(my_pow(p,count*B+1,MOD*(p-1))-1)/(p-1)%MOD;
		i+=count-1; 
	}
	printf("%lld\n",ans);
}


int main()
{
	doprime();
	LL A,B;
	while(scanf("%lld%lld",&A,&B)!=EOF)
	{
		solve(A,B);
	}
	return 0;
}

方法三： ((p^(eB+1)-1)/(p-1))%MOD，将除以p-1转换成乘以p-1的模的逆元。由于p-1并不一定与MOD互质，可以p%=MOD;