数学小结（持续更新）

数学小结

前言：这篇文章专门整理了各种有关数学的东西（包括原理和例题），主要是怕自己忘了。

整理的不齐全，而且写的都是比较基础的东西。

扩展欧几里得（exgcd）

问方程 $ax+by=c$ 的 $x$ 的最小非负整数解。

根据贝祖定理，若 $\gcd (a,b)|c$，则方程有解，否则无解。

那么我们先可以求出 $ax+by=\gcd (a,b)$ 的一组解，然后 $x$、$y$ 同时乘上 $\frac{c}{\gcd (a,b)}$，就可以得到 $ax+by=c$ 的解了。至于怎么求 $x$ 的最小非负整数解，我们先要将原方程约分，即 $a$、$b$、$c$ 都除以 $\gcd (a,b)$，然后再根据不定方程，让 $x$ 模上 $b$ 就好了。

现在的问题就是如何求出 $ax+by=\gcd (a,b)$ 的一组解。

考虑用类似求 $\gcd$ 的辗转相除法求解。

假设我们已经求出了 $a^{\prime }{x}^{\prime }+b^{\prime }{y}^{\prime }=\gcd (a^{\prime },b^{\prime })$ 的一组解，其中 $a^{\prime }=b$，$b^{\prime }=a\mathrm{m}\mathrm{o}\mathrm{d}b=a-\lfloor \frac{a}{b}\rfloor \times b$，现在考虑如何求 $x$、$y$ 的一组解。

显然，根据辗转相除法，$\gcd (a,b)=\gcd (b,a\mathrm{m}\mathrm{o}\mathrm{d}b)=\gcd (a^{\prime },b^{\prime })$。

那么有：

$$\begin{array}{rl}{a}^{\prime }{x}^{\prime }+b^{\prime }{y}^{\prime }& =\gcd (a^{\prime },b^{\prime })\\ b{x}^{\prime }+\left(a-\lfloor \frac{a}{b}\rfloor \times b\right)y^{\prime }& =\gcd (a,b)\\ a{y}^{\prime }+b\left(x^{\prime }-\lfloor \frac{a}{b}\rfloor \times y^{\prime }\right)& =\gcd (a,b)\end{array}$$

可以推得：

$$\{\begin{array}{l}x=y^{\prime }\\ y=x^{\prime }-\lfloor \frac{a}{b}\rfloor \times y^{\prime }\end{array}$$

这样就可以求出 $x$、$y$ 的一组解了。

考虑边界情况：当 $b=0$ 时，有 $ax=\gcd (a,b)=a$，得 $x=1$，$y=0$。

然后根据一开始说的方法就能求出 $x$ 的最小非负整数解了。

代码如下：

void exgcd(ll a,ll b,ll &x,ll &y)
{
	if(!b)
	{
		x=1,y=0;
		return;
	}
	exgcd(b,a%b,x,y);
	ll tmp=x;
	x=y;
	y=tmp-a/b*y;
}

求逆元

若 $a\times x\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$，则称 $x$ 为 $a$ 在模 $p$ 意义下的逆元，记作 $a^{-1}$，通常来算 $\frac{1}{a}\mathrm{m}\mathrm{o}\mathrm{d}p$ 的值。（保证 $\gcd (a,p)=1$）

1. 若 $p\in prime$，根据费马小定理，$a^{p-1}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$，即 $a\times a^{p-2}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$，那么答案即为 $a^{p-2}\mathrm{m}\mathrm{o}\mathrm{d}p$。

2. 若 $p�\in prime$，那么 $a\times x\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$ 可以看做是 $ax+py=1$ 的一组解。

代码就不放了。

中国剩余定理（crt）

问以下方程的最小非负整数解：

$$\{\begin{array}{l}x\equiv a_1(\mathrm{m}\mathrm{o}\mathrm{d}{p}_1)\\ x\equiv a_2(\mathrm{m}\mathrm{o}\mathrm{d}{p}_2)\\ \cdots \\ x\equiv a_k(\mathrm{m}\mathrm{o}\mathrm{d}{p}_k)\end{array}$$

其中 $p_1,p_2,\cdots ,p_k$ 为两两互质的数。

而这个定理貌似就是在构造一个特殊解。

设 $M=\prod _{i=1}^k{p}_i$，$m_i=\frac{M}{p_i}$，$t_i$ 为 $m_i$ 在模 $p_i$ 意义下的逆元，则有 $m_i\times t_i\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}{p}_i)$。

那么这个方程的一个特殊解就是：$x_0=\sum _{i=1}^k{a}_i{m}_i{t}_i$。

证明：

对于任意的 $i$，我们考虑 $x_0\equiv a_i(\mathrm{m}\mathrm{o}\mathrm{d}{p}_i)$ 是否成立。

将 $x_0$ 拆解成两部分之和：$x_0=\sum _{j=1,j�=i}^k{a}_j{m}_j{t}_j+a_i{m}_i{t}_i$。

显然，对于任意的 $j�=i$，$a_j{m}_j{t}_j\mathrm{m}\mathrm{o}\mathrm{d}{p}_i$ 肯定为 $0$，因为 $m_j$ 的定义就是除了 $p_j$ 之外的 $p_i$ 之积，得 $m_j\mathrm{m}\mathrm{o}\mathrm{d}{p}_i$ 为 $0$。

那么就可以得到 $\sum _{j=1,j�=i}^k{a}_j{m}_j{t}_j\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{p}_i)$。

然后因为 $m_i{t}_i\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}{p}_i)$，得 $a_i{m}_i{t}_i\equiv a_i(\mathrm{m}\mathrm{o}\mathrm{d}{p}_i)$。

所以 $x_0=\sum _{j=1,j�=i}^k{a}_j{m}_j{t}_j+a_i{m}_i{t}_i\equiv a_i(\mathrm{m}\mathrm{o}\mathrm{d}{p}_i)$。

那么 $x_0$ 就是 $x$ 的一个特殊解了。

证毕。

然后又因为 $x$ 一定是 $aM+b$ 的形式，所以 $x_{\min }=x_0\mathrm{m}\mathrm{o}\mathrm{d}M$。

代码就不放了。

扩展中国剩余定理（excrt）

问以下方程的最小非负整数解：

