poj 1845 Sumdiv 数论--等比数列和(逆元或者递归)

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

题解：先分解a的质因子，a=p1^t1*p2^t2*...*pk^tk(pi为质数)。再a^b=p1^(t1*b)*p2^(t2*b)*...pk^(tk*b)。选出所有的因子就是枚举所有的ti*b，求和可知sum=(1+p1+...p1^(t1*b))*(1+p2+...p2^(t2*b))*...*(1+pk+...+pk*(tk*b));

而求1+pi+pi^2+...pi^ci(ci=ti*b)有两种方法：

1.直接用二分递归求：举例来说，1+a+a^2+a^3+a^4=(1+a)*(1+a^2)+a^2;(1+a)=1*(1+a)。根据奇偶二分下去。只要只有一个数为止。

2.是写出通项公式求解:这是个等比数列，所以由等比公式可得1+pi+...+pi^ci=(pi^(ci+1)-1)/(pi-1)。接着就是快速幂和分数形式的取模了。分数形式的取模也有两种方法：

（1）逆元：A/B=A*B^(-1)，B^(-1)就是B的逆元；B*C%mod=1，则说C是B的逆元，求的方法有拓展欧几里德公式来求和定理求解。不题mod=9901是质数，所以可以直接定理求解(费马小定理：a^(mod-1)%mod=1,mod为素数),也可以用欧拉定理，也一样。。因而(A/B)%mod=A*B^(mod-2);

（2）变换模值：(A/B)%mod=(A%(mod*B))/B%mod。对B*mod取余，剩余的值必定是B的倍数，这种方法是用于mod和B小的时候，这题就刚刚好。

注意点：

（1）0^0=1,这题中0的无论多少次都是1.。。，不知是没有这类的数据还是原本就定义成这样。

（2）用通项公式求的时候需要注意的是A/B中，如果B是mod的倍数就不能用逆元。。why？因为逆元要求B*B^(-1)%mod=1,但B%mod=0恒成立，找不出模mod下的逆元。

逆元求分数取模代码：

耗时：16MS

#include <cstdio>
#include <cmath>
#include <iostream>
#include <algorithm>
using namespace std;
typedef __int64 LL;
const LL mod=9901;
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;
}
LL pow_mod(LL a,LL b,LL n)//修改后的求次方，避免了爆int64
{
    a=a%n;
    LL s=1;
    while(b)
    {
        if(b&1)
        {
            s=mul(s,a,n);
        }
        a=mul(a,a,n);
        b=b>>1;
    }
    return s;
}
int main()
{
    LL a,b;
    while(cin>>a>>b)
    {
        if(a<=1||b==0){cout<<1<<endl;continue;}
        LL ans=1,i,j,k,t,n,m;
        n=(LL)sqrt(a+0.5);
        for(i=2;i<=n;i++)
        {
            if(a%i==0)
            {
                //if((i-1)%mod==0)cout<<"*"<<i<<endl;
                //cout<<"*"<<endl;
                t=0;
                while(a%i==0){
                    a=a/i;
                    t++;
                }
                if((i-1)%mod==0)ans=ans*(pow_mod(i,t*b+1,mod*(i-1))/(i-1))%mod;
                else ans=ans*(pow_mod(i,t*b+1,mod)-1)*pow_mod(i-1,mod-2,mod)%mod;
            }
        }
        if(a>1)
        {
            //if((a-1)%mod==0)cout<<"*"<<a<<endl;
            if((a-1)%mod==0)ans=ans*(pow_mod(a,b+1,mod*(a-1))/(a-1))%mod;
            else ans=ans*(pow_mod(a,b+1,mod)-1)*pow_mod(a-1,mod-2,mod)%mod;
        }
        cout<<(ans+mod)%mod<<endl;
    }
    return 0;
}

 

递归二分代码：

耗时:32MS

#include <cstdio>
#include <cmath>
#include <iostream>
#include <algorithm>
using namespace std;
const int mod=9901;
int pow_mod(int a,int b)
{
	a=a%mod;
	int s=1;
	while(b)
	{
		if(b&1)
			s=(s*a)%mod;
		a=(a*a)%mod;
		b=b>>1;
	}
	return s;
}
int sum(int a,int b)//求1+a+a^2+...+a^b
{
	if(b==1)return 1;
	if(b&1)return (sum(a,b/2)*(1+pow_mod(a,b/2+1))+pow_mod(a,b/2))%mod;
	else return sum(a,b/2)*(1+pow_mod(a,b/2))%mod;
}
int main()
{
	int a,b;
	while(cin>>a>>b)
	{
	 	if(a<=1||b==0){cout<<1<<endl;continue;}
		int ans=1,i,j,k,t,n,m;
		n=(int)sqrt(a+0.5);
		for(i=2;i<=n;i++)
		{
			if(a%i==0)
			{
				t=0;
				while(a%i==0){
					a=a/i;
					t++;
				}
				ans=ans*sum(i,t*b+1)%mod;
			}
		}
		if(a>1)
			ans=ans*sum(a,b+1)%mod;
		cout<<(ans+mod)%mod<<endl;	
	}
	return 0;
}