三道线段树

2020-07-21

FDFZ JXD 作业

more >>

戳我看全文~ >>

那年那树那些事

2020-07-05

最小斯坦纳树

more >>

戳我看全文~ >>

多项式

2020-06-05

1 模板

快速傅里叶变换（FFT）

#include<bits/stdc++.h>
#define double long double 
#define pi acos(-1)
#define round(a) ((a)>0?((int)(a+0.5)):((int)(a-0.5)))
using namespace std;
typedef complex<double> cp;
const int N=(1<<20)+5;
int n,m,l,r[N]; cp a[N],b[N];

void fft(cp *a,int type) {
	for(int i=0;i<n;i++)
		if(i<r[i]) swap(a[i],a[r[i]]);
	for(int i=1;i<n;i*=2){
		cp d(cos(pi/i),type*sin(pi/i));
		for(int j=0;j<n;j+=i*2){
			cp b(1,0);
			for(int k=0;k<i;k++,b*=d){
				cp x=a[j+k], y=b*a[i+j+k];
				a[j+k]=x+y, a[i+j+k]=x-y;
			}
		}
	}
} 

int main() {
	scanf("%d",&n); m=n;
	for(int i=0;i<=n;i++) cin>>a[i];
	for(int i=0;i<=m;i++) cin>>b[i];
	m+=n;
	for(n=1;n<=m;n*=2) l++;
	for(int i=0;i<n;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
	fft(a,1), fft(b,1);
	for(int i=0;i<=n;i++) a[i]*=b[i];
	fft(a,-1);
	for(int i=0;i<=m;i++) printf("%d ",round(a[i].real()/n));
	return 0;
}

快速数论变换（NTT）

#include<bits/stdc++.h>
using namespace std;
const long long N=5e6+9,mod=998244353,g=3,ig=332748118;
typedef long long ll;
ll n,m,a[N],b[N],l,r[N];
ll ksm(ll a,ll b){
	if(b==0) return 1;
	else if(b==1) return a;
	else return ksm(a*a%mod,b/2)*(b%2?a:1)%mod;
}
void ntt(ll *f,ll type){
	for(ll i=0;i<n;i++)
		if(i<r[i]) swap(f[i],f[r[i]]);
	for(ll i=1;i<n;i*=2){
		ll wn=ksm((type>0)?g:ig,(mod-1)/(i*2));
		for(ll j=0;j<n;j+=i*2){
			for(ll k=0,w=1;k<i;k++,w=w*wn%mod){
				ll x=f[j+k],y=f[i+j+k]*w%mod;
				f[j+k]=(x+y)%mod,f[i+j+k]=(x-y+mod)%mod;
			}
		}
	}
	if(type<0){
		ll in=ksm(n,mod-2);
		for(ll i=0;i<n;i++) f[i]=f[i]*in%mod;
	}
}
int main(){
	scanf("%lld%lld",&n,&m);
	for(ll i=0;i<=n;i++) scanf("%lld",&a[i]);
	for(ll j=0;j<=m;j++) scanf("%lld",&b[j]);
	m+=n; for(n=1;n<=m;n*=2) l++;
	for (ll i=0;i<n;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
	ntt(a,1), ntt(b,1);
	for(ll i=0;i<n;i++) a[i]=a[i]*b[i]%mod;
	ntt(a,-1);
	ll in=ksm(n,mod-2);
	for(ll i=0;i<=m;i++) printf("%lld ",a[i]);
	return 0;
}

2. 二项式定理

$\tbinom{n}{k}$ 组合意义：从 $n$ 个不同元素的集合中选取 $k$ 个元素作为子集的方案数。

$\dbinom {n}{k}=\frac{n(n-1)...(n-k+1)}{k!}$

六边形性质（帕斯卡三角形中）

$\dbinom {r}{k}\dbinom {r+1}{k+2}\dbinom {r+2}{k+1}=\dbinom {r}{k+1}\dbinom {r+1}{k}\dbinom {r+2}{k+2}$

（1）对称恒等式

$\dbinom {n}{k}=\dbinom {n}{n-k}, n\in N,k\in Z$

（2）吸收恒等式

$\dbinom {r}{k}=\frac{r}{k}\dbinom {r-1}{k-1}$

推论：

$(r-k)\dbinom {r}{k}=r\dbinom {r-1}{k}$

（3）加法恒等式

$\dbinom {r}{k}=\dbinom {r-1}{k}+\dbinom {r-1}{k-1}$

二项式定理

$(x+y)^r=\sum_{k} x^{k}y^{r-k}$

范德蒙德卷积

$\sum_{k}{\dbinom{r}{m+k}\dbinom{s}{n-k}}=\dbinom{r+s}{m+n}\\ k-m代替k，n-m代替n后：\sum_{k}\dbinom{r}{k}\dbinom{s}{n-k}=\dbinom{r+s}{n}$

组合意义：$r$ 个 $0$ 和 $s$ 个 $1$ 中选取 $n$ 个数的方案数

3. 生成函数

普通型生成函数

定义

数列 $a_0,a_1,…$ 的生成函数，定义为以 $a_i$ 为系数的形式幂级数，即

$F(x)=\sum_{i=0}^{\infty} a_ix^i$

组合数

二项式定理：组合数 $C_{n}^{i}$ 的生成函数为 $F(x)=(1+x)^n$

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

例题1.1

现有 $n$ 个数，每个数为 $a_i$，可以取无限次，请问有多少方案可以让所有数的和为 $m$。

答案为

$F(x)=(1+x^{a_1}+x^{2a_1}+...)(1+x^{a_2}+x^{2a_2}+...)...(1+x^{a_n}+x^{2a_n}+...)\\ =\prod_{i=0}^{n}\sum_{k=0}^{\infty} x^{ka_i}$

的展开式中 $x^m$ 的系数。

其中，$x^m$ 相当于每一个因式中，选择一个数的指数相加为 $m$。所以得到的 $x^m$ 的项数即解，

同样，我们可以从 $dp$ 去想。设 $f(i,j)$ 为前 $i$ 个数中和为 $j$ 的情况。考虑决策，为总和做出 $+1$，$+a_i$，$+2a_i$，…的贡献，而这些就组成了这个因式的所有指数。

例题1.2

现有 $n$ 个数，每个数为 $a_i$，可以取 $b_i$ 次，请问有多少种方案可以让所有数的和为 $m$。

答案为

$F(x)=\prod_{i=0}^{n}\sum_{k=0}^{b_i} x^{ka_i}$

的展开式中 $x^m$ 的系数。

此题和例题1.1中唯一的不同是每个数都有了自己的最多取几个，所以决策变成了为总和做出 $+1$，$+a_i$，…，$+a_ib_i$ 的贡献，于是因式便不再是无穷级数。

例题1.3

KK P463 定制背包

根据题目列出生成函数

$F(x)=(1+x)(1+x^5+x^{10}+...)(1+x+x^2+x^3)(1+x^2+x^4+...)(1+x+x^2+x^3)(1+x^4+x^8+...)$

对于此类多项式，有等比数列求和&等比级数求和公式

$\sum_{i=0}^{n}x^{ik}=\frac{1-x^{(n+1)^k}}{1-x^k}\\ \sum_{i=0}^{\infty}x^{ik}=\frac{1}{1-x^{k}}$

于是我们可以将上述多项式化简为

$F(x)=(1+x)\frac{1}{1-x^5}\frac{1-x^5}{1-x}\frac{1}{1-x^2}\frac{1-x^4}{1-x}\frac{1}{1-x^4}\\ = \frac{1}{(1-x)^3}\\ = (\sum_{i=0}^{\infty}{x_i})^3$

看这个式子的组合意义，即将n拆成3个自然数的和。答案为 $C_{2}^{n+2}$ 。

例题2

（摘自 x义x的博客）求十进制下，$5$ 的数字个数为偶数的 $n$ 位数（不能有前导0）

先看 $dp$，$f_n$ 为 $5$ 的个数为偶数的 $n$ 位数个数，$g_n$ 为 $5$ 的个数为奇数的 $n$ 位数个数。

$f_n=9f_{n-1}+g_{n-1}\\ g_n=9g_{n-1}+f_{n-1}$

设 $f,g$ 的普通型生成函数为 $F(x),G(x)$，我们用错位相减法。怎么减呢？我们看 $F(x)=\sum_{i=0}^{\infty} f_ix^i$ ，那么为了减得我们已经知道的数（$f1,g1$ ），我们需要用 $F(x)-xF(x)$ 。我们看这题的转移方程，可以得到应该这样列式子

$F(x)-9xF(x)-xG(x)=8x\\ G(x)-9xG(x)-xF(x)=1x$

把 $F(x),G(x)$ 当做未知数去解方程。解得

$F(x)=\frac{8-71x}{(1-8x)(1-10x)}$

我们看到 $\frac{1}{1-8x},\frac{1}{1-10x}$ 都是可以用等比级数求和公式的。所以我们把上述式子做裂项得

$F(x)=\frac{1}{2}(\frac{7}{1-8x}+\frac{9}{1-10x})\\ =\frac{1}{2}\sum_{i=0}^{\infty} 7\times(8x)^i+9\times (10x)^i$

上述式子 $x^n$ 的系数显然是 $\frac{1}{2}(7\times 8^n+9\times 10^n)$。

例题3

UVA 12298 Super Poker II

有一副扑克，对于每一个合数都有四种花色的牌。

某些牌已经丢失，求出在剩下的牌中选出花色互不相同的四张且点数之和为 $n$ 的方案数。

列出生成函数

$F(x)=\prod_{i=1}^{4} \sum_{k=2}^{\infty}[k\not\in prime,k\not\in lost] k$

所以我们需要列出4个多项式，然后求出 $F(x)$，输出第 $a$ 项到第 $b$ 项的系数即可。

洛谷稍具体一些的题解：https://www.luogu.com.cn/blog/forever-captain/solution-uva12298

#include<bits/stdc++.h>
#define ll __int128
using namespace std;
const ll N=5e4+9,mod=4179340454199820289,g=3,ig=1393113484733273430;
ll n,m,l,r[N<<5],a[N<<5],b[N<<5],c[N<<5],d[N<<5];

ll ksm(ll a,ll b);
void ntt(ll *f,ll type);
bool vst[N];
void prime(ll n);

signed main() { 
    prime(50000);
    while(1) {
        long long ai,bi,ci;scanf("%lld%lld%lld",&ai,&bi,&ci);
        n=m=bi;
        memset(a,0,sizeof(b)), memset(b,0,sizeof(b)),
        memset(c,0,sizeof(c)), memset(d,0,sizeof(d));
        if(!ai&&!bi&&!ci) break;
        for(register ll i=0;i<=bi;i++) a[i]=b[i]=c[i]=d[i]=vst[i];
        for(register ll i=1;i<=ci;i++) {
            int g; char ch; cin>>g>>ch;
            if(ch=='S') a[g]=0;
            else if(ch=='H') b[g]=0;
            else if(ch=='C') c[g]=0;
            else if(ch=='D') d[g]=0;
        }
        //从这里到下面都是板子
        m<<=3; l=0;
        for(n=1;n<=m;n*=2) l++;
        for(ll i=0;i<=n;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
        ntt(a,1), ntt(b,1), ntt(c,1), ntt(d,1);
        for(ll i=0;i<=n;i++) a[i]=a[i]*b[i]%mod*c[i]%mod*d[i]%mod;
        ntt(a,-1);
        for(int i=ai;i<=bi;i++) printf("%lld\n",(long long)a[i]);
        puts("");
    } 
    return 0;
}

例题4

K阶差分和前缀和

给定一个长为 $n$ 的序列 a，求出其 $k$ 阶差分或前缀和。结果的每一项都需要对 $1004535809$ 取模。

$1\le n\le 10^5,\quad 0\le a_I\le 10^9,\quad 1\le k\le 10^{2333}$

对于生成函数 $F(x)$，$F(x)$ 的差分即 $(1-x)F(x)$ ，$F(x)$ 的前缀和即 $\frac{1}{1-x}F(x)$。

所以 $F(x)$ 的 $k$ 阶差分，前缀和分别为 $(1-x)^kF(x)$ 和 $(\frac{1}{1-x})^kF(x)$。

经过一系列我也不知道的玄学推导，得出

$K$ 阶差分为

$D(x)=A(x)F(x)\\ A(x)= \sum_{i=0}^{\infty}(-1)^i\dbinom {k}{i}x^i$

$K$ 阶前缀和为

$S(x)=A(x)F(x)\\ A(x)=\sum_{i=0}^{\infty}\dbinom{k+i-1}{i}x^i$

由于这是 NTT 模数，所以卷一卷就可以了。

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=1000009,mod=1004535809,ig=334845270,g=3;
int n,m,k,t,f[N],a[N],r[N]; string kk;
int ksm(int a,int b){
	if(b==0) return 1;
	else if(b==1) return a;
	else return ksm(a*a%mod,b/2)*(b%2?a:1)%mod;
}
void ntt(int *q,int type){
	for(int i=0;i<n;i++)
		if(i<r[i]) swap(q[i],q[r[i]]);
	for(int i=1;i<n;i*=2){
		int wn=ksm((type>0)?g:ig,(mod-1)/(i*2));
		for(int j=0;j<n;j+=i*2)
			for(int k=0,w=1;k<i;k++,w=w*wn%mod){
				int x=q[j+k],y=q[i+j+k]*w%mod;
				q[j+k]=(x+y)%mod,q[i+j+k]=(x-y+mod)%mod;
			}
	}
	if(type<0){
		int in=ksm(n,mod-2);
		for(int i=0;i<n;i++) q[i]=q[i]*in%mod;
	}
}
int orz(string s){
	int d=0;
    for(int i=0;i<s.size();i++)  d=(d*10+(s[i]-'0'))%mod;
    return d;
}
signed main(){
	cin>>n>>kk>>t; m=n;
	for(int i=1;i<=n;i++) scanf("%d",&f[i]);
	k=orz(kk);
	a[0]=1;
	for(int i=1;i<=n;i++)
		if(t==0) a[i]=a[i-1]*(k+i-1)%mod*ksm(i,mod-2)%mod;
		else a[i]=-1*a[i-1]*(k-i+1+mod)%mod*ksm(i,mod-2)%mod;
	int l=0; m+=n;
	for(n=1;n<=m;n*=2) l++;
	for(int i=0;i<n;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
	ntt(f,1), ntt(a,1);
	for(int i=0;i<n;i++) f[i]=f[i]*a[i]%mod;
	ntt(f,-1);
	for(int i=1;i<=m/2;i++) printf("%lld ",f[i]);
	return 0;
}

戳我看全文~ >>

莫比乌斯反演

2020-05-14

莫比乌斯反演

积性函数

定义

$\gcd (x,y)=1$ 且 $f(xy)=f(x)f(y)$，则 $f(n)$ 为积性函数。

性质

若 $f(x),g(x)$ 为积性函数，那么

$h(x)=f(x^p)\\h(x)=f^p(x)\\h(x)=f(x)g(x)\\h(x)=\sum_{d|x}f(d)g(\frac{x}{d})$

也为积性函数

例子

单位函数 $e(x)$
恒等函数 $id(x)$
常数函数 $1(x)$
除数函数 $\sigma_k(n)=\sum_{d|n}d^k\sigma_0(n)$
欧拉函数 $\varphi (n)=\sum_{i=1}^{n}[\gcd (i,n)=1]$
莫比乌斯函数 $\mu (n)$

所有积性函数都可以通过线性筛去筛出来。

莫比乌斯函数

定义式

$\mu(x)=\begin{cases}1& (x=1)\\0& (\exist d:d^2\mid n)\\(-1)^k & (x=p_1p_2p_3...p_k)\end{cases}$

计算方法

$O(n)$ 批量计算

线性筛筛质数，筛到 $k=x\times p$ 时

$\mu(k)=\begin{cases}0& (p\mid x,即p^2\mid k)\\± u_x& (p\nmid x)\end{cases}$

void mobius(){
	mu[1]=1;
	for(int i=2;i<=n;i++){
		if(!vst[i]) pr[++cnt]=i,mu[i]=-1;
		for(int j=1;j<=cnt&&i*pr[j]<=n;j++){
			vst[i*pr[j]]=1;
			if(i%pr[j]==0){mu[i*pr[j]]=0;break;}
            else mu[i*pr[j]]=-mu[i];
		}
	}
}

$O(\sqrt n)$ 单个计算

对 $1$ 到 $\sqrt{n}$ 做线性筛，由于大于 $\sqrt{n}$ 的质因数最多只有一个，所以最后判断一下即可。

性质

$\sum_{d|n} {\mu(d)}=\begin{cases}1& (n=1)\\0 (n>1)\end{cases}$

莫比乌斯反演

对于函数 $f(x)$，设 $g(x)=\sum_{x\mid d }f(d)$ ，d则有

$f(x)=\sum_{x\mid d}\mu (\frac{d}{x})g(d)$

也有，

$[n\mid m]\sum_{d\mid \frac{n}{m}}\mu(d)=[n=m]$

题目（收集自各大博客）

T1 求 $f(1)=\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=1]$

解法一

看到后面那个式子可以想到莫比乌斯反演中提到的第二个公式，$n=\gcd(i,j),m=1$，于是进行变形

$ys=\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{d\mid \gcd(i,j)}{\mu(d)}$

看到第三个求和的条件 $d\mid \gcd(i,j)$，可以转换为 $d\mid i,d\mid j$

$ys=\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{d\mid i,d\mid j}\mu(d)$

我们看到最难受的最后一个求和，想把它移到前面去，需要确定范围 $d\in [1,\min(n,m)]$

然后对于后面两个求和的范围，我们看到，既然 $d$ 已定，$d$ 一定是 $i,j$ 的约数，那么 $i,j$ 一定是 $d$ 的倍数

$ys=\sum_{d=1}^{\min(n,m)}\sum_{d\mid i}^{n}\sum_{d\mid j}^{m}\mu (d)$

把 $\mu$ 移到前面去

$ys=\sum_{d=1}^{\min(n,m)}\mu(d)\sum_{d\mid i}^{n}\sum_{d\mid j}^{m}1$

然后化简最后的两个求和

$ys=\sum_{d=1}^{\min(n,m)}\mu(d)\lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor$

可以 $O(n)$ 解决。

多组的话也可以线性筛去预处理，然后后面整除分块+前缀和解决。

解法二

设 $f(x)=\sum_{i=1}^{n}\sum_{j=1}^{m}[\gcd(i,j)=x],\quad g(x)=\sum_{x\mid d}f(d)$

$f(x)=\sum_{x\mid d}g(d)\mu(\frac{d}{x})\\ f(1)=\sum_{d=1}^{n}g(d)\mu(d)$

然后计算 $g(d)$

$g(x)=\sum_{i=1}^{n}\sum_{j=1}^{m} [x\mid \gcd(i,j)]\\ =\sum_{i=1}^{n/x}\sum_{j=1}^{m/x}1$

带入+整除分块即可

T2 【[POI2007]ZAP-Queries】求$f(x)=\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=x]$

这是一个比较基础的转变。如果 $i,j$ 都除以 $x$，那么 $\gcd$ 也除以 $x$。

$ys=\sum_{i=1}^{\lfloor\frac{n}{x}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{x}\rfloor} [\gcd(i,j)=1]\\ =\sum_{d=1}^{\min(\lfloor\frac{n}{x}\rfloor,\lfloor\frac{n}{x}\rfloor)}\mu(d)\lfloor\frac{n}{xd}\rfloor\lfloor\frac{m}{xd}\rfloor$

void xxs(int n){
	mu[1]=1;
	for(int i=2;i<=n;i++){
		if(!vst[i]) pr[++cnt]=i,mu[i]=-1;
		for(int j=1;j<=cnt&&i*pr[j]<=n;j++){
			vst[i*pr[j]]=1;
			if(i%pr[j]==0){mu[i*pr[j]]=0;break;}
            mu[i*pr[j]]=-mu[i];
		}
	}
	for(int i=1;i<=n;i++) s[i]=s[i-1]+mu[i];
}
int main(){
	int T; cin>>T; xxs(500000);
	while(T--){
		int n,m,d; scanf("%d%d%d",&n,&m,&d);
		n/=d,m/=d; long long ans=0;
		for(int l=1,r;l<=min(n,m);l=r+1){
			r=min(n/(n/l),m/(m/l));
			ans+=1ll*(n/l)*(m/l)*(s[r]-s[l-1]);
		}
		printf("%lld\n",ans);
	}
	return 0;
}

T3 【YY的GCD】求 $\sum_{i=1}^{n}\sum_{j=1}^{m}[\gcd(i,j)\in prime ]$

用莫比乌斯函数性质可得（或者定理提到的第二个东西）

$ys=\sum_{k=1}^{n}[k\in prime]\sum_{i=1}^{[\frac{n}{k}]}\sum_{j=1}^{[\frac{m}{k}]}[1=gcd(i,j)]\\ =\sum_{k=1}^{n}[k\in prime]\sum_{i=1}^{[\frac{n}{k}]}\sum_{j=1}^{[\frac{m}{k}]}\sum_{d\mid gcd(i,j)}{\mu(d)}\\ =\sum_{k=1}^{n}[k\in prime]\sum_{d=1}^{\min(\frac{n}{k},\frac{m}{k})}\mu(d)\sum_{i=1}^{[\frac{n}{kd}]}\sum_{j=1}^{[\frac{m}{kd}]} 1\\ =\sum_{k=1}^{n}[k\in prime]\sum_{d=1}^{\min(\frac{n}{k},\frac{m}{k})}\mu(d)\lfloor\frac{n}{kd}\rfloor\lfloor\frac{m}{kd}\rfloor$

然后会 TLE

显然罪魁祸首是这个 $k$，我们希望这个关于 $k$ 的东西可以放到后面然后预计算。

我们枚举的是 $k$ 和 $d$，我们考虑令 $t=kd$，然后移到前面去于是有

$ys=\sum_{t=1}^{\min(n,m)}\lfloor\frac{n}{t}\rfloor\lfloor\frac{m}{t}\rfloor\sum_{k\mid t,k\in prime}\mu(\frac{t}{k})$

后面的显然可以预计算。

#include<bits/stdc++.h>
const int N=2e7+9;
int pr[N],cnt,mu[N],f[N],s[N]; bool vst[N];
using namespace std;
void xxs(int n){
	mu[1]=1;
	for(int i=2;i<=n;i++){
		if(!vst[i]) pr[++cnt]=i,mu[i]=-1;
		for(int j=1;j<=cnt&&i*pr[j]<=n;j++){
			vst[i*pr[j]]=1;
			if(i%pr[j]==0){mu[i*pr[j]]=0;break;}
            mu[i*pr[j]]=-mu[i];
		}
	}
	for(int i=1;i<=cnt;i++)
		for(int j=1;pr[i]*j<=n;j++) f[j*pr[i]]+=mu[j];
	for(int i=1;i<=n;i++) s[i]=s[i-1]+f[i];
}
int main(){
	int T; cin>>T;
	xxs(10000000);
	while(T--){
		int n,m; scanf("%d%d",&n,&m);
		long long ans=0;
		for(int l=1,r;l<=min(n,m);l=r+1){
			r=min(n/(n/l),m/(m/l));
			ans+=1ll*(n/l)*(m/l)*(s[r]-s[l-1]);
		}
		printf("%lld\n",ans);
	}
	return 0;
}