$$\{\begin{array}{l}x\equiv a_1(\mathrm{m}\mathrm{o}\mathrm{d}{p}_1)\\ x\equiv a_2(\mathrm{m}\mathrm{o}\mathrm{d}{p}_2)\\ \cdots \\ x\equiv a_k(\mathrm{m}\mathrm{o}\mathrm{d}{p}_k)\end{array}$$

发现这类问题和普通中国剩余定理的区别就是少了这个条件：“其中 $p_1,p_2,\cdots ,p_k$ 为两两互质的数”。

这就导致了 $m_i$ 在模 $p_i$ 意义下的逆元可能根本不存在，所以不能那么做。

假设前 $i$ 个方程组的特殊解已经求出来为 $x_0$。

设 $M_i=\prod _{j=1}^i{p}_j$。

那么前 $i$ 个方程组的通解为 $x_i=x_0+M_i\times t$（$t\in Z$）。

我们现在就是要找到一个 $t$，使得 $x_0+M_i\times t\equiv a_i(\mathrm{m}\mathrm{o}\mathrm{d}{p}_i)$，即 $M_i\times t\equiv a_i-x_0(\mathrm{m}\mathrm{o}\mathrm{d}{p}_i)$。

这个用 exgcd 来做就可以了。

然后这么一直下去就能得到最后的解了。

代码如下：

//这里的m数组就是p数组

ll exgcd(ll a,ll b,ll &x,ll &y)
{
	if(!b)
	{
		x=1,y=0;
		return a;
	}
	ll gcd=exgcd(b,a%b,x,y);
	ll xx=x;
	x=y;
	y=xx-a/b*y;
	return gcd;
}

ll excrt()
{
	ll M=m[1],ans=aa[1];
	for(int i=2;i<=n;i++)
	{
		ll a=M,b=m[i],c=((aa[i]-ans)%b+b)%b,x,y;
		ll gcd=exgcd(a,b,x,y);
		if(c%gcd)
			return -1;
		ll mod=b/gcd;
		x=mul(x,c/gcd,mod);
		ans+=x*M;
		M*=mod;
		ans=(ans%M+M)%M;
	}
	return ans;
}

---

以上都应该是初一时应该学会的，结果没学好，又复习了一遍……

下面的就是最近才学的。

---

卢卡斯定理（Lucas）

原理

问 $\left(\genfrac{}{}{0ex}{}{n}{m}\right)\mathrm{m}\mathrm{o}\mathrm{d}p$ 的值，其中 $p\in prime$。

设 $n=ap+b$，$m=lp+r$。

那么 $\left(\genfrac{}{}{0ex}{}{n}{m}\right)$ 就相当于 $(x+1{)}^n$ 的展开式中 $x^m$ 的系数。

然后先证明一个东西：$(x+1{)}^p\equiv x^p+1(\mathrm{m}\mathrm{o}\mathrm{d}p)$。

证明比较简单，把 $(x+1{)}^p$ 展开，得 $(x+1{)}^p\equiv \sum _{i=0}^p(\genfrac{}{}{0ex}{}{p}{i})x^i(\mathrm{m}\mathrm{o}\mathrm{d}p)$，其中对于任意的 $0<i<p$，肯定有 $\left(\genfrac{}{}{0ex}{}{p}{i}\right)=\frac{p\times (p-1)\times \cdots \times (p-i+1)}{1\times 2\times \cdots \times i}$ 是 $p$ 的倍数，因为 $p$ 为质数，$1\sim i$ 均与 $p$ 互质。

证毕。

然后推式子：

$$\begin{array}{rl}& (x+1{)}^n\\ \equiv & (x+1{)}^{ap+b}\\ \equiv & (x+1{)}^{ap}\times (x+1{)}^b\\ \equiv & ((x+1{)}^p{)}^a\times (x+1{)}^b\\ \equiv & (x^p+1{)}^a\times (x+1{)}^b\\ \equiv & \left(\sum _{i=1}^a\left(\genfrac{}{}{0ex}{}{a}{i}\right)x^{ip}\right)\times \left(\sum _{j=1}^b\left(\genfrac{}{}{0ex}{}{b}{j}\right)x^j\right)(\mathrm{m}\mathrm{o}\mathrm{d}p)\end{array}$$

然后考虑如何找到这个式子中的第 $m$ 次项。

设 $m=lp+r$。

由于 $1\le j\le b\le p$，所以 $x^m=x^{lp+r}$ 只能是由第一个括号里面的 $x^{lp}$ 和第二个括号里面的 $x^r$ 拼凑而来，故 $x^m$ 的系数为 $\left(\genfrac{}{}{0ex}{}{a}{l}\right)\times \left(\genfrac{}{}{0ex}{}{b}{r}\right)$。

又因为 $x^m$ 本来的系数就是 $\left(\genfrac{}{}{0ex}{}{n}{m}\right)$，所以得：

$$\left(\genfrac{}{}{0ex}{}{n}{m}\right)=\left(\genfrac{}{}{0ex}{}{a}{l}\right)\times \left(\genfrac{}{}{0ex}{}{b}{r}\right)=\left(\genfrac{}{}{0ex}{}{\lfloor \frac{n}{p}\rfloor }{\lfloor \frac{m}{p}\rfloor }\right)\times \left(\genfrac{}{}{0ex}{}{n\mathrm{m}\mathrm{o}\mathrm{d}p}{m\mathrm{m}\mathrm{o}\mathrm{d}p}\right)$$

然后就可以递归求解了。

代码如下：

ll Lucas(ll n,ll m)
{
	if(!m) return 1;
	return Lucas(n/p,m/p)*C(n%p,m%p)%p;
}

至于怎么求 $\left(\genfrac{}{}{0ex}{}{n}{m}\right)(\mathrm{m}\mathrm{o}\mathrm{d}p)$（其中 $0\le m\le n<p$），这里介绍以下几种方法：

1. 通过杨辉三角 $O(p^2)$ 预处理 C 数组，询问 $O(1)$。

2. 用公式 $\left(\genfrac{}{}{0ex}{}{n}{m}\right)=\frac{n!}{(n-m)!m!}$ 来 $O(p)$ 求。

   代码如下：

    ll poww(ll a,ll b)
    {
    	ll ans=1;
    	while(b)
    	{
    		if(b&1) ans=(ans*a)%p;
    		a=(a*a)%p;
    		b>>=1;
    	}
    	return ans;
    }
   
    ll C(ll n,ll m)
    {
    	if(m>n) return 0;
    	m=min(m,n-m);
    	ll a=1,b=1;
    	for(int i=1;i<=m;i++)
    	{
    		a=(a*(n-i+1))%p;
    		b=(b*i)%p;
    	}
    	return a*poww(b,p-2)%p;
    }
   

3. $O(p)$ 预处理阶乘数组 fac和阶乘的逆元数组 inv，询问 $O(1)$。

例题

[ZJOI2010]排列计数

看到 $p_i>p_{\lfloor \frac{i}{2}\rfloor }$，联想到二叉树结构。

那么这就是一个结构固定的二叉树，需要满足子节点权值大于父节点权值。

类似一个小根堆。

设 $f_i$ 表示以 $i$ 为根的子树的方案数。显然这棵树里面有 $siz{e}_i$ 个权值。

那么 $i$ 节点的权值是这些权值里面最小的一个。

然后剩下的 $siz{e}_i-1$ 个权值分给左儿子 $siz{e}_{lc}$ 个，右儿子 $siz{e}_{rc}$ 个。

就有状态转移方程：$f_i=\left(\genfrac{}{}{0ex}{}{siz{e}_i-1}{siz{e}_{lc}}\right)\times f_{lc}\times f_{rc}$。

代码如下：

#include<bits/stdc++.h>

#define lc (u<<1)
#define rc (u<<1|1)
#define N 2000010
#define ll long long

using namespace std;

ll n,p,f[N],fac[N];
int size[N];

ll poww(ll a,ll b)
{
	ll ans=1;
	while(b)
	{
		if(b&1) ans=ans*a%p;
		a=a*a%p;
		b>>=1;
	}
	return ans;
}

ll C(ll n,ll m)
{
	if(n<m) return 0;
	if(!m) return 1;
	return fac[n]*poww(fac[n-m],p-2)%p*poww(fac[m],p-2)%p;
}

ll Lucas(ll n,ll m)
{
	if(n<m) return 0;
	if(!m) return 1;
	return C(n%p,m%p)*Lucas(n/p,m/p)%p;
}

void dfs(int u)
{
	if(u>n) return;
	dfs(lc),dfs(rc);
	size[u]=1+size[lc]+size[rc];
	f[u]=Lucas(size[u]-1,size[lc]);
	if(lc<=n) f[u]=(f[u]*f[lc])%p;
	if(rc<=n) f[u]=(f[u]*f[rc])%p;
}

int main()
{
	scanf("%lld%lld",&n,&p);
	fac[0]=1;
	for(int i=1;i<=n;i++)
		fac[i]=fac[i-1]*i%p;
	dfs(1);
	printf("%lld
",f[1]);
	return 0;
}

[SHOI2015]超能粒子炮·改

奇怪的推式子题。

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-vZFiEzQk-1604621224493)(http://192.168.102.138/JudgeOnline/upload/attachment/image/20180615/20180615104746_54736.png)]

代码如下：

#include<bits/stdc++.h>

#define N 2500
#define mod 2333
#define ll long long

using namespace std;

int t;
ll n,k,C[N][N],f[N][N];

void init()
{
	C[0][0]=1;
	for(int i=1;i<mod;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;
	}
	for(int i=0;i<mod;i++)
	{
		f[i][0]=1;
		for(int j=1;j<mod;j++)
			f[i][j]=(f[i][j-1]+C[i][j])%mod;
	}
}

ll Lucas(ll n,ll m)
{
	if(!m) return 1;
	if(n<m) return 0;
	return C[n%mod][m%mod]*Lucas(n/mod,m/mod)%mod;
}

ll F(ll n,ll k)
{
	if(k<0) return 0;
	if(!n||!k) return 1;
	if(n<mod&&k<mod) return f[n][k];
	return (f[n%mod][mod-1]*F(n/mod,k/mod-1)%mod+Lucas(n/mod,k/mod)*f[n%mod][k%mod]%mod)%mod;
}

int main()
{
	init();
	scanf("%d",&t);
	while(t--)
	{
		scanf("%lld%lld",&n,&k);
		printf("%lld
",F(n,k));
	}
	return 0;
}

[BJWC2018]上学路线

题解

Baby Step Giant Step（BSGS，北上广深）

原理

问 $a^x\equiv b(\mathrm{m}\mathrm{o}\mathrm{d}p)$ 的最小非负整数解，其中 $p\in prime$。

设 $t=\sqrt{p}$，$x=it-j$。（因为 $p\in prime$，所以 $0\le x\le p-1$）

$$\begin{array}{rl}{a}^x& \equiv b(\mathrm{m}\mathrm{o}\mathrm{d}p)\\ a^{it-j}& \equiv b(\mathrm{m}\mathrm{o}\mathrm{d}p)\\ a^{it}& \equiv b\times a^j(\mathrm{m}\mathrm{o}\mathrm{d}p)\end{array}$$

然后我们先可以 $O(\sqrt{p})$ 把 $b\times a_j$ 的所有值存进 hash 表里面。

然后我们再 $O(\sqrt{p})$ 枚举 $i$，在 hash 表里面找是否有 $a^{it}\equiv b\times a^j(\mathrm{m}\mathrm{o}\mathrm{d}p)$。

代码如下：

#include<bits/stdc++.h>

#define ll long long

using namespace std;

ll p,b,n,t,num,ans=-1;

map<ll,int>mp;

ll poww(ll a,ll b)
{
	ll ans=1;
	while(b)
	{
		if(b&1) ans=ans*a%p;
		a=a*a%p;
		b>>=1;
	}
	return ans;
}

int main()
{
	scanf("%lld%lld%lld",&p,&b,&n);
	num=t=sqrt(p);
	if(num*num<p) num++;
	for(int i=0;i<t;i++)
		mp[n*poww(b,i)%p]=i;
	for(int i=1;i<=num;i++)
	{
		ll tmp=poww(b,t*i);
		if(mp[tmp])
		{
			ans=t*i-mp[tmp];
			break;
		}
	}
	if(ans!=-1) printf("%lld
",ans);
	else puts("no solution");
	return 0;
}

例题

[SDOI2013]随机数生成器

用 吴potter 讲的方法求出通项公式，接下来就是一个裸的 BSGS 了。

代码如下：

#include<bits/stdc++.h>

#define int long long
#define ll long long

using namespace std;

int T,len,block,ans;
ll p,a,b,x1,t;

map<ll,int>mp;

ll poww(ll a,ll b)
{
	ll ans=1;
	while(b)
	{
		if(b&1) ans=ans*a%p;
		a=a*a%p;
		b>>=1;
	}
	return ans;
}

signed main()
{
	scanf("%lld",&T);
	while(T--)
	{
		scanf("%lld%lld%lld%lld%lld",&p,&a,&b,&x1,&t);
		if(x1==t)
		{
			puts("1");
			continue;
		}
		if(!a)
		{
			if(b==t) puts("2");
			else puts("-1");
			continue;
		}
		if(a==1)
		{
			if(!b)
			{
				if(x1==t) puts("1");
				else puts("-1");
				continue;
			}
			t=(t-x1+p)%p;
			t=t*poww(b,p-2)%p;
			printf("%lld
",t%p+1);
			continue;
		}
		t=(t+b*poww(a-1,p-2)%p)%p;
		ll tmp=(x1+b*poww(a-1,p-2)%p)%p;
		t=t*poww(tmp,p-2)%p;
		ans=-1;
		mp.clear();
		len=block=sqrt(p);
		while(len*block<p) block++;
		for(int i=1;i<=len;i++)
			mp[t*poww(a,i)%p]=i;
		for(int i=1;i<=block;i++)
		{
			ll tmp=poww(a,i*len);
			if(mp[tmp])
			{
				ans=i*len-mp[tmp];
				break;
			}
		}
		if(ans!=-1) printf("%lld
",ans%p+1);
		else puts("-1");
	}
	return 0;
}
/*
1
398098159 1 0 391559262 79657487
*/

矩阵树定理

原理

对于一个大小为 $n^2$ 的矩阵 $A$，定义其行列式 $\det (A)$ 为：（$\det (A)$ 也作 $|A|$）

$$\det (A)=\sum _P(-1{)}^{\mu (P)}\prod _{i=1}^{n}A(i,p_i)$$

其中 $P$ 是 $1\sim n$ 的一个排列，$\mu (P)$ 是排列 $P$ 的逆序对数。

也就是说，我们在矩阵的每一行选一个数（它们所在的列互不相同），然后乘起来。对于每一种选数方案，赋予它一个与逆序对有关的系数，再加起来。

矩阵树定理有以下的性质：

1. 交换矩阵的两行，行列式变号。

   也就等同于交换一个 $1\sim n$ 排列的两个数，整个序列改变的逆序对数为奇数。

   证明：

   

   设交换的这两个数位置是在 $i$、$j$ 位置，分别是 $x$、$y$。

   我们把序列分成三段，$A$ 段为 $1\sim (i-1)$，$B$ 段为 $(i+1)\sim (j-1)$，$C$ 段为 $(j+1)\sim n$。（如图）

   显然，交换 $x$ 和 $y$ 后，$x$ 和 $A$、$C$ 两段的相对位置没有改变，$y$ 和 $A$、$C$ 两段的相对位置也没有改变，所以由 $(A,x)$、$(A,y)$、$(x,C)$、$(y,C)$ 贡献的逆序对数没有改变。

   那么我们只需要考虑 $(x,B)$ 和 $(B,y)$ 两种情况改变的逆序对数就好了。

   设 $B$ 段中共有 $S$ 个数，其中比 $x$ 小的数有 $a$ 个，比 $y$ 大的数有 $b$ 个，那么这段中比 $x$ 大的数有 $S-a$ 个，比 $y$ 小的有 $S-b$ 个。

   可以知道交换前，这两种情况贡献的逆序对数共有 $a+b$ 个。

   交换后，这两种情况贡献的逆序对数共有 $(S-a)+(S-b)=2S-(a+b)$ 个。

   前后贡献的逆序对数之差为 $[2S-(a+b)]-(a+b)=2(S-a-b)$，为偶数。

   所以 $(x,B)$ 和 $(B,y)$ 两种情况改变的逆序对数为偶数对。

   注意到 $x$、$y$ 交换后，$(i,j)$ 可能从一开始不是逆序对变成了逆序对，也可能从一开始是逆序对变成了不是逆序对，变动对数为奇数。

   所以总的变动对数为奇数对。

   证毕。

2. 矩阵的某一行乘上 $k$，行列式也乘上 $k$。

3. $\left|\begin{array}{cc}a+a^{\prime }& b+b^{\prime }\\ c& d\end{array}\right|=\left|\begin{array}{cc}a& b\\ c& d\end{array}\right|+\left|\begin{array}{cc}{a}^{\prime }& b^{\prime }\\ c& d\end{array}\right|$

4. 若矩阵中某两行相同，则它的行列式为 $0$。

5. 用矩阵的一行加上另一行的倍数,行列式不变。

其中上面的第 $2\sim 4$ 点的证明在这篇题解有详细阐述，不会的可以直接去看。

发现第 $5$ 点操作很像高斯消元中的操作，所以可以考虑用类似高斯消元的方法把行列式变成一个右上三角矩阵，然后这个矩阵的行列式就是对角线上的乘积。

以上是关于行列式和如何求行列式（$O(n^3)$）。

接下来是矩阵树定理：

1. 无向图。设 $A$ 为邻接矩阵，$D$ 为度数矩阵（$D(i,i)$ 为节点 $i$ 的度数，其他的无值）。则基尔霍夫 （Kirchhoff）矩阵即为 : $K=D-A$。把 $K$ 去掉任意第 $i$ 行和第 $i$ 列后的矩阵 $K^{\prime }$ 的行列式即为答案。

2. 有根扩展。设根为 $rt$，则把 $K$ 去掉第 $rt$ 行和第 $rt$ 列后的矩阵 $K^{\prime }$ 的行列式即为答案。

3. 有向扩展。若为外向树，则度数矩阵统计的是入度；若为内向树，则度数矩阵统计的是出度。

4. 边权扩展。假设有一条边 $(u,v,w)$，那么可以看成是在 $(u,v)$ 之间，一共有 $w$ 条重边。

   那么我们把度数变成边权就好了。

   此时矩阵树定理求的就是 : 所有生成树边权乘积的总和。

以上扩展都可以参照这道题。

#include<bits/stdc++.h>

#define N 310
#define ll long long
#define mod 1000000007

using namespace std;

int n,m,t,tp;
ll ans=1,a[N][N];

ll poww(ll a,ll b)
{
	ll ans=1;
	while(b)
	{
		if(b&1) ans=ans*a%mod;
		a=a*a%mod;
		b>>=1;
	}
	return ans;
}

void work()//高斯消元
{
	for(int i=2;i<=n;i++)
	{
		if(!a[i][i])
		{
			for(int j=i+1;j<=n;j++)
			{
				if(a[j][i])
				{
					ans=-ans;
					swap(a[i],a[j]);
					break;
				}
			}
		}
		ll inv=poww(a[i][i],mod-2);
		for(int j=i+1;j<=n;j++)
		{
			ll tmp=a[j][i]*inv%mod;
			for(int k=i;k<=n;k++)
				a[j][k]=(a[j][k]-a[i][k]*tmp%mod+mod)%mod;
		}
	}
}

int main()
{
	scanf("%d%d%d",&n,&m,&t);
	for(int i=1;i<=m;i++)
	{
		int u,v,w;
		scanf("%d%d%d",&u,&v,&w);
		if(!t)
		{
			a[u][u]=(a[u][u]+w)%mod;
			a[v][v]=(a[v][v]+w)%mod;
			a[u][v]=(a[u][v]-w+mod)%mod;
			a[v][u]=(a[v][u]-w+mod)%mod;
		}
		else
		{
			a[v][v]=(a[v][v]+w)%mod;
			a[u][v]=(a[u][v]-w+mod)%mod;
		}
	}
	//以上是构造K矩阵
	work();
	for(int i=2;i<=n;i++)
		ans=ans*a[i][i]%mod;
	printf("%lld
",ans);
	return 0;
}

例题

[FJOI2007]轮状病毒

看题面显然是矩阵树定理。

可以使用高精度小数进行高斯消元，显然不够优雅。

发现题目给的图结构类似，不妨手推一下行列式长啥样。

发现得到的行列式如下：

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-AYqqOGyg-1604621224497)(http://192.168.102.138/JudgeOnline/upload/attachment/image/20180615/20180615085710_76849.png)]

再推一下行列式的性质，发现：$f_i=3f_{i-1}-f_{i-2}+2$（其中 $f_i$ 为 $n=i$ 的答案）。

然后用高精度递推就好了。

至于矩阵快速幂就不用了……

代码如下：

#include<bits/stdc++.h>

#define N 110

using namespace std;

struct Big_Number
{
	int len,a[50];
	Big_Number(){memset(a,0,sizeof(a));}
	Big_Number(int x)
	{
		len=0;
		memset(a,0,sizeof(a));
		while(x)
		{
			a[++len]=x%10;
			x/=10;
		}
	}
}f[N];

Big_Number operator + (Big_Number a,Big_Number b)
{
	Big_Number c;
	c.len=max(a.len,b.len);
	for(int i=1;i<=c.len;i++)
	{
		c.a[i]+=a.a[i]+b.a[i];
		if(c.a[i]>=10)
		{
			c.a[i]-=10;
			c.a[i+1]++;
		}
	}
	while(c.a[c.len+1]) c.len++;
	return c;
}

Big_Number operator - (Big_Number a,Big_Number b)
{
	Big_Number c;
	c.len=max(a.len,b.len);
	for(int i=1;i<=c.len;i++)
	{
		c.a[i]+=a.a[i]-b.a[i];
		if(c.a[i]<0)
		{
			c.a[i]+=10;
			c.a[i+1]--;
		}
	}
	while(!c.a[c.len]) c.len--;
	return c;
}

Big_Number operator * (Big_Number a,Big_Number b)
{
	Big_Number c;
	c.len=a.len+b.len-1;
	for(int i=1;i<=a.len;i++)
	{
		for(int j=1;j<=b.len;j++)
		{
			c.a[i+j-1]+=a.a[i]*b.a[j];
			c.a[i+j]+=c.a[i+j-1]/10;
			c.a[i+j-1]%=10;
		}
	}
	while(c.a[c.len+1]) c.len++;
	return c;
}

void print(Big_Number a)
{
	for(int i=a.len;i;i--)
		putchar(a.a[i]+'0');
	puts("");
}

int n;

int main()
{
	scanf("%d",&n);
	f[1]=Big_Number(1);
	f[2]=Big_Number(5);
	for(int i=3;i<=n;i++)
		f[i]=(f[i-1]*Big_Number(3))-f[i-2]+Big_Number(2);
	print(f[n]);
	return 0;
}

[HEOI2015]小 Z 的房间

题解

[SHOI2016]黑暗前的幻想乡

待填坑……

积性函数

原理

积性函数：对于一个函数 $f(x)$，若有 $f(ab)=f(a)f(b)$，其中 $\gcd (a,b)=1$，则称 $f(x)$ 为积性函数。

完全积性函数：对于一个函数 $f(x)$，若有 $f(ab)=f(a)f(b)$，则称 $f(x)$ 为完全积性函数。

常见的积性函数及其求法：

1. 欧拉函数（$\phi$）

   定义：$\phi (n)$ 是小于或等于 $n$ 的正整数中与 $n$ 互质的数的个数。

   求法：

   首先对于 $p\in prime$，肯定有 $\phi (p)=p-1$。

   然后对于 $p\in prime$，有 $\phi (p^k)=p^k-p^{k-1}$，因为除了 $p$ 的倍数以外的其他数都与 $p^k$ 互质。

   再加上它是个积性函数，就可以用线性筛求 $\phi$ 了。

2. 莫比乌斯函数（$\mu$）

   $$\mu (n)=\{\begin{array}{ll}1& n=1\\ (-1{)}^k& \text{n 无大于 1 的平方数因数，n 分解质因数为 n=p1p2⋯pk}\\ 0& \text{n 有大于 1 的平方数因数}\end{array}$$

   根据这个定义，我们就可以线性筛了。

常见的完全积性函数：

1. $ϵ(n)=[n=1]$。

2. $I(n)=1$。

3. $id(n)=n$。

例题

待填坑……

狄利克雷卷积

原理

设 $f,g$ 是两个积性函数，定义他们的狄利克雷卷积是：$(f\ast g)(n)=\sum _{d|n}f(d)g(\frac{n}{d})$。

它满足交换律、结合律。

然后对于任意的 $f$，有 $f\ast ϵ=f$。

常见的狄利克雷卷积有：

1. $\mu \ast I=ϵ$，即 $\sum _{d|n}\mu (d)=[n=1]$。

   证明：

   若 $n=1$，则易得原式等于 $1$。

   若 $n�=1$，将 $n$ 质因数分解得 $n=p_1^{a_1}{p}_2^{a_2}\cdots p_m^{a_m}$。

   设枚举到的 $d=p_1^{x_1}{p}_2^{x_2}\cdots p_m^{x_m}$

   按枚举到的 $d$ 分类讨论：

   1. 对于所有的 $x_i$，存在有至少一个 $x_i>1$。

      此时 $\mu (d)=0$。

   2. 对于所有的 $x_i$，都有 $x_i\le 1$。

      那么这部分的 $ans$ 为：

      $$ans=\sum _{i=0}^m\left(\genfrac{}{}{0ex}{}{m}{i}\right)(-1{)}^i$$

      其中第一个 $\sum _{i=0}^m$ 枚举所有的 $x_j$ 中（$1\le j\le m$），有 $i$ 个 $x_j$ 大于 $0$，即有 $i$ 个 $x_j$ 为 $1$。

      然后 $\left(\genfrac{}{}{0ex}{}{m}{i}\right)$ 意思是从所有的 $x_j$ 中（$1\le j\le m$），选出 $i$ 个 $x_j$ 为 $1$ 的方案数。

      然后 $(-1{)}^i$ 就是 $\mu (d)$。

      至于 $ans$ 为什么等于 $0$：

      1. 若 $m\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}2)$。

         则 $+\left(\genfrac{}{}{0ex}{}{m}{0}\right)$ 会和 $-\left(\genfrac{}{}{0ex}{}{m}{m}\right)$ 抵消，$-\left(\genfrac{}{}{0ex}{}{m}{1}\right)$ 会和 $+\left(\genfrac{}{}{0ex}{}{m}{m-1}\right)$ 抵消……以此类推，可得 $ans=0$。

      2. 若 $m\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}2)$。

         将每个 $\left(\genfrac{}{}{0ex}{}{m}{i}\right)$ 都画在杨辉三角里：

         

         $$\begin{array}{rl}ans=& \left(\genfrac{}{}{0ex}{}{m}{0}\right)-\left(\genfrac{}{}{0ex}{}{m}{1}\right)+\left(\genfrac{}{}{0ex}{}{m}{2}\right)-\cdots -\left(\genfrac{}{}{0ex}{}{m}{m-1}\right)+\left(\genfrac{}{}{0ex}{}{m}{m}\right)\\ =& \left(\genfrac{}{}{0ex}{}{m-1}{0}\right)-\left[\left(\genfrac{}{}{0ex}{}{m-1}{0}\right)+\left(\genfrac{}{}{0ex}{}{m-1}{1}\right)\right]+\left[\left(\genfrac{}{}{0ex}{}{m-1}{1}\right)+\left(\genfrac{}{}{0ex}{}{m-1}{2}\right)\right]-\cdots \\ & -\left[\left(\genfrac{}{}{0ex}{}{m-1}{m-2}\right)+\left(\genfrac{}{}{0ex}{}{m-1}{m-1}\right)\right]+\left(\genfrac{}{}{0ex}{}{m-1}{m-1}\right)\\ =& 0\end{array}$$

      所以无论如何，都有 $ans=0$。

   所以当 $n�=1$ 时，原式为 $0$。

   证毕。

2. $\phi \ast I=id$，即 $\sum _{d|n}\phi (d)=n$。

   证明：

   首先把式子转换一下：$\sum _{d|n}\phi (d)=\sum _{d|n}\phi (\frac{n}{d})$。

   考虑每一个 $\phi (\frac{n}{d})$ 代表着什么。

   可以把 $\phi (\frac{n}{d})$ 看成小于等于 $n$ 的所有正整数中，与 $n$ 的 $\gcd$ 为 $d$ 的个数。

   那么 $\sum _{d|n}\phi (\frac{n}{d})$ 就可以看成所有小于等于 $n$ 的正整数个数，得 $\sum _{d|n}\phi (\frac{n}{d})=n$。

   证毕。

3. $\mu \ast id=\phi$，即 $\sum _{d|n}\mu (d)\times d=\phi (n)$。

   证明：

   首先由（2）得 $\phi \ast I=id$。

   则 $\phi \ast I\ast \mu =id\ast \mu$。

   则 $\phi \ast ϵ=id\ast \mu$。

   即 $\phi =id\ast \mu$。

   证毕。

莫比乌斯反演

原理

定理：若 $g=f\ast I$，即 $g(n)=\sum _{d|n}f(d)$，则 $f=g\ast \mu$，即 $f(n)=\sum _{d|n}g(d)\mu (\frac{n}{d})$。

证明： $g=f\ast I$ 得 $g\ast \mu =f\ast I\ast \mu =f\ast ϵ=f$，证毕。

例题

其实很多莫比乌斯反演的题都可以不用莫反做。

[国家集训队]Crash的数字表格/【BZOJ2693】jzptab

题解

杜教筛

原理

求一个积性函数 $f$ 的前缀和，记 $S(n)=\sum _{i=1}^{n}f(i)$。

再找一个积性函数 $g$，先考虑它们的狄利克雷卷积的前缀和：

$$\begin{array}{rl}& \sum _{i=1}^n(f\ast g)(i)\\ =& \sum _{i=1}^n\sum _{d|i}f(d)g(\frac{i}{d})\\ =& \sum _{d=1}^{n}g(d)\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }f(i)\\ =& \sum _{d=1}^{n}g(d)S\left(\lfloor \frac{n}{d}\rfloor \right)\end{array}$$

考虑 $g(1)S(n)$ 是什么，发现：

$$\begin{array}{rl}g(1)S(n)& =\sum _{d=1}^{n}g(d)S\left(\lfloor \frac{n}{d}\rfloor \right)-\sum _{d=2}^{n}g(d)S\left(\lfloor \frac{n}{d}\rfloor \right)\\ & =\sum _{i=1}^n(f\ast g)(i)-\sum _{d=2}^{n}g(d)S\left(\lfloor \frac{n}{d}\rfloor \right)\end{array}$$

那么我们就可以找到或构造一个 $g$，使得 $\sum _{i=1}^n(f\ast g)(i)$ 和 $g$ 的前缀和可以快速算，然后数论分块递归求解。

然后可以预处理出前若干个答案的前缀和，并且加上记忆化，可以把时间复杂度优化成 $O(n^{\frac{2}{3}})$，其中要处理前 $n^{\frac{2}{3}}$ 个答案的前缀和。

代码如下：（【模板】杜教筛）

#include<bits/stdc++.h>

#define N 1664510
#define int long long
#define ll long long

using namespace std;

int t,n;
int cnt,prime[N+10],mu[N+10];
bool notprime[N+10];
ll phi[N+10];

map<int,ll>summu;
map<int,int>sumphi;

void init()
{
	mu[1]=phi[1]=1;
	for(int i=2;i<=N;i++)
	{
		if(!notprime[i])
		{
			prime[++cnt]=i;
			mu[i]=-1;
			phi[i]=i-1;
		}
		for(int j=1;j<=cnt&&i*prime[j]<=N;j++)
		{
			notprime[i*prime[j]]=1;
			if(!(i%prime[j]))
			{
				phi[i*prime[j]]=phi[i]*prime[j];
				break;
			}
			mu[i*prime[j]]=-mu[i];
			phi[i*prime[j]]=phi[i]*phi[prime[j]];
		}
	}
	for(int i=2;i<=N;i++)
	{
		mu[i]+=mu[i-1];
		phi[i]+=phi[i-1];
	}
}

ll getsummu(int n)
{
	if(n<=N) return mu[n];
	if(summu[n]) return summu[n];
	ll ans=1;
	for(int l=2,r;l<=n;l=r+1)
	{
		r=n/(n/l);
		ans-=(r-l+1)*getsummu(n/l);
	}
	return summu[n]=ans;
}

ll getsumphi(int n)
{
	if(n<=N) return phi[n];
	if(sumphi[n]) return sumphi[n];
	ll ans=1ll*n*(n+1)/2;
	for(int l=2,r;l<=n;l=r+1)
	{
		r=n/(n/l);
		ans-=(r-l+1)*getsumphi(n/l);
	}
	return sumphi[n]=ans;
}

signed main()
{
	init();
	scanf("%lld",&t);
	while(t--)
	{
		scanf("%lld",&n);
		printf("%lld %lld
",getsumphi(n),getsummu(n));
	}
	return 0;
}

例题

【BZOJ4916】神犇和蒟蒻

首先根据莫比乌斯函数的定义，第一问答案一定是 $1$。

然后考虑第二问：

首先 $\phi (i^2)=\phi (i)\times i$。

设 $f(n)=\phi (n)\times n$。

根据杜教筛的套路式：

$$g(1)S(n)=\sum _{i=1}^n(f\ast g)(i)-\sum _{d=2}^{n}g(d)S\left(\lfloor \frac{n}{d}\rfloor \right)$$

我们需要找到一个 $g$ 使得 $f\ast g$ 和 $g$ 自己的前缀和都可以快速算。

$$\sum _{i=1}^n(f\ast g)(i)=\sum _{i=1}^n\sum _{d|i}\phi (d)\times d\times g(\frac{n}{d})$$

发现把 $g$ 设为 $id$ 就好了，即 $g(n)=n$。

那么：

$$\begin{array}{rl}\sum _{i=1}^n(f\ast g)(i)& =\sum _{i=1}^n\sum _{d|i}\phi (d)\times d\times g(\frac{i}{d})\\ & =\sum _{i=1}^n\sum _{d|i}\phi (d)\times d\times \frac{i}{d}\\ & =\sum _{i=1}^n\sum _{d|i}\phi (d)\times i\\ & =\sum _{i=1}^{n}i\sum _{d|i}\phi (d)\\ & =\sum _{i=1}^n{i}^2\\ & =\frac{n(n+1)(2n+1)}{6}\end{array}$$

然后用杜教筛直接做就好了。

代码如下：

#include<bits/stdc++.h>
 
#define N 1000010
#define ll long long
#define mod 1000000007
 
using namespace std;
 
const int maxn=1000000;
 
int n;
int cnt,prime[N];
bool notprime[N];
ll phi[N];
 
map<int,ll>sum;
 
ll poww(ll a,ll b)
{
    ll ans=1;
    while(b)
    {
        if(b&1) ans=ans*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return ans;
}
 
void init()
{
    phi[1]=1;
    for(int i=2;i<=maxn;i++)
    {
        if(!notprime[i])
        {
            prime[++cnt]=i;
            phi[i]=i-1;
        }
        for(int j=1;j<=cnt&&i*prime[j]<=maxn;j++)
        {
            notprime[i*prime[j]]=1;
            if(!(i%prime[j]))
            {
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            }
            phi[i*prime[j]]=phi[i]*phi[prime[j]];
        }
    }
    for(int i=1;i<=maxn;i++)
        phi[i]=(phi[i-1]+i*phi[i]%mod)%mod;
}
 
ll dfs(int n)
{
    if(n<=maxn) return phi[n];
    if(sum[n]) return sum[n];
    ll ans=1ll*n*(n+1)%mod*(2*n+1)%mod*poww(6,mod-2)%mod;
    for(int l=2,r;l<=n;l=r+1)
    {
        r=n/(n/l);
        ans=(ans-1ll*(l+r)*(r-l+1)/2%mod*dfs(n/l)%mod+mod)%mod;
    }
    return sum[n]=ans;
}
 
int main()
{
    init();
    scanf("%d",&n);
    printf("1
%lld",dfs(n));
    return 0;
}

[SDOI2015]约数个数和/【BZOJ4176】Lucas的数论（前者是弱化版）

求 $\sum _{i=1}^n\sum _{j=1}^{m}d(ij)$，其中 $d(n)$ 为 $n$ 的因数个数。

先证明一个式子：$d(ij)=\sum _{x|i}\sum _{y|j}[\gcd (x,y)=1]$。

证明：

将 $ij$ 质因数分解，设 $ij=p_1^{c_1}{p}_2^{c_2}\cdots p_k^{c_k}$，$i=p_1^{a_1}{p}_2^{a_2}\cdots p_k^{a_k}$，$j=p_1^{b_1}{p}_2^{b_2}\cdots p_k^{b_k}$。那么易知对于任意的 $l$，均有 $a_l+b_l=c_l$。

我们将每一个质因数 $p_l$，我们分开考虑，因为质因数之间是互不影响的。

那么等式左边的贡献是 $p_l^{c_l}$ 的因数个数，即 $c_l+1$。

等式右边可以看成枚举 $p_l^{a_l}$ 的因子 $x$ 和 $p_l^{b_l}$ 的因子 $y$ 使得 $x,y$ 互质。那么只有 $3$ 种情况：$x=1$，$y=p_l^1,p_l^2,\cdots ,p_l^{b_l}$，共 $b_l$ 种情况；$y=1$，$x=p_l^1,p_l^2,\cdots ,p_l^{a_l}$，共 $a_l$ 种情况；$x=1$，$y=1$，共 $1$ 种情况。那么右边的总贡献为 $a_l+b_l+1=c_l+1$。

所以左右两边的贡献是相等的。

证毕。（其实证明比较感性）

那么就可以开始推式子了：（不妨设 $n<m$）

$$\begin{array}{rl}& \sum _{i=1}^n\sum _{j=1}^{m}d(ij)\\ =& \sum _{i=1}^n\sum _{j=1}^m\sum _{x|i}\sum _{y|j}[\gcd (x,y)=1]\\ =& \sum _{x=1}^n\sum _{y=1}^m[\gcd (x,y)=1]\sum _{x|i}\sum _{y|j}1\\ =& \sum _{x=1}^n\sum _{y=1}^m\lfloor \frac{n}{x}\rfloor \lfloor \frac{m}{y}\rfloor [\gcd (x,y)=1]\\ =& \sum _{x=1}^n\sum _{y=1}^m\lfloor \frac{n}{x}\rfloor \lfloor \frac{m}{y}\rfloor \sum _{d|\gcd (x,y)}\mu (d)\\ =& \sum _{d=1}^n\mu (d)\sum _{d|x}\lfloor \frac{n}{x}\rfloor \sum _{d|y}\lfloor \frac{m}{y}\rfloor \\ =& \sum _{d=1}^n\mu (d)\sum _{a=1}^{\lfloor \frac{n}{d}\rfloor }\lfloor \frac{n}{ad}\rfloor \sum _{b=1}^{\lfloor \frac{m}{d}\rfloor }\lfloor \frac{m}{bd}\rfloor \end{array}$$

式子就推到这了，弱化版和强化版的做法不一样。

弱化版：$T,n\le 5\times 10^4$。

$O(n)$ 预处理 $\mu$，$O(n\sqrt{n})$ 预处理 $g(n)=\sum _{i=1}^n\lfloor \frac{n}{i}\rfloor$。

然后 $O(T\sqrt{n})$ 询问。

代码如下：

#include<bits/stdc++.h>

#define N 50010
#define ll long long

using namespace std;

int t,n,m;
int cnt,prime[N],mu[N],sum[N];
bool notprime[N];
ll g[N];

void init()
{
	mu[1]=1;
	for(int i=2;i<=50000;i++)
	{
		if(!notprime[i])
		{
			prime[++cnt]=i;
			mu[i]=-1;
		}
		for(int j=1;j<=cnt&&i*prime[j]<=50000;j++)
		{
			notprime[i*prime[j]]=1;
			if(!(i%prime[j])) break;
			mu[i*prime[j]]=-mu[i];
		}
	}
	for(int i=1;i<=50000;i++)
		sum[i]=sum[i-1]+mu[i];
	for(int i=1;i<=50000;i++)
	{
		for(int l=1,r;l<=i;l=r+1)
		{
			r=i/(i/l);
			g[i]+=1ll*(r-l+1)*(i/l);
		}
	}
}

int main()
{
	init();
	scanf("%d",&t);
	while(t--)
	{
		scanf("%d%d",&n,&m);
		if(n>m) swap(n,m);
		ll ans=0;
		for(int l=1,r;l<=n;l=r+1)
		{
			r=min(n/(n/l),m/(m/l));
			ans+=1ll*(sum[r]-sum[l-1])*g[n/l]*g[m/l];
		}
		printf("%lld
",ans);
	}
	return 0;
}

强化版：$T=1,n\le 10^9$。

用杜教筛来优化 $\mu (n)$ 的前缀和。

#include<bits/stdc++.h>
 
#define N 1000010
#define ll long long
#define mod 1000000007
 
using namespace std;
 
int n;
int cnt,mu[N],sum[N],prime[N];
bool notprime[N];
 
map<int,int>mp;
 
void init()
{
    mu[1]=1;
    for(int i=2;i<=1e6;i++)
    {
        if(!notprime[i])
        {
            prime[++cnt]=i;
            mu[i]=-1;
        }
        for(int j=1;j<=cnt&&i*prime[j]<=1e6;j++)
        {
            notprime[i*prime[j]]=1;
            if(!(i%prime[j])) break;
            mu[i*prime[j]]=-mu[i];
        }
    }
    for(int i=1;i<=1e6;i++)
        sum[i]=sum[i-1]+mu[i];
}
 
ll calc(int n)
{
    ll ans=0;
    for(int l=1,r;l<=n;l=r+1)
    {
        r=n/(n/l);
        ans=(ans+1ll*(r-l+1)*(n/l))%mod;
    }
    return ans;
}
 
int summu(int n)
{
    if(n<=1e6) return sum[n];
    if(mp[n]) return mp[n];
    int ans=1;
    for(int l=2,r;l<=n;l=r+1)
    {
        r=n/(n/l);
        ans-=(r-l+1)*summu(n/l);
    }
    return mp[n]=ans;
}
 
int main()
{
    init();
    scanf("%d",&n);
    ll ans=0;
    for(int l=1,r;l<=n;l=r+1)
    {
        r=n/(n/l);
        ll tmp=calc(n/l);
        ans=(ans+(0ll+summu(r)-summu(l-1)+mod)%mod*tmp%mod*tmp%mod)%mod;
    }
    printf("%lld
",ans);
    return 0;
}

多项式

原根

定义：若 $g^{\phi (p)}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$，且 $\phi (p)$ 是 $g^x\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$ 的最小正整数解，则称 $g$ 为 $p$ 的原根。

若有 $p\in prime$，则 $\phi (p)=p-1$，那么 $g^{p-1}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$，那么有 $g^0\mathrm{m}\mathrm{o}\mathrm{d}p,g^1\mathrm{m}\mathrm{o}\mathrm{d}p,\cdots ,g^{p-2}\mathrm{m}\mathrm{o}\mathrm{d}p$ 与 $1\sim (p-1)$ 一一映射。（不一定 $g^i$ 与 $i$ 对应）

证明：

首先对于 $0\le i\le p-2$，肯定有 $g^i\mathrm{m}\mathrm{o}\mathrm{d}p\in [1,p-1]$。

因为如果 $g^i\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}p)$，那么对于 $j\ge i$，肯定有 $g^j\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}p)$，那么 $g^{\phi (p)}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$ 肯定不成立。

然后若存在 $g^i\equiv g^j(\mathrm{m}\mathrm{o}\mathrm{d}p)$（$0\le j<i\le p-2$），那么 $g^{i-j}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$，又 $i-j<\phi (p)$，则与原根的定义违背。

所以 $g^0\mathrm{m}\mathrm{o}\mathrm{d}p,g^1\mathrm{m}\mathrm{o}\mathrm{d}p,\cdots ,g^{p-2}\mathrm{m}\mathrm{o}\mathrm{d}p$ 与 $1\sim (p-1)$ 一一映射。

证毕。

多项式求逆

给定多项式 $F(x)$，求多项式 $G(x)$ 使得 $F(x)G(x)\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$，其中 $F(x)G(x)$ 的系数对一个数 $p$ 取模。（对于一个多项式 $F(x)$，$F(x)\mathrm{m}\mathrm{o}\mathrm{d}{x}^n$ 的意思是只看 $F(x)$ 中的第 $0$ 项到第 $n-1$ 项）

用递归求解。（以下的所有多项式都是在系数模了 $p$ 意义下的）

首先假设我们已经求出了 $F(x)H(x)\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\lceil \frac{n}{2}\rceil })$（记为 $1$ 式），现在要求 $F(x)G(x)\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$。

首先如果有 $F(x)G(x)\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$，那么一定有 $F(x)G(x)\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\lceil \frac{n}{2}\rceil })$（记为 $2$ 式），因为等号右边的那个 $1$ 只能在常数项出现，其他项的系数都是 $0$。

