poj 1845 Sumdiv(数论:欧拉函数+二分求等比数列前n项和+快速幂取模)

很凶残的一道题啊...

给定a,b,求a^b的所有因子之和

给定一个a先用欧拉函数求出唯一分解式f(a) = p1^e1*p2^e2*...*pk^ek

则变为求f(a)^b所有因子之和

可知所有因子之和为(1+p1+p1^2+...+p1^e1)*...*(1+pk+pk^2+...+pk^ek)

接着对每个等比序列用二分法求和

二分原理如下:

（1）若n为奇数,一共有偶数项，则：

1 + p + p^2 + p^3 +...+ p^n

= (1+p^(n/2+1)) + p * (1+p^(n/2+1)) +...+ p^(n/2) * (1+p^(n/2+1))

= (1 + p + p^2 +...+ p^(n/2)) * (1 + p^(n/2+1))

上式红色加粗的前半部分恰好就是原式的一半，那么只需要不断递归二分求和就可以了，后半部分为幂次式，将在下面第4点讲述计算方法。

（2）若n为偶数,一共有奇数项,则:

1 + p + p^2 + p^3 +...+ p^n

= (1+p^(n/2+1)) + p * (1+p^(n/2+1)) +...+ p^(n/2-1) * (1+p^(n/2+1)) + p^(n/2)

= (1 + p + p^2 +...+ p^(n/2-1)) * (1+p^(n/2+1)) + p^(n/2);

上式红色加粗的前半部分恰好就是原式的一半，依然递归求解

上面形式是经过把序列从中间分为两半,然后这两半从头到尾对应处理
如果看不懂上面的论述,可以自己用2^3和2^4带入试一下,这里就不再详细证明了

对于上面等式后半段的大数幂可以用快速幂取模来做

代码如下:

#include <cmath>
#include <cstdio>
#include <iostream>
#include <algorithm>
#define MAXN 25000
#define M 9901
#define LL long long
using namespace std;

struct Num {
LL p, e;
} a[MAXN];

LL euler_phi(LL n) {
LL m = (LL)sqrt(n+0.5);
LL cnt = 0, num;
for(int i=2; i<=m; ++i) {
if(n%i == 0) {
a[cnt].p = i;
num = 1;
n /= i;
while(n%i == 0) {
num++;
n /= i;
}
a[cnt++].e = num;
}
}
if(n > 1) {
a[cnt].p = n;
a[cnt++].e = 1;
}
return cnt;
}

LL Pow(LL a, LL b) {
LL ans = 1;
while(b) {
if(b & 1)
ans = ans*a%M;
a = a*a%M;
b >>= 1;
}
return ans%M;
}

LL get_sum(LL p, LL e) {
if(e == 0) return 1;
if(e & 1) {
return (get_sum(p, e/2)%M)*((1+Pow(p, e/2+1))%M)%M;
} else {
return ((get_sum(p, e/2-1)%M)*(1+(Pow(p, e/2+1))%M)+Pow(p, e/2)%M)%M;
}
}

int main(void) {
int n, m, cnt, ans;
cin >> n >> m;
ans = 1;
//n或m等于0的情况不要忘了考虑
if(!n) {
cout << 0 << endl;
return 0;
}
if(!m) {
cout << 1 << endl;
return 0;
}

cnt = euler_phi(n);
for(int i=0; i<cnt; ++i) {
//        cout << a[i].p << "\t" << a[i].e << endl;
ans = ans*get_sum(a[i].p, m*a[i].e)%M;
}
cout << (ans%M+M)%M << endl;
return 0;
}