AcWing 97. 约数之和—算数基本定理、递归

题目链接： AcWing 97. 约数之和
问题描述

分析
这道题的基础是约数的个数与和（基本算数定理求解）需要通过算数基本定理来解决问题，任何正整数$x$的质因子为$p_1,p_2,...,p_n$
那么$x$的任何约数都能由$(p_1^0+p_1^1,....,p_1^{a_1})(p_2^0+p_2^1,....,p_2^{a_2})\ast ...\ast (p_n^0+p_n^1,....,p_n^{a_n})$来表示，所以$x$的所有约数之和就为$(p_1^0+p_1^1,....,p_1^{a_1})(p_2^0+p_2^1,....,p_2^{a_2})\ast ...\ast (p_n^0+p_n^1,....,p_n^{a_n})$，这道题种可以先将$A$质因子分解，将质因子的指数再乘上$B$,但是有个问题,这里的$B$的范围较大，如果遍历算$(p_n^0+p_n^1,....,p_n^{a_n})$会超时，所有这里考虑用递归分治的方法来算$(p_n^0+p_n^1,....,p_n^{a_n})$。
其实之前考虑的是用等比数列的公式$\frac{p^n-1}{p-1}\%mod$和逆元来算，但是需要判断$p-1$是否是mod的倍数，这种情况就比较麻烦，如果题目比较友好，直接用公式+逆元来算会更好一些。

那么现在的问题就变成了如何快速求$(p_n^0+p_n^1,....,p_n^{a_n})$，这里的$a_n$很大，考虑用递归分治来做，
(1)若$a_n$为奇数,那么$p_n^0+p_n^1+....+p_n^{a_n}$一共有a_n+1，为偶数项，
前一半为 $p_n^0+p_n^1+...+p_n^{\frac{a_n-1}{2}}$，
后一半为 $p_n^{\frac{a_n+1}{2}}+...+p_n^{a_n}$，
可以发现后一半每一项都是前一项对应位置的$p_n^{\frac{a_n+1}{2}}$倍，因为$a_n$为奇数，我们把这个写成整除的形式就是$p_n^{\frac{a_n}{2}+1}$

get_sum(a,b)=get_sum(a,b/2)+a^(n/2+1)*get_sum(a,b/2)	
					=(ksm(a,b/2+1)+1)*get_sum(a,b/2)
get_sum(a,b)是求a^0+a^1+...+a^b
ksm是快速幂

(2)若$a_n$为奇数,那么可以求$p_n^{a_n}$+$(p_n^0+p_n^1,....,p_n^{a_n-1})$这里$(p_n^0+p_n^1,....,p_n^{a_n-1})$就变成了奇数项，可以按照(1)来求，$p_n^{a_n}$用快速幂来求

get_sum(a,b)=ksm(a,b)+get_sum(a,b-1)

代码如下:

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
#include<unordered_map>
using namespace std;
const int mod=9901;
int ksm(int a,int b){
    int res=1;
    a%=mod;//a可能比较大a*a会爆int
    while(b){
        if(b&1) res=(res*a)%mod;
        a=(a*a)%mod;
        b>>=1;
    }
    return res%mod;
}
int get_sum(int a,int b){
    if(b==0) return 1;
    if(b&1) return (ksm(a,b/2+1)+1)*get_sum(a,b/2)%mod;
    else return (ksm(a,b)+get_sum(a,b-1))%mod;
}
int main(){
    int a,b;
    cin>>a>>b;
    if(a==0) cout<<0;
    else if(b==0) cout<<1;
    else{
        unordered_map<int,int>mp;
        for(int i=2;i<=a/i;i++)
            while(a%i==0){
                mp[i]+=b;
                a/=i;
            }
        if(a>1) mp[a]+=b;
        int res=1;
        for(unordered_map<int,int>::iterator it=mp.begin();it!=mp.end();it++){
            res=(res*get_sum(it->first,it->second))%mod;

        }
        cout<<res;
    }
    return 0;
}