【poj1845】Sumdiv 递归+数学

这题很有意思。首先根据唯一分解定理我们可以得知a=p1^k1*p2^k2*k3......pn^kn，然后将所有的因数加起来再用一下乘法结合律就成了(p1^0+p1^1+p1^2......+p1^k1)......(pn^0+.......+pn^kn)这样一个式子，但是就成了求每一个的等比公式求和。但是不能直接使用公式，因为%并不适用于除法，所以：p^0+p^1+p^2......+p^k当k为奇数时{p^0+p^1+p^2+......+p^[(k-1)/2]}+p^[(k+1)/2](p^0+......p^[(k-1)/2])=(1+p^[(k+1)/2])*(p^0+......+p^(k-1)/2)最后%一下即可。偶数也是这样看程序吧。这样就能在log时间里求出等比公式之和。

#include<iostream>
#include<vector>
#include<cstdio>
#define ll long long
#define N 101000
#define mod 9901
using namespace std;
vector<pair<int,int> > vec;
ll ans=1;
int prime[N],is_prime[N],ToT,a,b;
void getprime() {
	prime[++ToT]=2;
	for (int i=3;i<N;i+=2) {
		if (!is_prime[i]) prime[++ToT]=i;
		for (int j=1;prime[j]*i<N;j++) {
			is_prime[prime[j]*i]=1;
			if (i%prime[j]==0) break; 
		}
	}
}
void divide(int x) {
	int QwQ;
	for (int i=1;i<=ToT;i++) {
		QwQ=0;
		for (;x%prime[i]==0;x/=prime[i]) QwQ++;
		if (QwQ) vec.push_back(make_pair(prime[i],QwQ));
	}
	if (x!=1) vec.push_back(make_pair(x,1));
}
ll q_pow(int x,int y) {
	ll ret=1,cur=x;
	for (;y;y>>=1) {
		if (y&1) ret=(ret*cur)%mod;
		cur=(cur*cur)%mod;
	}
	return ret;
}
ll solve(int x,int y) {
	if (!x) return 0;
	if (!y) return 1;
	if (y&1) return (1+q_pow(x,(y+1)/2))*solve(x,(y-1)/2)%mod;
	else return ((1+q_pow(x,y/2))*solve(x,y/2-1)%mod+q_pow(x,y))%mod;
}
int main() {
	getprime();
	scanf("%d%d",&a,&b);
	divide(a);
	for (int i=0;i<vec.size();i++) ans=(ans*solve(vec[i].first,b*vec[i].second))%mod;
	printf("%lld\n",ans);
	return 0;
}