bqlp的数学模板笔记
数学
数论
素数
素数计数函数：小于或等于\(x\)的素数的个数，用\(\pi(x)\)表示。随着\(x\)的增大，有这样的近似结果：\(\pi(x)\sim\frac{x}{ln(x)}\)
素数判定：
bool isPrime(a) {
if (a<2) return 0;
for (int i=2; i*i<=a; i++)
if (a%i==0) return 0;
return 1;
}
素数筛：
\(O(n\log\log n)\)
bool vis[maxn];
int pr[maxn], tot;
void init() {
int n=100000;
int m=sqrt(n+0.5);
memset(vis, 0, sizeof(vis));
for(int i=2; i<=m; i++) if(!vis[i])
for(int j=i*i; j<=n; j+=i) vis[j]=1;
for(int i=2; i<=n; i++) {
if(!vis[i]) pr[tot++]=i;
}
}
欧拉筛法，优化：取消合数的重复标记。\(O(n)\)
//phi, MAXN, vis, pri, cnt
void init() {
phi[1]=1;
for (int i=2; i<MAXN; i++) {
if(!vis[i]) {
phi[i]=i-1;
pri[cnt++]=i;
}
for(int j=0; j<cnt; j++) {
if (1ll*i*pri[j]>=MAXN) break;
vis[i*pri[j]]=1;
if (i%pri[j]) {
phi[i*pri[j]]=phi[i]*(pri[j]-1);
} else {
phi[i*pri[j]]=phi[i]*pri[j];
break;
}
}
}
}
求大数小区间\((a,b)\)间素数：
int cal(ll a, ll b){
memset(v, 0, sizeof(v));
for(int i=0; i<tot; i++){
if(pr[i]>b) break;
ll x=a/pr[i]+(a%pr[i]!=0);
if(x<=1) x=1;//x=2,就不包含本身
for(ll j=x*pr[i]; j<=b; j+=pr[i])
if(j>=a) v[j-a]=1;
}
int ans=0;
for(ll i=0; i<=b-a; i++){
if(!v[i]){
ans++;
}
}
return ans;
}
素性探测
检测\(n\)是否为素数。
Fermat 素性测试
使用费马小定理，在\([2,n-1]\)中不断取随机数\(a\)，检验\(a^{n-1}\equiv1\pmod n\)不成立时，一定不为素数。
对所有\(a\)成立时，可能为素数。
卡迈克尔数
对所有\([2,n-1]\)中的\(a\)都满足\(a^{n-1}\equiv1\pmod n\)，但\(n\)为合数。
eg：\(561=3\times11\times17\)
性质：若\(n\)为卡迈克尔数，则\(m=2^n-1\)也是卡迈克尔数，所以有无穷多个。
卡迈克尔数的形式一定不为\(p^e\)。
二次探测定理
原理：若\(p\)为奇素数，则\(x^2\equiv1\pmod p\)的解为\(x=1\)和\(x=p-1\)。
实现
分解\(n-1=2^t\times u\)，不断对\(u\)进行平方\(\dots\)
再极端的情况也只会把合数判成素数，不会把素数判成合数。
常用大素数：
18位素数：154590409516822759
19位素数：2305843009213693951 (梅森素数)
19位素数：4384957924686954497
一些伪素数（合数）：
开头中间结束各选了3个，
561
1105
1729
426821473
429553345
434330401
2851896013395343201
3589102252820237401
4954039956700380001
抄完模板后，以上用于验证模板抄错没。
#include<bits/stdc++.h>
using namespace std;
#define ll long long
ll mul(ll a,ll b,ll mod){
a%=mod;
b%=mod;
ll c=(long double)a*b/mod;
ll ans=a*b-c*mod;
return (ans%mod+mod)%mod;
}
ll pow_mod(ll x,ll n,ll mod){
ll res=1;
while(n){
if(n&1)
res=mul(res,x,mod);
x=mul(x,x,mod);
n>>=1;
}
return (res+mod)%mod;
}
bool Miller_Rabbin(ll a,ll n){
ll s=n-1,r=0;
while((s&1)==0){
s>>=1;r++;
}
ll k=pow_mod(a,s,n);
if(k==1)return true;
for(int i=0;i<r;i++,k=k*k%n){
if(k==n-1)return true;
}
return false;
}
bool isprime(ll n){
if(n==1) return false;
ll times=7;
ll prime[100]={2,3,5,7,11,233,331};
for(int i=0;i<times;i++){
if(n==prime[i])return true;
if(Miller_Rabbin(prime[i],n)==false)return false;
}
return true;
}
ll n;
int main(){
int T;
cin>>T;
while(T--){
scanf("%lld", &n);
if(isprime(n)) printf("yes\n");
else printf("no\n");
}
}
miller-rabin+Pollard_rho(kuangbin模板)
输入是\(<2^{63}\)
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define S 20//随即判定的次数，S越大出错率
ll mult_mod(ll a,ll b,ll c)//a*b%c
{
a%=c;
b%=c;
ll ret=0;
while(b)
{
if(b&1){ret+=a;ret%=c;}
a<<=1;
if(a>=c)a%=c;
b>>=1;
}
return ret;
}
ll pow_mod(ll x,ll n,ll mod)//x^n%mod
{
if(n==1)return x%mod;
x%=mod;
ll tmp=x;
ll ret=1;
while(n)
{
if(n&1) ret=mult_mod(ret,tmp,mod);
tmp=mult_mod(tmp,tmp,mod);
n>>=1;
}
return ret;
}
//一定是合数返回1，不一定返回0
bool check(ll a,ll n,ll x,ll t)
{
ll ret=pow_mod(a,x,n);
ll last=ret;
for(int i=1;i<=t;i++)
{
ret=mult_mod(ret,ret,n);
if(ret==1&&last!=1&&last!=n-1) return true;
last=ret;
}
if(ret!=1) return true;
return false;
}
//是素数或者伪素数为1， 非伪素数的合数为0
bool Miller_Rabin(ll n)
{
if(n<2)return false;
if(n==2)return true;
if((n&1)==0) return false;
ll x=n-1;
ll t=0;
while((x&1)==0){x>>=1;t++;}
for(int i=0;i<S;i++)
{
ll a=rand()%(n-1)+1;
if(check(a,n,x,t))
return false;
}
return true;
}
ll factor[100];//存因子，就算是2也不超过63个
//得到的factor是不排序的
int tol;
ll gcd(ll a,ll b)
{
if(a==0)return 1;
if(a<0) return gcd(-a,b);
while(b)
{
ll t=a%b;
a=b;
b=t;
}
return a;
}
//找到一个因子
ll Pollard_rho(ll x,ll c)
{
ll i=1,k=2;
ll x0=rand()%x;
ll y=x0;
while(1)
{
i++;
x0=(mult_mod(x0,x0,x)+c)%x;
ll d=gcd(y-x0,x);
if(d!=1&&d!=x) return d;
if(y==x0) return x;
if(i==k){y=x0;k+=k;}
}
}
//因式分解n
void findfac(ll n)
{
if(Miller_Rabin(n))
{
factor[tol++]=n;
return;
}
ll p=n;
while(p>=n)p=Pollard_rho(p,rand()%(n-1)+1);
//可以分为p,n/p，其中p可以不为素数
findfac(p);
findfac(n/p);
}
int main()
{
//srand(time(NULL));
ll n;
while(scanf("%lld",&n)!=EOF)
{
tol=0;
findfac(n);
for(int i=0;i<tol;i++)printf("%lld ",factor[i]);
printf("\n");
if(Miller_Rabin(n))printf("Yes\n");
else printf("No\n");
}
return 0;
}
卡迈克尔数
反素数
短小互推的
欧拉定理
\(gcd(a,m)=1\)，则\(a^{\varphi(m)}\equiv1(\mod m)\)
扩展欧拉定理
\(a^b\equiv\begin{cases}a^{b{\mod\varphi(p)}},\quad &gcd(a,p)=1 \\ a^b,\quad &gcd(a,p)\neq1,b<\varphi(p) \\ a^{b\mod\varphi(p)+\varphi(p)},\quad &gcd(a.p)\neq1,b\geq\varphi(p)\end{cases}\pmod p\)
费马小定理
\(p\)为素数, 且\(gcd(a,p)=1\), 则\(a^{p-1}\equiv1(modp)\)
裴蜀定理
\(a,\ b\)是不全为零的整数，则存在整数\(x,\ y\), 使得\(ax+by=gcd(a,b)\)
证明及求解见扩展欧几里得。
类欧几里得算法
\(a,b,c\)均为整数，求\(f(n,a,b,c)=\sum_{i=0}^n \lfloor\frac{ai+b}{c} \rfloor\\ g(n,a,b,c)=\sum_{i=0}^n i\lfloor\frac{ai+b}{c} \rfloor\\ h(n,a,b,c)=\sum_{i=0}^n \lfloor\frac{ai+b}{c} \rfloor^2\)
类似辗转相除法，\(O(\log n)\)
可单独算\(f\)
#define M 998244353//由题目给定
#define inv2 499122177
#define inv6 166374059
ll t,n,a,b,c;
struct node {
ll f, g, h;
};
node solve(ll a, ll b, ll c, ll n) {
node ans, tmp;
if(a==0) {
ans.f=(b/c)*(n+1)%M;
ans.g=(b/c)*n%M*(n+1)%M*inv2%M;
ans.h=(b/c)*(b/c)%M*(n+1)%M;
} else if(a>=c || b>=c) {
tmp=solve(a%c, b%c, c, n);
ans.f=(tmp.f+n*(n+1)%M*inv2%M*(a/c)%M+(n+1)*(b/c)%M)%M;
ans.g=(tmp.g+(a/c)*n%M*(n+1)%M*(2*n+1)%M*inv6%M+(b/c)*n%M*(n+1)%M*inv2%M)%M;
ans.h=(tmp.h+(a/c)*(a/c)%M*n%M*(n+1)%M*(2*n+1)%M*inv6%M+
(n+1)*(b/c)%M*(b/c)%M+2*(a/c)%M*tmp.g%M+2*(b/c)%M*tmp.f%M+
2*(a/c)%M*(b/c)%M*n%M*(n+1)%M*inv2%M)%M;
} else {
long long m=(a*n+b)/c;
tmp=solve(c, c-b-1, a, m-1);
ans.f=(n*(m%M)%M-tmp.f)%M;
ans.g=(n*(n+1)%M*(m%M)%M-tmp.f-tmp.h)%M*inv2%M;
ans.h=(n*(m%M)%M*((m+1)%M)%M-2*tmp.g-2*tmp.f-ans.f)%M;
}
return ans;
}
int main() {
int T;
cin>>T;
while(T--) {
scanf("%lld%lld%lld%lld",&n,&a,&b,&c);
node ans=solve(a,b,c,n);
printf("%lld %lld %lld\n",(ans.f+M)%M,(ans.g+M)%M,(ans.h+M)%M);
}
return 0;
}
线性同余方程
形如\(ax\equiv c\pmod b\)的方程。
原方程等价于\(ax+by=c\)，有整数解的充要条件为\(\gcd(a,b)\mid c\)。
特殊：\(\gcd(a,b)=1,\ x_0,y_0\)为方程一组解，通解为\(x=x_0+bt,\ y=y_0-at\)，\(t\)为任意整数。
int Exgcd(int a, int b, int &x, int &y) {
if (!b){
x=1, y=0;
return a;
}
int d=Exgcd(b, a%b, x, y);
int t=x;
x=y;
y=t-(a/b)*y;
return d;
}
bool liEu(int a, int b, int c, int& x, int& y){
int d=ex_gcd(a, b, x, y);
if(c%d!=0) return 0;
int k=c/d;
x*=k;
y*=k;
return 1;
}
ll x, y;
bool ans=liEu(a, b, c, x, y);
若\(gcd(a,b)=1\)，则\(ax=c\pmod b\)在\([0,b-1]\)上有唯一解。
若\(gcd(a,b)=d\)，则\(ax=c\pmod b\)在\([0,b/d-1]\)上有唯一解。
\(\therefore\)我的求\(ax\equiv c\pmod b\)：可以但没必要，这样算了两次，慢一些
ll x, y, ste=(a,b);
if(c%ste) return -1;
a/=ste, b/=ste, c/=ste;
exgcd(a, c, x, y);//求c=1
x=(x*c)%b;
x=(x%b+b)%b;//求最小非负x
中国剩余定理
中国剩余定理（Chinese Remainder Theorem, CRT ）
求解形如\(\begin{cases}x&\equiv &a_1\pmod {n_1}\\ x&\equiv &a_2\pmod{n_2}\\ x&\equiv &a_3\pmod{n_3}\\ &\vdots&\\ x&\equiv &a_n\pmod{n_k}\end{cases}\)的一元线性同余方程，其中\(n_1, n_2, n_3, \dots, n_k\)两两互质，求最小非负\(x\)。
求解：
\(n= \prod _{i=1} ^kn_i\)
第\(i\)个方程：\(m_i:m_i=\frac n{n_i},\\ m_i^{-1}:m_im_i^{-1}\equiv 1\pmod {n_i},\\ c_i: c_i=m_im_i^{-1}.\)
\(c_i\)不对\(n_i\)取模。
方程组唯一解：\(ans=\sum_{i=1}^ka_ic_i\pmod n\)
证明：不会证，不想证，懒得证。
//共r行，b[maxn]存n，a[maxn]存a
//辅助 Exgcd
ll b[11], a[11];//b=n, b>a
int r;
ll ans;
int main()
{
scanf("%d", &r);
for(int i=0; i<r; i++)//先读M，再读余数
scanf("%lld%lld", &b[i], &a[i]);
ll sb=1;
for(int i=0;i<r;i++)
sb*=b[i];
for(int i=0;i<r;i++)
{
ll m=sb/b[i], x, y;
Exgcd(m, b[i], x, y);
ans=((ans+m*x*a[i])%sb+sb)%sb;
}
printf("%lld\n",ans);
}
//防爆ll
ll mul(ll a,ll b, ll mod){
bool f=0;
if(a<0) f^=1, a=-a;
if(b<0) f^=1, b=-b;
a%=mod, b%=mod;
ll ans=0;
while(b){
if(b&1) ans=(ans+a)%mod;
b>>=1;
a=(a+a)%mod;
}
if(f) ans=-ans;
return ans;
}
ll Exgcd(ll a, ll b, ll &x, ll &y) {
if (!b) {
x=1, y=0;
return a;
}
ll d=Exgcd(b, a%b, x, y);
ll t=x;
x=y;
y=t-(a/b)*y;
return d;
}
//b[maxn]->n,a[maxn]->a
ll b[11], a[11];//b=n, b>a
int r;
ll zg()
{
ll ans=0;
ll sb=1;
for(int i=0;i<r;i++)
sb*=b[i];
for(int i=0;i<r;i++)
{
ll m=sb/b[i], x, y;
Exgcd(m, b[i], x, y);
ll t1=mul(m, x, sb);
ll t2=mul(t1, a[i], sb);
ans=((ans+t2)%sb+sb)%sb;
// ans=((ans+m*x*a[i])%sb+sb)%sb;
}
return ans;
}
应用：一些计数问题给出的模不是质数，但发现该模数没有平方因子，即由两两互质的的质数相乘得到，则可分别对素因子进行计算，再用CRT求解答案。
模数不互质解决方法 ？？
二次剩余
a为p的二次剩余：\(p\nmid a\)(a不是p的倍数)，且\(a\mod p\)为平方数。
a为p的非二次剩余：\(p\nmid a\)(a不是p的倍数)，且\(a\mod p\)不为平方数。
对二次剩余求解，即解方程：\(x^2\equiv a\pmod p\)
当p为奇素数时，使用Cipolla算法求解。（若\(p\)换为非素数，使用N次剩余）
满足\(x^2\equiv n\pmod p\)的\(n\)有\(\frac {p-1}2\)个(不含0)，非二次剩余也是有\(\frac {p-1}2\)个。
显然，若\(x^2\equiv n\pmod p\)，对\(x+kp,\ k\in {N^+}\)有相同\(n\)。
证：对不同\(n\)给出一对解\(x_1,\ x_2,\ x_1<p, x_2<p\)（若\(x_i>p\)先模\(p\)即可同\(n\)），则有\(x_1^2\equiv x_2^2\pmod p\)，则\(p\mid (x_1-x_2)(x_1+x_2)\)，因为\(\mid x_1-x_2\mid<p\)且\(p\nmid x_1-x_2\)，所以\(p\mid x_1+x_2\)，满足的\(x_i\)对共\(\frac {p-1}2\)种，由于为充要条件，所以\(x_1+x_2\neq p\)是，对应的\(n\)必不同。
故只需令\(x=1,2,3\dots \frac{p-1}2\)则可得除全部\(x_1,x_2\)与\(n\)的对应。
勒让德符号
\((\frac np)=\begin{cases} 1,&p\nmid n 且n是p的二次剩余\\ -1,&p\nmid n 且n不是p的二次剩余\\ 0,&p\mid n\end{cases}\)
欧拉判别准则
\((\frac np)\equiv n^{\frac{p-1}2}\pmod p\)
证明：由费马小定理\(n^{p-1}\equiv 1\pmod p\)则有\((n^{p-1}2-1)(n^{p-1}2+1)\equiv0\pmod p\)，两者中必有一个为\(p\)的倍数，所以\(n^{\frac{p-1}2}\mod p\)必为\(+1\)或\(-1\)，然后还没看懂。
Cipolla算法
解方程\(x^2\equiv n\pmod p\)：
找一个数\(a\)满足\(a^2-n\)为非二次剩余，随机数生成\(a\)，再用欧拉判别判断是否为非二次剩余。
自建一个类似复数域，令\(i^2=a^2-n\)，就可将所有数表示为\(a+bi\)形式，\(a,b\)是模意思上的数。
方程的解为\(x\equiv n^\frac12\equiv(a+i)^{\frac{p+1}2}\pmod p\)，\(i^2\)开方后有正负两个值，带入对应两个解。
特殊的：\(n=0\)，因为\(p\)为奇素数，所以只有一个解为0。其他，因为\(a\)随机选取，随机返回一个解\(x_1\)，另一个解为\(p-x_1\)。
特殊的：\(p=2\)，只有一个解为\(n\mod p\)。
复杂度为：先测试是否为二次剩余\(O(\log(\frac{p-1}2))\)。若为二次剩余再随机数生成一个非二次剩余，期望为\(O(1.5)\)。再建立复数域，计算复数的\(\frac{p-1}2\)次方，复杂度为\(O(\log(\frac{p-1}2))\)。故总复杂度，若为非二次剩余为\(O(\log(\frac{p-1}2))\)；若为二次剩余为\(O(\log(\frac{p-1}2)+1.5+\log(\frac{p-1}2))=O(2\log(\frac{p-1}2)+1.5)\)。
#include <bits/stdc++.h>
#define ll long long
using namespace std;
int t;
ll n, p;
ll w;
struct num { //建立一个复数域
ll x, y;//实部单位为1，虚部单位为i
};
num mul(num a, num b, ll p) { //复数乘法
num ans = {0, 0};
ans.x = ((a.x * b.x % p + a.y * b.y % p * w % p) % p + p) % p;
ans.y = ((a.x * b.y % p + a.y * b.x % p) % p + p) % p;
return ans;
}
ll binpow_real(ll a, ll b, ll p) { //实数快速幂
ll ans = 1;
while (b) {
if (b & 1) ans = ans * a % p;
a = a * a % p;
b >>= 1;
}
return ans % p;
}
ll binpow_imag(num a, ll b, ll p) { //复数快速幂
num ans = {1, 0};
while (b) {
if (b & 1) ans = mul(ans, a, p);
a = mul(a, a, p);
b >>= 1;
}
return ans.x % p;//最后虚部会等于0
}
ll cipolla(ll n, ll p) {//返回x_1, x_2=p-x_1, 返回大数或小数是随机的，
//因为随机获得a和虚部i^2=a^2-n
n %= p;
if (n == 0) return 0;
if (p == 2) return n;
if (binpow_real(n, (p - 1) / 2, p) == p - 1) return -1;
ll a;
while (1) { //生成随机数再检验找到满足非二次剩余的a
a = rand() % p;
w = ((a * a % p - n) % p + p) % p;
if (binpow_real(w, (p - 1) / 2, p) == p - 1) break;
}
num x = {a, 1};
return binpow_imag(x, (p + 1) / 2, p);
}
int main(){
srand(time(0));//不加也行
int T;
cin>>T;
while(T--){
scanf("%lld%lld", &n, &p);
if(n==0){
printf("0\n");
continue;
}
if(p==2){
printf("%lld\n", n%p);
continue;
}
//求p为奇素数,且n%p!=0,无解或解有两个
ll x1=cipolla(n, p);
if(x1==-1){
printf("Hola!\n");
continue;
}
ll x2=p-x1;
if(x1>x2) swap(x1, x2);
printf("%lld %lld\n", x1, x2);
}
return 0;
}
N次剩余
求解\(x^A\equiv B\pmod m\)，\(m\)可以为任意正整数，不限于素数。大体思路：
完全分解\(m\)，可得到\(\prod p_i^{a_i}\)，由CRT可求出\(x^A\equiv B\pmod{p_i^{a_i}}\)的解。
在每个方程中选出一个解，就可对应一个方程组（原方程）的解。
各方程的解数相乘，就是方程组解的个数。
要求所有的解，就要枚举每个方程选取每个解时，计算方程组解。（类似dfs）
下面求解\(x^A\equiv B\pmod M\)，其中\(M=p^a\)：
当\(B\equiv0\pmod M\)，则令\(x=p^t\)，\(t\times A\geq a\)即可。
当\(\gcd(B,M)=1\)，求\(M\)原根\(g\)（\(p^a\)形式的\(M\)一定有原根），\(\dots\)
当\(\gcd(B,M)\neq1\)，则\(\gcd(B,M)\)必为\(p^q\)，\(0<q<a\)的形式，则令\(x=p^t\)。要求\(p^{t\times A}\equiv p^q\pmod{p^a}\)，等价\(t\times A\equiv q\pmod a\)，求\(t\)，若有解则必\(A\mid q\)。\(\dots\)
#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
int A, B, mod;
//快速幂
int pow(int x, int y, int mod = 0, int ans = 1) {
if (mod) {
for (; y; y >>= 1, x = (LL) x * x % mod)
if (y & 1) ans = (LL) ans * x % mod;
} else {
for (; y; y >>= 1, x = x * x)
if (y & 1) ans = ans * x;
}
return ans;
}
//分解m，求x^A=B mod m，转化为x^A=B mod p_i^a_i，
//p_i^a_i为m的完全分解，从而转换为中国剩余定理可解决的问题
//tot为本质不同p_i个数
//prime为p_i从小到大，expo为p_i对应a_i，pk为对应p_i^a_i
//phi(i)为p_i^a_i的欧拉函数
struct factor {
int prime[20], expo[20], pk[20], tot;
void factor_integer(int n) {
tot = 0;
for (int i = 2; i * i <= n; ++i) if (n % i == 0) {
prime[tot] = i, expo[tot] = 0, pk[tot] = 1;
do ++expo[tot], pk[tot] *= i; while ((n /= i) % i == 0);
++tot;
}
if (n > 1) prime[tot] = n, expo[tot] = 1, pk[tot++] = n;
}
int phi(int id) const {
return pk[id] / prime[id] * (prime[id] - 1);
}
} mods, _p;
//求x对pk[i](M)的逆元
//x^-1=1*x^-1=x^phiM*x^-1=x^(phiM-1) mod M
int p_inverse(int x, int id) {
assert(x % mods.prime[id] != 0);
return pow(x, mods.phi(id) - 1, mods.pk[id]);
}
//求ax+by=gcd(a,b)
void exgcd(int a, int b, int &x, int &y) {
if (!b) x = 1, y = 0;
else exgcd(b, a % b, y, x), y -= a / b * x;
}
//求xret+mody=gcd(x, mod), if(ret>=2e31 && (mod&1)) ret++
int inverse(int x, int mod) {
assert(__gcd(x, mod) == 1);
int ret, tmp;
exgcd(x, mod, ret, tmp), ret %= mod;
return ret + ((ret >> 31) & mod);
}
vector<int> sol[20];
void solve_2(int id, int k) {
int mod = 1 << k;
if (k == 0) { sol[id].emplace_back(0); return; }
else {
solve_2(id, k - 1); vector<int> t;
for (int s : sol[id]) {
if (!((pow(s, A) ^ B) & (mod - 1)))
t.emplace_back(s);
if (!((pow(s | (1 << (k - 1)), A) ^ B) & (mod - 1)))
t.emplace_back(s | (1 << (k - 1)));
}
swap(sol[id], t);
}
}
// g^x = B (mod M) => g^iL = B*g^j (mod M) : iL - j
int BSGS(int B, int g, int mod) {
unordered_map<int, int> map;
int L = ceil(sqrt(mod)), t = 1;
for (int i = 1; i <= L; ++i) {
t = (LL) t * g % mod;
map[(LL) B * t % mod] = i;
}
int now = 1;
for (int i = 1; i <= L; ++i) {
now = (LL) now * t % mod;
if (map.count(now)) return i * L - map[now];
}
assert(0);
return -1;//不会走到这，防warning
}
int find_primitive_root(int id) {
int phi = mods.phi(id); _p.factor_integer(phi);
auto check = [&] (int g) {
for (int i = 0; i < _p.tot; ++i)
if (pow(g, phi / _p.prime[i], mods.pk[id]) == 1)
return 0;
return 1;
};
for (int g = 2; g < mods.pk[id]; ++g) if (check(g)) return g;
assert(0);
return -1;//不会走到这，防warning
}
// ax = b (mod M)
void division(int id, int a, int b, int mod) {
int M = mod, g = __gcd(__gcd(a, b), mod);
a /= g, b /= g, mod /= g;
if (__gcd(a, mod) > 1) return;
int t = (LL) b * inverse(a, mod) % mod;
for (; t < M; t += mod) sol[id].emplace_back(t);
}
void solve_p(int id, int B = ::B) {
int p = mods.prime[id], e = mods.expo[id], mod = mods.pk[id];
if (B % mod == 0) {
int q = pow(p, (e + A - 1) / A);
for (int t = 0; t < mods.pk[id]; t += q)
sol[id].emplace_back(t);
} else if (B % p != 0) {
int phi = mods.phi(id);
int g = find_primitive_root(id), z = BSGS(B, g, mod);
division(id, A, z, phi);
for (int &x : sol[id]) x = pow(g, x, mod);
} else {
int q = 0; while (B % p == 0) B /= p, ++q;
int pq = pow(p, q);
if (q % A != 0) return;
mods.expo[id] -= q, mods.pk[id] /= pq;
solve_p(id, B);
mods.expo[id] += q, mods.pk[id] *= pq;
if (!sol[id].size()) return;
int s = pow(p, q - q / A);
int t = pow(p, q / A);
int u = pow(p, e - q);
vector<int> res;
for (int y : sol[id]) {
for (int i = 0; i < s; ++i)
res.emplace_back((i * u + y) * t);
}
swap(sol[id], res);
}
}
vector<int> allans;
void dfs(int dep, int ans, int mod) {
if (dep == mods.tot) {allans.emplace_back(ans); return;}
int p = mods.pk[dep], k = p_inverse(mod % p, dep);
for (int a : sol[dep]) {
int nxt = (LL) (a - ans % p + p) * k % p * mod + ans;
dfs(dep + 1, nxt, mod * p);
}
}
//x^A=B mod mod
//m可以不为素数，输入顺序为A, mod, B
//输出第一行为不同解的个数，保证所有解的个数和小于等于1e6
//第二行为解，从小到大排序
void solve() {
cin >> A >> mod >> B, mods.factor_integer(mod);
allans.clear();
for (int i = mods.tot - 1; ~i; --i) {
sol[i].clear();
mods.prime[i] == 2 ? solve_2(i, mods.expo[i]) : solve_p(i);
if (!sol[i].size()) {return cout << 0 << '\n', void(0);}
}
dfs(0, 0, 1), sort(allans.begin(), allans.end());
cout << allans.size() << '\n';
for (int i : allans) cout << i << ' '; cout << '\n';
}
int main() {
//ios::sync_with_stdio(0), cin.tie(0);
int tc; cin >> tc; while (tc--) solve();
return 0;
}
卢卡斯定理
\(m\)为素数时，用于求解大组合数取模的问题。
模数为\(p\)，定理内容：\(C_n^m\bmod p=C_{\lfloor{\frac n p}\rfloor}^{\lfloor{\frac m p}\rfloor}\times C_{n\bmod p}^{m\bmod p}\bmod p\)
其中：\(m\bmod p,\ n\bmod p\ <p\)，可直接计算，
\(C_{\lfloor{\frac n p}\rfloor}^{\lfloor{\frac m p}\rfloor}\)继续使用卢卡斯定理求解，\(m=0\)时结束，返回\(1\)。
复杂度：\(O(f(p)+g(n)\log n)\)，\(f(p)\)为预处理的复杂度，\(g(n)\)为单次求组合数的复杂度。
ll Lucas(ll n, ll m, ll p) {
if(m==0) return 1;
return (C(n%p, m%p, p)*Lucas(n/p, m/p, p))%p;
}
exLucas定理
\(m\)不为素数。
完全分解：\(n=p_1^{k_1}p_2^{k_2}\cdots p_s^{k_s},\ \{n>1\}\)
结合中国剩余定理，列出\(i\)个方程\(C_n^m\equiv a_i\pmod{p_i^{q_i}}\)并求解。
然而，\(p_i^{q_i}\)也不一定为质数，以下为求解\(C_n^m\bmod p^t\)：
BSGS
BSGS（baby-step gaint-step），大步小步算法。求解离散对数，可以在\(O\sqrt p\)内求解\(a^x\equiv b\pmod p\)。注意区分N次剩余，这里求解指数，而N次剩余求解底数，且这里是\(p\)非\(m\)，要求\(\gcd(a,p)=1\)。得到\([0, p)\)中所有满足的\(x\)。
实际上， 就算\(p\)不是质数，只要\(\gcd(a,m)=1\)就可以用BSGS求解。下面我就写奇素数\(p\)了。
\(\because a^0\equiv1\pmod p\)且\(a^{\varphi(p)}\equiv1\pmod p\Rightarrow a^{p-1}\equiv 1\pmod p\)，出现循环，则\([0,p-2]\)为一个简化剩余系，其间出现所有（p-1个不同的）\(b\)。
因此，只需在\([0, p-2]\)中枚举\(x\)，暴力即可得到解。复杂度为\(p\)。
改变枚举方式，将一轮枚举变为两轮枚举：设\([0, p-2]\)内的两个数\(x_1\geq x_2\)，令\(a^{x_1+x_2}\equiv b\pmod p\)，则\(a^{x_1}\equiv b\times a^{-x_2}\pmod p\)。则原问题等价于是否存在\((x_1, x_2)\)对，使得上式成立。先枚举\(x_2\)存map（或者hash一下），再枚举\(x_1\)看左边值是否出现过，若出现过，则为一组解，枚举次数总和仍为\(p\)。
进一步约束\((x_1,x_2)\)的相对关系，减少枚举量：创建约束\(x=x_1\lceil\sqrt p\rceil-x_2\)，当\(x_2\geq\sqrt p\)，\(x_2-=\sqrt p,\ x_1+=1\)，使得\([0,p-2]\)的\(x\)可在\([0,\lceil\sqrt p\rceil]\)的\(x_1,x_2\)中表示出来。同上，依次枚举\(x_2, x_1\)，获得答案。时间复杂度\(2\sqrt p\)，若用map就加个\(\log(\sqrt p)\)。
当\(p\)换为\(m\)非质数时，仍要满足\(\gcd(a,m)=1\)。\(x\)在\([0,l]\)为一个循环节，计算所得\(b\)对应\(m\)的简化剩余系，故每个\(b\)与\(m\)互质。在以后的循环节，其它解依次加上\(l\)。以下模板为求在第一个循环节中的解。
kuangbin板子：
//hash
#define MOD 76543
int hs[MOD], head[MOD], nex[MOD], id[MOD], top;
void insert(int x,int y)
{
int k = x%MOD;
hs[top] = x, id[top] = y, nex[top] = head[k], head[k] = top++;
}
int find(int x)
{
int k = x%MOD;
for(int i = head[k]; i != -1; i = nex[i])
if(hs[i] == x)
return id[i];
return -1;
}
//a^x=b mod n, n可以不为素数，只要(a,n)=1即可
//无解返回-1，有解返回一个最小（在[0,l-1]中的）整数解
//l为a关于模m的阶，之后是循环，要全部解，依次之后加上l即可
int BSGS(int a,int b,int n)
{
memset(head,-1,sizeof(head));
top = 1;
if(b == 1)return 0;
int m = sqrt(n*1.0), j;
ll x = 1, p = 1;
for(int i = 0; i < m; ++i, p = p*a%n)insert(p*b%n,i);
for(ll i = m; ;i += m)
{
if( (j = find(x = x*p%n)) != -1 )return i-j;
if(i > n)break;
}
return -1;
}
int main(){
int a, b, p;
while(scanf("%d%d%d", &p, &a, &b)!=EOF){
int res=BSGS(a, b, p);
if(res==-1) printf("no solution\n");
else printf("%d\n", res);
}
return 0;
}
扩展BSGS
解决BSGS中\(\gcd(a,m)\not=1\)的情况。
解方程\(a^x\equiv b\pmod m\)，其中\(d=\gcd(a,m),\ d\not=1\)。
若\(d\not\mid b\)，则除\(x=0,\ b=1\)外无其他\((x,b)\)对满足。
若\(d\mid b\)，则\(a^x\equiv b\pmod m\)，令\(k\)为使得\(d_2=\gcd(a^k,m)\)最大前提下的最小整数，即不断让左侧乘\(a\)，直到\(d_2\)不再增大，则\(a^{x-k}\times \frac{a^k}{d_2}\equiv\frac{b}{d_2}\pmod{\frac m{d_2}}\)，
问题变为\(a^{x-k}\equiv\frac b{d_2}\times(\frac{a^k}{d_2})^{-1}\pmod {\frac m{d_2}}\)，\(\gcd(a,\frac m{d_2})=1\)。
若\(d_2\not\mid b\)，显然只有\(x=0,b=1\)一对解；
对\(d_2\mid p\)，只需在\([0,k]\)枚举\(x\)，后面部分直接BSGS，注意求得\(x\)要加上\(k\)。
以下采用将方程化为\(a^{x-k}\equiv\frac b{d_2}\times(\frac{a^k}{d_2})^{-1}\pmod {\frac m{d_2}}\)，\(\gcd(a,\frac m{d_2})=1\)的形式后，直接使用BSGS。所以，指数\(k\)最多\(31\)，加号连接所以前面部分忽略不计，直接考虑BSGS复杂度为\(k+2\sqrt{\frac m{d_2}}\)，不计\(k\)，则复杂度为\(2\sqrt{\frac m{d_2}}\)。也是只求出最小循环节的解，要别的继续加\(\text{ord}_{\frac m{d_2}}a\)。
kuangbin板子：
//hash用的，计算中不管数组大小
const int MAXN= 99991;
struct LINK{
ll data;
ll j;
ll next;
}HASH_LINK[1000000];
ll ad, head[MAXN];
ll Gcd(ll a, ll b){
return b ? Gcd(b, a % b) : a;
}
ll Ext_Gcd(ll a, ll b, ll &x, ll &y){
if(!b){
x = 1; y = 0;
return a;
}
ll r = Ext_Gcd(b, a % b, x, y);
ll t = x; x = y; y = t - a / b * y;
return r;
}
ll POWER(ll a, ll b, ll c){
ll ans = 1;
while(b){
if(b & 1) ans = ans * a % c;
a = a * a % c;
b >>= 1;
}
return ans;
}
void init(){
memset(head, -1, sizeof(head));
ad = 0;
}
ll Hash(ll a){
return a % MAXN;
}
void INSERT_HASH(ll i, ll buf){
ll hs = Hash(buf), tail;
for(tail = head[hs]; ~tail; tail = HASH_LINK[tail]. next)
if(buf == HASH_LINK[tail]. data) return;
HASH_LINK[ad]. data = buf;
HASH_LINK[ad]. j = i;
HASH_LINK[ad]. next = head[hs];
head[hs] = ad ++;
}
ll exBSGS(ll a, ll b, ll c){
ll i, buf, m, temp, g, D, x, y, n = 0;
for(i = 0, buf = 1; i < 100; i ++, buf = buf * a % c)
if(buf == b) return i;
D = 1;
while((g = Gcd(a, c)) != 1){
if(b % g) return -1; // g | b 不满足，则说明无解
b /= g;
c /= g;
D = D * a / g % c;
++ n;
}
init();
m = ceil(sqrt((long double) c));
for(i = 0, buf = 1; i <= m; buf = buf * a % c, i ++) INSERT_HASH(i, buf);
for(i = 0, temp = POWER(a, m, c), buf = D; i <= m; i ++, buf = temp * buf % c){
Ext_Gcd(buf, c, x, y);
x = ((x * b) % c + c) % c;
for(ll tail = head[Hash(x)]; ~tail; tail = HASH_LINK[tail].next)
if(HASH_LINK[tail]. data == x) return HASH_LINK[tail].j + n + i * m;
}
return -1;
}
//a^x=b mod m
int main(){
ll a, b, m;
while(scanf("%lld%lld%lld", &a, &m, &b)!=EOF && (a || b || m)){
ll res=exBSGS(a, b, m);
if(res==-1) printf("No Solution\n");
else printf("%lld\n", res);
}
}
原根
\(\gcd(a,m)=1\)，使\(a^l\equiv1\pmod m\)成立的最小\(l\)，当然\(l>0\)，称为\(a\)关于模\(m\)的阶，记作\(\text{ord}_ma\)。
作用：\(\text{ord}_ma=l\)，\(\because\ a^0\equiv1\pmod m,\ a^l\equiv1\pmod m\)，模数重复了，之后等价于模数依次乘\(a\)，故模数变化相同。因此，\(a^0,a^1,a^2,\dots,a^{l-1}\)是对模\(m\)的第一个循环节，其中没有相同的余数，此后一直循环。
即：若\(g\)是\(m\)的原根，则\(g^1,g^2,\dots ,g^{\varphi(m)}\)在模\(m\)意义下两两不同，恰好为\(1\)到\(p-1\)间与\(m\)互质的数，共\(\varphi(m)\)个。
\(\text{ord}_ma=l\)，则\(\text{ord}_ma^t=\frac l{(t,l)}\)，显然找最小正整数\(x\)满足\(l\mid tx\)，故\(x=\frac l{(t,l)}\)。
完全分解\(m=p_1^{a_1}p_2^{a_2}\dots p_k^{a_k}\)，显然\(\text{ord}_ma=[\text{ord}_{p_1}^{a_1},\text{ord}_{p_2}^{a_2},\dots,\text{ord}_{p_k}^{a_k}]\)
特殊的，当\(m\)为素数\(p\)：\(\text{ord}_pa=l\)，显然\(l\mid\varphi(p)\)，则对\(p\)阶为\(l\)的所有\(a\)本质不同\(a\mod p\)有\(\varphi(l)\)个。eg：\(\text{ord}_{17}2=8\)，\(\varphi(8)=4\)，当\(a=2,8,9,15\)时，有\(l=8\)。
特殊的：
\((a, m)=1\)，当\(\text{ord}_ma=\varphi(m)\)时，称\(a\)为\(m\)的一个原根。
\(a\)为\(m\)的一个原根充要条件为\(a,a^2,\dots ,a^{\varphi(m)}\)构成模\(m\)的一个既约剩余系。因为\(\gcd(a, m)=1\)，所以不会出现\(a^x\equiv b\pmod p\)使\(\gcd(b, m)!=1\)，又因为\(a^0\equiv a^{\varphi(m)-1}\equiv1\pmod m\)，形成循环，所以，当\((a,m)=1\)且\(\text{ord}_ma=\varphi(m)\)一定形成一个简化剩余系。
（简化剩余系也称既约剩余系或缩系，是\(m\)的完全剩余系中与\(m\)互素的数构成的子集）
判断原根存在：存在原根的\(m\)一定为：\(2,4,p^a,2p^a\)的形式，\(p\)为奇素数，\(a\)为正整数。(充要条件)
求一个原根：\((g,m)=1\)，设\(p_1,p_2,\dots,p_k\)为\(\varphi(m)\)的所有不同素因数，则\(g\)
为\(m\)原根的充要条件为：对所有素因数都有\(g^{\frac {\varphi(m)}{p_i}}\not\equiv 1\pmod m\)。
只需枚举最小原根，在\(O(m^{0.25})\)级别可枚举到。
求所有原根：先求出一个，则集合\(S=\{g^s\mid 1\leq s\leq\varphi(m),\gcd(s,\varphi(m))=1\}\)，为\(m\)的全部原根。即：若\(m\)存在原根则一定有\(\varphi(\varphi(m))\)个。
自己写的模板，虽然过了题，但是很慢，要是tle就提前打素数表，分解因式时大于2个判断，并提前退出，其他不变。
#define ll long long
ll qpm(ll a, ll b, ll m) {
a%=m;
ll res=1;
while(b>0) {
if(b&1) res=res*a%m;
a=a*a%m;
b>>=1;
}
return res;
}
//来自N次剩余
struct factor {
int m;
int prime[20], expo[20], pk[20], tot;
void factor_integer(int n) {
m=n;
tot = 0;
for (int i = 2; i * i <= n; ++i) if (n % i == 0) {
prime[tot] = i, expo[tot] = 0, pk[tot] = 1;
do ++expo[tot], pk[tot] *= i; while ((n /= i) % i == 0);
++tot;
}
if (n > 1) prime[tot] = n, expo[tot] = 1, pk[tot++] = n;
}
int phi() const {//与N次剩余此处修该
int res=m;
for(int i=0; i<tot; i++)
res=res/prime[i] * (prime[i] - 1);
return res;
}
}f, fp;
bool check(int a, int m){//O(m^0.25)
for(int i=0; i<fp.tot; i++){
if(qpm(a, f.phi()/fp.prime[i], m)==(ll)1)
return 0;
}
return 1;
}
//O(m^0.25)找到最小原根，再O(m)实验全部根
int rt(int m){
f.factor_integer(m);
if(m==2 || m==4 || (f.tot==1 && f.prime[0]!=2) || (f.tot==2
&& f.pk[0]==2));
else return -1;//不是2,4,p^2,2p^a，p是奇素数的形式
fp.factor_integer(f.phi());
for(int i=1; i<m; i++){
if(__gcd(i, m)!=1) continue;
if(check(i, m)) return i;
}
return -1;//不会到这一步
}
vector<int> v;
void art(int m){
v.clear();
int mrt=rt(m);
if(mrt==-1){
v.push_back(-1);
return;
}
for(int i=1; i<=f.phi(); i++){
if(__gcd(i, f.phi())!=1) continue;
v.push_back(qpm(mrt, i, m));
}
return;
}
void p(){
if(v[0]==-1){
printf("0\n\n");
return;
}
int si=v.size();
printf("%d\n", si);
sort(v.begin(), v.end());
for(int i=0; i<si; i++){
printf("%d ", v[i]);
}
printf("\n");
return;
}
int main(){
int m, T;
cin>>T;
while(T--){
scanf("%d", &m);
art(m);
p();
}
return 0;
}
杜教筛
用于求一些数论函数（积性函数）的前缀和，即对\(f\)，求\(S(n)=\sum_{i=0}^nf(i)\)。
想要快于\(O(n)\)的快速计算，要构造一个\(S(n)\)关于\(S(\lfloor\frac ni\rfloor)\)的递推式。
结论是：\(g(1)S(n)=\sum_{i=1}^n(f\ast g)(i)-\sum_{i=2}^ng(i)S(\lfloor\frac ni\rfloor)\)。
后面部分用数论分块可快速计算，只需要找到一个使前面部分快速计算的\(g\)即可。
常见数论函数前缀和对应的\(g\)函数：
（1）\(f=\mu,\ g=u(\equiv1),\ I=\mu\ast u,\)
\(\therefore S(n)=1-\sum_{i=2}^nS(\lfloor\frac ni\rfloor)\)
（2）\(f=\varphi,\ g=u,\ e=\varphi\ast u,\)
\(\therefore S(n)=\frac 12n(n+1)-\sum_{i=2}^nS(\lfloor\frac ni\rfloor)\)
乘法逆元
注意，一般写\(a>0,b>0, \gcd(a,b)=\gcd(a,b-a)\)，当\(ax+by=k\gcd(a,b)，y\)在原来基础上取个负数，\(gcd(a,b)=gcd(a,-b)=\gcd(a,a-b)\)。因为对\((a,a-b)\)进行一次辗转相除就可以回到\((a,a-b)\)，是相等的。所以不只有\((a,b)=(a,b-a)\)还有\((a,b)=(a,a-b)\)。
线性同余方程\(ax\equiv1\pmod b\)，则称\(x\)为\(a\mod b\)的逆元，记作\(a^{-1}\)。
费马小定理：\(p\)为质数且\(gcd(a,p)=1\Rightarrow a^{-1}=a^{p-2}\)，快速幂即可
扩展欧几里得：n
用于求\(ax+by=gcd(a, b)\)的一组解。
两边同除以\(gcd(a, b)\)后，等价于\(gcd(a, b)=1\)，求\(ax+by=1\)的一组解。
\(a=q_1b+r_1,\quad r_1=a\%b\)
\(b=q_2r_1+r_2,\quad r_2=b\%r_1\)
\(r_1=q_3r_2+r_3,\quad r_3=r_1\%r_2\)
\(\dots\)
\(r_{n-1}=q_{n+1}r_n,\quad r_{n+1}=0\)，结束。
\(gcd(a,b)=r_n\)；
则有\(Q_ka-P_kb=(-1)^{k-1}r_k,\ k=1,2,3,\dots,n\)，
其中\(P_0=1,\ P_1=q_1,\ P_k=q_kP_{k-1}+P_{k-2}\)
\(Q_0=0,\ Q_1=1,\ Q_k=q_kQ_{k-1}+Q_{k-2},\quad 2<=k<=n\)
计算：
我的想法：
使用元组\((x,\ y)\)表示\(ax+by\)的值，即：\(a=(1,0),\ b=(0,1)\)。有：
\(r_1=a-q_1b=(1,-q_1)\)
\(r_2=b-q_2r_1=(0, 1)-q_2(1,-q_1)=(-q_2，q_1q_2+1)\)
\(r_3=r_1-q_3r_2=(1,-q_1)-q_3(-q_2,q_1q_2+1)\)
\(\dots\)
\(r_n=r_{n-2}-q_nr_{n-1}=(x_n,y_n)\)
\(r_{n+1}=0=r_{n-1}-q_{n+1}r_n=(x_{n+1},y_{n+1})\)
\(\because r_n=gcd(a,b)\)
\(\therefore r_n=(x_n,y_n)\)为\(ax+by=gcd(a,b)\)的一个解。
exgcd：
\(ax_1+by_1=gcd(a, b)\)
\(bx_2+(a\mod b)y_2=gcd(b, a\mod b)=gcd(a, b)\)
\(ax_1+by_1=bx_2+(a\mod b)y_2\)
\(a\mod b=a-(\lfloor\frac{a}{b}\rfloor\times b)\)
\(ax_1+by_1=bx_2-(a-(\lfloor\frac{a}{b}\rfloor)\times b)y_2\)
\(x_1=y_2,\quad y_1=x_2-\lfloor\frac{a}{b}\rfloor y_2\)
ll Exgcd(ll a, ll b, ll &x, ll &y) {
if (!b) {
x=1, y=0;
return a;
}
ll d=Exgcd(b, a%b, x, y);
ll t=x;
x=y;
y=t-(a/b)*y;
return d;
}
ll Inv(ll a, ll m){
ll x, y;
ll d=exgcd(a, m, x, y);
if(d==1) return (x%m+m)%m;
return -1;
}
线性递推：1到n
\(O(n)\)
求\(a\)在膜\(m\)下的逆元：
\(1^{-1}\equiv1\pmod m\)
\(m=qa+r,\quad qa+r\equiv0\pmod m\)
两边同乘以\(a^{-1},\ r^{-1}\)：
\(qr^{-1}+a^{-1}\equiv0\pmod m\)
\(a^{-1}\equiv-qr^{-1}\pmod m\)
\(a^{-1}\equiv-(\frac{m}{a})\times(m\bmod a)^{-1}\pmod m\)
避免负号出现，原式加上\(m\)：
\(a^{-1}\equiv(m-\frac{m}{a})\times(m\bmod a)^{-1}\pmod m\)
inv[1]=1;
for(int i=2; i<=n; i++) inv[i]=(ll)(M-M/i)*inv[M%i]%M;
线性递推：1!到n!
\(O(n+\log p)\)
可用于求任意\(n\)个数的逆元，其中\(1\leq a_i<p\)，因为倒着递推，不必计算\(1\)到\(\min(a_i)-1\)的逆元。
使用费马小定理或扩展欧几里得计算最后一个数\(a_t\)的前缀积的逆元，
乘以最后一个数后得到\(a_{t-1}\)的前缀积的逆元，以此类推。
void init(){
fact[0]=1;
for(int i=1; i<maxn; i++){
fact[i]=fact[i-1]*i%M;
}
inv[maxn-1]=qpm(fact[maxn-1], M-2);
for (int i=maxn-2; i>=0; i--){
inv[i]=inv[i+1]*(i+1)%M;
}
}
线性递推：给定随机n个数
\(1\leq a_i\leq p\)的\(n\)个数存在\(a_1\)到\(a_n\)中，类似\(1!\)到\(n!\)，先求前缀积\(s_n\)，求\(s_n\)逆元\(invs_n\)，线性递推每个\(s_i\)逆元\(invs_i\)，最后一步\(inv_i=invs_i\times s_{i-1} \mod p\)。
//a[maxn]存给定数，inv[maxn]存对应逆元，n为a的大小，从1开始存
//辅助 s[maxn], invs[maxn], qpm
void init() {
s[0]=1;
for (int i=1; i<=n; i++) s[i]=s[i-1]*a[i]%M;
invs[n]=qpm(s[n], M-2);
for (int i=n; i>=1; i--) invs[i-1]=invs[i]*a[i]%M;
for (int i=1; i<=n; i++) inv[i]=invs[i]*s[i-1]%M;
}
常见数论函数
1）u(n), e(n), I(n)
\(u(n)\equiv1,\ n\geq1\)
\(e(n)=n,\ n\geq1\)
\(I(n)=\begin{cases}1,\quad n=1\\0,\quad n>1\end{cases}\)
2）d(n) 约数个数
\(n\)的所有正除数的个数，含1和\(n\)本身：\(d(n)=\sum_{d\mid n}1\)
完全分解：\(n=p_1^{k_1}p_2^{k_2}\cdots p_s^{k_s},\ \{n>1\}\)，\(d(1)=1\)
\(d(n)=\prod_{i=1}^s(k_i+1)\)
可乘函数可使用线性筛：
\(p\)存放质数，用d表示约数个数，num表示最小质因子出现次数。考虑\(n>1\)。
\(d(p)=2, num(p)=1\)
设\(d=\frac n p\)，\(p\)为\(n\)的最小质因子：
\(p\mid d:\quad num(n)=num(d)+1,\quad d(n)=\frac{d(d)\times (num(n)+1)}{num(d)+1}\)
考虑\(d,p\)互质，因为为可乘函数：
\(p\nmid d:\quad num(n)=1,\quad d(n)=d(d)\times d(p)\)
//数组 d, num, p
//bool v int tot, n
void pre(){
d[1]=1, tot=0;
for(int i=2; i<=n; i++){
if(!v[i]) p[++tot]=i, d[i]=2, num[i]=1;
for(int j=1; j<=tot && i<=n/p[j]; j++){
v[p[j]*i]=1;
if(i%p[j]==0){
num[i*p[j]]=num[i]+1;
d[i*p[j]]=d[i]/(num[i]+1)*(num[i*p[j]]+1);
break;
}
num[i*p[j]]=1;
d[i*p[j]]=d[i]*d[p[j]];
}
}
}
3）Omega(n) 可重复素因子个数
\(n\)的全部素因子的个数\(\Omega(n)\)：
\(\Omega(1)=0\)
\(\Omega(n)=\sum_{i=1}^sa_i\)
int b[maxn];
int cnt[maxn];
void pre() {
for(int i=2; i<maxn; i++) {
if(b[i]!=0) continue;
for(int j=i; j<maxn; j+=i) {
if(b[j]==0) b[j]=i;
}
}
cnt[2]=1;
for(int i=3; i<maxn; i++)
cnt[i]=cnt[i/b[i]]+1;
return;
}
4）omega(n) 不可重复素因子个数
\(n\)的不同素因子个数\(\omega(n)\)：
\(\omega(1)=0,\quad\omega(n)=s\)
bool vis[maxn];
int cnt[maxn];
void pre() {
for(int i=2; i<maxn; i++) {
if(vis[i]!=0) continue;
cnt[i]=1;
for(int j=i+i; j<maxn; j+=i) {
cnt[j]++, vis[j]=1;
}
}
return;
}
5）sigma_lambda(n) 约数幂和
\(n\)的正除数的幂和函数\(\sigma_\lambda(n),\ \lambda\)为实数：
\(\sigma_\lambda(n)=\sum_{d\mid n}d^\lambda\)
特殊的，约数和：\(\lambda=1\)时，\(\sigma(n)=\sum_{d\mid n}d\)
完全分解：\(n=p_1^{k_1}p_2^{k_2}\cdots p_s^{k_s},\ \{n>1\}\)，\(d(1)=1\)
\(\sigma(n)=\prod_{i=1}^s(1+p_i+p_i^2+\dots+p_i^{k_i})=\prod_{i=1}^s\sum_{j=0}^{k_i}p_i^j\)
又因：\(\sum_{j=0}^{k_i}p_i^j=\frac{p_i^{k_i+1}-1}{p_i-1}\)
\(\sigma(n)=\prod_{i=1}^s(1+p_i+p_i^2+\dots+p_i^{k_i})=\prod_{i=1}^s\frac{p_i^{k_i+1}-1}{p_i-1}\)
以下计算均采用前者推导。
可乘函数可使用线性筛：
\(p\)存放质数，\(f_i\)表示\(i\)的约数和，\(g_i\)表示\(i\)的最小质因子\(p\)的\(\sum_{j=0}^{k}p^j\)。考虑\(n>1\)。
\(g(p)=p+1,\quad f(p)=i+1\)
设\(d=\frac n p\)，\(p\)为\(n\)的最小质因子：
\(p\mid d:\quad g(n)=g(d)\times p+1,\quad f(n)=\frac{f(d)\times g(n)}{g(d)}\)
考虑\(d,p\)互质，因为为可乘函数：
\(p\nmid d:\quad g(n)=p+1,\quad f(n)=f(d)\times f(p)\)
void pre(){
g[1]=f[1]=1, tot=0;
for(int i=2; i<=n; i++){
if(!v[i]) p[++tot]=i, g[i]=i+1, f[i]=i+1;
for(int j=1; j<=tot && i<=n/p[j]; j++){
v[p[j]*i]=1;
if(i%p[j]==0){
g[i*p[j]]=g[i]*p[j]+1;
f[i*p[j]]=f[i]/g[i]*g[i*p[j]];
break;
}
g[i*p[j]]=p[j]+1;
f[i*p[j]]=f[i]*f[p[j]];
}
}
}
注：求\(\sigma(n)\)是否为奇数，根据第一个公式，\(p_i=2\)时，\(k_i\)为任意数可，对其他\(p_i,\ k_i\)必须为偶数，所以\(n\)为平方数或\(n/2\)为平方数为充要条件。求\(1\)到\(n\)中，\(\sigma(i)\)为奇数的个数，\(ans=sqrt(n)+sqrt(n/2)\)即可。
6）varphi(n) 互质个数
欧拉函数\(\varphi(n)\)：小于等于\(n\)，和\(n\)互质的个数。
1与任何正整数互质，\(\varphi(1)=1\)
求解：完全分解：\(n=p_1^{k_1}p_2^{k_2}\cdots p_s^{k_s}\)
推导表达式，由容斥原理：\(\varphi(n)=n+(-\frac{n}{p_1}-\frac{n}{p_2}-\dots-\frac{n}{p_s})+(\frac{n}{p_1p_2}+\frac{n}{p_1p_3}+\dots+\frac{n}{p_{s-1}p_s})\dots\varphi(n)=\sum_{S\subseteq\{p_1,p_2\dots p_s\}}(-1)^{\mid S\mid}\frac{n}{\prod_{p_i\in S} p_i}=n\prod_{i=1}^s(1-\frac{1}{p_i})\)
\(\varphi(n)=n\times\prod_{i=1}^{s}\frac{p_i-1}{p_i}\)
int phi(int n) {
int m=(int)sqrt(n+0.5);
int ans=n;
for (int i=2; i<=m; i++)
if (n%i==0) {
ans=ans/i*(i-1);
while(n%i==0) n/=i;
}
if (n>1) ans=ans/n*(n-1);
return ans;
}
性质：
（1）为可乘函数，\(gcd(a,b)=1\)：\(\varphi(a\times b)=\varphi(a)\times\varphi(b)\)。
（2）特殊的，\(n\)为奇数：\(\varphi(2n)=\varphi(n)\)，\(\because 2\times\frac{2-1}{2}=1\)
（3）\(n=\sum_{d\mid n}\varphi(d)\)，含\(1\)和\(n\)，证明：
法一：要证明\(e(n)=n\ (n\geq1)\)为\(\varphi(n)\)的莫比乌斯变换，
只需证明\(\varphi(n)\)为\(e(n)\)的反莫比乌斯变换即可。
\(e (n)\)的反莫比乌斯变换：\(\varphi(n)=\sum_{d\mid n}\mu(d)\frac n d\)
\(\varphi(n)\)从定义式正经变换：\(\varphi(n)=\sum_{1\leq d\leq n, (d,n)=1}1=\sum_{1\leq d\leq n}\sum_{l\mid (d, n)}\mu(l)=\sum_{l\mid n}\mu(l)\sum_{1\leq d\leq n,l\mid d}1=\sum_{l\mid n}\mu(l)\frac n l\)
可见上下两式相等，故\(\varphi(n)\)为\(e(n)\)的反莫比乌斯变换。
法二：\(gcd(a,n)=d\)，按\(d\)对\(a\)进行分类，令\(n=dk\)则\(gcd(\frac a d, k)=1\)。固定的\(d\)对应\(a\)的个数为
\(\varphi(k)=\varphi(\frac n d)\)，所以\(n=\sum_{d\mid n} \varphi(\frac n d)=\sum_{d\mid n} \varphi(d)\)。
（4）特殊的\(n=p^k\)：\(\varphi(n)=p^k-p^{k-1}\)
\(=(p-1)\times p^{k-1}=n\times\frac{p-1}{p}\)，由定义可得。
可用于求\(\sum_{i=1}^{n}\sum_{j=i+1}^n gcd(i, j)\)：
显然\(ans[n]=ans[n-1]+\sum_{i=1}^{n-1}gcd(i,n)\)，仅需证明此处\(\sum_{i=1}^{n-1}gcd(i,n)=\sum_{d\mid n}\varphi(\frac n d)\times d\)：
按\(d\)对\(a\)进行分类，令\(n=dk\)则\(gcd(\frac a d, k)=1\)。固定的\(d\)对应\(a\)的个数为
\(\varphi(k)=\varphi(\frac n d)\)，和为\(\varphi(\frac n d)\times d\)。\(f(n)=\sum_{d\mid n}\varphi(\frac n d)\times d=(\varphi\ast e)(n)\)，解题关键是\(f(n)\)的表达式也可用线性筛，所以不会超时。因为此题\(d<n\)而不是\(d\leq n\)，故内层循环从\(i*2\)开始，而不是从\(i\)开始。\((\varphi\ast e)(n)\)筛法：
for(int i=1; i<=n; i++){
for(int j=i+i; j<=n; j+=i){
f[j]+=(ll)i*(ll)phi[j/i];
}
}
筛法：
void phi_table(int n){
memset(phi, 0, sizeof(phi));
phi[1]=1;
for(int i=2; i<=n; i++)
if(!phi[i])
for(int j=i; j<=n; j+=i){
if(!phi[j]) phi[j]=j;
phi[j]=phi[j]/i*(i-1);
}
}
7）mu(n)
莫比乌斯函数\(\mu(n)\)：
\(\mu(n)=\begin{cases}1,\quad &n=1\\(-1)^s,\quad &n=p_1p_2\dots p_s,\ p_1<p_2<\dots<p_s\\0,\quad &其他\end{cases}\)
容斥原理的系数啊：利如，给定集合S，给出q个查询x，求S中与x不互质的数的个数。
注意题目中出现的其他说法：存在\(k(k>1)\)使\(k^2\mod x\)则\(f(x)=0\)，否则\(f(x)=1\)。
线性筛：
//mu p v
//tot
void pre() {
mu[1]=1;
for(int i=2; i<=maxn-13; i++) {//注意这里的上界
if(!v[i]) mu[i]=-1, p[++tot]=i;
for(int j=1; j<=tot && i*p[j]<=maxn-13; j++) {//改上界
v[i*p[j]]=1;
if(i%p[j]==0) {
mu[i*p[j]]=0;
break;
}
mu[i*p[j]]=-mu[i];
}
}
}
8）Lambda(n)
Mangoldt函数\(\Lambda(n)\)：
\(\Lambda(n)=\begin{cases}\log p,\quad &n=p^k,\ k\geq1\\0\quad &其他\end{cases}\)
9）lambda(n)
Liouville函数\(\lambda(n)\)：
\(\lambda(n)=(-1)^{\Omega(n)}\)
可乘函数（积性函数）
\(f(mn)=f(m)\times f(n),\quad\gcd(m,n)=1\)
欧拉函数\(\varphi(n)\)，莫比乌斯函数\(\mu(n)\)，正因子和\(\sigma(n)\)，正因子数\(d(n)\)。
均可采用线性筛法。
可乘函数的卷积仍是可乘函数，可用线性筛：
求\(f=g\ast h\)提前处理\(g,h\)后，线性筛法如下：
for(int i=1; i<=n; i++){
for(int j=i+i; j<=n; j+=i){
f[j]+=g[i]*h[j/i];
}
}
完全（绝对）可乘函数
\(f(mn)=f(m)f(n)\)
\(u(n),\ e(n),\ I(n)\)
Dirichlet 乘积
数论分块
引理1：\(\forall a,b,c\in\mathbb{Z},\lfloor\frac{a}{bc}\rfloor=\lfloor\frac{\lfloor\frac{a}{b}\rfloor}{c}\rfloor\)
引理2：\(\forall n\in N,\vert\{\lfloor\frac{n}{d}\rfloor\mid d\in N\}\vert\leq\lfloor2\sqrt{n}\rfloor\)，\(\vert V\vert\)表示集合\(V\)的元素个数
解决类似计算\(\sum_{i=1}^n \lfloor\frac ki\rfloor\)的问题，其中\(n,\ k\)为常量：
对每个\(i(i\leq k)\)要找到相应最大\(j(i\leq j\leq k)\)，使得\(\lfloor\frac ki\rfloor=\lfloor\frac kj\rfloor\)，即可按\(j\)分块计算，显然，\(j=\lfloor\frac k{\lfloor\frac ki\rfloor}\rfloor\)。所以，想办法让求和表达式中出现\(\lfloor\frac ki\rfloor\)，且其他部分的求和可合并后\(O(1)\)直接计算即可。
例：求\(\sum_{i=1}^n (k\mod i)\)
\(=\sum_{i=1}^n k-i\lfloor\frac ki\rfloor\)
ll n, k;
int main() {
scanf("%lld%lld", &n, &k);
ll ans=n*k;
for(ll l=1, r; l<=n; l=r+1){
if(k/l!=0) r=min(k/(k/l), n);
else r=n;
ans-=(k/l)*(l+r)*(r-l+1)/2;
//(k/l)为题干中k/i，后部分为区间求和
}
printf("%lld\n", ans);
return 0;
}
标准计算(只计算到i=n，而不是无穷大)，若要乘以某个函数，提前处理其前缀和。（二维见莫比乌斯反演），若想计算到无穷大，令\(n==k\)即可。计算\(\sum_{i=1}^n\lfloor\frac ki\rfloor\)，若\(i\)固定从\(s\)开始，而不是\(1\)，把\(l\)初始值设为\(s\)即可：
ll fk(ll n, ll k){//n为i从1到n，k为k/i
ll ans=0;
for(ll l=1, r; l<=n; l=r+1){
if(k/l!=0) r=min(k/(k/l), n);
else r=n;//l大于k时
ans+=k/l*(r-l+1);//此处写带k/i的表达式
// ans+=(s[r]-s[l-1])*(k/l);
}
return ans;
}
常见卷积
设\(f(n),\ h(n)\)是两个数论函数，
\(h(n)=\sum_{d\mid n}f(d)g(\frac{n}{d})\)为\(f(n),\ g(n)\)的Dirichlet 乘积或卷积，记作：
\(h=f\ast g,\quad (f\ast g)(n)=\sum_{ab=n}f(a)g(b),\)
满足结合律和交换律。
\(f\ast g\ast h=\sum_{abc=n}f(a)g(b)h(c)\)
\(I(n)=\begin{cases}1,\quad n=1\\0,\quad n>1\end{cases}\)为卷积单位元素。
常见卷积：
\(I=\mu\ast u=\sum_{d\mid n}\mu(d)\)
\(d=u\ast u=\sum_{d\mid n}1\)
\(\sigma=e\ast u=\sum_{d\mid n}d\)
\(e=\varphi\ast u=\sum_{d\mid n}\varphi(d)\)
\(\varphi=\mu\ast e=\sum_{d\mid n}d\cdot\mu(\frac nd)\)
互为Dirichlet逆：\(I=\mu\ast u\)
若\(F=f\ast u\)，则称\(F\)为\(f\)的莫比乌斯变换，即：
\(F(n)=\sum_{d\mid n}f(d)\)，等于\(f\ast \mu^{-1}\)，即再卷积一个莫比乌斯函数可变为原函数。
还原：\(f=F\ast\mu,\ f(n)=\sum_{d\mid n}F(d)\mu(\frac n d)\)
常见F为f的莫比乌斯变换有：
\(d(n),\quad u(n)\)， 正除数个数\(=\sum_{d\mid n}1\)
\(I(n),\quad \mu(n)\)，\([n=1]=\sum_{d\mid n}\mu(d)\)
\(e(n),\quad \varphi(n)\)，\(n=\sum_{d\mid n}\varphi(d)\)
\(\sigma(n),\quad e(n)\)， 正除数和\(=\sum_{d\mid n}d\)
\(\log n,\quad\Lambda(n)\)，\(\log n=\sum_{d\mid n}\Lambda(d)\)
莫比乌斯反演
带莫比乌斯函数希望计算不是\(O(n)\)，先分解质因数，dfs质因数选与不选，虽然可以卡，但总的来说1e8可以开个根。
很难求出原函数的值，而容易求出其倍数和或约数和 ，则可通过莫比乌斯反演简化运算，求得原函数的值。
\(F=f\ast u\Rightarrow F\ast\mu=f\)
\(F(n)=\sum_{d\mid n}f(d)\)
\(f(n)=\sum_{d\mid n}F(d)\mu(\frac{n}{d})\)
因为\([n=1]=\sum_{d\mid n}\mu(d)\)，将\(n\)换为想计数的变量，即可求得该变量等于1的数量，例：
\([\gcd(i,j)=1]\iff\sum_{d\mid\gcd(i,j)}\mu(d)\)
数论分块的二维形式，求\(\sum \lfloor\frac ni\rfloor\lfloor\frac mi\rfloor\)，\(i,j\)均到无穷大：
int sol(int n, int m) {
int ans=0;
for (int l=1, r; l<=min(n, m); l=r+1) {
r=min(n/(n/l), m/(m/l));
ans+=(mu[r]-mu[l-1])*(n/l)*(m/l);
}
return ans;
}
若\(i,j\)不是到无穷大时，参考数论分块代码。
一些变换：
\(\sum_{1\leq i\leq n\\\gcd(i,n)=1} \frac 1i=\sum_{d\mid n}\mu(d)h_{\frac nd}\frac 1d\)，其中\(h\)为调和级数，可直接线性预处理。
FT
快速傅里叶变换
FFT，求多项式卷积。\(O(n\log n)\)
若求正常定义的卷积\(h=f\ast g,\ (f\ast g)(n)=\sum_{a+b=n}f(a)g(b)\)，值为\(h\)的\(n\)次项系数。
输入\(f,\ h\)从低次到高次多项式的系数表达，求\(g=f*h\)。\(f,\ h\)做\(DFT\)转换为点标记的表示方法\(F,\ H\)，对应点相乘后，得到\(G=F\cdot H\)的离散方式的各点，再对\(G\)做\(IFT\)转换为系数表达式。
蝴蝶变换：把位置\(i\)的元素直接移到\(i\)进行\(\log_2len\)次递归后的位置，为\(i\)的二进制数的翻转，eg：\(11010\Rightarrow01011\)，而取消递归。
\(IFT\)： 把多项式\(g\)的离散傅里叶变换\(G\)结果作为另一个多项式的系数，取单位根的倒数即\(\omega_n^i\ \{i=0,-1,\dots,-(n-1)\}\)作为\(x\)代入\(G\)，得到的\(G(x)\)再除以\(n\)，就是\(g\)的各项系数。
注：数组大小乘2。
struct Complex {
double x, y;
Complex(double _x = 0.0, double _y = 0.0) {
x = _x;
y = _y;
}
Complex operator-(const Complex &b) const {
return Complex(x - b.x, y - b.y);
}
Complex operator+(const Complex &b) const {
return Complex(x + b.x, y + b.y);
}
Complex operator*(const Complex &b) const {
return Complex(x * b.x - y * b.y, x * b.y + y * b.x);
}
};
// 蝴蝶变换 len必须为2的幂，i与二进制反转位置交换
void change(Complex y[], int len) {
int k;
for (int i = 1, j = len / 2; i < len - 1; i++) {
if (i < j) swap(y[i], y[j]);
k = len / 2;
while (j >= k) {
j = j - k;
k = k / 2;
}
if (j < k) j += k;
}
}
// len必须为2的幂
// on=-1逆运算是带入共轭复数再除n
void fft(Complex y[], int len, int on) {
change(y, len);
for (int h = 2; h <= len; h <<= 1) {
Complex wn(cos(2 * PI / h), sin(on * 2 * PI / h));
for (int j = 0; j < len; j += h) {
Complex w(1, 0);
for (int k = j; k < j + h / 2; k++) {
Complex u = y[k];
Complex t = w * y[k + h / 2];
y[k] = u + t;
y[k + h / 2] = u - t;
w = w * wn;
}
}
}
if (on == -1) {
for (int i = 0; i < len; i++) {
y[i].x /= len;
}
}
}
Complex x1[maxn<<2], x2[maxn<<2];
//len为最高次+1，x[i]表示x^i对应的系数
int main() {
int len1, len2, len=1;
cin>>len1>>len2;
len1++, len2++;
double t;
while(len<len1*2 || len<len2*2) len<<=1;
for(int i=0; i<len1; i++){
scanf("%lf", &t);
x1[i]=Complex(t, 0);
}
for(int i=len1; i<len; i++) x1[i]=Complex(0, 0);
for(int i=0; i<len2; i++){
scanf("%lf", &t);
x2[i]=Complex(t, 0);
}
for(int i=len2; i<len; i++) x2[i]=Complex(0, 0);
fft(x1, len, 1);
fft(x2, len, 1);
for(int i=0; i<len; i++) x1[i]=x1[i]*x2[i];
fft(x1, len, -1);
int m=len1+len2-1;
for(int i=0; i<m; i++)
printf("%d ", (int)(x1[i].x+0.5));
printf("\n");
return 0;
}
快速沃尔什变换
\(C_k=\sum_{i|j=k}A_i\times B_j\)
FWT，将多项式卷积中的\(*\)换为其他位运算符。
不取模的话，开long long，把逆元改为除2即可。
void FWT_or(int *a, int N, int opt)
{
for(int i=1;i<N;i<<=1)
for(int p=i<<1,j=0;j<N;j+=p)
for(int k=0;k<i;++k)
if(opt==1)a[i+j+k]=(a[j+k]+a[i+j+k])%MOD;
else a[i+j+k]=(a[i+j+k]+MOD-a[j+k])%MOD;
}
void FWT_and(int *a, int N, int opt)
{
for(int i=1;i<N;i<<=1)
for(int p=i<<1,j=0;j<N;j+=p)
for(int k=0;k<i;++k)
if(opt==1)a[j+k]=(a[j+k]+a[i+j+k])%MOD;
else a[j+k]=(a[j+k]+MOD-a[i+j+k])%MOD;
}
void FWT_xor(int *a, int N, int opt)
{
for(int i=1;i<N;i<<=1)
for(int p=i<<1,j=0;j<N;j+=p)
for(int k=0;k<i;++k){
int X=a[j+k],Y=a[i+j+k];
a[j+k]=(X+Y)%MOD;a[i+j+k]=(X+MOD-Y)%MOD;
if(opt==-1) a[j+k]=1ll*a[j+k]*inv2%MOD,a[i+j+k]=1ll*a[i+j+k]*inv2%MOD;
}
}
快速幂
求\(x\)逆元时，可以先将在分母处\(x%M\)，再计算\(qpm(x,M-2)\)，这在单独计算\(n!\)逆元时比较方便。
ll qpm(ll a, ll b) {
ll res=1;
while(b>0) {
if(b&1) res=res*a%M;
a=a*a%M;
b>>=1;
}
return res;
}
快速乘
求a>1e9 && b>1e9但模M时，可防爆ll：
ll mul(ll a,ll b, ll mod){
bool f=0;
if(a<0) f^=1, a=-a;
if(b<0) f^=1, b=-b;
a%=mod, b%=mod;
ll ans=0;
while(b){
if(b&1) ans=(ans+a)%mod;
b>>=1;
a=(a+a)%mod;
}
if(f) ans=-ans;
return ans;
}
常见数列
斐波那契数
\(F_0=0\)，\(F_1=1\)，\(F_n=F_{n-1}+F_{n-2}\)
前几项：0，1，1，2，3，5，8，13，21，34，55，89
性质：1.\(F_{(n-1)}F_{(n+1)}-F_n^2=(-1)^n\)
2.\(F_{(n+k)}=F_kF_{(n+1)}+F_{k-1}F_n\) \(\Rightarrow\) \(F_{2n}=F_n(F_{n+1}+F_{n-1})\)
3.\(\Rightarrow\)
4.\(gcd\)性质：\((F_m,F_n)=F(m,n)\)
5.\(F_n=F_{n-1}+F_{n-3}+F_{n-5}\dots\)也是类似斐波那契数(\(F_0=1\)，其余不变)
eg题目：给出\(n\)个格子，在一些格子中放灯（\(r\ge1\)，即至少照亮本身）使\(n\)被灯光精确覆盖的方案数：\(f_1=f_2=1,f_3=2\)；存在完全分类枚举最左边的那个灯的格子数（\(2r-1\)）即可则令一端方案数为\(f_{n-(2r-1)}\)，\(f_n=f_{n-1}+f_{n-3}+\dots\)，因为\(n\)为偶数最后两位为\(f_1+f_3,f_1=f_2\)，\(2,3\)连续则可以按\(f_n=f_{n-1}+f_{n-2}\)递推到\(f_n\)，\(n\)为奇数最后两位为\(f_0+f_2\)，但这里\(n==0\)不分配灯也是一种可行的方案数（即只有一个最左边的灯），所以\(f_0=1\neq0\)，\(f_0=f_1\)，\(1,2\)连续，则可继续递推到\(f_n\)。
卡特兰数
\(H_0=1\)，\(H_1=1\)
前几项：1，1，2，5，14，42，132
描述：
1.将\(n\)分为\(n_1,n_2\)，将\(n_1\)分为\(n_{11},n_{12}\)，...分法数，例如多边形中划分三角形划分种数，\(H_n=\begin{cases} \sum_{i=1}^{n}H_{i-1}H_{n-i},n\geq2,n\in{N_+}\\1,n=0,1 \end{cases}\)
2.\(num(a)==num(b)==n\)，当对任意\(p\)有\(\sum_{i=1}^{p}a\geq\sum_{i=1}^{p}b\)即将\(n\)个\(+1,-1\)填入长\(2n\)的数组，\(\forall s_i\ge0\)的排列数，例如在栈中存取数的方案数。按这样分类是等价1的：按第一次前缀和\(s=0\)时位置分类，只有一次遇到\(s=0\)即\(pos=2n\)，\(H_{n-1}=H_{n-1}H_0\)，遇到两次及以上\(s=0\)，前部分至少保留1位其余方案数\(H_{i-1}\)，后部分不用保留方案数\(H_{n-i}\)，所以总方案数\(\sum_{i=1}^nH_{i-1}H_{n-i}\)。
\(H_n=\frac{C_{2n}^n}{n+1}=C_{2n}^n-C_{2n}^{n-1}\)，\(a_1\)与\(a_2\)无区别，或认为\(a_i\)间顺序已定，公式等价2：按第一次前缀和\(s=-1\)时位置分类，前部分\(+1\)与后部分\(-1\)个数和为\(n-1\)，前部分\(-1\)与后部分\(+1\)个数和为\(n+1\)，所以排列数为\(C_{2n}^{n-1}\)。或用二维平面要过\(y=x+1\)，第一次经过后反转，最终到\((x-1, y+1)\)更为明显。
一般使用递推式：\(H_n=\frac{H_{n-1}\times(4n-2)}{n+1}\)，可由上述推导验证。
Harmonic(调和级数)
\(h_n=\frac 11+\frac 12+\frac 13+\dots+\frac 1n\)
存;储：
\(n<2^{31}\)，且多次查询。不能开那么大的数组，正常求法是分段打表，仅存100整数倍函数值，再递推。可以整数取逆元，也可算小数。
不太正常的做法，当\(n\)较大时，只能算小数：
欧拉常数\(\gamma=\lim_{n\to+\infty}[(\sum_{k=1}^n\frac 1k)-\ln(n)]=\int_1^\infty(\frac 1{\lfloor x\rfloor}-\frac 1x)dx\)
\(h_n\)为发散无穷级数，但\(S_n=(\sum_{k=1}^n\frac 1k)-\ln(n)\)存在极限\(\gamma\)。
\(\gamma= 0.57721566490153286060651209 \dots\)
\(\therefore \lim_{n\to+\infty}h_n \approx \gamma+\ln(n)\)
另外一般的\(h_n=\gamma+\ln(n)+\epsilon_n\)，其中\(\epsilon_n\approx\frac1{2n}\)
常见常数
e
\(e=1+\sum \frac 1{i!}= 2.718281828459\dots\)
pi
\(\arctan=x-\frac{x^3}3+\frac{x^5}5-\dots\)，带入\(x=1\)得到
莱布尼兹级数：\(\frac{\pi}4=\sum_{i=0}^\infty\frac{(-1)^k}{2k+1}\)
收敛很慢，基本上没用。
快速收敛公式：
拉马努金公式：\(\frac1\pi=\frac{2\sqrt2}{99^2}\sum_{k=0}^\infty\frac{(4k!)(1103+26390k)}{(k!)^4396^{4k}}\)
每计算一项得到8个有效数字。其变体，现有最快的
Chudnovsky公式：\(\frac1\pi=\frac1{53360\sqrt{640320}}\sum_{k=0}^\infty(-1)^k\frac{(6k!)}{(k!)^3(3k)!}\times\frac{13591409+545140134k}{640320^{3k}}\)
每计算一项得到14个有效数字。
特殊的：
BBP公式：\(\pi=\sum_{k=0}^\infty[\frac1{16^k}(\frac4{8k+1}-\frac2{8k+4}-\frac1{8k+5}-\frac1{8k+6})]\)计算一项可得到1.25个有效数字，虽然慢，但可在\(O(n)\)计算\(\pi\)的特定位(16进制)。分成四部分计算，两边同乘以\(16^{n-1}\)整数部分不要，小数部分相加，最后\(ans=res-(int)res\)，若小于零则加1。计算完前\(n-1\)项后，可适当添加一些常数加项防止没进位，加22位都过题了。
泰勒公式
\(f(x)=\sum_{i=0}^\infty\frac {f^i(x_0)}{i!}(x-x_0)^i\)
麦克劳林：令\(x_0=0,\ f(x)=\sum_{i=0}^\infty\frac {f^i(0)}{i!}x^i\)
常见：
\(e^x=\sum_{i=0}^\infty \frac 1{i!}x^i\) 1
\(\ln(1+x)=\sum_{i=0}^\infty \frac{(-1)^i}{i+1}x^{i+1}\) 2
\(\frac 1{1-x}=\sum_{i=0}^\infty x^i\) 3
\(\frac 1{e^n-1}=\sum_{i=1}^\infty e^{-in}\) 3
组合数学
排列组合
\(C_n^k+C_n^{k-1}=C_{n+1}^k\)分为第一个选与不选
\((x_1+x_2+\cdots+x_t)^n=\sum(_{n_1n_2\cdots n_t}^{n})x_1^{n_1}x_2^{n_2}\cdots x_t^{n_t}\)，其中含满足\(n_1+n_2+\cdots+n_t=n\)的所有非负整数，\((_{n_1n_2\cdots n_t}^{n})=C_{n}^{n_1}\times C_{n-n_1}^{n_2}\times C_{n-n_1-n_2}^{n_3}\times\dots\times C_{n_t}{n_t}=\frac{n!}{n_1!\times n_2!\times\dots\times n_t!}\)意思是从\(n\)中选出\(n_1,n_2,\cdots,n_t\)共\(t\)类\(n\)项，有多少种选法。
http://codeforces.com/contest/1614/problem/C因为对称\(2^{n-1}=\C_n^1+\C_n^3+\dots\)，所以计算一个01串的所有（不连续子串的异或）的和是一定的。eg：100:1+0+0+1+1+0=4，110=1+1+0+0+1+1+0=4，111=1+1+1+0+0+1=4
设a个1，b个0，则和为\(\sum_{i=1}^{i+=2,i<=a}\C_a^i*2^b=2^{a-1}*2^b=2^n-1\)，注意，全0和为0
ll qpm(ll a, ll b) {
ll res=1;
while(b>0) {
if(b&1) res=res*a%M;
a=a*a%M;
b>>=1;
}
return res;
}
ll fac[maxn];
void init(){
fac[0]=1;
for(int i=1; i<maxn; i++) fac[i]=fac[i-1]*i%M;
}
ll C(ll n, ll m){
return fac[n]*qpm(fac[m], M-2)%M*qpm(fac[n-m], M-2)%M;
}
其它数学
Cantor
排列组合的字典序，从0开始，结果经常加1。
康托展开预处理后，可在\(o(1)\)查询较小表格、字符串等映射。
int cantor(string s){
int ans=0;
for(int i=0; i<9; i++){
int k=0;
for(int j=i+1; j<9; j++)
if(s[i]>s[j]) k++;
ans+=fac[9-i-1]*k;
}
return ans;
}
string decantor(int ans){
int num=0;
string s;
set<int> st;
for(int i=0; i<9; i++) st.insert(i);
while(num<9){
auto it=st.begin();
for(int i=0; i<ans/fac[9-num-1]; i++) it++;
s+=*it+'0';
st.erase(it);
ans%=fac[9-num-1];
num++;
}
return s;
}
数位dp
例题模板
//不要连续的62和4
int dp[20][2], num[20];
ll dfs(int len, int sta, int lim) {//sta表示前一个是6，lim表示前缀是上限，可能一些题目有xy，限制与前缀无关只xy当前位互相限制，就没有sta
if(len<1) return 1;
if(!lim && dp[len][sta]!=-1) return dp[len][sta];
int ans=0, mi=lim ? num[len] : 9;
for(int i=0; i<=mi; i++) {
if(i==4 || (sta && i==2)) continue;
ans+=dfs(len-1, i==6, lim && i==mi);
}
if(!lim) dp[len][sta]=ans;
return ans;
}
int jj(int x) {
int cnt=0;
while(x) {
num[++cnt]=x%10;
x/=10;
}
return dfs(cnt, 0, 1);
}
int main() {
memset(dp, -1, sizeof(dp));
int n, m;
while(cin>>n>>m && n && m) {
cout<<jj(m)-jj(n-1)<<endl;
}
return 0;
}
https://codeforces.com/contest/1245/problem/F
给出\(0\le l\le r\le10^9\)求\(l\le a,b\le r\)满足\(a+b=a\ xor\ b\)，即\(a\&b=0\)。
其中一直错的一点：把\(x_1,x_2\)写为int，会使一些\(mi\)从\(0\)变为\(1\)，比如\(1>>32\)会等于\(1\)而不是\(0\)，注意这一点。所以要不把\(len\)从31开始，要不把\(x\)改为\(ll\)。
#include<bits/stdc++.h>
using namespace std;
#define ll long long
ll dp[43][2][2];
ll dfs(int len, int lim1, int lim2, ll x1, ll x2) {
if(len<0) return 1;
if(dp[len][lim1][lim2]!=-1) return dp[len][lim1][lim2];
dp[len][lim1][lim2]=0;
int mi1=lim1?((x1>>len)&1):1, mi2=lim2?((x2>>len)&1):1;
for(int i=0; i<=mi1; i++) {
for(int j=0; j<=mi2; j++){
if(i&j) continue;
dp[len][lim1][lim2]+=dfs(len-1, lim1 && (i==mi1),
lim2 && (j==mi2), x1, x2);
}
}
return dp[len][lim1][lim2];
}
ll cal(int x1, int x2){
if(x1<0 || x2<0) return 0;
memset(dp, -1, sizeof(dp));
return dfs(32, 1, 1, x1, x2);
}
int main() {
int T, a, b;
scanf("%d", &T);
while(T--){
scanf("%d%d", &a, &b);
printf("%lld\n", cal(b, b)-cal(a-1, b)*2+cal(a-1, a-1));
}
return 0;
}
计算基础
\(1^2+2^2+3^2+\cdots+n^2=\frac{n\times(n+1)\times(2n+1)}{6}\)
\(1^3+2^3+3^3+\cdots+n^3=(1+2+3+\cdots+n)^2=(\frac{n\times(n+1)}{2})^2\)
每个\(k\)次项都可以由\(k+1\)项多项式表示，在一些题目推导出结果为一些求和多项式后，可对每次项进行化简再合并，使\(O(n)\)到\(O(1)\)。以下附出前10次：
求区间\([i, j]\)的两两组合乘积的和：
即\([1,3]=2+3+6\)：
法一：\(dp[i][j+1]=dp[i][j]+a[j]*(sum[j-1]-sum[i-1])\)
\(dp[i+1][j]=dp[i][j]-a[i]*(sum[j]-sum[i])\)
法二：\(2w[i][j]=(sum[j]+sum[i-1])^2-(squaresum[j]-squaresum[i-1])\)
求\(n^k\)的最高\(m\)位：
令\(10^w=n^k\)，\(w\)为小数，\(x=(ll)w,\ y=w-x\)，则\(10^x=100\dots0\)，\(10^y=\frac{n^k}{10^x}\)，相当于右移\(x\)位，最后留一位在小数点前。所以，\((ll)10^y\)位\(n^k\)的最高位，\((ll)10^{y+m-1}\)为\(n^k\)的最高\(m\)位。
关键在于把\(k\)次方在\(\log10\)中变成乘\(k\)：\(w=\log10(n^k)=k\times\log10(n)\)
\(10^y\)使用\(pow(10, y)\)即可。
注：c++中\(\log(n)\)是以\(e\)为底，而不是\(10\)，\(\log10(n)\)才是以十为底。且结果可为小数，此处利用可为小数，关键在于我不用关注如何计算\(10\)的小数次方。输出变(ll)前一定要加eps，不然会减一。
求\(lcm(i,j)=n,\ i\leq j\)的\((i,j)\)对数：
完全分解：\(n=p_1^{k_1}p_2^{k_2}\cdots p_s^{k_s},\ \{n>1\}\)。
仅看\(p_i\)，\((i,j)\)组合数有\(2*(k_i+1)-1\)，最后再加上一次\(i==j\)的情况，所以\(ans=(\prod_{i=1}^s(2\times k_i+1)+1)/2\)
求\(n!\)中\(p\)的个数：
即\(n!\)中可有多少个\(p\)相乘。
完全分解：\(n!=\prod_{p\leq n}p^{\sum[\frac n{p^r}]}\)。证明见数论笔记本。
对于素因子\(p\)的指数：\(h=[\frac n p]+[\frac n{p^2}]+\dots=\sum_{r=1}^{\infty}[\frac n {p^r}]\)
ll check(ll m, ll p){
ll ans=0;
while(m){
ans+=m/p;
m/=p;
}
return ans;
}
求\(ax+by=c\)且\(x\geq0,b\geq0\)的\((x,y)\)对数，\(\gcd(a,b)=1\)：
对固定\(a,b\)，变动\(c\)，外层循环枚举\(x\)，内层循环枚举\(y\)，给\(vis[ax+by]\)赋值1，即求出\(a*b\)(即\(lcm(a,b)\)内的)次数。对于\(c>a*b\)，在查询\(c\)时，\(cnt=c/(ab)+vis [c\%ab]\)。
求\(ax+(a+1)y+(a+2)z==b\)，\(x,y,z\)为非负整数的\((x,y,z)\)对数：
只选\(a,a+1\)等于\(b\)时\(ans++\)，且\(y>=2\)时可将\(2(x+1)=x+(x+2)\)，所以贡献为\(y/2+1\),同理计算\(a+1,a+2\)等于\(b\)时，贡献相同。
其中，将\(2(a+1)=a+(a+2)\)的变化，不会产生重复的\((x,y,z)\)，为什么呢，因为相当于按把\(x,z\)全部尽可能移向\(y\)分类，此时的\((x,z)\)不同，则最终不同。
其中有重复计算\((k,y,k)\)，若存在\((a+1)\mid b\)可计算后减去重复计算的\(\frac{\frac b{a+1}}2+1\)。或者直接把\(x,y\)的初始值设为\(0,1\)或\(1,0\)。
//s[i]存可变动a使ax+(a+1)y+(a+2)z=i的(x,y,z)对数。
//若要计算特定a，去掉最外层循环即可
void pre2() {
for(ll x=1; x<=n/3-1; x++) {//枚举a
ll d=n;
for(ll i=0; x*i<=d; i++) {
for(ll j=0; x*i+(x+1)*j<=d; j++) {
s[x*i+(x+1)*j+3*x+3]+=j/2+1;
}
}
for(ll i=1; (x+2)*i<=d; i++) {
for(ll j=0; (x+2)*i+(x+1)*j<=d; j++) {
s[(x+2)*i+(x+1)*j+3*x+3]+=j/2+1;
}
}
}
}
求有限\(\pm2^i\)表示数：
有\(2^i\)若干个，能否用\(\sum\pm2^i\)表示\(n\)。直接指数从大到小贪心。
for(int i=25; i>=0; i--){//题目i最大25
while(cnt[i]--){
if(w>=0) w-=a[i];
else w+=a[i];
}
}
if(w) cout<<"No\n";
else cout<<"Yes\n";
四边形不等式
若二元函数满足当\(a\le b\le c\le d\),有\(w(a,d)+w(b,c)\ge w(a,c)+w(b,d)\)，则称二元函数满足四变形不等式。
判定： 若\(a<b\)，\(w(a,b+1)+w(a+1,b)≥w(a,b)+w(a+1,b+1)\)则二元函数\(w\)满足四变形不等式。
应用： 方程\(f_i=min_{0\le j<i}\{f_j+w(j,i)\}\)中\(w\)满足四边形不等式，则决策点单调递增。
多项式
拉格朗日插值
\(f(x)=\sum_{i=1}^ny_i\prod_{i\neq j}\frac{x-x_j}{x_i-x_j},\ O(n^2)\)
特殊的，当\(x_i\)为等差数列时，可以预处理分母部分为阶乘逆元，复杂度降为\(O(nlog(n))\)
高精度
平平给的模板
typedef long long LL;
struct BigInteger {
static const int BASE = 100000000;
static const int WIDTH = 8;
vector<LL> s;
BigInteger(long long num = 0) {
*this = num;
}
BigInteger operator = (long long num) {
s.clear();
do {
s.push_back(num % BASE);
num /= BASE;
} while(num > 0);
return *this;
}
BigInteger operator = (const string & str) {
s.clear();
int x, len = (str.size() - 1) / WIDTH + 1;
for(int i = 0 ; i < len ; i++ ) {
int ed = str.size() - i * WIDTH;
int sta = max(0,ed - WIDTH);
sscanf(str.substr(sta,ed - sta).c_str(),"%d",&x);
s.push_back(x);
}
strip();
return *this;
}
bool operator < (const BigInteger & b) const {
if(s.size() != b.s.size()) return s.size() < b.s.size();
for(int i = s.size() - 1 ; i >= 0 ; i-- )
if(s[i] != b.s[i]) return s[i] < b.s[i];
return false;
}
bool operator > (const BigInteger & b) const {
return b < *this;
}
bool operator <= (const BigInteger & b) const {
return !(b < *this);
}
bool operator == (const BigInteger & b) const {
return !(b < *this) && !(*this < b);
}
bool operator != (const BigInteger & b) const {
return !(*this == b);
}
BigInteger operator * (const BigInteger & b) const {
///printf("*1\n");
BigInteger c;
c.s.resize(s.size() + b.s.size() + 1,0);
for(int i = 0 ; i < (int)s.size() ; i++ )
for(int j = 0 ; j < (int)b.s.size() ; j++ )
c.s[i + j] += s[i] * b.s[j];
for(int i = 0 ; i < (int)s.size() + (int)b.s.size() ; i++ ) if(c.s[i] > BASE) {
c.s[i + 1] += c.s[i] / BASE;
c.s[i] %= BASE;
}
c.strip();
///printf("*2\n");
return c;
}
int mult(const BigInteger & d,const BigInteger & b) const {
int L = 0,R = BASE - 1;
while(L < R) {
int m = (L + R) / 2;
if((2 * m) != L + R) ++m;
BigInteger p = b * m;
if(p == d) return m;
else if(p < d) L = m;
else R = m - 1;
}
return L;
}
BigInteger operator / (const BigInteger & b) {
///printf("\1\n");
if(b == 0) return -1;
BigInteger c,d;
c.s.resize(s.size());
for(int i = s.size() - 1 ; i >= 0 ; i-- ) {
d.s.insert(d.s.begin(),s[i]);
d.strip();
int cnt = mult(d,b);
d = d - b * cnt;
c.s[i] = cnt;
}
///printf("\2\n");
c.strip();
return c;
}
BigInteger operator - (const BigInteger & b) {
BigInteger c;
c.s.resize(s.size(),0);
for(int i = 0 ; i < (int)c.s.size() ; i++ )
c.s[i] = s[i] - (i < (int)b.s.size() ? b.s[i] : 0);
int i = c.s.size() - 1;
for( ; i >= 0 ; i-- ) if(c.s[i]) break;
for(i-- ; i >= 0 ; i-- ) if(c.s[i] < 0) {
int j = i + 1,k = 0;
while(j < (int)c.s.size() && !c.s[j]) {
j++;
k++;
}
c.s[j]--;
c.s[i] += BASE;
for( ; k > 0 ; k-- ) c.s[i + k] = BASE - 1;
}
c.strip();
return c;
}
void strip() {//去掉首位0
while(!s.empty() && !s.back()) s.pop_back();
if(s.empty()) s.push_back(0);
return ;
}
BigInteger operator + (const BigInteger & b) const {
BigInteger c;
c.s.clear();
for(int i = 0 ,g = 0 ; ; i++ ) {
if(g == 0 && i >= (int)s.size() && i >= (int)b.s.size()) break;
int x = g;
if(i < (int)s.size()) x += s[i];
if(i < (int)b.s.size()) x += b.s[i];
c.s.push_back(x % BASE);
g = x / BASE;
}
return c;
}
};
istream& operator >> (istream &in,BigInteger & x) {
string s;
if(!(in >> s)) return in;
x = s;
return in;
}
ostream& operator << (ostream &out ,const BigInteger & x) {
out << x.s.back();
for(int i = x.s.size() - 2 ; i >= 0 ; i-- ) {
char buf[20];
sprintf(buf,"%08lld",x.s[i]);
for(int j = 0 ; j < (int)strlen(buf) ; j++ ) out << buf[j];
}
return out;
}
const int N = 300;
BigInteger f[N],p[N][N],k;
int n,m,m_comple,ans[N];
int vis[N],sub[N],comple[N];
BigInteger frac(int x);
void dfs(BigInteger accum,int sub_i,int taken);
BigInteger P(int x,int y);
int main() {
BigInteger a = BigInteger((long long)1e18);
BigInteger b = BigInteger((long long)1e18);
a=a*b;
cout<<a;
return 0;
}
BigInteger frac(int x) {
if(x <= 1) return 1;
if(f[x] != 0) return f[x];
f[x] = frac(x - 1) * x;
/// printf("f[%d] = ",x);
///cout << f[x] << endl;
return f[x];
}
void dfs(BigInteger accum,int sub_i,int taken) {
if(taken == n) {
printf("%d",ans[0]);
for(int i = 1 ; i < n ; i++ ) printf(" %d",ans[i]);
printf("\n");
return ;
}
BigInteger pp;
int flag = 1,i = 0;
while(true) {
if(vis[comple[i]]) {
++i;
continue;
}
if(i >= m_comple) {
ans[taken] = sub[sub_i];
dfs(accum,sub_i + 1,taken + 1);
return ;
}
if(sub_i < m && flag && sub[sub_i] < comple[i]) {
///printf("11111\n");
flag = 0;
pp = P(n - taken - 1,m - sub_i - 1);
///cout << k << " 1 " << accum << " 1 " << pp << endl;
///printf("11111\n");
if(k <= accum + pp) {
ans[taken] = sub[sub_i];
dfs(accum,sub_i + 1,taken + 1);
return ;
} else accum = accum + pp;
} else {
///printf("222222\n");
pp = P(n - taken - 1,m - sub_i);
///cout << k << " 2 " << accum << " 2 " << pp << endl;
if(k <= accum + pp) {
ans[taken] = comple[i];
vis[comple[i]] = 1;
dfs(accum,sub_i,taken + 1);
return ;
} else {
accum = accum + pp;
++i;
}
}
}
return ;
}
BigInteger P(int x,int y) {
if(x <= 0) x = 1;
if(y <= 0) y = 1;
if(p[x][y] != 0) return p[x][y];
///printf("x = %d;y = %d\n",x,y);
///cout << frac(x) << "--------" << frac(y);
return p[x][y] = frac(x) / frac(y);
}
概率
期望逆推的时候，为什么期望逆推，我是这样想的：
对投骰子的事件，从\(A\)可以从\(B_i(1\le i\le6),B_i\le A\)，走来，其中\(B_i\)是原因，\(A\)是结果。
先看几个基本公式：
条件概率公式：\(P(A|B_i)=\frac{P(AB_i)}{P(B_i)}\)
全概率公式：\(P(A)=\sum P(B_i)P(A|B_i)=\sum P(AB_i)\)
全期望公式：\(E(A)=E(E(A|B))=\sum P(B_i)E(A|B_i)\)
贝叶斯定理：\(P(B_x|A)=\frac{P(B_x)P(A|B_x)}{\sum P(B_i)P(A|B_i)}\)
不论到\(P(A)\)还是\(E(A)\)，这里\(P(B_i)\)都是不均匀的，要是硬核地用全期望公式计算，需要存\(P(i)\)和\(E(i)\)。这里的\(E(i)\)并不表示当走到了\(i\)期望走了多少步（不是真正数学上的数学期望\(E(A)\)），因为求\(E(A)\)的前提是\(A\)发生，而代码中的\(E(A)\)还有\(1-P(A)\)时候根本没有\(A\)，这里的\(EP\)相除就得到数学期望。这里的\(E\)都是对\(P\)加权了的，这样存再加上逆推会很好解决回到开始点的问题（1有0.5回到1，0.5来到2）。这样的\(P(B_i)\)不可计算，因为对\(A\)，每个\(P(AB_i)\)不仅跟\(P(A|B_i)\)有关（题目明确描述转移），还跟\(P(B_i)\)有关（需要自己计算），怎么才可以避开自己计算\(P(B_i)\)呢?
现在想要快速计算\(E(A)\)的方法，先替换上面\(E(A)\)的含义，让他变成真正的数学期望，令\(E(A)=\frac{E(A)}{P(A)}\)。这样可以规范化\(\sum P(B_i)=1\)，同时把\(B_i\)替换为\(A\)到达的点，现在\(A\)为\(B_i\)的原因。这样，当\(A\)存在时，每一个\(P(AB_i)\)中的\(P(B_i)\)是等分了。比如每个数0.5概率回到1，0.5概率+1，当前在2，0.5到3，0.5到1，那\(E(2)=0.5E(1)+0.5E(3)+1\)，即使1出现次数比3多很多次，但结果根本不关注这个，因为那些1是别的点要去的，现在我就是2，\(P(2)=1\)就是全部，我只要考虑自己的下一步怎么走，\(1,3,4,5,\dots\)到1都是别人的生活啦。
逆推让问题变成了先知道结果\(B_i\)，然后求原因\(A\)，而全期望公式中\(B_i\)可以是就是\(A\)的结果。如果困惑为什么\(B_i\)是\(A\)的结果为什么还可以有\(E(A|B_i)\)的写法，那就看看贝叶斯公式，他求的不就是先验概率，知道的是原因\(B_i\)发生后结果\(A\)的情况，想求结果\(A\)发生的前提下，原因\(B_i\)发生的概率。所有公式中，\(A,B\)只是两个字母，数学公式才不知道哪个字母是什么事情，自然也不分原因和结果，之前一直想不明白是我一直默认两个字母的因果关系了，而这里只需要把全概率公式中\(B_i\)当作结果就好了。之前的期望dp我陷入了把\(B_i\)当原因的定式，以至于我连最基础的事件递推都困惑了，直到看到了贝叶斯公式才警醒了我，先看到的也可以是结果，就像体会到\(\max(\max())\)的递归，他们的逻辑都是在顺序上很严密的。现在才发现以前学概率论课的时候完全没体会到原因和结果之间的奥妙，只是单纯的推导了他就以为学懂了他，还没找到他的意义就把他看低，原来研究他的意义不只是已知原因想要知道结果，已知结果向原因推到在数学上原来都是适用的，甚至在期望dp的应用中体现出更优秀的效率，不过还好现在再开始了解他，还来得及快乐。
绝对值
快速求带变量的纯加法，eg:\(|1-x|+|3-7x|+|-5+2x|+\dots\)，一次计算遍历一遍，当\(x\)变化时如何快速计算？
对这种形式显然当\(b>0\)：\(|a+bx|=\begin{cases}a+bx,\quad &x\geq-\frac a b\\-a-bx,\quad &x<-\frac a b\end{cases}\)
可以预处理时让b永远为正数（\(=|-a-bx|\)），对\(-\frac a b\)进行排序即可。
浙公网安备 33010602011771号