POJ 1845

此題需要注意的一個細節時，若MOD|P或MOD|(P-1)，此時不能應用費馬小定理求逆元的方法。

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
#define LL __int64
using namespace std;

const LL MOD=9901;

using namespace std;
LL prime[7100],np;
bool isprime[7100];

void prim(){
	memset(isprime,true,sizeof(isprime));
	np=0;
	isprime[1]=false;
	for(LL i=2;i<(LL)7100;i++){
		if(isprime[i]){
			prime[np++]=i;
			for(LL j=i*i;j<(LL)7100;j+=i)
			isprime[j]=false;
		}
	}
}

LL quick(LL a,LL b,LL m){
	LL ans=1;
	a%=m;
	while(b){
		if(b&1){
			ans=(ans*a)%m;
		//	cout<<ans<<endl;
		}
		b>>=1;
		a=(a*a)%m;
	}
	return ans;
}

int main(){
	prim();
	LL a,b;
	while(scanf("%I64d%I64d",&a,&b)!=EOF){
		LL ans=1; LL cnum;
		for(int i=0;i<np&&prime[i]<=a;i++){
			cnum=0;
			if(a%prime[i]==0){
				if(prime[i]%MOD==0) continue;
				else if((prime[i]-1)%MOD==0){
					while(a%prime[i]==0){
						cnum++;
						a/=prime[i];
					}
					cnum=cnum*b;
					ans=(ans*((cnum+1)%MOD))%MOD;
				}
				else{
					while(a%prime[i]==0){
						cnum++;
						a/=prime[i];
					}
					cnum=cnum*b+1;
				//	cout<<cnum<<endl;
					LL p=quick(prime[i],cnum,MOD);
					p=((p-1)%MOD+MOD)%MOD;
					LL q=quick(prime[i]-1,MOD-2,MOD);
					ans=(ans*((p*q)%MOD))%MOD;
				}
			}
		}
	//	cout<<ans<<endl;
		if(a>1){
			if(a%MOD==0);
			else if((a-1)%MOD==0){
				cnum=1;
				cnum=cnum*b;
				ans=(ans*((cnum+1)%MOD))%MOD;
			}
			else{
				cnum=1;
				cnum=cnum*b+1;
				LL p=quick(a,cnum,MOD);
				p=((p-1)%MOD+MOD)%MOD;
				LL q=quick(a-1,MOD-2,MOD);
				ans=(ans*((p*q)%MOD))%MOD;
			}
		}
		printf("%I64d\n",ans);
	}
	return 0;
}