由 $1$ 式减 $2$ 式得 $F(x)\left(H(x)-G(x)\right)\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\lceil \frac{n}{2}\rceil })$。

由于 $F(x)�\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\lceil \frac{n}{2}\rceil })$，所以 $H(x)-G(x)\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\lceil \frac{n}{2}\rceil })$。

然后等式两边平方，因为 $H(x)-G(x)\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\frac{n}{2}})$，即 $H(x)-G(x)$ 的第 $0$ 项到第 $\lceil \frac{n}{2}\rceil -1$ 项系数都是 $0$，那么 $(H(x)-G(x){)}^2$ 肯定满足第 $0$ 项到第 $n-1$ 项系数都是 $0$，所以有：

$$H(x{)}^2+G(x{)}^2-2G(x)H(x)\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$$

两边同时乘 $F(x)$ 得：

$$\begin{array}{rl}F(x)H(x{)}^2+F(x)G(x{)}^2-2F(x)G(x)H(x)& \equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)\\ F(x)H(x{)}^2+G(x)-2H(x)& \equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)\end{array}$$

上面第二步的化简是用 $F(x)G(x)\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$ 化简的，注意不能用 $F(x)H(x)\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\lceil \frac{n}{2}\rceil })$，因为两个式子 $\mathrm{m}\mathrm{o}\mathrm{d}$ 的东西不一样。

