【POJ 1845】Sumdiv

题目描述

求 $A^B$ 的所有约数之和 $\mod 9901$。$1≤A,B≤5\times10^7$。

算法分析

将 $A$ 分解质因数得 $A=p_0^{c_0}p_1^{c_1}...P_n^{c_n}$，则 $A^B=p_0^{c_0B}p_1^{c_1B}...P_n^{c_nB}$，它的约数集合就是 $\{p_0^{k_0}p_1^{k_1}...P_n^{k_n}|0≤k_i≤c_iB\}$。

根据乘法原理，$A^B$ 的约数之和为 $(p_0^0+p_0^1+...+p_0^{c_0B})\times(p_1^0+p_1^1+...+p_1^{c_1B})\times...\times(p_n^0+p_n^1+...+p_n^{c_nB})$。

上面乘法运算的每个因数都是一个等比数列，但由于公式涉及除法，而取模运算对除法不满足分配率，因此不能用等比数列求和公式求解。

我们考虑分治，将 $x^0+x^1+...+x^n$ 拆成 $x^0+x^1+...+x^{\lfloor\frac{n-1}{2}\rfloor}+(x^0+x^1+...+x^{\lfloor\frac{n-1}{2}\rfloor})\times x^{\lfloor\frac{n-1}{2}\rfloor+1}$，若元素个数是奇数，则还要再乘上 $x^n$，这样就把问题规模缩小了一半，配合快速幂，可以在 $O(log_2n)$ 的时间复杂度求解。

实现时要提高分解质因数的效率，枚举到 $\sqrt{n}$ 然后特判比枚举到 $n$ 要快，还要注意乘法运算可能溢出，使用 $64$ 位整数储存中间结果。

代码实现

#include <cstdio>
#include <vector>
const int mod=9901;
std::vector<int> prime;
std::vector<int> times;
void divide(int num) {
    for(int i=2,end=num;i*i<=end;++i) {
        if(num%i==0) {
            prime.push_back(i);
            int cnt=0;
            while(num%i==0) {
                ++cnt;
                num/=i;
            }
            times.push_back(cnt);
        }
    }
    if(num>1) {
        prime.push_back(num);
        times.push_back(1);
    }
}
int pow(int n,int k) {
    int ans=1;
    while(k) {
        if(k&1) ans=ans*n%mod;
        n=(long long int)n*n%mod;k>>=1;
    }
    return ans;
}
int powsum(int n,int k) {
    if(k==0) return 1;
    int temp=powsum(n,(k-1)/2);
    int ans=(temp+temp*pow(n,(k-1)/2+1)%mod)%mod;
    if(!(k&1)) ans=(ans+pow(n,k))%mod;
    return ans;
}
int main() {
    int a,b;scanf("%d%d",&a,&b);
    divide(a);
    for(int i=0,end=times.size();i<end;++i) times[i]*=b;
    int ans=1;
    for(int i=0,end=prime.size();i<end;++i) ans=ans*powsum(prime[i],times[i])%mod;
    printf("%d\n",ans);
    return 0;
}