先說題意：輸入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;
}