那么就可以得到：$G(x)=2H(x)-F(x)H(x{)}^2$。

那么我们就能通过 $\lceil \frac{n}{2}\rceil$ 的答案求出 $n$ 时的答案了，所以一直向下递归。

最后当递归到 $1$ 的时候，要求 $F(x)G(x)\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}{x}^1)$，也就是 $F(x)$ 的常数项和 $G(x)$ 的常数项乘起来模 $p$ 为 $1$，那么只用求 $F(x)$ 常数项的逆元就好了。

总时间复杂度为：（设 $n=2^m$，以下的 $\log$ 均以 $2$ 为底）

$$\sum _{i=0}^{m-1}{2}^i\log (2^i)<\sum _{i=0}^{m-1}{2}^i\log (n)=\left(2^{m-1+1}-1\right)\log (n)=(n-1)\log (n)<n\log n$$

代码如下：

#include<bits/stdc++.h>

#define N 200010
#define ll long long

using namespace std;

const int G=3,mod=998244353;

int n,rev[N<<1];
ll a[N<<1],b[N<<1],c[N<<1];

ll poww(ll a,ll b)
{
	ll ans=1;
	while(b)
	{
		if(b&1) ans=ans*a%mod;
		a=a*a%mod;
		b>>=1;
	}
	return ans;
}

