POJ——1845 Sumdiv （尤拉篩+快速冪+遞迴二分）
Consider two natural numbers A and B. Let S be the sum of all natural divisors of A^B. Determine S modulo 9901 (the rest of the division of S by 9901).
Input
The only line contains the two natural numbers A and B, (0 <= A,B <= 50000000)separated by blanks.
Output
The only line of the output will contain S modulo 9901.
Sample Input
2 3Sample Output
15Hint
2^3 = 8.
The natural divisors of 8 are: 1,2,4,8. Their sum is 15.
15 modulo 9901 is 15 (that should be output)
題意：求a的b次方的約數和。
題解：
這道題涉及到的知識：
1. 整數的唯一分解定理：
任意正整數都有且只有一種方式寫出其素因子的乘積表示式。
A=(p1^k1)*(p2^k2)*(p3^k3)*....*(pn^kn) 其中pi均為素數
2.約數和公式：
對於已經分解的整數A=(p1^k1)*(p2^k2)*(p3^k3)*....*(pn^kn)
有A的所有因子之和為：
S = (1+p1+p1^2+p1^3+...p1^k1) * (1+p2+p2^2+p2^3+….p2^k2) * (1+p3+ p3^3+…+ p3^k3) * .... * (1+pn+pn^2+pn^3+...pn^kn)
所以A的B次方=(p1^ (k1*B) )*(p2^ (k2*B) )*(p3^ (k3*B) )*....*(pn^( kn*B) )
S = (1+p1+p1^2+p1^3+...p1^(k1*B)) * (1+p2+p2^2+p2^3+….p2^(k2*B)) * (1+p3+ p3^3+…+ p3^ (k3*B)) * .... * (1+pn+pn^2+pn^3+...pn^( kn*B))
公式記住就好了！！！
3.同餘模公式：
(a+b)%m=(a%m+b%m)%m
(a*b)%m=(a%m*b%m)%m
4.二分求等比數列和，也可以逆元。
：
因為這個題目有1這一項所以第4個知識點反著用。
其他注意看程式碼註釋！
#include <iostream>
using namespace std;
typedef long long ll;
const int mod=9901;
const int MAX = 1e6+100;
int pr[MAX];
bool vis[MAX];
struct hh{
int ji;
int sum;
}w[MAX];
ll quick(ll a,ll b,ll c){//快速冪
ll ans=1;
while(b){
if(b&1) ans=(ans*a)%c;
b>>=1;
a=(a*a)%c;
}
return ans%c;
}
void getpr(){ //尤拉篩素數 時間複雜度O(n)
int cnt=0;
vis[0]=vis[1]=true;
for (int i = 2; i < MAX;i++){
if(!vis[i]) pr[cnt++]=i;
for (int j = 0; j < cnt&& i*pr[j]< MAX;j++){
vis[i*pr[j]]=true;
if(i%pr[j]==0) break;
}
}
}
ll jisuan(ll a,ll b ){//遞迴二分求等比數列 1+p+p^2+p^3+…+p^n
if(b==0) return 1;
else if(b&1) return jisuan(a,b/2)*(1+quick(a,b/2+1,mod))%mod;
else return (jisuan(a,b/2-1)*(1+quick(a,b/2+1,mod))+quick(a,b/2,mod))%mod;
}
int main(){
ll a,b;
getpr();
cin >> a >> b;
int ans=0;
for (int i = 0; pr[i]*pr[i] <= a;i++){//分解質因數
if(a%pr[i]==0){
w[ans].ji=pr[i];
while(a%pr[i]==0){
a/=pr[i];
w[ans].sum++;
}
ans++;
}
}
if(a!=1){//如果本身是質數加特判
w[ans].ji=a;
w[ans].sum++;
ans++;
}
ll cnt=1;
for (int i = 0; i < ans;i++){
cnt=cnt*jisuan(w[i].ji,w[i].sum*b)%mod;
}
cout << cnt << endl;
return 0;
}