本文已参与「新人创作礼」活动，一起开启掘金创作之路。

@TOC

Sumdiv

传送门

问题描述

求$A^B$的所有约数之和$mod\ 9901(0\leq A,B\leq 5*10^7)$

问题分析

把$A$质因数分解，表示为$p_1^{k_1}*p_2^{c_2}*\cdots*p_n^{c_n}$（其中$p_i$都是质数）

那么$A^B=(p_1^{c_1}*p_2^{c_2}*\cdots*p_n^{c_n})^B=p_1^{B*c_1}*p_2^{B*c_2}*\cdots*p_n^{B*c_n}$。

则$A^B$的所有约数为集合$\{p_1^{k_1}*p_2^{k_2}*\cdots*p_n^{k_n}\}$，其中$0\leq k_i\leq B*c_i(1\leq i\leq n)$

根据乘法分配律，$A^B$的所有约数的和就是$(1+p_1+p_1^2+\cdots+p_1^{B*c_1})*(1+p_2+p_2^2+\cdots+p_2^{B*c_2})*\cdots*(1+p_n+p_n^2+\cdots+p_n^{B*c_n})$

有关质因数分解和约数的内容将在后面详细介绍。

上述式子中每个括号中的内容都是等比数列，如果使用等比数列求和公式，就需要用到除法。但是答案要对$9901$取模，取模运算只对加、减、乘有分配率，不能直接对分子、分母分别取模后再做除法。我们可用换一种思路，使用分治法进行等比数列求和。

问题：使用分治法求$sum(p,c)=1+p+p^2+\cdots+p^c=\ ?$

若$c$为奇数： ​ $sum(p,c)=(1+p+\cdots+p^{\frac{c-1}{2}})+(p^\frac{c+1}{2}+\cdots+p^c)$

​ $=(1+p+\cdots+p^{\frac{c-1}{2}})+p^\frac{c+1}2*(1+p+\cdots+p^{\frac{c-1}{2}})$

​ $=(1+p^\frac{c+1}2)*sum(p,\frac{c-1}2)$

若$c$为偶数，同理有：

​ $sum(p,c)=(1+p^\frac{c}{2})*sum(p,\frac{c}2-1)+p^c$

每次分治（递归）之后，问题规模均会缩小一半，配合快速幂即可在$O(log\ c)$的时间内求出等比数列。

代码如下：

#include <bits/stdc++.h>
using namespace std;
#define mem(a) memset(a, 0, sizeof(a))
#define dbg(x) cout << #x << " = " << x << endl
#define fi(i, l, r) for (int i = l; i < r; i++)
#define cd(a) scanf("%d", &a)
typedef long long ll;
#define mod 9901
int power(int a, int b)
{
    int ans=1%mod;
    while(b)
    {
        if(b&1)ans=(ll)ans*a%mod;
        a=(ll)a*a%mod;
        b>>=1;
    }
    return ans;
}
int sum(int p, int c)
{
    if(c==0)return 1;
    if(c%2)return ((1+power(p,(c+1)/2))*sum(p,(c-1)/2))%mod;
    else return ((1+power(p,c/2))*sum(p,c/2-1)+power(p,c))%mod;
}
unordered_map<int,int>prime;
void divid(int n)
{
    for(int i=2;i*i<=n;i++)
        while(n%i==0)
            prime[i]++,n/=i;
    if(n>1) prime[n]++;
}
int main()
{
    int a,b;
    cin>>a>>b;
    divid(a);
    int ans=1;
    for(unordered_map<int,int>::iterator it=prime.begin();it!=prime.end();it++)
    {
        ans=ans*sum(it->first,b*it->second)%mod;
    }
    if(a==0)ans=0; // 注意a,b不同时为0，但是a可能是0！
    cout<<ans<<endl;
    return 0;
}

同步发文于我的CSDN,原创不易，转载请附上原文链接哦~
Tisfy：letmefly.blog.csdn.net/article/det…