void NTT(ll *a,int limit,int opt)
{
	for(int i=0;i<limit;i++)
		if(i<rev[i])
			swap(a[i],a[rev[i]]);
	for(int mid=1;mid<limit;mid<<=1)
	{
		int len=mid<<1;
		ll gn=poww(G,(mod-1)/len);
		for(int i=0;i<limit;i+=len)
		{
			ll g=1;
			for(int j=0;j<mid;j++,g=g*gn%mod)
			{
				ll x=a[i+j],y=a[i+mid+j];
				a[i+j]=(x+g*y)%mod,a[i+mid+j]=((x-g*y)%mod+mod)%mod;
			}
		}
	}
	if(opt==-1)
	{
		reverse(a+1,a+limit);
		ll tmp=poww(limit,mod-2);
		for(int i=0;i<limit;i++)
			a[i]=a[i]*tmp%mod;
	}
}

void solve(int n)
{
	if(n==1)
	{
		b[0]=poww(a[0],mod-2);
		return;
	}
	solve(ceil(1.0*n/2.0));
	int bit=0,limit=1;
	while(limit<(n<<1)) 
		bit++,limit<<=1;
	for(int i=0;i<limit;i++)
		rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
	for(int i=0;i<n;i++)
		c[i]=a[i];
	for(int i=n;i<limit;i++)
		c[i]=0;
	NTT(b,limit,1),NTT(c,limit,1);
	for(int i=0;i<limit;i++)
		b[i]=((2*b[i]-c[i]*b[i]%mod*b[i]%mod)%mod+mod)%mod;
	NTT(b,limit,-1);
	for(int i=n;i<limit;i++) 
		b[i]=0;
}

