[ 题解 ] [ 数学 ] [ JZOJ5809 ] 数羊

题面

牧羊人 A 和牧羊人 B 总是很无聊，所以他们要玩一个游戏。A 有 $a$ 只羊，B 有 $b$ 只羊。他们想要知道 $a^b$ 的因子和是多少。这就很为难两个牧羊人了，由于答案太大，你能不能告诉我答案取模 $9901$ 的数。

Example In #1

2 3

Example Out #1

15

对于 $100$ 的数据，$0\le a,b\le 50000000$。

题解

举例：

$$\begin{array}{rl}S& =(1+2+4)(1+3+9)\\ & =(2+4+3+9)+(2\times 3)+(2\times 9)+(4\times 3)+(4\times 9)\\ & =7\times 13=91\end{array}$$

可以看出 $S$ 即为 $2^3\times 3^2=72$ 的因数和。
推广到一般：

$$\begin{array}{r}S=(1+p_1^{1b}+p_1^{2b}+\cdots +p_1^{k_{1}b})\times \cdots \times (1+p_m^{1b}+p_m^{2b}+\cdots +p_m^{k_{m}b})\end{array}$$

等比数列求和公式：

$$S_n=\frac{a_1\times (1-q^n)}{1-q}=\frac{a_1-a_{n}q}{1-q}=\frac{a_{n}q-a_1}{1-q}$$

$n$ 为项数，$S_n$ 为和，$q$ 为公比，$a_1$ 为首项，$a_n$ 为末项。
此题中：

$$\begin{array}{rl}& n=m(m为质因数个数)\\ & q=p_i\\ & a_1=1\\ & a_n=p_i^{k_{i}b}\end{array}$$

代入式中，可见：

$$\begin{array}{rl}S& =\prod _{i=1}^m\frac{p_i^{k_{i}b+1}-1}{p_i-1}\end{array}$$

其中，除法取模用费马小定理求乘法逆元。

$$\begin{array}{rl}S& =\prod _{i=1}^m(p_i^{k_{i}b+1}-1)\times (p_i-1{)}^{(mod-2)}\end{array}$$

#include <iostream>
#include <utility>
#include <vector>
#include <cmath>

const int MOD = 9901;
using i64 = long long;

i64 q_pow(i64 a, int b)
{
    i64 res = 1;
    while (b > 0)
    {
        if (b % 2 == 1)
            res = res * a % MOD;
        a = a * a % MOD;
        b /= 2;
    }

    return res;
}

bool is_prime(int n)
{
    for (int i = 2; i <= std::sqrt(n); i++)
        if (n % i == 0)
            return false;
    return true;
}

int main()
{
    int a, b;
    std::cin >> a >> b;

    std::vector<std::pair<int, int>> fac;
    for (int i = 2; i <= std::sqrt(a); i++)
    {
        if (is_prime(i) && a % i == 0)
        {
            fac.emplace_back(i, 0);
            while (a % i == 0)
            {
                a /= i;
                fac.back().second++;
            }
            fac.back().second *= b;
        }
    }
    if (a != 1)
        fac.emplace_back(a, b);

    i64 ans = 1;
    for (auto f : fac)
    {
        i64 res = (q_pow(f.first, f.second + 1) - 1) % MOD;
        res = res * q_pow(f.first - 1, MOD - 2) % MOD;
        ans = ans * res % MOD;
    }

    std::cout << ans;
    return 0;
}