數學小結

前言：這篇文章專門整理了各種有關數學的東西（包括原理和例題），主要是怕自己忘了。

整理的不齊全，而且寫的都是比較基礎的東西。

擴充套件歐幾里得（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\notin 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\ne i}^k{a}_j{m}_j{t}_j+a_i{m}_i{t}_i$。

顯然，對於任意的 $j\ne 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\ne 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\ne 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\n",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\n",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\n",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\n",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\n",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\n",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\ne 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\ne 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\n",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\n%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\n",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\n",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)\not\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 = { a i ^ } A=\{\widehat{a_i}\} A={ai​ ​}， B = { b i } B=\{b_i\} B={bi​}，$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 , ⋯   , a s } S=\{a_1,a_2,\cdots,a_{s}\} S={a1​,a2​,⋯,as​}，它們對應的集合是 { g b 1 , g b 2 , ⋯   , g b s } \{g^{b_1},g^{b_2},\cdots,g^{b_s}\} {gb1​,gb2​,⋯,gbs​}，設集合 T = { b 1 , b 2 , ⋯   , b s } T=\{b_1,b_2,\cdots,b_s\} T={b1​,b2​,⋯,bs​}。

原來的問題是：從 $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\n",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\ne i}^n\frac{x-x_j}{x_i-x_j}$$

注意到當 $x=x_i$ 時，${\ell }_i(x)=1$；當 $x=x_j$（$j\ne 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\n",(ans+mod)%mod);
	return 0;
}