AK F.*ing leetcode 流浪计划之费马小定理与组合数取模

欢迎关注更多精彩
关注我，学习常用算法与数据结构，一题多解，降维打击。

费马小定理与证明

参考 https://zhuanlan.zhihu.com/p/594859227

费马小定理：如果p是一个质数，而正整数a不是p的倍数，那么a^(p-1)≡1（mod p）。

证明：

引理1 S={1a, 2a, 3a, …, (p-1)a}，S里所有数都不是p的倍数

由于 1到p-1都小于p, a也不能被p整除，说明a, 1到p-1 这些数都没有因子p。

由此引理1得证明。

引理2 S中所有元素对p取模不为0，且正好为集合 N = {1,2,3,…, p-1}

由引理1可知，对于任意小于p的两个不相同整数, i, j, ia%p != ja%p !=0。

假设 ia%p = ja%p，则 i%p * a%p = j%p * a%p, 得出i=j，与条件矛盾。

所以，对于S%p所有数字不相同，总个数为p-1个，由鸽巢原理可知，引理2正确。

根据引理2，可得到 $\displaystyle \prod_S \%p = \prod_N \%p$

两边整理得到 $a^{p-1}(p-1)!\%p=(p-1)!\%p$

由于gcd(p, (p-1)!)=1, 得到 a^(p-1)%p=1。

逆元

定义

已知整数a,p , a与p互质，如果存在一个数c<p使得 a*c%p=1, 则称c为a在模p下的逆元，记c=a^-1 %p;
通俗理解就是求(1/a)%p。

用费马小定理求逆元

a^(p-1)≡1（mod p）

可知 a*a^(p-2)%p=1

从定义可知 a在模p下的逆元为 a^(p-2)%p, 一般p是比较大的质数，可以使用快速幂求解。

typedef long long lld;lld MOD = 998244353;const int N = 1e6 + 10;// a的逆元为 fastmi(a, p-2)
lld fastmi(lld a, lld n) {lld ans = 1;while (n) {if (n & 1)ans = ans*a%MOD;a *= a;a %= MOD;n >>= 1;}return ans;
}

组合数取模

已经知p是一个比较大的质数（超过10的7次），求组合数C(n, m) , 0<=m<=n。

打表法

当n<=1000时，可以使用扬辉三角打表法。
用一个下三角矩阵存储组合数结果。
利用公式C[i][j] =( C[i - 1][j-1] + C[i - 1][j])%p，求解;
在这里插入图片描述


lld C[1010][1010];void initC() {C[0][0] = 1;C[1][0] = C[1][1] = 1;for (int i = 2; i < 1010; i++) {C[i][0] = C[i][i] = 1;for (int j = 1; j < i; ++j) C[i][j] =( C[i - 1][j-1] + C[i - 1][j])%MOD;}
}

逆元法

$C_n^m=\frac{A_n^m}{m!} = \frac{n!}{m!(n-m)!}$

从上式可以看出当m不是很大时，m!的逆元是可以在O(n)log§时间内求出来的。

当n不大时也可以直接迭代得到n!%p。

如果n-m不大也可以直接求得 $A_n^m%p$


lld inv[N];
// 计算n! 的逆元, inv[n]=(1/n!)%p = (n!)^(p-2)%p
void initInv() {lld sum = 1;inv[0] = 1;for (lld i = 1; i < N; ++i) {sum *= i;sum %= MOD;inv[i] = fastmi(sum, MOD - 2);}
}// 计算C(n, m) = n!/(n-m)!%p * inv[m]%p
lld C(lld n, lld m) {lld sum = 1;for (int i = n; i > (n - m); --i) sum = sum * i % MOD;return sum * inv[m] % MOD;
}

题目一

https://codeforces.com/contest/1967/problem/C

题目大意

题目基础相关知识：树状数组

定义lowbit(x) 是由x二进制最低位的1和后面的0组成的数字，例如，lowbit(12)=4, lowbit(8)=8.
取K值
arr是一个长度为n数组。

假如一个长度为n的数组sum 满足 $sum_i=\left \{\displaystyle \sum_{j=i-lowbit(i)+1}^{i}arr[j]\right \}mod\ 998244353$ ，则称sum为arr的二叉索引树。定义sum=f(arr)。

下图可以解释上述定义，红线为累加给上级关系。

树状数组示意图

给定一个数组a, 整数k，

