【題解】POJ 1845 Sumdiv
【題目大意】
給定\(A\)和\(B\)，求\(A^B\)的所有約數之和，對\(9901\)取模。
（對於全部數據，\(0<= A <= B <=50,000,000\)）
【樣例輸入】
2 3【樣例輸出】
15【算法關鍵詞】
- 數論
- 綜合模板
二分，乘法逆元
【題解】
不管什麽題首先思考的肯定是暴力解法。起碼可以騙分啊，當然，如果能一眼標算，那再好不過了。
這道題暴力做法就不說了，其實仔細思考也不會真的打暴力吧。。。
看見約數，首先想到的應該就是數論的有關知識。那麽這道題實際上就是在求\(A^B\)的約數之和，問題的難點就在於求解答案的時間優化。
首先，可以思考唯一分解定理，即任意一個自然數都可分解且只能分解成以下形式：
其中，\(p\)為質因數，\(k\)為自然數。
那如何實現呢？
首先通過線性篩把題目範圍內的質數篩出來存在一個數組裏，然後枚舉數組裏的質數是否能被\(A\)整除，即枚舉\(p\)，如果能，就以\(p\)為底枚舉\(k\)，並存入相對應的數組。
於是就可分解為：
\[
A^B=p_1^{k_1*B}*p_2^{k_2*B}*p_3^{k_3*B}*...*p_m^{k_m*B}
\]
設\(a\)所有約數和為\(ans(a)\)
對於兩個互質的數\(a\),\(b\),必定有\(ans(ab)=ans(a)*ans(b)\)
令\(a\)約數為\(x_1\),\(x_2\)...\(x_m\)，\(b\)約數為\(y_1\),\(y_2\),...\(y_n\)
由\(a\),\(b\)互質,\(a_1->a_m\)與\(b_2->b_n\)無公共元素
那麽
\[
ans(ab)=\sum_{i=1,j=1}^{i<=m,j<=n}x_i*y_j\\quad=\sum_{i=1}^{m}x_i*\sum_{j=1}^{n}y_j\=ans(a)*ans(b)
\]
得證 這不是積性函數的性質嗎
同理，\(p_i^{k_i*B}\)與\(p_j^{k_j*B}\)互質
所以:
\[
ans(A^B)=\prod_{i=1}^msum_i
\]
其中:
\[ sum_i=\sum_{j=0}^{k_i*B}p_i^j \]
即:
\[ ans(A^B)=\prod_{i=1}^{m}\sum_{j=0}^{k_i*B}p_i^j \]
（理解不透徹的可以把這個拆開看看）
公式推到這裏，其實就可以寫代碼了，在這裏，還要註意\(sum_i\)的計算。
\[ sum_i=p_i^0+p_i^1+p_i^2+...+p_i^{k_i} \]
這裏有兩種方法。
二分+遞歸思想。
若\(k_i\)為奇數:
\[ sum(p, k) = sum(p, k/2)*(1+p^{n/2+1}) \]
若\(k_i\)為偶數:
\[ sum(p, k) = sum(p, k/2-1)*(1+p^{n/2+1})+p^{n/2} \]
證明的話可以把上面的式子拆開來看，這裏不再贅述。- 等比數列通項公式。
\[ sum(p, k) = \frac{p^{k+1}-1}{p-1} \]
當然，如果在模的意義下這樣做就需要乘法逆元的相關知識，這裏就只給出法一的代碼。法二的代碼可以自己嘗試下。其實是因為我懶
#include<cmath>
#include<ctime>
#include<queue>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
#define debug() puts("FBI WARNING!")
#define ll long long
#define R register
using namespace std;
/*常量申明*/
const int MAX = 500000;
const int BIG_SEVEN = 777777;
const int MEIZI = 100000;
const int mod = 9901;
/*函數申明*/
void ola(); void test(); void divide(); void cls();void solve();
ll binary_pow(ll, ll); ll sum(ll, ll); ll read();
/*變量申明*/
ll a, b, s, prim[BIG_SEVEN], p[MEIZI], k[MEIZI];
bool vis[MAX + 5];
int cnt, tot;
int main(){
ola();
while (~scanf("%lld %lld", &a, &b)) {
cls();
divide();
solve();
//test();
}
return 0;
}
inline ll read(){ // 快讀 (貌似沒用上)
ll f = 1, x = 0;char ch;
do { ch = getchar(); if (ch == '-') f = -1; } while (ch < '0'||ch>'9');
do {x = x*10+ch-'0'; ch = getchar(); } while (ch >= '0' && ch <= '9');
return f*x;
}
inline void ola() { // 歐拉篩
vis[1] = 1;
for (R int i = 2;i < MAX; ++i) {
if (!vis[i]) {
prim[cnt++] = i;
}
for (R int j = 0;j < cnt ; ++j) {
if (i*prim[j] > MAX) break;
vis[i*prim[j]] = 1;
if (i%prim[j] == 0) break;
}
}
}
inline void cls() {tot = 0, cnt = 0;} // 多組數據初始化
inline void divide() { // 唯一分解定理
memset(k, 0, sizeof(k));ll rose = a;
for (R int i = 0; prim[i] <= rose/prim[i]; ++i) {
if (rose%prim[i] == 0) {
p[tot] = prim[i];
while (rose%prim[i] == 0) {
// debug();
k[tot]++;
rose /= prim[i];
}
tot++;
}
}
if (rose != 1) {
p[tot] = rose;
k[tot++] = 1;
}
}
inline void solve() { // 累加答案
ll ans = 1;
for (R int i = 0; i < tot; ++i) {
ans *= (sum(p[i], b*k[i]) %mod);
ans %= mod;
}
printf("%lld\n", ans);
}
inline void test() { // 測試函數
/*check ola()*/
for (R int i = 0;i < cnt; ++i) {
cout << prim[i] << " ";
}
printf("%.2lf\n",(double)clock()/CLOCKS_PER_SEC);
/*check divide()*/
for (R int i = 0;i < tot; ++i) {
printf("%lld^%lld\n", p[i], k[i]);
}
cout << tot;
}
inline ll binary_pow(ll a, ll b) { // 二進制實現的快速模冪
ll ans = 1, tmp = a%mod;
while (b) {
if (b&1) {ans *= tmp; ans %= mod;}
b >>= 1;
tmp *= tmp; tmp %= mod;
}
return ans;
}
inline ll sum(ll p, ll k) { // 約數和二分實現（p^0+p^1+...+p^k）
if (p == 0) return 0;
if (k == 0) return 1;
if (k & 1) return ((1+binary_pow(p, k/2+1)) %mod * sum(p, k/2) %mod) %mod; // 奇數
else return ((1+binary_pow(p, k/2+1)) %mod*sum(p, k/2-1) + binary_pow(p, k/2) %mod) %mod;
}三玖天下第一！
【題解】POJ 1845 Sumdiv