Sumdiv（经典数学问题，快速幂，唯一分解定理，约数和定理，递归求等比，同余

• 题目描述
  有两个自然数a和b(a,b≤50000000)
  求a的b次方的所有约数之和模9901
• 输入格式
  一行，包含由空格分隔的两个自然数a和b
• 输出格式
  一行，a的b次方的约数和模9901
• 样例输入
  2 3
• 样例输出
  15
• 样例解释
  8的约数是1,2,4,8， 它们的总和是15
  15模9901是15

题意：求A^B的所有约数之和。
用到的知识:
1，唯一分解定理:
A=(p1ⁿ¹p2ⁿ²p3ⁿ³…pk^nk); 把底数A分解 其中p_i为素数
2，约数和定理：就是把A所有的约数
由上面我们，我们可以得出A=(p1ⁿ¹p2ⁿ²p3ⁿ³…pk^nk)；
sum=(1+p1¹+p1²+……p1ⁿ¹) * (1+p2¹+p2²+……p2ⁿ²) * ……(1+pk¹+pk²……+pk^nk);
3,同余模公式:
(a+b)%m=(a%m+b%m);
(ab)%m=(a%mb%m);

解法：

1.素因子分解

A首先对第一个素数2不断取模，A%2==0时 ，记录2出现的次数+1，A/=2；（这样会节省一部分时间）
当A%2!=0时，则A对下一个连续素数3不断取模…

以此类推，直到A==1为止。

注意特殊判定，当A本身就是素数时，无法分解，它自己就是其本身的素数分解式。

2.求a^b的所有约数之和

用递归二分求等比数列1+pi+pi^2+pi^3+...+pi^n：

（1）若n为奇数,一共有偶数项，则：
1 + p + p^2 + p^3 +…+ p^n
= (1+p^(n/2+1)) + p * (1+p^(n/2+1)) +…+ p^(n/2) * (1+p^(n/2+1))
= (1 + p + p^2 +…+ p^(n/2)) * (1 + p^(n/2+1))

上式红色加粗的前半部分恰好就是原式的一半，那么只需要不断递归二分求和就可以了，后半部分为幂次式，将在下面第4点讲述计算方法。
（2）若n为偶数,一共有奇数项,则:
1 + p + p^2 + p^3 +…+ p^n
= (1+p^(n/2+1)) + p * (1+p^(n/2+1)) +…+ p^(n/2-1) * (1+p^(n/2+1)) + p^(n/2)
= (1 + p + p^2 +…+ p^(n/2-1)) * (1+p^(n/2+1)) + p^(n/2);
上式红色加粗的前半部分恰好就是原式的一半，依然递归求解

#include <map>
#include <queue>
#include <string>
#include<iostream>
#include<stdio.h>
#include<string.h>
#include <algorithm>
#include <math.h>
typedef long long ll;
using namespace std;
const int maxn=10000;
#define
const int mod=9901;
const int MOD=9901;

int n,m;
int dis[maxn],vis[maxn];


/*void prime(){
for(int i=2;i<maxn;i++)
{
if(dis[i]==0) vis[++cnt]=i;
for(int j=1;j<=cnt&&i*vis[j]<maxn;j++)
{
dis[i*vis[j]]=1;
if(i%vis[j]==0) break;
}
}
}*/

ll qpow(ll a, ll b){//快速幂
ll sum=1;
while(b){
if(b&1) sum=sum*a%MOD;
b>>=1;
a=a*a%MOD;
}
return sum;
}

ll f(ll a, ll b){
if(b==0){
return 1;
}else if(b%2){
return (f(a,b/2)*(1+qpow(a,b/2+1)))%MOD;
}else {
return (f(a,b/2-1)*(1+qpow(a,b/2+1))+qpow(a,b/2))%MOD;
}
}


int main() {
scanf("%d%d",&n,&m);
int cnt=0;
for(int i=2;i*i<=n;)//唯一分解定理
{
if(n%i==0){
dis[cnt]=i;
vis[cnt]=0;
while(n%i==0){
vis[cnt]++;
n/=i;
}
cnt++;
}
if(i==2) i++;
else i+=2;
}

if(n!=1){
dis[cnt]=n;
vis[cnt++]=1;
}


int sum=1;
for(int i=0;i<cnt;i++){
sum=(sum*(f(dis[i],vis[i]*m)%MOD))%MOD;
}

cout<<sum<<endl;
return 0;
}