洛谷题单指南-数学基础问题-P1593 因子和

原题链接：https://www.luogu.com.cn/problem/P1593

题意解读：计算a^b的因子和。

解题思路：

1、如何计算因子和

我们知道，对于一个整数x，分解质因数后x = p1^a1*p2^a2*...*pn^an

其所有因子和可以表示为(p1⁰+p1¹+p1²+...+p1^a1)*(p2⁰+p2¹+p2²+...+p2^a2)*...*(pn⁰+pn¹+pn²+...+pn^an)

这个式子拆开就是因子p1^a*p2^b*...*pn^c的所有可能组合之和，a可以取0~a1，b可以取0~a2，c可以取0~an

接下来，对于一个数a^b的因子和

将a分解质因数p1^a1*p2^a2*...*pn^an, a^b = p1^a1*b * p2^a2*b * ... * pn^an*b

对于每一个质因数pi，计算pi⁰+pi¹+pi²+...+pi^ai*b, 再把结果相乘。

2、如何计算p⁰+p¹+p²+...+p^x

可以采用分治法

设f(p, x)表示计算p⁰+p¹+p²+...+p^x

当x == 0时，f(p, x) = 1

当x是奇数时，p⁰+p¹+p²+...+p^x一共有偶数项，f(p, x) = f(p, x/2) + f(p, x / 2) * p^x/2+1

例如x=5时，f(p, x/2) = p⁰+p¹+p², f(p, x/2) * p^x/2+1 = p³+...+p⁵, 所以f(p, x) = f(p, x/2) + f(p, x / 2) * p^x/2+1 = f(p, x/2) * (1 + p^x/2+1)

当x是偶数时，p0+p1+p2+...+px一共有奇数项，f(p, x) = f(p, x/2-1) + p^x/2+ f(p,x/2-1) * p^x/2+1

例如x=6时，f(p, x/2-1) = p⁰+p¹+p²表示左半部分, p^x/2=p3表示中间项，f(p,x/2-1) * p^x/2+1=p⁴+p⁵+p⁶表示右半部分，所以f(p, x) = f(p, x/2-1) + p^x/2+ f(p,x/2-1) * p^x/2+1 = f(p, x/2-1) * (1 + p^x/2+1 ) + p^x/2

递归即可实现上述过程。

//计算：(p^0 + p^1 + p^2 + ... + p^x 
LL f(LL p, LL x)
{
    if(x == 0) return 1;
    if(x % 2 == 1) return f(p, x / 2) * (1 + ksm(p, x / 2 + 1, MOD)) % MOD; //ksm(a, k, M)是计算a^k % M
    else return (f(p, x / 2 - 1) * (1 + ksm(p, x / 2 + 1, MOD)) % MOD + ksm(p, x / 2, MOD)) % MOD;
}

3、如何快速计算a^k % M

在上述分治过程中，涉及计算p^x/2，p^x/2+1，如果采用普通的循环相乘求幂，时间复杂度过高，可以采用快速幂算法。

快速幂：a^k % M

普通算法要把a乘k次，本题k在10^7范围，过大。

快速幂的思想是把k当做二进制，例如计算3^12，

12用二进制表示为1100，也就是2^3+2^2 = 8 + 4，所以3^12 = 3^(8+4) = 3^8 * 3^4

如何用程序模拟上述过程呢？

另初始x = 3，设答案res = 1

只需要从低到高遍历每一个二进制位，如果为1，则res = res * x ，再把x = x * x 

这样只需要遍历4次，而不是12次。快速幂的代码为（每一步加上了取模）：

//计算a^k % p
int ksm(int a, int k, int M)
{
    int res = 1;
    while(k) 
    {
        if(k & 1) res = 1ll * res * a % M; //如果k的最二进制低位为1
        k = k >> 1; //k右移一位，考察下一个二进制位
        a = 1ll * a * a % M; //a每次翻倍
    }
}

100分代码：

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;

const int MOD = 9901;
long long a, b, ans = 1;

//快速幂：a^k % M
LL ksm(LL a, LL k, LL M)
{
    LL res = 1;
    while(k)
    {
        if(k & 1) res = res * a % M;
        k = k >> 1;
        a = a * a % M;
    }
    return res;
}

//计算：(p^0 + p^1 + p^2 + ... + p^x) % MOD
LL f(LL p, LL x)
{
    if(x == 0) return 1;
    if(x % 2 == 1) return f(p, x / 2) * (1 + ksm(p, x / 2 + 1, MOD)) % MOD;
    else return (f(p, x / 2 - 1) * (1 + ksm(p, x / 2 + 1, MOD)) % MOD + ksm(p, x / 2, MOD)) % MOD;
}

int main()
{
    cin >> a >> b;
    //分解质因数
    for(LL i = 2; i * i <= a; i++)
    {
        int cnt = 0;
        if(a % i == 0)
        {
            while(a % i == 0)
            {
                a /= i;
                cnt++;
            }
        }
        //cout << i << "^" << cnt * b << endl;
        //质因数i，指数cnt*b
        ans = ans * f(i, cnt * b) %  MOD;
    }
    if(a > 1) ans = ans * f(a, b) % MOD; //质因数a，指数b
    cout << ans;

    return 0;
}