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;
}
}