Link
https://www.luogu.org/problemnew/show/P4714
Description
给定
n
,
k
n,k
n,k ，计算
∑
d
k
∣
n
∑
d
k
−
1
∣
d
k
⋯
∑
d
4
∣
d
5
∑
d
3
∣
d
4
∑
d
2
∣
d
3
∑
d
1
∣
d
2
1
\sum\limits_{d_k|n}\sum\limits_{d_{k-1}|d_k}\cdots\sum\limits_{d_4|d_5}\sum\limits_{d_3|d_4}\sum\limits_{d_2|d_3}\sum\limits_{d_1|d_2}1
dk∣n∑dk−1∣dk∑⋯d4∣d5∑d3∣d4∑d2∣d3∑d1∣d2∑1
Analysis
「你们竞赛就学这个啊？太简单了吧。」
这个式子很丑啊 你又不好化要是好化那直接a了
这种情况下，你要不考虑从实际意义分析，要不考虑从小情况开始
很明显每多一层
∑
\sum
∑ 就会有一些东西被重复计数
是 什 么 呢 ？
举例吧。
2
2
3
3
5
2^23^35
22335
1
+
2
+
3
+
5
+
2
2
+
23
+
25
+
3
2
+
35
+
2
2
3
+
2
2
5
+
3
3
+
3
2
5
+
2
2
3
2
+
2
2
35
+
3
3
5
+
2
2
3
2
5
1+2+3+5+2^2+23+25+3^2+35+2^23+2^25+3^3+3^25+2^23^2+2^235+3^35+2^23^25
1+2+3+5+22+23+25+32+35+223+225+33+325+2232+2235+335+22325
推起来感觉很生成函数
(
x
2
+
x
+
1
)
(
x
3
+
x
2
+
x
+
1
)
(
x
+
1
)
(x^2+x+1)(x^3+x^2+x+1)(x+1)
(x2+x+1)(x3+x2+x+1)(x+1) ……（虽然实际上算的并不是这个——
3
⋅
4
⋅
1
3\cdot4\cdot1
3⋅4⋅1
多套一层
∑
\sum
∑
(
3
+
2
+
1
)
(
4
+
3
+
2
+
1
)
⋅
1
(3+2+1)(4+3+2+1)\cdot1
(3+2+1)(4+3+2+1)⋅1
(
3
+
2
+
1
+
2
+
1
)
(
4
+
3
+
2
+
1
+
3
+
2
+
1
+
2
+
1
+
1
)
⋅
1
(3+2+1+2+1)(4+3+2+1+3+2+1+2+1+1)\cdot1
(3+2+1+2+1)(4+3+2+1+3+2+1+2+1+1)⋅1
⋯
\cdots
⋯
我们考虑
f
0
(
4
)
f_0(4)
f0(4) → 4 + 3 + 2 + 1
f
1
(
4
)
f_1(4)
f1(4) → (4+3+2+1) + (3+2+1) + (2+1) + 1
¿
f
k
+
1
(
x
)
=
∑
i
=
1
x
f
k
(
i
)
f_{k+1}(x)=\sum\limits_{i=1}^{x}f_k(i)
fk+1(x)=i=1∑xfk(i)
发现是一种
k
k
k 次前缀和的形式 ← 这是一个经典问题
所以一个
p
i
a
i
{p_i}^{a_i}
piai被求和
k
k
k 次之后产生的贡献是一个
k
k
k 次前缀和
k
k
k 次前缀和的话，任意某个位置上的结果都是组合数
如果忘记具体的形式，画转移矩阵弄幂或者找规律打表都可以直接搞出答案
（貌似也可以直接用矩阵加速（反正都逃不过大数分解。。。
如果你没推过
好吧，怎么推？上面那个式子看着很可以搞事情啊
要求的东西是
f
k
(
a
i
)
f_k(a_i)
fk(ai)
考虑
f
−
2
(
4
)
f_{-2}(4)
f−2(4) → 1
f
−
1
(
4
)
f_{-1}(4)
f−1(4) → 4
这种情况下我们可以看作有一堆点从
f
−
2
f_{-2}
f−2 出发
求到达
f
k
(
a
i
)
f_k(a_i)
fk(ai) 的方案数。
我们假设
f
k
(
a
i
)
f_k(a_i)
fk(ai) 在某个诡异的（？）坐标系里的坐标
(
k
,
a
i
)
(k,a_i)
(k,ai)
x
=
k
,
y
=
a
i
x=k,y=a_i
x=k,y=ai
那么。
每次
x
x
x 必定 +1
每次
y
y
y 必定不变小
所以。相当于求
∑
j
Δ
y
j
=
a
i
\sum_j\Delta y_j=a_i
∑jΔyj=ai
经典问题，用插板就可以解决。
更简单地
n
=
∏
p
i
a
i
n=\prod{p_i}^{a_i}
n=∏piai
I
k
(
n
)
=
∏
(
a
i
+
k
+
1
a
i
)
=
∏
(
a
i
+
k
+
1
k
+
1
)
=
∏
(
k
+
2
)
⋯
(
a
i
+
k
+
1
)
1
⋯
a
i
I^k(n)=\prod\binom{a_i+k+1}{a_i}=\prod{a_i+k+1\choose k+1}=\prod\frac{(k+2)\cdots(a_i+k+1)}{1\cdots a_i}
Ik(n)=∏(aiai+k+1)=∏(k+1ai+k+1)=∏1⋯ai(k+2)⋯(ai+k+1)
啥意思？
k
=
0
k=0
k=0 的时候很显然，就是在
a
i
+
1
a_i+1
ai+1 个
p
i
x
{p_i}^x
pix 里面选一个啊。
k
=
−
1
k=-1
k=−1 就更不用说了。
k
>
0
k>0
k>0 呢？相当于是在
a
i
+
1
a_i+1
ai+1 个空里插
k
+
1
k+1
k+1 个板（空代表0到n次幂，插一个空就选了这个幂）
你问我怎么是插板？首先插板插下去的是无序的，所以可以一律当成是从（幂）大到小插进去的
这样的话就相当于（可重复地）插这么些个板板，一种插板方案对应了展开之后的一个因子
以上。
举个例子吧。
1
[
]
p
[
]
p
2
[
]
p
3
[
]
1[\quad]p[\quad]p^2[\quad]p^3[\quad]
1[]p[]p2[]p3[]
你随便插一下
1
[
]
p
[
  
∣
    
]
p
2
[
  
∣
∣
  
]
p
3
[
]
1[\quad]p[\;|\;\;]p^2[\;||\;]p^3[\quad]
1[]p[∣]p2[∣∣]p3[]
就相当于是：
第一次展开
p
3
p^3
p3 的因子
p
2
p^2
p2
第二次展开
p
2
p^2
p2 的因子
p
2
p^2
p2
第三次展开
p
2
p^2
p2 的因子
p
p
p
这就是一种方案。
代码就比较精污了，你不得不上Pollard-ρ+Miller-Rabbin来分解N
回忆一下啊
Miller-Rabbin ①
a
ϕ
(
p
)
≡
1
(
m
o
d
p
)
a^{\phi(p)}\equiv1\pmod{p}
aϕ(p)≡1(modp) 于是做费马素性检测
a
p
−
1
≡
1
a^{p-1}\equiv1
ap−1≡1
② 二次探测
a
2
≡
1
a^2\equiv1
a2≡1 或者
a
2
≡
−
1
a^2\equiv-1
a2≡−1
选取前 ⑨ 个素数（~23） 作为费马素性检测的
a
a
a 即可。
Pollard-ρ：忘了
利用
x
2
+
α
x^2+\alpha
x2+α 生成伪随机数列
a
n
a_n
an
利用循环节碰撞最小素因子
p
p
p （没办法搞
2
2
2 所以得单独做）
客观存在数列
a
n
m
o
d
  
