Number theory to find the sum of factors of a ^ B (A and B are 1e12 levels)

1. First of all, consider a Quality factor decomposition , Screen out 1e6 Internal a The prime factor and power of , if a And the quality factor （ Not screened into 1 Words ） There is one greater than 1e6 Large qualitative factor of （ At most one ）（ Fill in and find all the prime factors ）

2. take b Put in the power exponent

3. Consider a generating function

a^b The sum of the factors =(p1^0+p1^1.....+p1^k1)(p2^0+.....+p2^k2).....(pm^0+...pm^km)

Each term of this generating function is a different factor （ here p by a The quality factor of ）

4. The last is to deal with the sum of a series of equal ratio sequences , The first method is equal ratio summation , But because of mod Too small Possible q by 1 The situation of , Now I want to go back to summation （ individual ）, The rest are normal

The second method is direct recursion log Find the continuous sum of the sequence of equal ratios （ The template is as follows ）

（ Be careful qmul Prevent multiplication overflow （ Greater than 1e9 Easy to ））

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxx=1000006;
const int mod=9901;
int pri[maxx];
bool ispri[maxx];
int cnt;
long long cc[maxx];
void prime()
{
	cnt=1;
	memset(ispri,1,sizeof(ispri));
	ispri[0]=ispri[1]=0;
	for(int i=2;i<=maxx;++i)
	{
		if(ispri[i])
		pri[cnt++]=i;
		for(int j=1;j<cnt&&pri[j]*i<maxx;++j)
		{
			ispri[pri[j]*i]=0;
			if(i%pri[j]==0)
			break;
		}
	}
}
ll qmul(ll a,ll b)
{
	ll ans=0;
	while(b)
	{
		if(b&1)
		ans=(ans+a)%mod;
		a=(a+a)%mod;
		b>>=1;
	}
	return ans;
}
ll qpow(ll a,ll b)
{
	ll ans=1;
	while(b)
	{
		if(b&1)
		ans=qmul(ans,a);
		a=qmul(a,a);
		b>>=1;
	}
	return ans;
}
ll m(ll x)
{
	return qpow(x,mod-2);
}
ll work(ll a,ll b)
{
	if(b==0ll) return 1ll;
	if(b==1ll) return (1+a)%mod;
	if(b%2==0) return ((1ll+qpow(a,b/2))%mod*work(a,(b-1)/2)%mod*a%mod+1ll)%mod;
	return (1ll+qpow(a,(b+1ll)/2))%mod*work(a,b/2)%mod;
}
int main ()
{
	freopen("spring.in","r",stdin);
	freopen("spring.out","w",stdout);
	ll a,b;
	cin>>a>>b;
	prime();
	--cnt;
	long long k;
	int maxx;
	for(int i=1;i<=cnt;++i)
	{
		if(a%pri[i]==0)
		{
			maxx=i;
			cc[i]=1;
			a/=pri[i];
		    while(1)
		  {
			if(a%pri[i]==0)
			{
				cc[i]++;
			a/=pri[i];
			}
			else break;
		  }
		}
		cc[i]*=b;
	}
	ll ans=1;
	for(int i=1;i<=maxx;++i)
	{
		if(cc[i])
		ans=(ans*work(pri[i],cc[i]))%mod;
	}
	if(a>1)
	ans=(ans*work(a,b))%mod;
	cout<<ans;
	return 0;
}