洛谷P1593 因子和 题解

题目传送门。

看到 \(a,b \le 5 \times 10^7\)，这显然是一道数学题。

每一个数都可以分解成若干个素数的乘积，\(a^b\) 的因子也是一样的。

我们很容易想到对 \(a\) 进行质因数分解，设 \(a = \prod_{i = 1}^k p_i^{e_i}\)，显然有 \(a^b = \prod_{i = 1}^k p_i^{e_i \times b}\)，而 \(a^b\) 的因子和就是 \(\prod_{i = 1}^k \sum_{j = 0}^{e_i \times b} p_i^j\)，直观一点就是 \((1+p_1+p_1^2+\dots+p_1^{e_1 \times b}) \times (1+p_2+p_2^2+\dots+p_2^{e_2 \times b}) \times \dots \times (1+p_k+p_k^2+\dots+p_k^{e_k \times b})\)，此时如果直接求是可以被卡掉的。

如何优化？显然 \(1,p_i,p_i^2,\dots,p_i^{e_i \times b}\) 是一个等比数列，那么根据等比数列求和公式（对于该公式的推导请看这里），我们有 \(\sum_{j = 0}^{e_i \times b} p_i^j = \frac{p_i^{e_i \times b+1}-1}{p_i-1}\)，所以原式变为 \(\prod_{i = 1}^k \frac{p_i^{e_i \times b+1}-1}{p_i-1}\)，此时只需要维护 \(p_i-1\) 在模 \(9901\) 意义下的逆元就行了，因为 \(9901\) 是素数（可以写程序验证），所以可以用费马小定理求逆元。

如果写程序的时候为了加快效率而不开 long long 的话，那么一定要记得 \(9901\) 这个模数很小，求快速幂的时候对与初始的底数一定要记得取模。

当 \(p_i-1 \equiv 0 \pmod {9901}\)，那么也就是说 \(p_i-1\) 在模 \(9901\) 的意义下不存在逆元，此时我们不能再用等比数列公式进行求和了，但是我们可以尝试直接求这个等比数列在模 \(9901\) 意义下的和。因为 \(p_i-1 \equiv 0 \pmod {9901}\) 等价于 \(p_i \equiv 1 \pmod {9901}\)，也就是说存在 \(k \in \mathbb{N}^+\) 使得 \(p_i = k \times 9901+1\)，也就是说此时 \(p_i^n\) 可以转化为 \((k \times 9901+1)^n\)，用二项式定理进行展开：

\[(k \times 9901+1)^n = \sum_{i = 0}^n \binom{n}{i} (k \times 9901)^i 1^{n-i} \]

显然只有在 \(i = 0\) 时 \(\binom{n}{i} (k \times 9901)^i 1^{n-i} \equiv 1 \pmod {9901}\)，其它情况都是 \(\binom{n}{i} (k \times 9901)^i 1^{n-i} \equiv 0 \pmod {9901}\)，也就是说 \((k \times 9901+1)^n \equiv 1 \pmod {9901}\)，可得 \(\forall n \in \mathbb{N},\ p_i^n \equiv 1 \pmod {9901}\)，即：

\[\sum_{j = 0}^{e_i \times b} p_i^j \equiv e_i \times b+1 \pmod {9901} \]

于是就知道了 \(p_i-1 \equiv 0 \pmod {9901}\) 时该等比数列在模 \(9901\) 意义下的和，是 \((e_i \times b+1) \bmod 9901\)。

代码：

#include<bits/stdc++.h>
using namespace std;
const int mod = 9901;
const int N = 1e5+5;
struct node
{
	int p;
	int e;
}s[N];
int main()
{
	ios_base::sync_with_stdio(0);
	cin.tie(0);
	cout.tie(0);
	int a,b,num = 0;
	cin >> a >> b;
	for(int i = 2;i*i<=a;++i)
	{
		int cnt = 0;
		while(a%i == 0)
		{
			++cnt;
			a/=i;
		}
		if(cnt)
		{
			s[++num] = {i,cnt*b};
		}
	}
	if(a>1)
	{
		s[++num] = {a,b};
	}
	int sum = 1;
	for(int i = 1;i<=num;++i)
	{
		int ni = 1,x = (s[i].p-1)%mod,y = mod-2;
        if(!x)
        {
            sum = sum*((s[i].e+1)%mod)%mod;
            continue;
        }
		while(y)
		{
			if(y&1)
			{
				ni = ni*x%mod;
			}
			x = x*x%mod;
			y>>=1;
		}
		int ans = 1;
		x = s[i].p%mod,y = s[i].e+1;
		while(y)
		{
			if(y&1)
			{
				ans = ans*x%mod;
			}
			x = x*x%mod;
			y>>=1;
		}
		sum = sum*(ans-1+mod)%mod*ni%mod;
	}
	cout << sum;
	return 0;
}