p
a_n\mod p
anmodp 由生日悖论，循环节期望出现位置
O
(
N
)
O(\sqrt{N})
O(N)
具体：用两个指针
x
x
x 和
2
x
2x
2x 在
a
∞
a_\infty
a∞ 上走
等一个
g
c
d
(
a
b
s
(
a
x
−
a
2
x
)
,
P
)
≠
1
gcd(abs(a_x-a_{2x}),P)\ne1
gcd(abs(ax−a2x),P)̸=1 且
≠
P
\ne P
̸=P 此时找到
P
P
P 的一个（有用哒）因子
按步批量处理gcd加速。
实际上并不需要严格分解成多个素数幂的乘积。尽量偷懒嗷
另外还可以先用
1
0
6
10^6
106 以内的素数筛
n
n
n 然后剩下的一定是
1
  
o
r
  
p
  
o
r
  
p
q
  
o
r
  
p
2
1\;or\;p\;or\;pq\;or\;p^2
1orporpqorp2
然后……搞一搞？？玄学
话说我一开始写的超假的miller rabbin居然只挂了一个点。。。我这个pollardrho好像也有一点假
我的 码力 怎么 这么 差 啊 啊 啊
下面是 仍然可能有锅的代码
这个选手的代码可以拿去当wc题让人找错 论写一道题出几十个bug是什么感受
提醒一下，long double乘法取模最后要%p+p%p
不要问我为什么，虽然我觉着理论上它不会出负数。。。不过既然这是事实
我就煤办法腊（）
#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<cctype>
#include<ctime>
using namespace std;
const long long p = 998244353;
const int MAXN = 100;
int Id[MAXN];
long long n, k, Pi[MAXN], Ai[MAXN], Pk[MAXN], Ak[MAXN];
int Prime[9] = {2, 3, 5, 7, 11, 13, 17, 19, 23};
long long gcd(const long long& a, const long long& b)
{
return !b?a:gcd(b,a%b);
}
long long Mul(long long a, long long b, const long long& p)
{
static long long t;
t = (long double) a * b / p;
return ((a * b - t * p) % p + p) % p;
}
long long qpowm(long long a, long long b, const long long& p)
{
if (b <= 0) return 1;
long long ret = 1;
while (b)
{
if (b & 1) ret = Mul(ret, a, p);
a = Mul(a, a, p);
b >>= 1;
}
return ret;
}
long long PowList[100];
bool Miller(const long long& n)
{
if (n == 2 || n == 3 || n == 5 || n == 7 || n == 11 || n == 13 || n == 17 || n == 19 || n == 23) return 1;
if (n < 23) return 0;
register bool flag;
register long long sub = n - 1;
while (!(sub&1)) sub >>= 1;
for (register long long t, gkd, i = 0; i < 9; ++i)
{
gkd = sub; PowList[0] = 0; t = qpowm(Prime[i], gkd, n);
while (gkd <= n - 1)
{
PowList[++PowList[0]] = t;
t = Mul(t, t, n);
gkd <<= 1;
}
if (PowList[PowList[0]] != 1) return 0;
for (register int j = PowList[0] - 1; j >= 1; --j)
{
if (PowList[j] == n - 1) break;
if (PowList[j] != 1) return 0;
}
}
return 1;
}
long long Pollard_Find(const long long& n)
{
register long long alph = 1, x, y, LastX, LastY, Step = pow(n, 0.1), Prod;
while (true)
{
x = 2; y = 2;
while (true)
{
LastX = x; LastY = y; Prod = 1;
for (register long long i = 1; i <= Step; ++i)
{
x = Mul(x, x, n); x += alph; if (x >= n) x -= n;
y = Mul(y, y, n); y += alph; if (y >= n) y -= n;
y = Mul(y, y, n); y += alph; if (y >= n) y -= n;
Prod = Mul(Prod, abs(x-y), n);
}
Prod = gcd(Prod, n);
if (Prod == 1) continue;
if (Prod != n) return Prod;
x = LastX; y = LastY;
for (register long long i = 1; i <= Step; ++i)
{
x = Mul(x, x, n); x += alph; if (x >= n) x -= n;
y = Mul(y, y, n); y += alph; if (y >= n) y -= n;
y = Mul(y, y, n); y += alph; if (y >= n) y -= n;
Prod = gcd(n, abs(x-y));
if (Prod == n) break;
if (Prod != 1) return Prod;
}
break;
}
++alph;
}
}
void Pollard_Rho(const long long& n)
{
if (n <= 1) return;
if (Miller(n)) { Pi[++Pi[0]] = n, ++Ai[Pi[0]]; return;}
long long k = Pollard_Find(n);
Pollard_Rho(k); Pollard_Rho(n/k);
}
inline bool cmp(const int& a, const int& b)
{
return Pi[a] < Pi[b];
}
int main()
{
scanf("%lld%lld", &n, &k);
if (!(n&1))
{
Pi[++Pi[0]] = 2, Ai[1] = 1; n >>= 1;
while (!(n&1)) n>>=1, ++Ai[Pi[0]];
}
Pollard_Rho(n);
for (register int i = 1; i <= Pi[0]; ++i)
{
Id[i] = i;
}
sort(Id + 1, Id + 1 + Pi[0], cmp);
if (Pi[0]) Pk[++Pk[0]]=Pi[Id[1]], Ak[1]=Ai[Id[1]];
for (register int i = 2; i <= Pi[0]; ++i)
{
if (Pi[Id[i]] == Pi[Id[i-1]]) ++Ak[Pk[0]];
else Pk[++Pk[0]] = Pi[Id[i]], Ak[Pk[0]] = Ai[Id[i]];
}
long long Ans = 1;
for (register int i = 1; i <= Pk[0]; ++i)
{
for (register long long x = 1; x <= Ak[i]; ++x)
{
Ans = Mul(Ans, k + 1 + x, p);
Ans = Mul(Ans, qpowm(x, p - 2, p), p);
}
}
printf("%lld", Ans);
return 0;
}
328
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