int main()
{
	scanf("%d",&n);
	for(int i=0;i<n;i++)
		scanf("%lld",&a[i]);
	solve(n);
	for(int i=0;i<n;i++)
		printf("%lld ",b[i]);
	return 0;
}

多项式开根

跟多项式求逆差不多，用递归求解。

假设我们现在要求 $G(x)$ 使得 $G(x{)}^2\equiv F(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$，已经求出了 $H(x{)}^2\equiv F(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\lceil \frac{n}{2}\rceil })$。（记为 $1$ 式）

显然若 $G(x{)}^2\equiv F(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$，则 $G(x)\equiv \sqrt{F(x)}\equiv H(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\lceil \frac{n}{2}\rceil })$。

$$\begin{array}{rl}G(x)& \equiv H(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\lceil \frac{n}{2}\rceil })\\ G(x)-H(x)& \equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\lceil \frac{n}{2}\rceil })\\ (G(x)-H(x){)}^2& \equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)\\ G(x{)}^2+H(x{)}^2-2G(x)H(x)& \equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)\\ F(x)+H(x{)}^2-2G(x)H(x)& \equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)\\ G(x)& \equiv \frac{F(x)+H(x{)}^2}{2H(x)}(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)\end{array}$$

然后就能递归求逆做了。

