[数论+二分求等比数列]POJ 1845 Sumdiv

题目链接 ：http://poj.org/problem?id=1845

题目大意 ：求A^B 所有因子的和% 9901

思路 ：A = p1^a1 * p2^a2 * ... ... * pn^a3 那么A的所有因子之和为 sigma ( pi^0 + pi^1 + pi^2 + ... ... + pi^ai )

A^B 的所有因子之和为 sigma ( pi^0 + pi^1 + pi^2 + ... ... + pi^（ai*B ))

对于 ( pi^0 + pi^1 + pi^2 + ... ... + pi^（ai*B )) 这是一个等比数列，可以用二分求解，时间复杂度为logN

Code：

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cmath>
using namespace std;
#define foru(i, a, b) for (int i=a; i<=b; i++)
#define ford(i, a, b) for (int i=a; i>=b; i--)
#define ll __int64
#define M 9901
#define N 10001

bool f
;
int pri
, m = 0 ;
void getpri(){
foru(i, 2, N-1){
if (! f[i]) {m++; pri[m] = i;}
foru(j, 1, m){
if (i * pri[j] > N-1) break;
f[i * pri[j]] = 1;
if (i % pri[j] == 0) break;
}
}
}

ll g(ll a, ll b){
ll s = 1;
while (b){
if ((b&1) == 1) s = (s*a)%M;
b >>= 1;
a = (a*a)%M;
}
return s;
}

ll x, y;
ll getsum(ll p, ll n){
if (! n) return 1ll;
if (n&1) return (getsum(p, n/2) * (1 + g(p, n/2+1)))%M;
return (getsum(p, n/2-1) * (1 + g(p, n/2+1)) + g(p, n/2))%M;
}

ll get(){
ll sum = 1ll;
foru(i, 1, m){
if (x == 1) break;
ll t = 0ll;
while (x%pri[i] == 0){
t ++;
x /= pri[i];
if (x == 1) break;
}
if (t == 0) continue;
sum = (sum * getsum(pri[i], t*y))%M;
}
if (x > 1) sum = (sum * getsum(x, y))%M;
while (sum < 0) sum += M;
return sum;
}

int main(){
//freopen("K.txt", "r", stdin);
getpri();
while (scanf("%I64d %I64d", &x, &y)!= EOF ){
if (! x && ! y) printf("%I64d\n", 1ll);
else if (! x) printf("%I64d\n", 0ll);
else if (! y)printf("%I64d\n", 1ll);
else printf("%I64d\n", get());
}

return 0;
}

Tips：

二分求等比数列pi^0 + pi^1 + pi^2 + ... ... + pi^n

1°若n是偶数，那么等比数列一共有奇数项，（1 + pi^1 + pi^2 + ... ... + pi^(n/2-1）) * (1 +pi^(n/2+1） ) + pi^(n/2） 最高次 (n/2-1) +（n/2+1) = n, 加上中间漏掉的n/2项

2°若n是奇数，那么等比数列一共有偶数项，（1 + pi^1 + pi^2 + ... ... + pi^(n/2）) * (1 +pi^(n/2+1） ) 最高次 (n/2） +（n/2+1) = n

ll getsum(ll p, ll n){
if (! n) return 1ll;
if (n&1) return (getsum(p, n/2) * (1 + g(p, n/2+1)))%M;
return (getsum(p, n/2-1) * (1 + g(p, n/2+1)) + g(p, n/2))%M;
}