一本通 1633：【例 3】Sumdiv

题目传送门

---

今天早上考试考了这道题

题意：求$ A^{B} $所有约数之和%9901的结果。

思路：[暴力]快速幂+线性判约数再求和，30分。

[正解]看到求约数之和，很自然想到唯一分解定理，对于正整数N， N = $ a_1^{b_1}a_2\dots a_n^{b_n} $ 而言，它的约数之和为 $ (1+a_1+a_1^2+\dots +a_1^{{b_1})(1+a_2+a_2}2+\dots + a_2^{b_2})\dots (1+a_n+a_n^2+\dots +a_n^{b_n}) $
而 $ A^B $ 只需将一直加到 $ a_i^{b_{i}* B} $ 即可。直接暴力求上述多项式的解似乎有50分。

接着考虑优化，不难看出对于每一个多项式，都是一个等比数列的和，即$ 1+a_i+a_i^2+\dots +a_i^{b_i} = \frac{a_i^{bi+1}-1}{a_i-1} \(， 所以 \) 1+a_i+a_i^2+\dots +a_i^{b_i} \equiv \frac{a_i^{bi+1}-1}{a_i-1} $ (mod p)，那么这里我们只需要求出 $ a_i^{b_i+1}-1 $ 以及 $ a_i-1 $ 在%p下的逆元即可。逆元这里可以用费马小定理求。（我也不晓得为什么递推求会WA掉几个点）

对了，如果 $ a_i-1 $ 恰好是p的倍数怎么办呢？这种情况下逆元是不存在的。但是由于 $ a_i-1\equiv 0 $ (mod p), 那么 $ a_i\equiv 1 $ (mod p),这样一来
$ 1+a_i+a_i^2+\dots +a_i^{b_i}\equiv 1+1+1+\dots + 1 $ (mod p), 即 $ (b_i+1) $ %p.

---

Code:

#include <iostream>
#include <cstdio>
#include <cstring>
#include <queue>
#include <vector>
#include <algorithm>
#include <cmath>
using namespace std;
//Mystery_Sky
//
#define M 5000101
#define INF 0x3f3f3f3f
#define ll long long
#define Mod 9901
inline ll read()
{
	ll x=0, f=1;char c=getchar();
	while(c<'0'||c>'9') {if(c=='-')f=-1;c=getchar();}
	while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();}
	return x*f;
}
ll A, B;
ll prime[M], tot, check[M], inv[M], cnt;
ll a[M], b[M];

inline void get_inv()
{
	inv[1] = 1;
	int p = Mod;
	for(int i = 2; i <= a[cnt]; i++) {
		inv[i] = ((p-p/i) * inv[p%i]) % p;
		//printf("%d ", inv[i]);
	}
}

inline ll quick_pow(ll x, ll k)
{
	ll ret = 1;
	while(k) {
		if(k & 1) ret = (ret * x) % Mod;
		k >>= 1;
		x = (x * x) % Mod;
	}
	return ret % Mod;
}

int main() {
	A = read(), B = read();
	ll n = A;
	for(int i = 2; i * i <= n; i++) {
		if(n % i == 0) {
			a[++cnt] = i;
			while(n % i == 0) {
				n /= i;
				b[cnt]++;
			}
		}
	}
	if(n > 1) a[++cnt] = n, b[cnt] = 1;
	ll ans = 1;
	for(int i = 1; i <= cnt; i++) {
		if((a[i]-1) % M == 0) {
			ans = ((ll) (B * b[i] + 1) % Mod * ans % Mod) % Mod;
			continue;
		}
		else {
			ll x = quick_pow(a[i], (ll)B * b[i]+1);
			x = (x - 1 + Mod) % Mod;
			ll y = quick_pow(a[i]-1, Mod-2);
			ans = ((ll) x * y % Mod * ans % Mod) % Mod;
		}
	}
	printf("%lld\n", (ans + Mod) % Mod);
	return 0;
}