$f^k(a)=\left \{\begin{array}{l}f(a) \ \ \ k=1\\f(f^{k-1}(a))\ \ \ k>1\end{array} \right .$

问题：给定结果sum,k 求最初始的数组arr。
n 为arr长度, 1≤n≤2⋅10⁵,1≤k≤10⁹

问题解析

既然是累加关系，那么可以尝试找到累加系数关系。之后就可以使用高斯消元法求解每个未知数。

设c[i]为sumi的累加系数。

$sum_i=\left \{\displaystyle \sum_{j=i-lowbit(i)+1}^{i}arr[j]*c[i][j]\right \}mod\ 998244353$

假设已经知道系数，我们可以从低到高枚举i, 然后把用sumi把相关联的上级减掉，剩下的sum数组就是答案。

下面来寻找系数与k的关系

手动模拟一下：

假设现在数组长度为8。
当k=1时，
系数c如下
在这里插入图片描述
横坐标为arr[i], 纵坐标为sum[i]。

k=2, 在上面基础上再模拟一遍
在这里插入图片描述
k=3，在上面基础上再模拟一遍

上面我们只看1，2，4，8行（从树上看是连续上级关系），

从上表观察发现每次K增加1，相应系数就是把k-1对应的系数全部加起来。
例如k=3时，c[8]1 = (c[8][1] + c[4][1]+c[2][1]+c[1][1])(k=2)

经过简单的分析可以得到以下一个表。

在这里插入图片描述

将上表命名为mat, 横坐标为k-1, 纵坐标为当前数字与目标数字的层数相。
可以看出下行就是上一行的前缀和。

至此，我们已经找到了k与系数的关系，但是k太大了，直接计算或者打表都不可能。

对上表进行旋转，可以得到扬辉三角，

在这里插入图片描述
也就是说mat[i][j] 与组合数C(n, m)存在一个关系。
通过线性变换可以得到 mat[i][j] = C(i+j, i)

实现细节

根据上面的推导

sum[1] = arr[1];
sum[2] = c[1][1]*arr[1] + arr[2];
sum[3] = arr[3];
sum[4] = c[4][1]*arr[1] + c[4][2]*arr[2] + c[4][3]*arr[3] + arr[4];
sum[5] = arr[5];
sum[6] = c[6][5]*arr[5] + arr[6];
sum[7] = arr[7];
sum[8] = c[8][1]*arr[1] + c[8][2]*arr[2] + c[8][3]*arr[3] +c[8][4]*arr[4] + c[8][5]*arr[5] + c[8][6]*arr[6] + c[8][7]*arr[7] + arr[8];
…

d = dep[i]-dep[j], 在二叉树上的深度
$c[i][j] = mat[d][k-1]=C(d+k-1, d) = inv[d]*A_{d+k-1}^d$

可从i=1开始用c[][]*sum[i] 去减掉所有上级树上加成项。
相当于d每次会加1，
$mat[d+1][k-1]=C(d+1+k-1, d+1) = inv[d+1]*A_{d+1+k-1}^{d+1}$
inv已经求好了。

$A_{d+1+k-1}^{d+1}=A_{d+k-1}^d * (d+1+k-1)$
根据上述递推式，可以一边递推一边减掉加成项。
这样可以平均O(1)求出系数c.

代码

#include <iostream>
#include <stdio.h>
#include <queue>
#include <string.h>
#include <stack>
#include <vector>
#include <map>
#include <algorithm>
#include <cmath>using namespace std;typedef long long lld;lld MOD = 998244353;const int N = 1e6 + 10;lld fastmi(lld a, lld n) {lld ans = 1;while (n) {if (n & 1)ans = ans * a % MOD;a *= a;a %= MOD;n >>= 1;}return ans;
}lld inv[N];void initInv() {lld sum = 1;inv[0] = 1;for (lld i = 1; i < N; ++i) {sum *= i;sum %= MOD;inv[i] = fastmi(sum, MOD - 2);// if (i < 10)cout <<sum<<", "<< inv[i] << ", " << (inv[i]*sum%MOD) << endl;}
}int lowbit(int x) {return x & (-x);
}lld arr[N];lld C[1010][1010];void initC() {C[0][0] = 1;C[1][0] = C[1][1] = 1;for (int i = 2; i < 1010; i++) {C[i][0] = C[i][i] = 1;for (int j = 1; j < i; ++j) C[i][j] = (C[i - 1][j - 1] + C[i - 1][j]) % MOD;}
}void solve() {initInv();initC();int t;cin >> t;while (t--) {lld n, k;cin >> n >> k;for (int i = 1; i <= n; ++i)cin >> arr[i];for (int i = 1; i <= n; ++i) {lld isum = 1;for (int j = i + lowbit(i), d = 1; j <= n; j += lowbit(j), d++) {isum = isum * ((k - 1 + d) % MOD) % MOD;//cout << j<<", "<< d<<", "<< isum * inv[d] % MOD<<endl;arr[j] -= isum * arr[i] % MOD * inv[d] % MOD;arr[j] %= MOD;//cout << "mat" << k - 1 << ", " << d << endl;//cout << "C" << k - 1 +d<< ", " << d<< "=" << isum * inv[d] % MOD<< endl;/*arr[j] -= C[k-1+d][d] * arr[i] % MOD;arr[j] %= MOD;*///cout << "C" << k - 1 + d << ", " << d << "=" << C[k - 1 + d][d] << endl;}}for (int i = 1; i <= n; ++i) {if (arr[i] < 0)arr[i] += MOD;cout << arr[i] << " ";}cout << endl;}
}int main() {solve();return 0;
}
/*
2
8 100
1 2 1 4 1 2 1 8
6 7
1 4 3 17 5 162
8 1
1 2 1 4 1 2 1 8
6 2
1 4 3 17 5 16*/lld CC(lld n, lld m) {lld sum = 1;for (int i = n; i > (n - m); --i) sum = sum * i % MOD;return sum * inv[m] % MOD;
}