洛谷 4714.约数个数和

这题看起来很毒瘤,然而也确实很毒瘤……

设答案为 fK(N)f_K(N)
考虑证明对于任意的 KKfKf_K 都是积性函数。

首先,有 fK(N)=dNfK1(d)f_K(N) = \sum\limits_{d | N} f_{K - 1}(d)
fK=fK11f_K = f_{K - 1} * \mathbb 1
易证两个积性函数的狄利克雷卷积仍然是积性函数,这里略过。
于是根据 f0=11f_0 = \mathbb 1 * \mathbb 1,证毕。

现在考虑 fKf_K 的值。
由于 fKf_K 是积性函数,设 NN 的质因子分解形式为 i=1mpici\prod\limits_{i = 1}^m {p_i}^{c_i}fK(N)=i=1mfK(pici)f_K(N) = \prod\limits_{i = 1}^m f_K({p_i}^{c_i})

考虑 pp 为素数时 fK(pc)f_K(p^c) 的值。
显然有 f0(pc)=c+1f_0(p^c) = c + 1
假设把 ccpp 排成一横排,在其中插一个板(两头都可以插),这个板之前的 pp 乘起来就是 pcp^c 的一个约数。
所以 f0(pc)=(c+11)f_0(p^c) = \dbinom{c + 1}{1}
我们猜测 f1(pc)=(c+22)f_1(p^c) = \dbinom{c + 2} 2
相当于先插一个板,然后对这个板之前的这个因子求 f0f_0 的值,所以相当于插两个板。
所以 fK(pc)=(c+K+1K+1)=(c+K+1c)f_K(p^c) = \dbinom{c + K + 1}{K + 1} = \dbinom{c + K + 1} c

由于 cc 最多只有 5959,所以二项式系数可以直接暴力计算。

质因子分解时,必须先用筛法筛出较小的素数然后试除,然后对剩下的数使用 Pollard-Rho,因为这是一个随机算法。
注意快速乘。

代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
#include <cstdio>
#include <cmath>
#include <algorithm>
#include <tr1/unordered_map>
using namespace std;
using namespace tr1;
const int CNT = 1e7;
const long long mod = 998244353;
const int MX = 59;
unordered_map<long long,int> vis;
int mx;
long long inv[MX + 5];
inline long long fmul(long long a,long long b,long long mod)
{
return (a * b - (long long)((long double)a / mod * b) * mod + mod) % mod;
}
inline long long fpow(long long a,long long b,long long mod)
{
long long ret = 1;
for(;b;b >>= 1)
(b & 1) && (ret = fmul(ret,a,mod)),a = fmul(a,a,mod);
return ret;
}
inline bool witness(long long x,long long p)
{
if(fpow(x,p - 1,p) ^ 1)
return 0;
long long k = p - 1;
while(!(k & 1))
{
k >>= 1;
long long t = fpow(x,k,p);
if((t ^ 1) && (t ^ p - 1))
return 0;
if(t == p - 1)
return 1;
}
return 1;
}
inline bool Miller_Rabin(long long n)
{
if(n == 1)
return 0;
const int base[] = {2,3,7,61,24251,19260817};
for(register int i = 0;i < 6;++i)
if(n == base[i])
return 1;
else if(!witness(base[i],n))
return 0;
return 1;
}
void Pollard_Rho(long long n)
{
if(Miller_Rabin(n))
{
mx = max(mx,++vis[n]);
return ;
}
const long long c = 1e9 + 9;
const int k = 200;
long long t1 = 1e9 + 7;
long long t2 = (fmul(t1,t1,n) + c) % n;
long long p = 1;
int i = 0;
while(t1 ^ t2)
{
++i,p = fmul(p,abs(t1 - t2),n);
if(!p)
{
long long d = __gcd(abs(t1 - t2),n);
Pollard_Rho(d),Pollard_Rho(n / d);
return ;
}
if(!(i % k))
{
long long d = __gcd(p,n);
p = 1;
if((d ^ 1) && (d ^ n))
{
Pollard_Rho(d),Pollard_Rho(n / d);
return ;
}
}
t1 = (fmul(t1,t1,n) + c) % n,t2 = (fmul(t2,t2,n) + c) % n,t2 = (fmul(t2,t2,n) + c) % n;
}
long long d = __gcd(p,n);
if((d ^ 1) && (d ^ n))
Pollard_Rho(d),Pollard_Rho(n / d);
}
long long n,k;
long long ans = 1;
int prime[CNT + 5],cnt,bz[CNT + 5];
int main()
{
for(register int i = 2;i <= CNT;++i)
{
if(!bz[i])
prime[++cnt] = i;
for(register int j = 1;j <= cnt && i * prime[j] <= CNT;++j)
{
bz[i * prime[j]] = 1;
if(!(i % prime[j]))
break;
}
}
scanf("%lld%lld",&n,&k);
for(register int i = 1;i <= cnt && prime[i] <= n;++i)
for(;!(n % prime[i]);mx = max(mx,++vis[prime[i]]),n /= prime[i]);
if(n > 1)
Pollard_Rho(n);
inv[1] = 1;
for(register int i = 2;i <= mx;++i)
inv[i] = fmul((mod - mod / i),inv[mod % i],mod);
for(unordered_map<long long,int>::iterator it = vis.begin();it != vis.end();++it)
for(register int i = 1;i <= it->second;++i)
ans = fmul(fmul(ans,it->second + k - i + 2,mod),inv[i],mod);
printf("%lld\n",ans);
}

arknights