POJ-1845-Sumdiv 等比数列求和/数学/（二分法/逆元法/变换取模法）

题意：输入a，b，求a^b的所有因子之和

http://poj.org/problem?id=1845

分解a的质因数a=p1^t1*p2^t1........

每个质因数对sum的贡献： 当除去质因数p1时的因数和为sum，当计入p1时，因子和变成sum*p1^0+sum*p1^1+sum*p1^2......+sum*p1^t1

也就是所有的sum=【1+p1+p1^2+p1^3+...+p1^t1】*【p2.....】【p3...】

然后由于是a^b,所以最后是

sum=sum=【1+p1+p1^2+p1^3+...+p1^（t1*b）】*【p2.....】【p3...】

显然就是求关于a的所有质因数的一个 等比数列之和前n项和.

求这个等比数列有几种方法：

1、直接二分递归求：  

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

当n=奇数，会有偶数项，我们把1和p^(n/2+1)放在一起 p和p^(n/2+1+1)放在一起，以此类推

最后提公因式，得到【1+p^(n/2+1)】*【1+p+p^2+p^3+....+p^(n/2-1)】右边恰为原式的一半

当n=偶数，同理得到【1+p^(n/2)】*【1+p+p^2+....p^(n/2)】

因此二分递归则可以得到答案

2、逆元，用公比求和公式  1+pi+...+pi^n=(pi^(n+1)-1)/(pi-1)

  由于涉及到除法，且mod=9901为素数，所以可以用费马小定理求逆元,只是要注意mod比较小，

当【prim[i]-1】%mod==0（分母是mod的倍数）时，逆元不存在，不过此时恰好公比为1啦，答案就是Sn=n

3、变换取模，(A/B)%mod=(A%(mod*B))/B%mod,此题的数有点大.5e7....所以最大的质因子可能是 1e7左右。然后 本公式中MOD=mod*b≈MOD*1e7=1e10>1e9  那么计算过程中，两个1e10相乘就溢出了。。。这里错了好久才找到。。如果一定要用可以用高精度模拟一下这里的乘法

1二分法：

#include <cstdio>
#include <cmath>
#include <cstring>
#include <string>
#include <algorithm>
#include <iostream>
#include <queue>
#include <map>
#include <set>
#include <vector>
#include<stack>
using namespace std; 
int ok;
__int64 times[105]; 
int prim[105];
void ff(int x )
{
	for (int i=2;i*i<=x;i++)
	{
		if (x%i==0)
		{
			prim[++ok]=i;
			while(x%i==0) 
			{
				x=x/i;
				times[ok]++;
			}
		}
	}
	if (x!=1) 
	{
		prim[++ok]=x;
		times[ok]=1;
	}
	
}
typedef __int64 LL;  
 
__int64 pow_m(__int64 a,__int64 b,__int64 c)
{
	__int64 ans=1; 
	__int64 tmp=a%c;
	
	while(b!=0)
	{ 	
		if (b&1) 
		ans=ans*tmp%c;    
 		tmp=tmp*tmp%c; 
		b=b>>1;
	} 
	return ans;	 
}
__int64 mod=9901;
__int64 sum(__int64 p,__int64 n)  //递归二分求 (1 + p + p^2 + p^3 +...+ p^n)%mod  
{                          //奇数二分式 (1 + p + p^2 +...+ p^(n/2)) * (1 + p^(n/2+1))  
	if(n==0)               //偶数二分式 (1 + p + p^2 +...+ p^(n/2-1)) * (1+p^(n/2+1)) + p^(n/2)  
		return 1;  
	if(n%2)  //n为奇数,  
		return (sum(p,n/2)*(1+pow_m(p,n/2+1,mod)))%mod;  
	else     //n为偶数  
		return (sum(p,n/2-1)*(1+pow_m(p,n/2+1,mod))+pow_m(p,n/2,mod))%mod;  
}  
int main()
{ 
	__int64 m,n;
	memset(prim,0,sizeof(prim));
	scanf("%I64d%I64d",&n,&m);
	__int64 ans=1;
	if (n<=1||!m)
	{
		printf("1\n");return 0;
	} 
	ok=0;
	ff(n);
	int i;
	for (i=1;i<=ok;i++)
	{ 
		//求 (1 + p + p^2 + p^3 +...+ p^（t[i]*m）)%mod  
		 
		ans=ans*sum(prim[i],times[i]*m)%mod;  

		}
	while(ans<0) ans+=mod;
	printf("%I64d\n",ans%mod);
	return 0;
} 