时间复杂度 $O(n\log n)$。

代码如下：

#include<bits/stdc++.h>

#define N 400010
#define ll long long
#define mod 998244353

using namespace std;

int n,rev[N];
ll inv2,f[N],g[N],ff[N],gg[N],fff[N];

ll poww(ll a,ll b)
{
	ll ans=1;
	while(b)
	{
		if(b&1) ans=ans*a%mod;
		a=a*a%mod;
		b>>=1;
	}
	return ans;
}

void NTT(int limit,int opt,ll *a)
{
	for(int i=0;i<limit;i++)
		if(i<rev[i])
			swap(a[i],a[rev[i]]);
	for(int mid=1;mid<limit;mid<<=1)
	{
		int len=mid<<1;
		ll gn=poww(3,(mod-1)/len);
		if(opt==-1) gn=poww(gn,mod-2);
		for(int i=0;i<limit;i+=len)
		{
			ll omega=1;
			for(int j=0;j<mid;j++,omega=omega*gn%mod)
			{
				ll x=a[i+j],y=a[i+mid+j];
				a[i+j]=(x+omega*y)%mod,a[i+mid+j]=((x-omega*y)%mod+mod)%mod;
			}
		}
	}
	if(opt==-1)
	{
		ll tmp=poww(limit,mod-2);
		for(int i=0;i<limit;i++)
			a[i]=a[i]*tmp%mod;
	}
}

void Inv(int n,ll *f,ll *g)
{
	g[0]=poww(f[0],mod-2);
	int bit=0,limit,now;
	for(now=1;now<(n<<1);now<<=1)
	{
		bit++,limit=now<<1;
		for(int i=0;i<now;i++)
			ff[i]=f[i];
		for(int i=0;i<limit;i++)
			rev[i]=(rev[i>>1]>>1)|((i&1)?now:0);
		NTT(limit,1,ff),NTT(limit,1,g);
		for(int i=0;i<limit;i++)
			g[i]=((2*g[i]-ff[i]*g[i]%mod*g[i]%mod)%mod+mod)%mod;
		NTT(limit,-1,g);
		for(int i=now;i<limit;i++)
			ff[i]=g[i]=0;
	}
	for(int i=0;i<now;i++) ff[i]=0;
	for(int i=n;i<now;i++) g[i]=0;
}

void Sqrt(int n,ll *f,ll *g)
{
	g[0]=sqrt(f[0]);
	int bit=0,limit,now;
	for(now=1;now<(n<<1);now<<=1)
	{
		bit++,limit=now<<1;
		Inv(now,g,gg);
		for(int i=0;i<now;i++)
			fff[i]=f[i];
		for(int i=0;i<limit;i++)
			rev[i]=(rev[i>>1]>>1)|((i&1)?now:0);
		NTT(limit,1,fff),NTT(limit,1,g),NTT(limit,1,gg);
		for(int i=0;i<limit;i++)
			g[i]=inv2*gg[i]%mod*((fff[i]+g[i]*g[i])%mod)%mod;
		NTT(limit,-1,g);
		for(int i=0;i<limit;i++)
			gg[i]=0;
		for(int i=now;i<limit;i++)
			fff[i]=g[i]=0;
	}
	for(int i=0;i<now;i++) fff[i]=0;
	for(int i=n;i<now;i++) g[i]=0;
}

int main()
{
	inv2=poww(2,mod-2);
	scanf("%d",&n);
	for(int i=0;i<n;i++)
		scanf("%lld",&f[i]);
	Sqrt(n,f,g);
	for(int i=0;i<n;i++)
		printf("%lld ",g[i]);
	return 0;
}

多项式求导

设 $f(x)=\sum _{i\ge 0}{a}_i{x}^i$，则：

$$\begin{array}{rl}{f}^{\prime }(x)=& \underset{d\to 0}{\lim }\sum _{i\ge 0}\frac{a_i(x+d{)}^i-a_i{x}^i}{d}\\ =& \underset{d\to 0}{\lim }\sum _{i\ge 0}{a}_i\frac{(x+d{)}^i-x^i}{d}\end{array}$$

我们现在只看 $\frac{(x+d{)}^i-x^i}{d}$ 怎么求：

$$\begin{array}{rl}& \frac{(x+d{)}^i-x^i}{d}\\ =& \frac{\sum _{j=0}^i\left(\genfrac{}{}{0ex}{}{i}{j}\right)x^j{d}^{i-j}-x^i}{d}\\ =& \frac{\sum _{j=0}^{i-1}\left(\genfrac{}{}{0ex}{}{i}{j}\right)x^j{d}^{i-j}+x^i{d}^0-x^i}{d}\\ =& \frac{\sum _{j=0}^{i-1}\left(\genfrac{}{}{0ex}{}{i}{j}\right)x^j{d}^{i-j}}{d}\\ =& \sum _{j=0}^{i-1}\left(\genfrac{}{}{0ex}{}{i}{j}\right)x^j{d}^{i-j-1}\end{array}$$

则 $j<i-1$ 时，$i-j-1>0$，又由于 $d$ 趋近于 $0$，则 $d^{i-j-1}=0$。

所以：

