D

https://atcoder.jp/contests/arc130/tasks/arc130_d

给一个树,对节点赋值1~n,两两互不相同,

限制:树上的点要么小于它直接相邻所有点的值,要么大于它直接相邻所有点的值

求方案数.

题解

显然任意选一个为根, 大于和方案和小于的方案对称性,可以只算一个然后乘2

也可以,看成二分图,但是是求满足拓扑序的赋值方案数


注意到大于和小于关系,只需要颠倒就能实现,所以这里我们只讨论 根节点大于所有直接相邻的子节点

我们考虑树上DP

状态设计为 ans[i][j] = 节点i在它和它的所有子树节点中,排序从小到大刚好为j的方案数

对于一个选定的节点,有子树T1,T2,T3,T4,…

设它们的根节点为r1,r2,r3,r4,…

设它们的节点个数为s1,s2,s3,s4,…

那么, 我们考虑维护合并这些子树.

因为 我们这里考虑的是 当前节点要大于所有子树的根节点

所以 合并的时候维护 dp[i] = 已合并的子树中,根节点最大值恰好为i的方案数


以T1和T2为例

分成两种, T1的根更大,和T2的根更大(可以对称处理)

下面讨论 T1根更大

设r1在T1 中的位置为i(1-index)

合并后的序列为(i-1个T1中的点 和 j个T2中的点) r1 ((s1-i)个T1中的点 和 (s2-j)个T2中的点)

注意在两个子树分别的内部,顺序已经固定, 所以只有交叉的顺序不定

方案为 $C((i-1)+j,i-1) \cdot C((s1-i)+(s2-j),s1-i) \cdot ans[r1][i]$

再考虑, T2中的r2的位置,因为T1根更大,所以r2只能在前j个中

$dp[i+j] += C(i+j-1,i-1) \cdot C(s1+s2-i-j,s1-i) \cdot ans[r1][i] \cdot (ans[r2][1]+\cdots +ans[r2][j])$

右边的和可以前缀和O(1),

组合数可以预处理递推O(1),

所以

我们通过 循环一遍T1的节点数,再一遍T2的节点数数,就能得到dp, 每个+= 是O(1)

这里总代价n2,而我错误估计以为n3就没写,哭了

乍看上去,所有节点处理一次,每个节点是子树之间合并,个数的乘积是n3复杂度

实际上,考虑贡献, 一次合并子树的两个子树节点个数乘积大小的代价, 相当于,两个子树之间的节点建立一个配对关系.又因为任何两个节点至多建立两次(讨论根更大)关系,所以总代价为$O(n^2)$

代码

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
#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++)

int n; // <=3000
vector<int> p2[3010]; // 树连接关系
// [node i][从小到大第j个] = 方案数
vector<vector<ll> > ans(3010, vector<ll>(1, 0));

ll qC[3010][3010]; // 预处理组合数

void mergeTable(const vector<ll>& arr0,
const vector<ll>& arr1,
vector<ll>& res) {

vector<ll> presum(arr1.size(), arr1[0]);
rep(i, 1, arr1.size()) { presum[i] = (arr1[i] + presum[i - 1]) % MOD; } // 前缀和
// arr0 根更大 , root(arr0) > root(arr1)
// 虽然每次代价是 子树大小之和 x arr0.size(),看起来是三次方,
// 但实际上,从贡献角度思考 总代价n2
rep(i, 0, arr0.size()) {
rep(j, 0, arr1.size()) {
// arr0: i 个 左侧, arr0.size()-i-1 个 右侧
// arr1: j+1 个 左侧, arr1.size()-j-1 个右侧
// res: i+j+1 个左侧, arr0.size()+arr1.size()-i-j-2 个右侧
// presum = presumof arr1
// 就是表达式 加上了MOD
(res[i + j + 1] +=
(((arr0[i] * qC[i + j + 1][i]) % MOD *
qC[arr0.size() + arr1.size() - i - j - 2][arr0.size() - i - 1]) %
MOD * presum[j]) %
MOD) %= MOD;
}
}
}

void dfs(int idx, int fa/*父节点*/, int g/*大于小于关系*/) {
vector<ll> dp = {}; // [最大的子树 根节点 位置] = 方案数
for (auto item : p2[idx]) {
if (item == fa)
continue;
dfs(item, idx, g ^ 1);
if (g) { // 逆向
// 只用一次也不需要恢复了
reverse(ans[item].begin(), ans[item].end());
}
if (!dp.size()) {
dp = ans[item];
} else {
vector<ll> res(dp.size() + ans[item].size(), 0);
mergeTable(dp, ans[item], res); // dp 根更大
mergeTable(ans[item], dp, res); // ans[item] 根更大
dp = res;
}
}
ans[idx] = vector<ll>(dp.size() + 1, 0);
if (dp.size() == 0) {
// 叶子节点
ans[idx][0] = 1;
} else {
ll cnt = 0;
rep(i, 0, dp.size() + 1) { // 实际上也是前缀和, 当前节点 大于 所有子树根节点
ans[idx][i] = cnt;
(cnt += dp[i]) %= MOD;
}
}
if (g) {
reverse(ans[idx].begin(), ans[idx].end());
}
}