2 逆元法

#include <cstdio>
#include <cmath>
#include <cstring>
#include <string>
#include <algorithm>
#include <iostream>
#include <queue>
#include <map>
#include <set>
#include <vector>
#include<stack>
using namespace std; 
int ok;
__int64 times[105]; 
int prim[105];
void ff(int x )
{
	for (int i=2;i*i<=x;i++)
	{
		if (x%i==0)
		{
			prim[++ok]=i;
			while(x%i==0) 
			{
				x=x/i;
				times[ok]++;
			}
		}
	}
	if (x!=1) 
	{
		prim[++ok]=x;
		times[ok]=1;
	}
	
}
typedef __int64 LL;  
 
__int64 pow_m(__int64 a,__int64 b,__int64 c)
{
	__int64 ans=1; 
	__int64 tmp=a%c;
	
	while(b!=0)
	{ 	
		if (b&1) 
		ans=ans*tmp%c;     

 	tmp=tmp*tmp%c; 
		b=b>>1;
	} 
	return ans;	 
}
__int64 mod=9901; 
int main()
{ 
	__int64 m,n;
	memset(prim,0,sizeof(prim));
	scanf("%I64d%I64d",&n,&m);
	__int64 ans=1;
	if (n<=1||!m)
	{
		printf("1\n");return 0;
	} 
	ok=0;
	ff(n);	//因数分解
	int i;
	for (i=1;i<=ok;i++)
	{ 
		//求 (1 + p + p^2 + p^3 +...+ p^（t[i]*m）)%mod   
		if((prim[i]-1)%mod==0)ans=ans*(times[i]*m+1)%mod;	//公比为1的情况
        else ans=ans*(pow_m(prim[i],times[i]*m+1,mod)-1)*pow_m(prim[i]-1,mod-2,mod)%mod;  //逆元
		}
	while(ans<0) ans+=mod;
	printf("%I64d\n",ans%mod);
	return 0;
} 

3、变换取模法：

#include <cstdio>
#include <cmath>
#include <cstring>
#include <string>
#include <algorithm>
#include <iostream>
#include <queue>
#include <map>
#include <set>
#include <vector>
#include<stack>
using namespace std; 
int ok;
__int64 times[105]; 
int prim[105];
void ff(int x )
{
	for (int i=2;i*i<=x;i++)
	{
		if (x%i==0)
		{
			prim[++ok]=i;
			while(x%i==0) 
			{
				x=x/i;
				times[ok]++;
			}
		}
	}
	if (x!=1) 
	{
		prim[++ok]=x;
		times[ok]=1;
	}
	
}
typedef __int64 LL;  

LL mul(LL a,LL b,LL n)//大数乘法，直接相乘会爆int64，需要逐位相乘  
{  
    LL s=0;  
    while(b)  
    {  
        if(b&1)  
            s=(s+a)%n;  
        a=(a*2)%n;  
        b=b>>1;  
    }  
    return s;  
} 
__int64 pow_m(__int64 a,__int64 b,__int64 c)
{
	__int64 ans=1; 
	__int64 tmp=a%c;
	
	while(b!=0)
	{ 	
		if (b&1) 
	 	ans=mul(ans,tmp,c); 

tmp=mul(tmp,tmp,c);
		b=b>>1;
	} 
	return ans;	 
}
__int64 mod=9901; 
int main()
{ 
	__int64 m,n;
	memset(prim,0,sizeof(prim));
	scanf("%I64d%I64d",&n,&m);
	__int64 ans=1;
	if (n<=1||!m)
	{
		printf("1\n");return 0;
	} 
	ok=0;
	ff(n);
	int i;
	for (i=1;i<=ok;i++)
	{ 
		//求 (1 + p + p^2 + p^3 +...+ p^（t[i]*m）)%mod   
			__int64 tmp=pow_m(prim[i],times[i]*m+1,mod*(prim[i]-1))-1;//快速幂过程可能溢出
		tmp=tmp/(prim[i]-1)%mod;
		ans=ans*tmp%mod;   
		}
	while(ans<0) ans+=mod;
	printf("%I64d\n",ans%mod);
	return 0;
}