POJ 1845 Sumdiv(因子分解+快速幂+二分求和)

题意：给你A，B，让求A^B所有的因子和模上9901

思路：A可以拆成素因子的乘积: A = p1^x1 * p2^x2 *...* pn^xn

那么A^B = p1^(B*x1) * p2^(B*x2) *...* pn^(B*xn)

那么A^B所有的素因子和就是

(p1^0 + p1^1 + p1^2 + ... + p1^(B*x1) ) * (p2^0 + p2^1 + ... + p2^(B*x2) ) * ... * (pn^0 + pn^1 + ... + pn^(B*xn))

可以看出每一个括号内都是等比数列，但是不要用等比数列公式，因为有除法(刚开始我用除法，然后求了模的逆元，wa到爽死)，因为不一定满足乘法逆元所需要的条件，除数与模数可能不互素(除数可能是模数的多少倍)。既然不能用公式，那么就要借助于二分了。比如如下式子求和:A1+A2+A3+A4 = A1+A2+A2(A1+A2)。通过这个式子发现，只要求出来A2就行了，然后只要计算一次A1+A2，就可以省掉一半的计算量。那么同理A1+A2也可以继续往下分。

现在推广到一般式。A1+A2+...+An

1) n为偶数: A1+A2+...+An = A1+A2+ ...+A(n/2)+ A(n/2)(A1+A2+...+A(n/2))

2) n为奇数: A1+A2+...+An = A1+A2+ ...+A(n/2)+ A(n/2)(A1+A2+...+A(n/2)) + An

推出来这些就可以递归求解了。

注：找素因子时，打个素数表，只需要打到sqrt(n)就行了，因为只可能在sqrt(n)里面，如果有比sqrt(n)大的两个素因子，乘积自然就大于n了，所以只需要sqrt(n)就可以了。因为就算有一个大素数和一个小素数相乘得来，那么在约掉小素数的时候，只剩大素数了，这会就直接跑到循坏外判断了。

#include <cstdio>
#include <iostream>
#include <cstring>
#include <algorithm>
using namespace std;

typedef long long ll;
const int maxn = 10010;
const ll mod = 9901LL;
ll A, B;
struct Factor {
ll fac;
ll cnt;
}factor[maxn];
int tot;
bool prime[maxn + 100];
int pr[maxn];//素数表
int pr_cnt;
void init_prime()
{
memset(prime, true, sizeof(prime));
prime[0] = prime[1] = false;
for (int i = 2; i * i <= maxn; i++)
if (prime[i])
for (int j = i + i; j < maxn; j += i)
prime[j] = false;
pr_cnt = 0;
for (int i = 2; i <= maxn; i++)
if (prime[i])
pr[pr_cnt++] = i;
}
void init()//找到所有的素因子
{
tot = 0;
memset(factor, 0, sizeof(factor));
for (int i = 0; i < pr_cnt && pr[i] <= A; i++)
{
if (A % pr[i] == 0)
{
factor[tot].fac = pr[i];
while (A % pr[i] == 0)
{
factor[tot].cnt++;
A /= pr[i];
}
factor[tot].cnt *= B;
tot++;
}
}
if (A > 1)
{
factor[tot].fac = A;
factor[tot++].cnt = B;
}
}
ll quickpow(ll a, ll b, ll mod)
{
ll ans = 1LL;
while (b)
{
if (b & 1) ans = ans * a % mod;
a = a * a % mod;
b >>= 1;
}
return ans % mod;
}
ll binary_pow(ll a, ll b, ll mod)//计算等比数列的和
{
if (b == 0) return 1LL;
if (b == 1) return a;
ll ans = 0;
if (b & 1)
{
ans = quickpow(a, b, mod);
ans = (ans + (quickpow(a, b / 2, mod) + 1LL) % mod * binary_pow(a, b / 2, mod)) % mod;
}
else
ans = (quickpow(a, b  / 2, mod) + 1LL) % mod * binary_pow(a, b / 2, mod) % mod;
return ans;
}
void solve()
{
if (B == 0)
{
puts("1");
return;
}
if (A == 0)
{
puts("0");
return;
}
init();
ll ans = 1LL;
for (int i = 0; i < tot; i++)
{
ans = ans * (binary_pow(factor[i].fac, factor[i].cnt, mod) + 1LL) % mod;
}
cout << ans << endl;
}
int main()
{
init_prime();
while (cin >> A >> B)
{
solve();
}
return 0;
}