$$\sum _{j=0}^{i-1}\left(\genfrac{}{}{0ex}{}{i}{j}\right)x^j{d}^{i-j-1}=\left(\genfrac{}{}{0ex}{}{i}{i-1}\right)x^{i-1}{d}^0=i{x}^{i-1}$$

那么：

$$\begin{array}{rl}{f}^{\prime }(x)=& \underset{d\to 0}{\lim }\sum _{i\ge 0}{a}_i\frac{(x+d{)}^i-x^i}{d}\\ =& \sum _{i\ge 1}{a}_{i}i{x}^{i-1}\end{array}$$

同样的，若知道 $f^{\prime }(x)=\sum _{i\ge 0}{a}_i{x}^i$，那么 $f(x)=\sum _{i\ge 1}\frac{a_{i-1}}{i}{x}^i$。

一些特殊形式的求法

一、求 $\sum _{k=1}^n{x}^k\sum _{i-j=k}{a}_i{b}_j$

平常的 FFT 是求 $\sum _{k=1}^n{x}^k\sum _{i+j=k}{a}_i{b}_j$，那么我现在让你求 $\sum _{k=1}^n{x}^k\sum _{i-j=k}{a}_i{b}_j$。

设 $S(k)=\sum _{i-j=k}{a}_i{b}_j$，那么题目是要求 $\sum _{k=1}^n{x}^{k}S(k)$。

设 $\widehat{a_i}=a_{n-i}$，那么：

$$\begin{array}{rl}S(k)=& \sum _{i-j=k}{a}_i{b}_j\\ =& \sum _{(n-i)-j=k}\widehat{a_i}{b}_j\\ =& \sum _{i+j=n-k}\widehat{a_i}{b}_j\end{array}$$

那么：

$$S(n-k)=\sum _{i+j=n-(n-k)}\widehat{a_i}{b}_j=\sum _{i+j=k}\widehat{a_i}{b}_j$$

那么就可以设多项式 $A=\{\widehat{a_i}\}$，$B=\{b_i\}$，$C=A\times B$，则 $S(n-k)=[x^k]C$。

所以将 $C$ 翻转一下就是我们要求的了。

例题

[SDOI2015]序列统计

首先 $x>1$ 所以数列里不可能取 $0$。

然后又考虑到原根的性质：$g^0\mathrm{m}\mathrm{o}\mathrm{d}m,g^1\mathrm{m}\mathrm{o}\mathrm{d}m,\cdots ,g^{m-2}\mathrm{m}\mathrm{o}\mathrm{d}m$ 与 $1\sim (m-1)$ 一一映射，所以考虑用 $g^0\mathrm{m}\mathrm{o}\mathrm{d}m,g^1\mathrm{m}\mathrm{o}\mathrm{d}m,\cdots ,g^{m-2}\mathrm{m}\mathrm{o}\mathrm{d}m$ 代替 $1\sim (m-1)$。

假设 $x$ 对应的是 $g^y$，题目给定的 $S=\{a_1,a_2,\cdots ,a_s\}$，它们对应的集合是 $\{g^{b_1},g^{b_2},\cdots ,g^{b_s}\}$，设集合 $T=\{b_1,b_2,\cdots ,b_s\}$。

原来的问题是：从 $S$ 中选 $n$ 个可以重复的数排成一个序列，使得这 $n$ 个数的乘积为 $x$ 的方案数。

那么现在的问题是：从 $T$ 中选 $n$ 个可以重复的数排成一个序列，使得这 $n$ 个数的加和为 $y$ 的方案数。

那么这就可以用 NTT 来搞了。

具体过程就是设 $c_i$ 表示 $i$ 这个数有没有在 $T$ 中出现，设 $c_i$ 的生成函数为 $F(x)$，那么 $F(x{)}^n$ 就能表示这个加和的过程了，用快速幂实现。

#include<bits/stdc++.h>
 
#define M 8010
#define ll long long
 
using namespace std;
 
const int mod=1004535809;
 
int n,m,x,size;
int G,mp[M];
int limit,bit,rev[M<<2];
ll f[M<<2],ff[M<<2],g[M<<2];
bool vis[M];
 
ll poww(ll a,ll b,ll mod)
{
    ll ans=1;
    while(b)
    {
        if(b&1) ans=ans*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return ans;
}
 
void init()
{
    for(int i=2;;i++)
    {
    	memset(vis,0,sizeof(vis));
    	int tmp=1;
    	bool flag=true;
        for(int j=0;j<=m-2;j++)
        {
        	if(!tmp||vis[tmp])
        	{
        		flag=false;
        		break;
			}
			vis[tmp]=1;
        	tmp=tmp*i%m;
		}
		if(flag)
		{
			G=i;
			tmp=1;
			for(int j=0;j<=m-2;j++)
			{
				mp[tmp]=j;
				tmp=tmp*G%m;
			}
			break;
		}
    }
}
 
void NTT(ll *a,int opt)
{
    for(int i=0;i<limit;i++)
        if(i<rev[i])
            swap(a[i],a[rev[i]]);
    for(int mid=1;mid<limit;mid<<=1)
    {
        int len=mid<<1;
        ll gn=poww(3,(mod-1)/len,mod);//这里不用G而是用3是因为这里的原根是针对模数的，而不是m（一直50分没调出来就是在这qaq）
        if(opt==-1) gn=poww(gn,mod-2,mod);
        for(int i=0;i<limit;i+=len)
        {
            ll omega=1;
            for(int j=0;j<mid;j++,omega=omega*gn%mod)
            {
                ll x=a[i+j],y=omega*a[i+mid+j]%mod;
                a[i+j]=(x+y)%mod,a[i+mid+j]=(x-y+mod)%mod;
            }
        }
    }
    if(opt==-1)
    {
        ll tmp=poww(limit,mod-2,mod);
        for(int i=0;i<limit;i++)
            a[i]=a[i]*tmp%mod;
    }
}
 
void mul(ll *f,ll *g)
{
    NTT(f,1),NTT(g,1);
    for(int i=0;i<limit;i++)
        f[i]=f[i]*g[i]%mod;
    NTT(f,-1),NTT(g,-1);
    for(int i=0;i<m;i++)
        f[i]=(f[i]+f[i+m])%mod;
    for(int i=m;i<limit;i++)
        f[i]=g[i]=0;
}
 
void polypow(ll *f,ll *g,int b)
{
    limit=1,bit=0;
    while(limit<(m<<1))
        limit<<=1,bit++;
    for(int i=0;i<limit;i++) 
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
    while(b)
    {
        if(b&1) mul(g,f);
        for(int i=0;i<limit;i++)
			ff[i]=f[i];
        mul(f,ff);
        b>>=1;
    }
}

int main()
{
    scanf("%d%d%d%d",&n,&m,&x,&size);
    init();
    for(int i=1;i<=size;i++)
    {
        int x;
        scanf("%d",&x);
        if(!x) continue;
        f[mp[x]]++;
    }
    m--;
    g[0]=1;
    polypow(f,g,n);
    printf("%lld
",g[mp[x]]);
    return 0;
}

拉格朗日插值

原理

给出 $n$ 个点 $(x_i,y_i)$，可知这 $n$ 个点唯一确定一个多项式 $f(x)$。

现在询问 $k$，求 $f(k)$。

用高斯消元可以做到预处理 $O(n^3)$，每次询问 $O(n)$。

但是现在只有一次询问，且 $n$ 很大。

所以我们考虑用询问更快的方法求解。

拉格朗日插值是通过构造的方法，在 $O(n^2)$ 的时间内进行询问。

设拉格朗日基本多项式为：

$${\ell }_i(x)=\prod _{j=1,j�=i}^n\frac{x-x_j}{x_i-x_j}$$

注意到当 $x=x_i$ 时，${\ell }_i(x)=1$；当 $x=x_j$（$j�=i$）时，${\ell }_i(x)=0$。

那么就能把这个多项式 $f(x)$ 构造出来了：

$$f(x)=\sum _{i=1}^n{y}_i{\ell }_i(x)$$

显然，这个多项式满足对于任意的 $i$，$f(x_i)=y_i$。

那么这个多项式就是我们要找的多项式。

#include<bits/stdc++.h>

#define N 2010
#define ll long long
#define mod 998244353

using namespace std;

int n;
ll k,x[N],y[N];

ll poww(ll a,ll b)
{
	ll ans=1;
	while(b)
	{
		if(b&1) ans=ans*a%mod;
		a=a*a%mod;
		b>>=1;
	}
	return ans;
}

int main()
{
	scanf("%d%lld",&n,&k);
	for(int i=1;i<=n;i++)
		scanf("%lld%lld",&x[i],&y[i]);
	ll ans=0;
	for(int i=1;i<=n;i++)
	{
		ll sum=1;
		for(int j=1;j<=n;j++)
		{
			if(i==j) continue;
			sum=sum*(k-x[j])%mod*poww(x[i]-x[j],mod-2)%mod;
		}
		ans=(ans+y[i]*sum)%mod;
	}
	printf("%lld
",(ans+mod)%mod);
	return 0;
}