poj 1845 所有因子和

传送门

题意：求A^B的所有因子和。

刚开始想到这就是一个等比数列的乘积嘛，不过还是要优化一下的。

一个数有唯一分解定理，可以分成若干质数相乘，

若对一个数n进行素数分解，n=p1^a1*p2^a2*p3^a3*...*pk^ak
那么n的所有正因子之和sum=(1+p1+...+p1^a1)*(1+p2+...+p2^a2)*...*(1+pk+...+pk^ak)

然后可以用等比数列求和公式(pk^(ak+1)-1)/(pk-1)求每项的和，再累乘。
用等比数列求1+pk+...+pk^ak时候要注意几点：

1.这里有除法，所以模的时候要将除以分母转化成乘以分母的逆元
a = (b/c) ==> a%m = b*c^(m-2)%m ( m为素数 )
证明：
b = a * c
根据费马小定理 a^(p-1)= 1 %p;(p是素数且a不能整除p)
所以 c^(m-1)%m=1%m
因此 a % m = a*1%m = a * c^(m-1)%m = a*c*c^(m-2)%m = b*c^(m-2)%m;

2.等比求和公式要注意，分母pk-1不能为0。
当pk%mod=1的时候， (1+pk+...+pk^ak)%mod=ak+1

3.当pk%mod=0的时候，(1+pk+...+pk^ak)=1，可以直接pass即可

#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <map>
#include <algorithm>
using namespace std;

#define pi acos(-1)
#define endl '\n'
#define srand() srand(time(0));
#define me(x,y) memset(x,y,sizeof(x));
#define foreach(it,a) for(__typeof((a).begin()) it=(a).begin();it!=(a).end();it++)
#define close() ios::sync_with_stdio(0); cin.tie(0);
#define FOR(x,n,i) for(int i=x;i<=n;i++)
#define FOr(x,n,i) for(int i=x;i<n;i++)
#define W while
#define sgn(x) ((x) < 0 ? -1 : (x) > 0)
#define bug printf("***********\n");
typedef long long LL;
const int INF=0x3f3f3f3f;
const LL LINF=1e18+7;
const int dx[]= {-1,0,1,0,1,-1,-1,1};
const int dy[]= {0,1,0,-1,-1,1,-1,1};
const int maxn=2005;
const int maxx=1e5+100;
const double EPS=1e-7;
//const int mod=998244353;
template<class T>inline T min(T a,T b,T c)
{
	return min(min(a,b),c);
}
template<class T>inline T max(T a,T b,T c)
{
	return max(max(a,b),c);
}
template<class T>inline T min(T a,T b,T c,T d)
{
	return min(min(a,b),min(c,d));
}
template<class T>inline T max(T a,T b,T c,T d)
{
	return max(max(a,b),max(c,d));
}
inline LL Scan()
{
	LL Res=0,ch,Flag=0;
	if((ch=getchar())=='-')Flag=1;
	else if(ch>='0' && ch<='9')Res=ch-'0';
	while((ch=getchar())>='0'&&ch<='9')Res=Res*10+ch-'0';
	return Flag ? -Res : Res;
}

LL a,b,mod=9901;
LL ans=1;
LL tot,pr[maxx],vis[maxx];  
void init()//素数表  O(nsqrt(n))  
{  
    for(int i=2;i<maxx;i++)  
    {  
        if(!vis[i])pr[tot++]=i;  
        for(int j=0;j<tot&&i*pr[j]<maxx;j++)  
        {  
            vis[i*pr[j]]=1;  
            if(i%pr[j]==0)break;  
        }  
    }  
}  
LL _pow(LL a,LL n)  
{  
    LL ret=1;  
    while(n)  
    {  
        if(n&1)  ret=ret*a%mod;  
        a=a*a%mod;  
        n>>=1;  
    }  
    return ret;  
}  
LL inv(LL x)  
{  
    return _pow(x,mod-2);  
}  
void solve()
{
	for(int i=0;i<tot;i++)
	{
		int c=0;
		if(a%pr[i]==0)
		{
			while(a%pr[i]==0)
			{
				c++; a/=pr[i];
			//	cout<<c<<" "<<pr[i]<<endl;
			}
			if(pr[i]%mod==0) continue;
			else if(pr[i]%mod==1) 
			{
				ans=(ans*(c*b+1))%mod;
			}
			else 
			{
				LL ret=(_pow(pr[i],(c*b+1))-1+mod)%mod;
				ans=(ans*ret*inv(pr[i]-1))%mod;
			}
		}
	}
	if(a>1)
	{
		if(a%mod==0)
			return ;
		else if(a%mod==1)
		{
			ans=(ans*(b+1))%mod;
		}
		else 
		{
			LL ret=(_pow(a,(b+1))-1+mod)%mod;
			ans=(ans*ret*inv(a-1))%mod;
		}
	}
}
int main()
{
	init();
	while(~scanf("%lld%lld",&a,&b))
	{
		ans=1;
		solve();
		cout<<ans<<endl;
	}
}