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