Lagrange inversion theorem 拉格朗日反演

解决的问题

给定F(x), 找G(x) 使得 G(F(x)) = x, 也就是找反函数

这里有一点要用到,但是没有证明的是 $G(F(x)) = x$ 则 $F(G(x)) = x$, 只能说在$F(x)$的值域内$F(G(F(x)) = F(x)$即$F(G(x)) = x$可以证明, 但在值域外不知道如何证明…

结论

如果$F,G$满足

  • 0次项系数均为0: $[x^0]F(x)=[x^0]G(x)=0$
  • 1次项系数均非0: $[x^1]F(x) \neq 0,[x^1]G(x) \neq 0$

那么 $\lbrack x^n \rbrack G(x)^k = \frac{k}{n} \lbrack x^{n-k} \rbrack \left(\frac{x}{F(x)}\right)^n.$

$\lbrack x^n \rbrack G(x) = \frac{1}{n} \lbrack x^{n-1} \rbrack \left(\frac{x}{F(x)}\right)^n.$

证明

辅助lemma: 任何 $F(0)=0$(0次系数为0), $F’(0)\neq 0$ (?存在非0次的非0系数?????), 有对于$\forall k \in \mathbb{Z}$

$\lbrack x^{-1} \rbrack F’(x) F(x)^k = \lbrack k = -1\rbrack$

即$k=-1$时$-1$次系数为$1$,否则$-1$次系数为$0$

证明lemma:

对于$k\neq -1$, 显然求导法则$\displaystyle F’(x) F(x)^k = \frac{\left ( F(x)^{k+1} \right)’}{k+1}$, 所以$\le 0$的x的幂次的系数都是0

对于$k = -1$, $F(x) = \sum_{i>0} a_i x^i$

$\displaystyle \frac{F’(x)}{F(x)} = \frac{a_1+2a_2x+3a_3x^2+\cdots}{x(a_1+a_2x+a_3 x^2+\cdots}= x^{-1} \frac{1 + 2\frac{a_2}{a_1}x + \cdots}{1 + \frac{a_2}{a_1}x + \cdots}$

也就是右侧这个分式除完以后是$1+k_1x+k_2x^2+\cdots$的样子, 因此 -1 次方的系数是 1, lemma 证毕.

这里 注意到 像$\frac{1}{1-x}$这种 对应的 $x^{-1}$的系数是0,因为它虽然写在分母上,但实际是泰勒展开到正,而$\frac{1}{x}$ 是不能在0点泰勒展开成正的?而去查 1/x的泰勒展开,都是在点1的展开


因此$G(F(x))^k = x^k,k\ge 1,k\in \mathbb{Z}$, 的$G$ 满足条件

$(\sum_{i} ([x^i]G^k(x))F(x)^i)’\cdot F’(x) = kx^{k-1}$ ( 同时对$x$求导, 以及按位展开式 $G^k(x)=\sum_{i}([x^i]G^k(x))x^i$

展开$\sum_i i(\lbrack x^i\rbrack G^k(x) ) F(x)^{i-1} F’(x) = kx^{k-1}$, (基本的求导法则 $(ax^i)’ = iax^{i-1}$, 注意到前两项都是常系数 而非生成函数

$\sum_{i} i (\lbrack x^i \rbrack G^k(x)) F^{i-1-n}(x) F’(x) = kx^{k-1}F^{-n}(x).$ ( 同乘上$F^{-n}$, 这里想用lemma, 也就是$[x^{-1}]F^?(x)F’(x)$, 同时最终目标是$[x^?]G(x)$

$\lbrack x^{-1}\rbrack \left(\sum_{i} i (\lbrack x^i \rbrack G^k(x)) F^{i-1-n}(x) F’(x)\right) = \lbrack x^{-1}\rbrack \left(kx^{k-1}F^{-n}(x)\right).$ (提取$-1$次项目的系数

因为左侧$i (\lbrack x^i \rbrack G (x))$ 整个都是系数, 以及上面的lemma, 左侧只有$i-1-n=-1$即$i=n$时, $[x^{-1}]F^{i-1-n}(x)F’(x)$生成系数的内容才为$1$,其它则是$0$,

$n[x^n]G^k(x) = [x^{-1}]kx^{k-1}F^{-n}(x)$

变形一下,就有了最初要证明的$\displaystyle \lbrack x^n \rbrack G(x) = \frac{k}{n} \lbrack x^{n-k} \rbrack \left(\frac{x}{F(x)}\right)^n.$


这玩意 百科上说建立了函数方程和幂级数之间的联系 !!??

使用示例: 卡特兰数Catalan number

$c_0 = 1$

$c_{n+1} = \sum_{i=0}^n c_{i} \cdot c_{n-i}$

令$C$为以卡特兰数 1,1,2,5,14为系数的生成方程

令$F(x) = C(x) - 1$, 保证0次项系数为0, 1次系数非0

那么$F(x) = x(F(x)+1)^2$ , 根据卡特兰数本身推导的定义的到的

令$G(x) = \frac{x}{(x+1)^2}$

那么有$G(F(x)) = \frac{F(x)}{(F(x)+1)^2} = x$

那么$c_n$ 就有了

因此

$\begin{aligned} [x^n]F(x) &= \frac{1}{n} \lbrack x^{n-1} \rbrack \left(\frac{x}{G(x)}\right)^n \\ &= \frac{1}{n} \lbrack x^{n-1} \rbrack (x + 1)^{2n} \\ &= \frac{1}{n} \binom{2n}{n-1} = \frac{1}{n+1} \binom{2n}{n}, \end{aligned}$

相关

https://atcoder.jp/contests/abc222/editorial/2765

abc222ex

拉格朗日反演

超理论坛 拉格朗日反演

洛谷 拉格朗日反演

D

https://atcoder.jp/contests/arc127/tasks/arc127_d

题目大意,(atcoder 的这个题没有“背景”,直接去看原题,就是题意

给等长数组A,B, 对于两两坐标i,j, 计算$min(A_i \oplus A_j, B_i \oplus B_j )$

求所有min的和

范围

$n \leq 250000 , A_i,B_i < 2^{18}$

先做一些简单的分析,N 有点大,$N^2$的话肯定超时, 那么基本范围是$N , N log(N)$ 左右的

$2^{18} = 262144$


min(x,y)+max(x,y) = x+y


一些 $\oplus$ 的常见知识, $(A \oplus B) = (A + B) \bmod 2$, 其中采取 每位不进位允许超过2的任意值, 这个好处是能计数

例如 5 = 101,7=111, 9 = 1001

$(0,1,0,1)+(0,1,1,1)+(1,0,0,1) = (1,2,1,3) $

$(1,2,1,3) \bmod 2 = (1,0,1,1) $

所以$5 \oplus 7 \oplus 9 = 11$, 同时我们知道最高位1个1,2个0


a < b 那么 $a \oplus b$ 的最高位1出现在b

反过来

$a \oplus b$ 最高位出现在b , 那么 a < b

(归纳法易证


$a \oplus (b+c+d+e+\cdots) =$ 后面部分通过上面按位不进位加和统计的每一位 0 个数 或1个数,在看与a对应位是否相等

通过前缀和或者扫描记录当前, 长度为n的数组的两两$\oplus$的和 的时间复杂度O(N) 就能完成


综上,我们可以直接算出 min, 也可以去算max然后拿总和减去min

算法

把所有数看成18位,不足高位补零

我们 关心$(A_i \oplus A_j) \oplus (B_i \oplus B_j) $ 的最高位的bit从哪里来,这样我们就知道哪个最大了

$(A_i \oplus A_j) \oplus (B_i \oplus B_j) = (A_i \oplus B_i) \oplus (A_j \oplus B_j) $


令$C_i = A_i \oplus B_i$

w = 18

把$C_i$分类成两组, 一组是 第$w$位为$1$的,另一组是$w$位为$0$的

每组内,循环这个方法并且w-1


对于组之间的,

$C_i$的$w$位为1, $C_j$的$w$位为0 , 对于高于w位的,因为通过上面的分治,两边保持一致,也就是$\oplus$以后都是0了不需要考虑

对于$w$位为1 里面我们还可以分类为 $A_i$ 的w位是0或1两类

$G_{x,y}$ w位是x,A的w位是y

$G_{1,0},G_{0,0}$, 高位1来自B的$\oplus$

$G_{1,0}, G_{0,1}$, 高位1来自A的$\oplus$

$G_{1,1}, G_{0,0}$, 高位1来自A的$\oplus$

$G_{1,1}, G_{0,1}$, 高位1来自B的$\oplus$

也就是分组的时间复杂度就是$2^L$

下面考虑

$G_{i_0,i_1,i_2,\cdots} $ 与 $G_{j_0,j_1,j_2,\cdots}$ 中 如果高位来自A的$\oplus$

那么贡献为

$ \sum B_{i_?} \oplus B_{j_?}$

也就是,上面想求总和提到过的方法, 通过按位计数,拿着一块求和的单个复杂度为$O(L\cdot size(G))$, 总复杂度为$O(L\cdot N)$


然后注意处理$C_i$ 一致的组,这样的组里面,两两的异或值相同, 也需要上面的按位计数的求和

Code

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
#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
#define MOD 1000000007
#define rep(i,a,n) for (ll i=a;i<n;i++)
#define per(i,a,n) for (ll i=n-1;i>=a;i--)
#define pb push_back

int A[250010];
int B[250010];

const int L = 18;

// 两两xor和 A[i0[?]]
ll xorSum(vector<int> & i0){
if(i0.size() == 0)return 0;
vector<int> bits = vector<int>(20,0);
ll ans = 0;
rep(i,0,i0.size()){
if(i > 0){
rep(b,0,L){
ans += ((A[i0[i]]>>b) % 2)?
((ll)1 << b) * (i - bits[b]):
((ll)1 << b) * (bits[b]);
}
}
rep(b,0,L){
bits[b]+=(A[i0[i]]>>b)%2;
}
}
return ans;
}

// 任意左Arr[i0[?]] xor 任意右Arr[i1[?]] 的和
ll xorSum(vector<int> & i0,vector<int> & i1, int * Arr){
if(i0.size() == 0 || i1.size() == 0)return 0;
vector<int> bits = vector<int>(20,0);
rep(i,0,i0.size()){
rep(b,0,L){
bits[b]+=(Arr[i0[i]]>>b)%2;
}
}
ll ans = 0;
rep(i,0,i1.size()){
rep(b,0,L){
ans += ((Arr[i1[i]]>>b) % 2)?
((ll)1 << b) * (i0.size() - bits[b]):
((ll)1 << b) * (bits[b]);
}
}
return ans;
}

// pwr 以上的位, idxs里面 A[i]^B[i] 两两相等
ll f(int pwr,vector<int> & idxs){
if(pwr < 0){
// (Ai xor Aj) xor (Bi xor Bj) = (Ai xor Bi) xor (Aj xor Bj) = 0
return xorSum(idxs);
}
// groups
vector<int> gC[2] = {{},{}}; // w 位的 C
vector<int> gCA[2][2] = {{{},{}},{{},{}}};// w 位 C 和 A的 零一
for(auto i:idxs){
int C = A[i]^B[i];
gC[(C>>pwr)%2].pb(i);
gCA[(C>>pwr)%2][(A[i]>>pwr)%2].pb(i);
}
ll ans = 0;
rep(i,0,2){
if(gC[i].size() < 2) continue;
ans += f(pwr-1,gC[i]);
}
// bit
rep(bA0,0,2){
rep(bA1,0,2){
ans += (bA0 ^ bA1)?
xorSum(gCA[0][bA0],gCA[1][bA1], B):
xorSum(gCA[0][bA0],gCA[1][bA1], A);
}
}
return ans;
}

int main(){
int n;
cin>>n;
rep(i,0,n){
scanf("%d",A+i);
}
rep(i,0,n){
scanf("%d",B+i);
}
vector<int>idxs;
rep(i,0,n){
idxs.pb(i);
}
printf("%lld",f(L,idxs));
return 0;
}

Ref

https://atcoder.jp/contests/arc127/editorial/2697

题目

https://projecteuler.net/problem=233

题意

$f(N) = $ 和 $(0,0)(0,N),(N,0),(N,N)$ 共圆的整点的个数

例如 $f(10000) = 36$

求 $1 \leq i \leq 10^{11}$ 使得$f(i)=420$ 的所有$i$的和

题解

我目前想到的

对称性

首先 圆上的点 (x,y) 满足

$(x-\frac{N}{2})^2+(y-\frac{N}{2})^2 = \frac{N^2}{2}$

根据对称性

我们只用考虑 $\frac{N}{2} < x < y $ 的个数,然后乘8 加4(对角线的点)

$420 = 52 \cdot 8 +4$

所以这里我们其实 需要求满足上述两个条件的整点 方案为52的即可

奇偶

若 N 为奇数

$(2x-N)^2+(2y-N)^2 = 2N^2$, 要解的方程是$a^2+b^2=2c^2, (0 < a < c < b )$且$a,b,c$均为奇数的形式

若 N 为偶数,令 $N=2k$

$(x-k)^2+(y-k)^2 = 2k^2$, 要解的方程是$a^2+b^2=2c^2$ 的形式

考虑 到$\bmod 4$ 的情况

$ 0^2 = 0 \bmod 4 $

$ 1^2 = 1 \bmod 4 $

$ 2^2 = 0 \bmod 4 $

$ 3^2 = 1 \bmod 4 $

所以 $(a,b,c)$ 同奇偶

因此如果k仍然是偶数,则继续除以2,反复迭代直到c为奇数

最终我们只用考虑 奇数的形式

也得到结论如果N是偶数,那么$f(N) = f(\frac{N}{2})$

公约数

若 $gcd(a,b) = k > 1$,那么显然 $gcd(a,c) = k$

所以$(a,b,c)$同时除以$k$依然成立

同理$(a,c),(b,c)$ 只用考虑 $(a,b,c)$ 两两互质,非互质的是通过互质的倍数得到的

令 $g(c) = $ 使得 $a^2+b^2=2c^2 ,(0 < a < c < b),(a,b,c)$ 两两互质 且 全为奇数(废话,可以不加这条因为是等式的结论) 的方案数

$c$是奇数时 $ f(c) = \sum_{c = 0 \bmod i}g(i)$

变形 与 直角三角形

$a^2+b^2=2c^2, (0 < a < c < b )$且$a,b,c$均为奇数的形式

思路1

$b^2 - 2c^2 = -a^2$ 然后 佩尔方程类似的去解

如果能得到不同a的初始解,剩余的解很容易得到,但是初始解还是只能枚举

思路2

令$d = \frac{a+b}{2}$,$e = \frac{b-a}{2}$

则$d,e$为整数且奇偶互异,

带入 $(d-e)^2+(d+e)^2=2c^2$

即$d^2+e^2=c^2$

这和 直角三角形的基础解是一一对应的,(直角三角形基础解的讲解 见 PE075 的题解)

所以 $k(a,b,c)$ 与 $k(d,e,c)$ 一一对应

所以 问题变成了

c 为直角三角形边 恰好 52 种的方案数

计算方案数 与 方案数性质

根据 直角三角形基础解的性质

若 满足$c = k (m^2+n^2), (m>n>0) , (m+n)=1 \bmod 2$

则 c 是一个直角三角形的斜边

问题是 数量级在$10^{11}$ 虽然 m 和n 都不到$10^6$了,但是基于时间和空间 依然 无法暴力

费马平方和定理

奇质数能表示为两个平方数之和的充分必要条件是该质数被4除余1

https://yexiaorain.github.io/Math/Fermat_s_two_square_theorem/

回到题目

我们之前的结论是:

若 满足$c = k (m^2+n^2), (m > n > 0) , (m+n)=1 \bmod 2$

则 $c$ 是一个直角三角形的斜边, 求有52个解的所有的$c$的和

现在有的理论基础是:

若$a,b$ 互质,则$a^2+b^2$所有质因子都是$4n+1$的形式, 且这些质因子自身唯一分解

反过来$4n+1$的质数的乘积 能表示成 $a^2+b^2$,但是可能$a,b$不互质, 例如$5^4 = 5^2(5^2) = 5^2(3^2+4^2)=15^2+20^2$

所以上面的一个数$c$质因数分解,则 模4余3的质因子和2全在k中,模4余1的因子可以在k中也可以不在k中

$k(m^2+n^2) = c = 2^{c_1} p_1^{a_1} \cdots p_r ^{a_r} q_1^{b_1} \cdots q_s ^{b_s}$

$k$是$2^{c_1}q_1^{b_1} \cdots q_s ^{b_s}$ 的倍数

$p_1^{a_1} \cdots p_r ^{a_r}$是$(m^2+n^2)$ 的倍数

表示方法来源的唯一性, 一个仅由4n+1素数乘积构成的数的平方数分解,全部恰好由其任意一个素数的表示法与剩余部分的表示法生成

我们下面只关心$X = m^2+n^2 = p_1^{d_1} \cdots p_r ^{d_r}$的部分, 这里幂次用d区别了上面的a,表示确定选择的部分,仅仅是上面a的部分约数即可

任意拆一个质因子

$X = (p_1^{d_1} \cdots p_i^{d_i-1}\cdots p_r ^{d_r})\cdot p_i$

已知对于左侧任意一个表示法$a^2+b^2$, 和 右侧的一个唯一的表示法$c^2+d^2$

由费马平方和定理第一条, 可以对$X = (a^2 + b^2)(c^2 + d^2) = (ac+bd)^2+(ad-bc)^2 = (ac-bd)^2+(ad+bc)^2$

也就是 对于给定$(a,b), p $ 能产生考虑正负的16个点,这8个点一组,每组内的点按坐标轴和对角线对撑

$((a,b), p) => (\pm(ac+bd),\pm(ad-bc)) or (\pm(ac-bd),\pm(ad+bc))$ 或交换横纵坐标


反过来,考虑$X = e^2+f^2 = $ 的一个任意的表示法

$X = e^2+f^2 = (p_1^{d_1} \cdots p_i^{d_i-1}\cdots p_r ^{d_r})\cdot (c^2+d^2)$

由上面费马平方和定理第二条

$(p_1^{d_1} \cdots p_i^{d_i-1}\cdots p_r ^{d_r}) = (\frac{ec+fd}{p})^2+(\frac{ed-fc}{p})^2$ 能表示成整数的平方和 或 $= (\frac{ec-fd}{p})^2+(\frac{ed+fc}{p})^2$能表示成平方和的倍数

所以,至少有一个拆解

如果是第一个

$((\frac{ec+fd}{p})^2+(\frac{ed-fc}{p})^2)(c^2+d^2) = (c\frac{ec+fd}{p}+d\frac{ed-fc}{p})^2+(c\frac{ed-fc}{p}-d\frac{ec+fd}{p})^2 = (\frac{e(c^2+d^2)}p)^2+(\frac{-f(c^2+d^2)}p)^2 = e^2+f^2$

第二个同理可证,说明了刚好为逆运算

也就说明了,对于每一个$X$的表示,其所有的可行拆解对应逆运算能得到$X$对应的表示,那么在考虑所有正负的情况下,通过 费马定理第一条,能不漏的计算出$X$的所有表示


然后整理一下值,我们发现$(\pm a,\pm b), (\pm b,\pm a)$ 所产生的16个点不会因为$a,b$的正负和顺序受到影响

所以我们仅去考虑$0 < a < b$, 也就是x正向轴,和y=x正向 之间区域的点作为a,b ,同时产生的点也仅考虑在这个区域的点。

那么有,根据费马定理的第一步,每组$((a,b),p)$ 产生两个新的点

点的计数(TODO

这样保证了点的不漏,每多一个质数,原来的一个点就能产生新的两个点,数量显然是$2^{d_1+d_2+\cdots+d_r-1}$

接下来考虑不重复计数

$((a,b), p = c^2+d^2) => ((ac+bd),(ad-bc)) or ((ac-bd),(ad+bc))$

若 $ac+bd = ad+bc $, 即$(a-b)(c-d) = 0$ 不会成立

若 $ac+bd = | ac - bd | $, 即$bd=0$或$ac =0$ 不会成立

所以 对于一个$((a,b),p)$ 产生的两个点必不相同


考虑不同的$a,b$, $a_1^2+b_1^2 = a_2^2+b_2^2$

$((a_1,b_1),p) => ((a_1c+b_1d),(a_1d-b_1c)) or ((a_1c-b_1d),(a_1d+b_1c))$

$((a_2,b_2),p) => ((a_2c+b_2d),(a_2d-b_2c)) or ((a_2c-b_2d),(a_2d+b_2c))$

若 $a_1c+b_1d = a_2c+b_2d$,

?????????????????????

目前看 相关资料,表示Dummit & Foote, 3rd ed., p 291有相关证明

点的计数

引理定义们

除法: 对于所有$a,b\in \mathbb{Z},b \neq 0$,有唯一$q,r\in \mathbb{Z}, 0\leq r < |b|$ 使得 $a=bq+r$

整除: 对于$a,b\in Z$ $a$是$b$的除数,当存在一个整数$x$, 使得$ax = b$, 写作$a \mid b$,不整除写作$a \nmid b$

恒成立的整除: $1 \mid a,-1 \mid a, a \mid a, -a \mid a$

质数: 仅有上述恒成立的整除关系,且为正,非1

互质: $a,b$ 互质,当它们公共的除数,仅有正负1

最大公约数: $gcd(a,b)$ 所有a和b共有的约数中最大的一个

合数: $a = p_1\cdot p_2 \cdots p_k$, 唯一表示

Ref

Brahmagupta-Fibonacci Identity

Fermat’s theorem on sums of two squares

Fermat Sum of Two Squares Calculator

Fermat’s 4n+1 Theorem

SUM OF TWO SQUARES

F

https://codeforces.com/contest/1556/problem/F

比赛n个队伍$(n\leq 14)$

第i个队伍有个$a_i(1\leq a_i \leq 10^6)$

i队和j队打,分别的获胜概率为$\frac{a_i}{a_i+a_j}$ 和 $\frac{a_j}{a_i+a_j}$

任意两个队之间打且仅打一场

如果 i直接战胜j 或者 间接战胜 i 战胜 ? 战胜 ? 战胜 j ,都算战胜j, (也就意味着 i战胜了j,j也同时战胜了i)

如果一个队战胜了其它所有队,那么它是winnner

求winner个数的期望值($\mod 10^9+7$ 意义下)

题解

我的初步思路

抽象成数据结构

n个点,任意两点之间有且仅有一条有向边,有向边的方向通过上述概率决定,求能走到所有点的点的个数期望值

边数 最大 $\frac{14 \cdot 13}{2} = 91$

如果一点是winner,那么进入该点的都是winner

如果一系列点构成环,那么它们要么同时是 winnner,要么同时不是winnner

单winnner 存在一个点入度为0,其它点入度全大于0,剩余点的边无论怎么排,都是单winner

双winner 不可能, 因为 两点互相能到,任意两点之间直连是仅有一条的单向边,路径上必定有其它的点,那么其它的点也是winner,必定多余2个winner

三winner 同双winner的互相可达原理,和可达路径上的点也是winner的原理,三个winner 构成环,且3个winner以外的边和它们连接全是出度,没有入度,同样剩余点的边无论怎么排都是3winner

同理

m个winner,意味着,这m个两两可达,且对m个以外的边全是出度

期望来了,

每次选m个点,让它们两两可达的概率和所有其它边是出度的概率 * m 再求和

其它边都是出度的概率好求

n=14 枚举点的话,也只有 $2^{14} = 16384$

那么问题变成了,给定m个点,求其内部两两可达的概率

和上面类似的,如果 A -> B, 那么如果B能到 其它所有,无论A的剩余边怎么连都是A也能到其它所有,且显然 B -> A还有一条非直连的简单路径

但这只是保证了A,B两个点,并不是两两可达

剩余

首先官方题解的拆分是一样的, 也是聚集到如何计算 P(winnners)

这里容斥 + 切

也就是 在winnners中选取集合subs, 集合中两两可达

subs 中所有向 $ winnners - subs$ 都是胜利的方案,这样保证了这些都是两两可达的。 这里需要证明能覆盖所有非两两可达,因为如果不是所有两两可达,那么对其中局部两两可达的进行缩点,缩点后,整个拓扑有序的,所以能找到一条切割方案, 让一侧战胜另一侧,源点两两可达(因为所有点之间都有有向路径,不可能存在多余一个的0入度的点,所以至多一个无入度的点,那么这个点缩点前就是subs),所以上面的是充要的

下面证明,不同subs选取的互斥性

对于同一个winnners, subs0 和 subs1 的点选取不同,则至少有一个点,在且仅在其中一个集合中,那么subs0 和 subs1 中的能到subs的点集已经不同了,那么subs0和subs1 不会重复

所以 这里说是容斥,但是实际上,是选取的互斥,又全覆盖的划分

P(集合) = 集合中两两可达的概率

G(集合A,集合B) = 集合A 任何人,战胜 集合B 所有人的概率 (切的概率)

F(集合) = 仅集合内全是winners(集合外全不是winnners)的概率

ALL 所有人

$|集合| = 集合大小$

$答案 = sum( F(winners) * |winnners|)$ 概率乘值,期望公式

$F(winners) = P(winnners) \cdot G(winnners, ALL - winnners)$ 上面的结论,winnners 打败所有非winners

$P(winners) = 1 - sum ( P(sub) \cdot G(sub, winnners - sub) )$ 上面的互斥拆分

G 没啥好说,就是切上

$G(X,Y) = \prod_{\forall x\in X, \forall y \in Y} \frac{a_x}{a_x+a_y}$

公式递推起来似乎可以做了

然而,状态看起来有 $(2^n)^2$ 个,还不包括状态等,超界了

$G(X,Y) = \prod_{\forall x \in X} (\prod_{\forall y \in Y \frac{a_x}{ax+a_y}})$

于是对于每个x,我们可以 $O(2^n)$ 计算,所有的x对于不同的Y,可以 $O(n 2^n)$计算出

那么对于一个给定$G$可以 O(n) 计算出

没太懂,这里$G_{sidefrom,sideto}$ 是怎么回事

代码

这里我想的有问题, 没有正确估计到复杂度,

虽然是枚举 i去算集合,sub 是i的子集,但是 复杂度并不是$O(2^n \cdot 2^n \cdot n)$ , 分别是 i,sub,和G(sub,i-sub);

下面次数也是大概,没有那么准确,比如空和全没有排除

sub 的个数 $\sum_{i=1}^n C_i^n \cdot 2^i$

G内部计算次数 $\sum_{i=1}^n C_i^n \cdot \sum_{j=1}^{i} C_j^i \cdot j = 3^{n-1} \cdot n$, n = 14 时表达式的值为 $3^{13} \cdot 14 = 22320522$

证明见wolframalpha 或下方

当然通过代码cnt++也可以统计

1
2
3
4
n = 14
i = 16383
cross = 4766585
cross ret = 22320522
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
#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
#define MOD 1000000007
#define rep(i,a,n) for (ll i=a;i<n;i++)
#define per(i,a,n) for (ll i=n-1;i>=a;i--)
#define pb push_back
// 末尾零的个数
#define ctz __builtin_ctz
#define popc __builtin_popcount
const int N = 16, N2 = 16400; // 2**14=16384

int n;
int a[N];
ll G[N][N]; // G[i][j] = a[i]/(a[i]+a[j])
ll cr[N][N2]; // [点][bit mask state]
ll prob[N2];

ll PowerMOD(ll a, int n, ll c = 1) {
for (; n; n >>= 1, a = a * a % MOD)
if (n & 1) c = c * a % MOD;
return c;
}

// 点集A战胜点集B,O(n)
int cross(int A, int B) {
ll ret = 1;
// 每次去掉A的最后一个1
for (; A; A &= A - 1) {
ret = ret * cr[ctz(A)][B] % MOD;
}
return ret;
}

int main() {
ll ans = 0;
cin >> n;
int ALL = ~(-1 << n);
rep(i,0,n){
scanf("%d",a+i);
}

// G[i][j] = a[i]/(a[i]+a[j])
// 没有特殊处理 i == j
rep(i,0,n){
rep(j,0,n){
G[i][j] = PowerMOD(a[i] + a[j], MOD - 2, a[i]);
}
}
// O(n * 2^n) , 每个点i,战胜点集合j(bit mask)
// 没有特殊处理 i \in j 之后计算不会使用即可
rep(i,0,n){
cr[i][0] = 1;
rep(j,1,ALL+1){
cr[i][j] = cr[i][j & (j - 1)] * G[i][ctz(j)] % MOD;
}
}
// ans = sum{ |popc(i)| * prob[i] * cross(i,j) }
// winner数 winner内部互连 winnner战胜非winner
//
// prob[i] = 1 - sum { prob(sub) * cross(sub, i - sub) }

// 次数
// n = 14
// i = 16383
// cross = 4766585
// cross ret = 22320522
rep(i,1,ALL+1){
prob[i] = 1;
// 这里对 i 的 子集全枚举
for (int j = i & (i - 1); j; j = (j - 1) & i){
(prob[i] -= prob[j] * cross(j, i - j)) %= MOD;
}
(ans += prob[i] * cross(i, ALL - i) % MOD * popc(i)) %= MOD;
}
printf("%lld\n", (ans+MOD) % MOD );
return 0;
}

$\sum_{i=1}^n C_i^n \cdot \sum_{j=1}^{i} C_j^i \cdot j $

= $\sum_{i=1}^n \frac{n!}{i!(n-i)!} \cdot \sum_{j=1}^{i} \frac{i!}{j!(i-j)!} \cdot j $

= $\sum_{i=1}^n \frac{n!}{i!(n-i)!} \cdot \sum_{j=1}^{i} \frac{(i-1)!}{(j-1)!(i-j)!} \cdot i $

= $n \cdot \sum_{i=1}^n \frac{(n-1)!}{(i-1)!(n-i)!} \cdot \sum_{j=1}^{i} \frac{(i-1)!}{(j-1)!(i-j)!} $

= $n \cdot \sum_{i=0}^{n-1} C_{i}^{n-1} \cdot \sum_{j=0}^{i-1} C_{j}^{i-1} $

= $n \cdot \sum_{i=0}^{n-1} C_{i}^{n-1} \cdot (1+1)^{i-1}$

= $n \cdot \sum_{i=0}^{n-1} C_{i}^{n-1} \cdot 2^{i-1}$

= $n \cdot (1+2)^{n-1}$

= $n \cdot 3^{n-1}$

总结

  1. 看到14,12,20这类数,能一下向bit想,是习惯
  2. 互斥拆分还需要多练习增加经验
  3. 高效内置库,末尾0个数,所有1个数#define ctz __builtin_ctz,#define popc __builtin_popcount
  4. 乘 除 模 同优先级能省掉部分括号 加快代码编写
  5. bit 子集枚举 for (int j = i & (i - 1); j; j = (j - 1) & i){
  6. clang overflow 提示 开关 https://clang.llvm.org/docs/UndefinedBehaviorSanitizer.html

ref

官方题解

yhx-12243 code

B

题意

n <= 4e6, 1e8<mod<1e9, mod is prime, 时间6s,内存128MB

$f(1) = 1$

$f(x) = \sum_{y=1}^{x-1}f(y) + \sum_{z=2}^{x}f(\lfloor \frac{x}{z} \rfloor)$

给定n和mod, 求f(x)%mod

先看时间复杂度

虽然题目给了6s, 但是 这里n已经是4e6了

显然左边可以前缀和

$ f(x) = presumf(x-1) + \sum_{z=2}^{x}f(\lfloor \frac{x}{z} \rfloor)$

左边 O(1)了

右边,注意到$\lfloor \frac{x}{z} \rfloor$ 的取值个数是$O(\sqrt{n})$ 的

这一块即使用了连续段的优化依然整个是$O(n^{1.5})$的复杂度

剩余我没想到的题解部分

考虑 S(x+1) 和 S(x) 的差别

对于加法,多一个 S(x+1-1 = x),

对于除法,多一个 S((x+1)/(x+1) = 1)

对于 i > 1, 可能 S(x+1) 中是 S(i) 而 S(x)中是 S(i-1)

举例

1
2
3
4
5
6
10:的除序列
5,3,2,2,1,1,1,1,1
11:的除序列
5,3,2,2,1,1,1,1,1,1(多)
12:的除序列
6(变),4(变),3(变),2,2(变),1,1,1,1,1,1(多)

$\lfloor \frac{x+1}{i} \rfloor \neq \lfloor \frac{x}{i} \rfloor $ , 原始值变化为$\frac{1}{i} \leq 1$

$\lfloor \frac{x+1}{i} \rfloor = 1 + \lfloor \frac{x}{i} \rfloor $

$\frac{1}{i}$ 的增加在跨越一个整数时必定是整数

$ x+1 = ki $

所以 $i$ 是 x的约数时才会变

约数看成倍数, 计算到x时 -> diff[kx] += f[x]-f[x-1]

直接 变成递推,时间$O(n \log n)$, 空间$O(n)$

Code

on codeforces

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
#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
#define rep(i,a,n) for (int i=a;i<n;i++)

int n;

ll MOD;

// 计算之前的是结果,之后的是差异
ll f[4000010];


int main(){
cin>>n>>MOD;
f[1] = 1;
rep(i,2,n+1){
if(i == 2){
// 特殊 没有第一个 f[i-1]
f[i] = 2;
}else{
(f[i] += f[i-1] + (f[i-1] + f[1]))%=MOD;
}
rep(j,2,n+1){
ll ij=i*j;
if(ij > n)break;
(f[ij] += f[i]-f[i-1])%=MOD;
}
}
printf("%lld\n",(f[n]+MOD)%MOD);
return 0;
}

总结

关键点在 直接换个思路,换成计算差异,然后这里的差异能从 整数部分不同推导到 刚好是约数时会变,从而就能完成效率优化

D

题意

t <= 1e5 组测试

每组

0 <= m < n <= 2e5

所有组m保证 sum{m} <= 2e5, 不对n有相应保证

每次是长度为n的数组(不提供具体值),告诉你m个操作,每次操作x,y (y < x) ,意思是把位置在x移动到y(一次插入排序的操作,也就是 y 是 最大的插入位置)

例如

1 2 3 2

如果把 最后一个2 ,可以插入在2前,也可以是3前,这里只能选择3前, 所以给的(x y)是 (4 3)

问原始序列所有值在[1~n]之间的合法方案数, MOD 998244353

题解

首先 n 不限制,那么可能会有离散化相关的内容

假设序列开始是[a0 a1 a2 …] 那么根据m次操作,有唯一的移动结果,因为每次都指定的x和y

对于移动以后,必定 非严格单调递增

那么问题变成了两个,找到可能一样的值的组(通过能确定的不一样的值做切分),这些组如何得到答案

换个计数法

令c表示,在最终序列中,相邻值“一定不同”的个数(也就是有c个小于号(来源交换推导出的大小关系)),而最终序列中剩余的相邻的位置,可能相同可能不同

那么把这c个一定不同的位置-1,变成了 所有相邻位置可能相同可能小于, 范围由n变为了n-c

也就是 n个数,范围在[1~ n-c]中,非严格单调递增的方案数

也就是 n 个相同的球,放入 n-c 个不同盒子中(可以空盒子)

也就是 n+n-c 个相同的球,放入 n-c 个不同盒子中(不能空盒子)

挡板法,有 2n-c-1 个间隙

显然方案数 C(2n-1-c,n)

也就有了,如果能的到c就能得到答案

计算必定小于个数

离线+倒序

维护一个位置集合S,初始含有1到n

按照插入的倒序,(知道xi是单调递减的)

对于(xi,yi),

1
2
3
4
5
6
7
// 操作前 (反向操作会让 [xi~n]这部分之后再也不会使用)
[1....yi-1][yi...xi-1][xi][xi+1...n]
// 操作后
[1....yi-1][xi][yi...xi-1][xi+1...n]
// 因为我们是反向操作,所以 关注的是这两个
! !
p q

p = S中第yi小的

q = S中第yi+1小的

标记q是 前面相邻被拆入过值的,然后删除p (因为xi是单调递减,xi 到n,在倒序处理过程中,都不再回用到,和把p变成n+1一样效果

而找第k小和移除一个数,这种可以通过线段树,树状数组来完成

效率

注意到不要做成 O(n log n), 要是O(m log n), 一个办法是,首次建立足够大的树状数组/线段树。 每轮询问结束后,再回滚,这样建树 O( max n ), 操作都是 O ( m log n )

这样 算法+效率问题都解决了

代码

狗都不手写平衡树

倒序+树状数组+二分法查找第k大+组合数

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
#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
#define MOD 998244353
#define rep(i,a,n) for (ll i=a;i<n;i++)
#define per(i,a,n) for (ll i=n-1;i>=a;i--)
#define pb push_back
const double pi = acos(-1.0);

const int N = 400000;
ll fac[N+10];

ll pwr(ll v,ll mi){
ll r =1;
while(mi){
if(mi%2) (r*=v)%=MOD;
(v*=v)%=MOD;
mi/=2;
}
return r;
}

ll inv(ll v){
return pwr(v,MOD-2);
}

ll C(ll n,ll m){
return (((fac[n]*inv(fac[m]))%MOD)*inv(fac[n-m]))%MOD;
}

int n,m;

ll ta [2*N+10];

ll lowbit(ll v){
return v & -v;
}

void tinit(){
rep(i,1,2*N+1){
ta[i] = lowbit(i);
}
}

void tadd(int i,int v){
for(;i < 2*N+1;i+=lowbit(i)){
ta[i] += v;
}
}

ll tpre(int i){
ll s = 0;
for(;i!=0;i-=lowbit(i)){
s+=ta[i];
}
return s;
}

ll getk(int k){
ll l = 0;
ll r = N+1;
while(l+1 < r){
int mid = (l+r)/2;
ll p = tpre(mid) ;
if(p>=k)r=mid;
else l=mid;
}
return r;
}



ll getc(vector<int>&ys){
set<int> used ;
vector<int> rm;
per(i,0,m){
int y = ys[i];
used.insert(getk(y+1));
int rmi = getk(y);
rm.pb(rmi);
tadd(rmi,-1);
}
int r = used.size();
per(i,0,m){
tadd(rm[i],1);
}

return r;
}

int main(){
fac[0]=1;
rep(i,1,N+1){
fac[i] = (fac[i-1]*i)%MOD;
}
tinit();
int t;
cin>>t;
while(t--){
vector<int> ys;
scanf("%d %d",&n,&m);
rep(i,0,m){
int x,y;
scanf("%d %d",&x,&y);
ys.pb(y);
}
ll c= getc(ys);
printf("%lld\n",C(2*n-c-1,n));

}
return 0;
}

Ext

pbds

对于 题解提到的 Policy-Based Data Structures 见下面洛谷日报,总之是在std以外,c++的一个扩展库里提供了对应的平衡树实现,可以不用自己去写平衡树

Rope

1
2
#include<ext/rope>
using namespace __gnu_cxx;
1
rope<变量类型>变量名称;

简单来说 rope是个超级string 内部实现平衡树,也就是可以不用上面的所有知识,直接正着做XD, 本意可能是用作巨大文本编辑时使用

crope = rope<char>

https://github.com/gcc-mirror/gcc/blob/master/libstdc%2B%2B-v3/include/ext/rope

复杂度 看https://www.geeksforgeeks.org/stl-ropes-in-c/ 里说 操作都是log,所以总的复杂度应该是 m log n

Ref

官方题解

洛谷日报39期 pbds

gnu online doc pbds

E

题意

无向连通图

保证每点出度大于2

无自环和重边

点<=1000

边<=2000

每个点有 1<= a[i],b[i] <= 1e9

初始值x, 从点1出发, 每次首次经过一个点,需要满足 当前x > a[i],首次经过该点后会 x+=b[i]

可以重复走边走点

额外限制,走完一条边以后不能立刻重走该边, 也就是不能 a->b->a

求能够经过所有点的最小的x

每次100组测试,每次 测试内容的 sum{n}<=1000 sum{m}<=2000

题解

显然,x单调递增,也就是如果我们能超过最大点,则所有点都能走,那么x的下界是 max{最大a点的a-sum{其它b},1点相邻点的最小a}+1

上界 直接最大值+1

所以如果给定一个值,我们能校验其合法性的话, 那么可以二分 [既然是2分了,其实不考虑上下界 直接[1~1e9+1]开始二分也行]

定义set表示已经访问过的集合, 我们找寻两种增广路径

  1. 到set以外的一些点的简单路径,并返回到set中, set中的点a->out->out->out…->set中的点b,我们把这些out的点加入到set中。注意因为首先我们知道set中的点任意两点至少有一条简单路径,现在能通过外部建立一条简单路径,那么必定可以仅在set中任意的行走,走到任意的点, 对于点a点b是同一个点,就更显然,我们有一个环,能在set中走环上任意点,然后即可走到set上任意点

  2. 不一定返回原来的set中,但是能在外部走成环 set中的点->out->out->out->out中出现过的点, 同上加入set中,我们又可达任意的点

以上两种都解决了不能立刻原路返回的限制

上面两种是必要的方法,但是充分性需要证明一下,也就是没有其它的方案

如果从set中走出,无法达成上面两种,则所有出走都是死路,因为连通性和度保证不会因为边的路径限制而是死路,只会因为 不能立刻回头+a[i] 的限制而成死路,因为每走一步,既不成环也不回到set,那么out没走过的必然-1,数量离散有下界,单调递减,必有界,

证明如果全是死路,则当前set无法完成,如果全是死路,则所有路径均不会消耗所有点,且无法扩展,出走后无法行走,如果消耗完所有点,那么说明x比所有a都打,必定可以成环矛盾

再证 每次增广的顺序不会影响,初始set,

对于一个可行的x和可行的增广顺序,显然x单调递增,如果有不同的方案导致死路,那么死路前的set,可以按照合法的增广顺序操作,每一次操作的x值都大于等于 原合法顺序的值,必定可行,所以增广的顺序和初始x值是否可行没有关系,

同时 一个 set的 最终x = sum{b}+初始x

那么问题就变成

对于给定的x,通过上述两种方案增广直到全覆盖,或者全死路为止

O(效率 log(max(a)))

那么 v->out->out->u(out) ,记录到u的时x的值 , 那么如果有另一条路径 到达u,两条路径中大x的值向小x值走动,那么可以和它构成环, 所以dfs即可

代码

暴力bfs过了,没缩点 时间复杂度, 每次 bfs( O(m+n) ), 增广次数O(n), 所以 总时间复杂度 O(n(m+n)log(max(a)))

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
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
#define t3 1000+10
#define t4 10000+10
#define t5 100000+10
#define t6 1000000+10
#define MOD 1000000007
#define rep(i,a,n) for (ll i=a;i<n;i++)
#define per(i,a,n) for (ll i=n-1;i>=a;i--)
#define pb push_back
const double pi = acos(-1.0);

int n,m;

ll a[1010];
ll b[1010];

vector<int> p2[1010];


vector<ll> get_circle(ll x ,vector<ll> & vis, vector<ll> bfs){
vector<ll> cur = vector(1010,(ll)0);
vector<ll> fa = vector(1010,(ll)0);
int st = 0;
while(st < bfs.size()){
int p = bfs[st];
for(auto item:p2[p]){
if(item == fa[p])continue; // 不能立即返回
// item set 内
if(vis[item]){
// p set 内
if(vis[p]){
continue;
}
vector<ll> ret ;
while(p && !vis[p]){
ret.pb(p);
p = fa[p];
}
assert(ret.size());
return ret;
}
// item非set内
if(cur[item] != 0){ // 访问过
// p->item <- old
// 这个比较有必要吗,因为 cur[x] > a[x] 是肯定的, 所以最小值一定在a中出现,必定至少成立一个不等式
// if(cur[item] > a[p] // old->item -> p -> ...
// || max(cur[p],x) > a[item]) { // p->item -> old
vector<ll> ret ;
while(p && !vis[p]){
ret.pb(p);
p = fa[p];
}
p = item;
while(p && !vis[p]){
ret.pb(p);
p = fa[p];
}
assert(ret.size());
return ret;
// }
}else{
// 未访问过 且可达
if( max(cur[p],x) > a[item] ){
bfs.pb(item);
cur[item] = max(cur[p],x) + b[item];
fa[item] = p;
}
}
}
st++;
}
return {};
}

bool test(ll startx){
vector<ll> vis = vector(1010,(ll)0);
vector<ll> s; // set 中的点
s.pb(1);
vis[1] = 1;
int cnt = 1;
while(true){
auto r = get_circle(startx, vis,s);

if(!r.size()){
return false;
}
for(auto item:r){
if(!vis[item]){
// TODO 缩点
vis[item] = 1;
s.pb(item);
cnt++;
startx += b[item];
}
}
if(cnt == n){
return true;
}
}
}


int main(){
int t;
cin>>t;
while(t--){
scanf("%d %d",&n,&m);
rep(i,2,n+1){
scanf("%lld",a+i);
}
rep(i,2,n+1){
scanf("%lld",b+i);
}
rep(i,0,m){
int u,v;
scanf("%d %d",&u,&v);
p2[u].pb(v);
p2[v].pb(u);
}
ll l = 0;
ll r = 1'000'000'001;
while(l+1 < r){
ll mid = (l+r)/2;
if(test(mid)){
// printf("%lld ok\n",mid);
r = mid;
}else{
// printf("%lld not ok\n",mid);
l = mid;
}
}
printf("%lld\n",r);
// clear;
rep(i,1,n+1){
p2[i] = {};
}
}
return 0;
}

0%