int main() {
rep(i, 0, 3001) {
rep(j, 0, i + 1) {
qC[i][j] = j == 0 ? 1 : (qC[i - 1][j] + qC[i - 1][j - 1]) % MOD;
}
}

scanf("%d", &n);
rep(i, 0, n - 1) {
int u, v;
scanf("%d %d", &u, &v);
p2[u].push_back(v);
p2[v].push_back(u);
}
dfs(1, 0, 0);
ll res = 0;
rep(i, 0, n) { (res += ans[1][i]) %= MOD; }
printf("%lld\n", (2 * res) % MOD); // 大于小于对称性 乘2
return 0;
}

submissions

官方题解

总结

  1. 可以当作知识点的: 子树全部 节点数相乘的代价,总代价是$O(n^2)$
  2. 复杂度也可以在纸上不要脑补

E

https://atcoder.jp/contests/arc130/tasks/arc130_e

一个数组,每次选最小值中的一个+1,记录下标。

形成了一个下标数组。

现在给你下标数组,求满足下标数组的字典序最小的原数组

题解

题目样例为例

1 1 4 4 2 1

如果 我们能给他们分组

1 | 1 4 | 4 2 1

那么,很容易知道合法的最小结果

1 (操作前大小为1)| 1 4 (操作前大小为2) | 4 2 1 (操作前大小为3)

我们只关心

1 | 1 4 (操作后所有为3)

所以原数组为

3-2 3-0 3-0 3-1 = 1 3 3 2


所以核心问题变成去找这个层之间的分割线

见下面代码中的注释

last(pos) = 这个位置值相同的值的上一个位置

cur[value] = 位置,辅助计算last的数组

cnt 当前层的个数

mx 上一层至少的结束的位置/当前层至少的起始位置

f(pos)= pos是有效的层级分割值,或者-1

代码

jiangly 大佬的代码 加的注释的

https://atcoder.jp/contests/arc130/submissions/27614814

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

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

// 若按照 最后一次出现的顺序排序1~n
// 最优结果的起始最小值一定是1
// 那么所有最优结果 在所有操作后 一定满足 AN-1<=A1<=A2<=A3<=...<=AN

int main() {
int N, K;
scanf("%d %d", &N, &K);

vector<int> i(K);
rep(k, 0, K) {
scanf("%d", &i[k]);
i[k]--; // (0-index)
}

vector<int> last(K); // 这个位置上同一个数上一次出现的位置 last[pos] = oldpos
vector<int> cur(N, -1); // 最后一次出现的位置 cur[value] = pos
rep(k, 0, K) {
last[k] = cur[i[k]];
cur[i[k]] = k;
}

vector<int> f(K + 1, 0); // 相当于找有效分割线 = 分割线的层级,从0开始
int cnt = 0;
int mx = -1;
rep(k, 0, K) {
mx = max(mx, last[k]); // 更新最近上一组数 的最大下标
if (last[k] == -1) { // 当前值首次出现
cnt++; // 属于组的个数+1
}
// -1 无效值
// 最大位置 和 当前位置之间个数小于 cnt 则 无效
// 上一个位置无效则 无效,
// 因为如果当前是分割线,那么当前到上一个分割线的距离一定恰好是cnt
f[k + 1] =
(mx >= k + 1 - cnt || f[k + 1 - cnt] == -1) ? -1 : 1 + f[k + 1 - cnt];
}

int c = -1; // 分割最小层数
int len = -1; //
rep(k, mx + 1,
K + 1) { // 对于最后一组分割,最后一组里各不相同,且组长度 <= cnt
// 上一个分割线有效
// 尚无合法的,或当前层数小于等于 上一个方案
if (f[k] >= 0 && (c == -1 || c >= f[k])) {
c = f[k]; // 更新层数
len = k; // 最后一组的起点位置
}
}
// 没有有效的最后一组
if (c == -1) {
cout << "-1\n";
return 0;
}

// 从最后的结果反推初始值,
// 这里排除了最后一组,所以在最后一组之前,所有值是相等的
vector<int> A(N, c + 1); // 最后的结果
rep(k, 0, len) { A[i[k]]--; }
rep(i, 0, N) { printf("%d%c", A[i], " \n"[i == N - 1]); }

return 0;
}

总结

跳表什么都会,但是能把细节分析出来,感觉是熟练度而不是什么缺失的知识点。属于这次学会了,下次可能还是不会的题。

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

0%