关于数论

整除性理论

*P1835 素数密度（区间筛）

给定两个整数\(L,R(1<=L<=R<=2^{31},R-L<=10^6)\)，求闭区间 \([L,R]\) 中的质数个数。

推论 1：若\(a\)为合数，则\(a\)能被表示成\(a=pq\)，其中 \(p,q>1\)。易证\(p,q\)中一定有一个不超过\(\sqrt a\)（想象一下若\(p,q\)都超过\(\sqrt a\)则有\(a<pq\)）。更严格地，若\(a\)为合数，则一定存在质数 \(p|a\)，且\(p\leq \sqrt a\)。

乍看这题的数据范围大的吓人，只有\(R-L\)的范围比较正常。但若是有推论 1和埃氏筛的铺垫，便不难得到一个正确的突破口：理论上只要利用不超过 $\sqrt {2147483647}\leq 50000 $ 的质数就可以筛除题目范围内的所有合数。

事实上，再联系到\(R-L\leq1000000\)和筛法流程，一个区间筛的算法大致成型：预先求出 50000 以内的所有质数（发现只有不到6000个），再对每个质数进行一遍\([L,R]\)上的合数筛除。最后统计一下没被筛过的就是质数了。

因为只需要利用 a[L..R] 的数，所以代码中将这段区间平移成了 a[0…R-L]，是一个节省空间的好方法。

#include<bits/stdc++.h>
using namespace std;
#define maxn 1000010
int n,N,x,pri[10000];
typedef long long LL;
bool a[maxn];
int Eratosthenes_Sieve(int n,int pri[])
{
    for(int i=2;i*i<=n;i++)
    	if(a[i]==0)
    		for(int j=i<<1;j<=n;j+=i)
    			a[j]=1;
    int cnt=0;
    for(int i=2;i<=n;i++)
        if(!a[i])pri[cnt++]=i;
    return cnt;
}
int main()
{
    int cnt=Eratosthenes_Sieve(50000,pri);
    LL L,R;
    scanf("%lld %lld",&L,&R);
    memset(a,0,sizeof(a));
    for(int i=0;i<cnt;i++)
        for(LL j=max(2ll,(L-1)/pri[i]+1)*pri[i]/*找到[L,R]上第一个开始筛的数*/;j<=R;j+=pri[i])
            a[j-L]=1;
    int ans=0;
    for(LL i=L;i<=R;i++)
        if(a[i-L]==0)ans++;
    printf("%d",ans);
    return 0;
}

*P3601 签到题（区间筛plus）

我们定义一个函数：qiandao(x)为小于等于x的数中与x不互质的数的个数。

这题作为签到题，给出l和r，要求求\(\sum_{i=l}^r qiandao(i)~mod~666623333\)。

对于100%的数据，\(1 \leq l \leq r \leq 10^{12}\)，\(r-l \leq 10^6\)。

和上面一题一样，我们希望利用埃氏筛筛出\(\varphi(l..r)\)。

地球人都知道一个数\(x\)的质因子\(p\)对\(\varphi(x)\)的贡献是\(*=\frac{p-1}{p}\)

只需要特判大于平方根的因子即可。

#include <bits/stdc++.h>
using namespace std;
#define maxn 1000010
#define mod 666623333
int n, N, pri[maxn];
typedef long long LL;
bool a[maxn];
LL phi[maxn], K[maxn];
int Eratosthenes_Sieve(int n, int pri[])
{
    for (int i = 2; i * i <= n; i++)
        if (a[i] == 0)
            for (int j = i << 1; j <= n; j += i)
                a[j] = 1;
    int cnt = 0;
    for (int i = 2; i <= n; i++)
        if (!a[i])
            pri[cnt++] = i;
    return cnt;
}
int main()
{
    int cnt = Eratosthenes_Sieve(1000000, pri);
    LL L, R, ans = 0;
    scanf("%lld %lld", &L, &R);
    for (LL i = L; i <= R; i++)
        phi[i - L] = K[i - L] = i;
    for (int i = 0; i < cnt; i++)
        for (LL p = pri[i], j = max(2ll, (L - 1) / p + 1) * p; j <= R; j += p)
        {
            LL x = j - L;
            phi[x] = phi[x] / p * (p - 1);
            for (; K[x] % p == 0;)
                K[x] /= p;
        }
    for (LL i = L; i <= R; (ans += i - phi[i - L]) %= mod, i++)
        if (K[i - L] != 1)
            phi[i - L] = phi[i - L] / K[i - L] * (K[i - L] - 1);
    printf("%lld", ans);
    return 0;
}

*UOJ#48 [UR#3]核聚变反应强度 （枚举）

著名核物理专家 Picks 提出了核聚变特征值这一重要概念。

核聚变特征值分别为 \(x\) 和 \(y\) 的两个原子进行核聚变，能产生数值为 \(sgcd(x,y)\) 的核聚变反应强度。

其中， \(sgcd(x,y)\) 表示 \(x\) 和 \(y\) 的次大公约数，即能同时整除 \(x,y\) 的正整数中第二大的数。如果次大公约数不存在则说明无法核聚变， 此时 \(sgcd(x,y)=?1\)。

现在有 \(n\) 个原子，核聚变特征值分别为 \(a_1,a_2,…,a_n\)。然后 Picks 又从兜里掏出一个核聚变特征值为 \(a_1\) 的原子，你需要计算出这个原子与其它 \(n\) 个原子分别进行核聚变反应时的核聚变反应强度，即 \(sgcd(a_1,a_1),sgcd(a_1,a_2),…,sgcd(a_1,a_n)\)。

对于 \(100\%\) 的数据，\(n\leq10^5,a_i\leq10^{12}\)

性质思考：

• \(x=\prod p_i^{\alpha_i},y=\prod p_i^{\beta_i},d=\prod p_i^{min(\alpha_i,\beta_i)}\Rightarrow sgcd|gcd\)
• \(sgcd=gcd/min\{x:x|gcd\}\)
• \(O(n\sqrt{gcd})\)？
• 这个n看起来丢不掉，能不能去sqrt？
• 为啥题目一直用\(a_1\)？
• 我会了！

从小到大用\(a_1\)质因子试除gcd，注意long long。

#include<bits/stdc++.h>
#define maxn 50
using namespace std;
typedef long long LL;
#define G c=getchar()
inline LL read()
{
	LL x=0,f=1;char G;
	while(c>57||c<48){if(c=='-')f=-1;G;}
	while(c>47&&c<58){x=x*10+c-48;G;}
	return x*f;
}
LL p[maxn],tot;
LL x,y,w;
LL gcd(LL x,LL y){return y?gcd(y,x%y):x;}
int main()
{
	int n=read();
	y=x=read();
	for(LL i=2;i*i<=x;i++)
		if(x%i==0)
			for(p[tot++]=i;x%i==0;)x/=i;
	if(x-1)p[tot++]=x;
	printf("%lld ",tot?y/p[0]:-1);
	for(int i=2,j;i<=n;i++)
	{
		x=read();w=gcd(x,y);
		int f=0;
		for(j=0;j<tot;j++)
			if(w%p[j]==0){f=1;break;}
		printf("%lld ",f?w/p[j]:-1);
	}
	return 0;
}

不定方程

*P4167 [Violet]樱花（签到题）

求不定方程

\[\frac{1}{x} + \frac{1}{y} = \frac{1}{n!} \]

的正整数解(x,y)的数目。

对于 \(100\%\) 的数据，保证\(1 \le n\le 10^6\)。

变形\((x-n!)(y-n!)=(n!)^2\)

由一一对应关系可以知道解的数目等于\((n!)^2\)的约数个数

\(n!\)好大，但是明显质因子不超过n，用约数个数公式即可

#include<bits/stdc++.h>
#define maxn 1000010
#define mod 1000000007
using namespace std;
int prime[maxn],vis[maxn],tot,n;
long long ans=1,tmp,x,y;
int main()
{
	for(int i=2;i<=maxn;i++)
	{
		if(!vis[i])prime[tot++]=i;
		for(int j=0;j<tot&&prime[j]*i<maxn;j++)
		{
			vis[prime[j]*i]=1;
			if(i%prime[j]==0)break;
		}
	}
	scanf("%d",&n);
	for(int i=0;prime[i]<=n;i++)
	{
		tmp=0;x=n;y=prime[i];
		for(;x;x/=prime[i])
			tmp+=x/y;
		ans=(ans*(tmp<<1|1)%mod)%mod;
	}
	printf("%d",ans);
}

*P2508 [HAOI2008]圆上的整点（勾股方程结论）

求一个给定的圆(\(x^2+y^2=r^2\))，在圆周上有多少个点的坐标是整数。

对于 \(100\%\) 的数据，\(n<=2000 000 000\)

勾股方程的所有解：

\(x=d\frac{v^2-u^2}{2},y=duv,r=d\frac{v^2+u^2}2\)

所以我们只要轻易地枚举\(2r\)的因子\(d\)，对于每个\(d\)我们用\(O(\sqrt{\frac{r}{d}})\)的时间枚举\(u\)，代入\(r\)的计算式得出\(v^2\)计算\(v^2\)是否为完全平方数及\(u\)与\(v\)是否互质

这样可以枚举出一个象限内的整点个数 然后输出\((ans+1)?4\)即可

复杂度O(轻松跑得过)

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
LL r, ans;
LL gcd(LL x, LL y) { return y ? gcd(y, x % y) : x; }
bool check(LL u, LL V)
{
    int v = LL(sqrt(V));
    if (V == v * v)
        return gcd(u, v) == 1;
    return 0;
}
LL calc(LL x)
{
    LL tmp = 0;
    for (LL u = 1; u * u * 2 < x; u++)
        tmp += check(u, x - u * u);
    return tmp;
}
int main()
{
    scanf("%lld", &r);
    for (LL d = 1; d * d <= 2 * r; d++)
        if (2 * r % d == 0)
            ans += calc(2 * r / d) + (d * d == 2 * r ? 0 : calc(d));
    printf("%lld\n", ans * 4 + 4);
    return 0;
}

同余

*P1082 同余方程（扩展欧几里得）

求关于\(x\)的同余方程 \(a x \equiv 1 \pmod {b}\) 的最小正整数解。

对于 \(100\%\) 的数据，\(2 ≤a, b≤ 2,000,000,000\)。

#include<bits/stdc++.h>
void gcd(int a,int b,int& d,int& x,int &y)
{
	if (b)gcd(b,a%b,d,y,x),y-=x*(a/b);
	else d=a,x=1,y=0;
}
int main()
{
	int a,b,x,y,d;
	scanf("%d %d",&a,&b);
	gcd(a,b,d,x,y);
	int aa=b/d;
	printf("%d",(x%aa+aa)%aa);
	return 0;
}

*P2312 解方程（大整数取模）

已知多项式方程：

\[a_0+a_1x+a_2x^2+\cdots+a_nx^n=0 \]

求这个方程在 \([1,m]\) 内的整数解（\(n\) 和 \(m\) 均为正整数）。

对于 \(100\%\) 的数据：\(0<n\le 100,|a_i|\le 10^{10000},a_n≠0,m<10^6\) 。

#include <bits/stdc++.h>
using namespace std;
#define maxn 1000010
typedef long long LL;
LL a[110][4], mod[4] = {10007, 11003, 12007, 13001};
int b[maxn], n, m;
string qwq;
bool check(LL x)
{
    for (int k = 0; k < 4; k++)
    {
        LL tmp = a[n][k];
        for (int i = n - 1; i >= 0; i--)
            tmp = (tmp * x + a[i][k]) % mod[k];
        if (tmp)
        {
            for (int i = x; i <= m; i += mod[k])
                b[i] = 1;
            return 0;
        }
    }
    return 1;
}
int main()
{
    scanf("%d%d", &n, &m);
    for (int i = 0; i <= n; i++)
    {
        cin >> qwq;
        for (int k = 0; k < 4; k++)
        {
            LL flg = 1, tmp = 0;
            for (int j = 0; j < qwq.length(); j++)
                if (qwq[j] == '-')
                    flg = -1;
                else
                    tmp = (tmp * 10 + qwq[j] - 48) % mod[k];
            a[i][k] = tmp * flg;
        }
    }
    int ans = 0;
    for (int x = 1; x <= m; x++)
        if (!b[x] && check(x))
            ans++;
    printf("%d\n", ans);
    for (int i = 1; i <= m; i++)
        if (!b[i])
            printf("%d\n", i);
}

*Vijos 1915 解方程加强版（根号技巧）

对于 \(100\%\) 的数据，\(0 < n ≤ 100， ∣a_i∣ ≤ 10^{10000} ，a_n ≠ 0，m ≤ 100000000 = 10^8\)。

\(O(n^3+n\sqrt m)\)的算法，我们选取素数\(p_0,p_1\)，满足\(p_0^2\geq m,p_1^2>m\)，然后在\(\pmod {p_i}\)意义下算出来\(0<=x<p_i\)的所有根，这一部分的时间复杂度是\(O(n\sqrt{m})\)的。

可以说明，找到的两组根都不超过\(n\)个，我们记为\(t_1,t_2,...,t_u\)和\(s_1,s_2,...,s_v\)。

那么我们断言原问题的可行根\(x\)一定可以通过某一个\(t_i\)和某一个\(s_j\)组合得到，再利用\(O(n)\)的方法判定是不是根。

这一部分的时间复杂度是\(O(n^3)\)

#include <bits/stdc++.h>
#define maxn 110
#define pri 4
using namespace std;
typedef long long LL;
int prime[4] = {10007, 10009, 12007, 13001}, a[110][4];
int s[110][2], la[2], anss[110], ans, n, m;
string qwq;
void gcd(LL a, LL b, LL &d, LL &x, LL &y) //ax+by=d
{
    if (!b)
        d = a, x = 1, y = 0;
    else
        gcd(b, a % b, d, y, x), y -= x * (a / b);
}
LL CRT(LL a1, LL p1, LL a2, LL p2)
{
    LL s = a2 - a1, k1, k2, pp = p1 * p2, d;
    gcd(p1, -p2, d, k1, k2);
    s = (s % pp + pp) % pp;
    k1 = (k1 % pp + pp) % pp;
    k1 = k1 * s % pp;
    return (k1 * p1 % pp + a1 % pp) % pp;
}
bool check(LL x)
{
    if (x > m)
        return 0;
    for (int k = 0; k < 4; k++)
    {
        LL tmp = a[n][k];
        for (int i = n - 1; i >= 0; i--)
            tmp = (tmp * x + a[i][k]) % prime[k];
        if (tmp)
            return 0;
    }
    return 1;
}
int main()
{
    scanf("%d%d\n", &n, &m);
    for (int i = 0; i <= n; i++)
    {
        cin >> qwq;
        for (int k = 0; k < pri; k++)
        {
            LL flg = 1, tmp = 0;
            for (int j = 0; j < qwq.length(); j++)
                if (qwq[j] == '-')
                    flg = -1;
                else
                    tmp = (tmp * 10 + qwq[j] - 48) % prime[k];
            a[i][k] = tmp * flg;
        }
    }
    for (int A = 0; A < 2; A++)
        for (int i = 1; i <= prime[A]; i++)
        {
            LL tmp = a[n][A];
            for (int j = n - 1; j >= 0; j--)
                tmp = (tmp * i + a[j][A]) % prime[A];
            if (!tmp)
                s[++la[A]][A] = i;
        }
    int ans = 0;
    for (int I = 1; I <= la[0]; I++)
        for (int J = 1, w; J <= la[1]; J++)
            if (check(w = CRT(s[I][0], prime[0], s[J][1], prime[1])))
                anss[++ans] = w;
    if (!ans)
        puts("0");
    else
    {
        sort(anss + 1, anss + 1 + ans);
        printf("%d\n", ans);
        for (int i = 1; i <= ans; i++)
            printf("%d\n", anss[i]);
    }
    return 0;
}

数论杂题

*P3951 小凯的疑惑（小学奥数）

小凯手中有两种面值的金币\(a\) 和 \(b\)，两种面值均为正整数且彼此互素。每种金币小凯都有 无数个。在不找零的情况下，仅凭这两种金币，有些物品他是无法准确支付的。现在小 凯想知道在无法准确支付的物品中，最贵的价值是多少金币？注意：输入数据保证存在 小凯无法准确支付的商品。

对于$ 100%$的数据：$1 \le a,b \le 10^9 $。

#include <bits/stdc++.h>
using namespace std;
int main()
{
	long long a,b;
	cin>>a>>b;
	cout<<a*b-a-b;
}

*P5020 货币系统（猜结论）

在网友的国度中共有 \(n\) 种不同面额的货币，第 \(i\) 种货币的面额为 \(a[i]\)，你可以假设每一种货币都有无穷多张。为了方便，我们把货币种数为 \(n\)、面额数组为 \(a[1..n]\) 的货币系统记作 \((n,a)\)。

在一个完善的货币系统中，每一个非负整数的金额 \(x\) 都应该可以被表示出，即对每一个非负整数 \(x\)，都存在 \(n\) 个非负整数 \(t[i]\) 满足 \(a[i] \times t[i]\) 的和为 \(x\)。然而， 在网友的国度中，货币系统可能是不完善的，即可能存在金额 \(x\) 不能被该货币系统表示出。例如在货币系统 \(n=3\), \(a=[2,5,9]\) 中，金额 \(1,3\) 就无法被表示出来。

两个货币系统 \((n,a)\) 和 \((m,b)\) 是等价的，当且仅当对于任意非负整数 \(x\)，它要么均可以被两个货币系统表出，要么不能被其中任何一个表出。

现在网友们打算简化一下货币系统。他们希望找到一个货币系统 \((m,b)\)，满足 \((m,b)\) 与原来的货币系统 \((n,a)\) 等价，且 \(m\) 尽可能的小。他们希望你来协助完成这个艰巨的任务：找到最小的 \(m\)。

对于 \(100\%\) 的数据，满足 \(1 ≤ T ≤ 20, 1\leq n\leq100,1\leq a[i] \leq25000\)。

#include <bits/stdc++.h>
using namespace std;
#define maxv 25010
#define maxn 110
int T, n, a[maxn];
bitset<maxv> v;
int main()
{
    for (scanf("%d", &T); T; T--)
    {
        scanf("%d", &n);
        for (int i = 1; i <= n; i++)
            scanf("%d", a + i);
        sort(a + 1, a + 1 + n);
        v.reset();
        v[0] = 1;
        int tmp = 0;
        for (int i = 1; i <= n; i++)
            if (!v[a[i]])
            {
                tmp++;
                for (int w = a[i]; w <= a[n]; w <<= 1)
                    v |= v << w;
            }
        printf("%d\n", tmp);
    }
}

*POI2003 Sums (Stage III - day 1)（环上背包）

题意：给定\(n\)个数\(a_i\)，\(T\)次询问，每次问一个数\(b_i\)能否表示成这\(n\)个数之和，一数可以任意次使用。

对于$ 100%$的数据：\(1 \le n \le 5000, a_i\leq50000,T\leq10000,b_i\leq10^9\)。

\(mn\)表示\(n\)个数中最小的数，那么只要考虑\(\%mn\)意义下每个被表示数中最小的数，用f[i]表示。
是个环状背包，所以需要在环上跑两遍。

方程：f[k]=min(f[k],f[k-a[i]]+a[i])

边界：f[0]=0;f[1..mn-1]=inf

检验：f[x%mn]<=x?等价于可以拼出。

#include<bits/stdc++.h>
#define maxn 50010
using namespace std;
typedef long long LL;
#define G c=getchar()
inline int read()
{
	int x=0,f=1;char G;
	while(c>57||c<48){if(c=='-')f=-1;G;}
	while(c>47&&c<58)x=x*10+c-48,G;
	return x*f;
}
int m,i,w,j,st,p,q,n,f[maxn],a[maxn/10],mn=maxn;
bool vis[maxn],cnt;
int main()
{
	n=read();
	for(i=1;i<=n;i++)a[i]=read(),mn=min(mn,a[i]);
	for(i=1;i<mn;i++)f[i]=2000000000;
	for(i=1,w,j;i<=n;i++)
		for(w=a[i]%mn,j=1;j<=2;j++,cnt^=1)
			for(st=0;vis[st]==cnt&&st<mn;st++)
				for(p=st,q=st+w>=mn?st+w-mn:st+w;vis[p]==cnt;vis[p]=cnt^1,p=q,q=q+w>=mn?q+w-mn:q+w)
					if(f[p]+a[i]<f[q])f[q]=f[p]+a[i];
	m=read();
	for(;m--;)
	{
		i=read();
		if(i>=f[i%mn])puts("TAK");
		else puts("NIE");
	}
	return 0;
}

可以当单源最短路做：建出带权有向图，点的标号从0到a[0]-1，i号点向(i+a[j])%a[0]号点连边，边权为a[j]。

一遍Dijkstra求出单源最短路即可完成预处理。

#include<bits/stdc++.h>
using namespace std;
#define maxn 50010
#define G c=getchar()
inline int read()
{
	int x=0,f=1;char G;
	while(c>57||c<48){if(c=='-')f=-1;G;}
	while(c>47&&c<58)x=x*10+c-48,G;
	return x*f;
}
int a[maxn],f[maxn],mn;
typedef pair<int, int> P;
priority_queue<P, vector<P>, greater<P>> Q;
int main()
{
    int n = read();
    for (int i = 0;i<n;i++)
        a[i] = read();
    mn = a[0];
    for (int i = 1; i < mn; i++)
        f[i] = 200000000;
    Q.push(P(0, 0));
    while(!Q.empty())
    {
        P t = Q.top();
        Q.pop();
        if(f[t.second]<t.first)
            continue;
        for (int u = t.second, i = 1; i < n;i++)
            if(f[u]+a[i]<f[(u+a[i])%mn])
                Q.push(P(f[(u + a[i]) % mn] = f[u] + a[i], (u + a[i]) % mn));
    }
    for (int x, m = read(); m--; puts(f[x % mn] <= x ? "TAK" : "NIE"))
        x = read();
}

组合计数

*P2822 组合数问题（签到题）

组合数 \(C_n^m\) 表示的是从 \(n\) 个物品中选出 \(m\) 个物品的方案数。举个例子，从 \((1,2,3)\) 三个物品中选择两个物品可以有 \((1,2),(1,3),(2,3)\) 这三种选择方法。根据组合数的定义，我们可以给出计算组合数 \(C_n^m\) 的一般公式：

\[C_n^m=\frac{n!}{m!(n-m)!} \]

其中 \(n!=1\times2\times\cdots\times n\)；特别地，定义 \(0!=1\)。

小葱想知道如果给定 \(n,m\) 和 \(k\)，对于所有的 \(0\leq i\leq n,0\leq j\leq \min \left ( i, m \right )\) 有多少对 \((i,j)\) 满足 \(C_i^j\) 是 \(k\) 的倍数。

对于 \(100\%\) 的数据，\(n\leq2000,m\leq2000,k\leq21,t\leq10^4\)

离线+前缀和

#include<bits/stdc++.h>
#define maxn 2010
using namespace std;
#define G c=getchar()
inline int read()
{
    int x=0,f=1;char G;
    for(;'0'>c||c>'9';G)if(c=='-')f=-1;
    for(;'0'<=c&&c<='9';G)x=10*x+c-48;
    return x*f;
}
int s[maxn][maxn],c[maxn][maxn],D=-1,d[2],g[2],K,T;
int main()
{
    T=read(),K=read();
    for(int i=0;i<=2000;i++)
        s[i][0]=c[i][0]=1;
    for (int i = 1; i <= 2000; i++)
        for (int j = 1; j <= i; j++)
            c[i][j] = (c[i - 1][j - 1] + c[i - 1][j]) % K;
    for(int i=1;i<=2000;i++)
        for(int j=1;j<=2000;j++)
            s[i][j]=s[i-1][j]+s[i][j-1]-s[i-1][j-1]+(c[i][j]==0&&i>=j);
    for(int n,m;T--;)
        n=read(),m=read(),printf("%d\n",s[n][m]);
    return 0;
}

*P4071 [SDOI2016]排列计数（简单组合计数）

求有多少种长度为 n 的序列 A，满足以下条件：

1 ~ n 这 n 个数在序列中各出现了一次

若第 i 个数 A[i] 的值为 i，则称 i 是稳定的。序列恰好有 m 个数是稳定的

满足条件的序列可能很多，序列数对 \(10^9+7\) 取模。

对于 \(100\%\) 的数据，$ T = 500000 \(，\) n \leq 1000000 \(，\) m \leq 1000000 $。

枚举稳定的数，剩余的一定是错排

\[ans=\binom{n}{m}d(n-m) \]

#include <bits/stdc++.h>
#define mod 1000000007
#define maxn 1000010
#define LL long long
using namespace std;
LL fac[maxn], f[maxn];
int t, n, m;
LL Pow(LL a, LL b)
{
    LL ret = 1;
    for (; b; b >>= 1, a = a * a % mod)
        if (b & 1)
            ret = ret * a % mod;
    return ret;
}
int main()
{
    fac[0] = 1;
    for (int i = 1; i <= 1000000; i++)
        fac[i] = fac[i - 1] * i % mod;
    f[0] = 1, f[1] = 0;
    for (int i = 2; i <= 1000000; i++)
        f[i] = (f[i - 1] + f[i - 2]) % mod * (i - 1) % mod;
    for (scanf("%d", &t); t--;)
    {
        scanf("%d%d", &n, &m);
        printf("%lld\n", fac[n] * Pow(fac[m], mod - 2) % mod * Pow(fac[n - m], mod - 2) % mod * f[n - m] % mod);
    }
    return 0;
}

*[HackerRank] Kundu and Tree（快乐容斥）

树包含N个点和N-1条边。树的边有2中颜色红色('r')和黑色（'b')。给出这N-1条边的颜色，求有多少节点的三元组(a,b,c)满足：节点a到节点b、节点b到节点c、节点c到节点a的路径上，每条路径都至少有一条边是红色的。

注意(a,b,c), (b,a,c)以及所有其他排列被认为是相同的三元组。输出结果对1000000007取余的结果。
\(1\leq N\leq50000\)

快乐容斥：既然求至少有一条红色的边，那么找到完全没有红色的边减一减就行了。

#include<bits/stdc++.h>
#define maxn 50010
#define maxm 100010
#define mod 1000000007
using namespace std;
typedef long long LL;
#define G c=getchar()
inline int read()
{
	int x=0,f=1;char G;
	while(c>57||c<48){if(c=='-')f=-1;G;}
	while(c>47&&c<58)x=x*10+c-48,G;
	return x*f;
}
#define AE(u,v) to[Si]=v,nxt[Si]=idx[u],idx[u]=Si++
int idx[maxn],nxt[maxm],to[maxm],Si;
int vis[maxn];LL n,x=0;
void dfs(int u)
{
	x++;
	vis[u]=1;
	for(int i=idx[u];i+1;i=nxt[i])if(!vis[to[i]])dfs(to[i]);
}
int main()
{
	memset(idx,-1,sizeof(idx));
	n=read();
	for(int i=1;i<n;i++)
	{
		int u,v;char c;
		scanf("%d %d %c",&u,&v,&c);
		if(c=='b')AE(u,v),AE(v,u);
	}
	LL sum=0;
	for(int i=1;i<=n;i++)
	{
		x=0;
		if(!vis[i])dfs(i);
		sum+=x*(x-1)/2*(n-x)+x*(x-1)*(x-2)/6;
	}
	sum=n*(n-1)*(n-2)/6-sum;sum=sum%mod+mod;
	printf("%lld",sum%mod);
	return 0;
}

*51nod P1486 大大走格子（容斥计数）

有一个\(h\)行\(w\)列的棋盘，里面有\(n\)格子是不能走的，现在要求从左上角走到右下角的方案数。（只能向右和向下走）

\(1 ≤ h, w ≤ 10^5, 1 ≤ n ≤ 2000\)

dp[i]表示从第i个点出发走到右下角的方案数

先考虑一个点：地球人都会的样子

两个点：画图

多个点：推广

#include<bits/stdc++.h>
#define maxn 200010
#define mod 1000000007
using namespace std;
typedef long long LL;
#define G c=getchar()
inline int read()
{
	int x=0,f=1;char G;
	while(c>57||c<48){if(c=='-')f=-1;G;}
	while(c>47&&c<58)x=x*10+c-48,G;
	return x*f;
}
LL inv[maxn],fac[maxn],dp[2010];
LL Pow(LL x,LL y)
{
	LL r=1;
	for(;y;y>>=1,x=x*x%mod)if(y&1)r=r*x%mod;
	return r;
}
int C(int n,int m)//C(n,m)
{
	return fac[n+m]*inv[n]%mod*inv[m]%mod;
}
struct poi
{
	int x,y;
}a[2010];
inline bool cmp(const poi &a,const poi &b){return a.x==b.x?a.y>b.y:a.x>b.x;}
int main()
{
	fac[0]=1;inv[0]=Pow(1,mod-2);
	for(int i=1;i<maxn;i++)fac[i]=fac[i-1]*i%mod,inv[i]=Pow(fac[i],mod-2);
	int n=read(),m=read(),k=read();
	for(int i=1;i<=k;i++)
		a[i].x=read(),a[i].y=read();
	a[0].x=a[0].y=1;
	sort(a,a+1+k,cmp);
	for(int i=0;i<=k;i++)
	{
		dp[i]=C(n-a[i].x,m-a[i].y);
		for(int j=0;j<i;j++)
			if(a[i].x<=a[j].x&&a[i].y<=a[j].y)
			{
				(dp[i]-=dp[j]*C(a[j].x-a[i].x,a[j].y-a[i].y)%mod)%=mod;
				if(dp[i]<0)dp[i]+=mod;
			}
	}
	printf("%lld",dp[k]);
	return 0;
}

P5023 填数游戏（还不会，没看题解_(:з」∠)_）

小 D 特别喜欢玩游戏。这一天，他在玩一款填数游戏。

这个填数游戏的棋盘是一个\(n \times m\)的矩形表格。玩家需要在表格的每个格子中填入 一个数字（数字 \(0\)或者数字 \(1\)），填数时需要满足一些限制。

下面我们来具体描述这些限制。

为了方便描述，我们先给出一些定义：

• 我们用每个格子的行列坐标来表示一个格子，即（行坐标，列坐标）。（注意： 行列坐标均从 \(0\) 开始编号）

• 合法路径 \(P\)：一条路径是合法的当且仅当：

  1. 这条路径从矩形表格的左上角的格子\((0,0)\)出发，到矩形的右下角格子 \((n - 1,m - 1)\)结束；
  2. 在这条路径中，每次只能从当前的格子移动到右边与它相邻的格子，或者 从当前格子移动到下面与它相邻的格子。
     例如：在下面这个矩形中，只有两条路径是合法的，它们分别是\(P_1\)： \((0,0)\) → \((0,1)\) → \((1,1)\)和\(P_2\)：\((0,0)\) → \((1,0)\) → \((1,1)\)。

对于一条合法的路径 \(P\)，我们可以用一个字符串\(w(P)\)来表示，该字符串的长度为\(n + m - 2\)，其中只包含字符“\(R\)”或者字符“\(D\)”， 第 \(i\) 个字符记录了路径$ P$ 中第 \(i\) 步的移动 方法，“ \(R\)”表示移动到当前格子右边与它相邻的格子，“ \(D\)”表示移动到当前格子下面 与它相邻的格子。例如，上图中对于路径\(P_1\)，有\(w(P_1) =\) "\(RD\)"；而对于另一条路径\(P_2\)， 有\(w(P_2) =\) "\(DR\)"。

同时，将每条合法路径 \(P\) 经过的每个格子上填入的数字依次连接后，会得到一个长 度为\(n + m - 1\)的 \(01\) 字符串，记为 \(s(P)\)。例如，如果我们在格子\((0,0)\)和\((1,0)\)上填入数字 \(0\)，在格子\((0,1)\)和\((1,1)\)上填入数字 1（见上图红色数字）。那么对于路径\(P_1\)，我们可以得 到\(s(P_1) =\) "\(011\)",对于路径\(P_2\)，有\(s(P_2) =\) "\(001\)"。

游戏要求小 D 找到一种填数字$ 0\(、\) 1$ 的方法，使得对于两条路径\(P_1\)， \(P_2\)，如果\(w(P_1) > w(P_2)\)，那么必须\(s(P_1) ≤ s(P_2)\)。我们说字符串 \(a\) 比字符串 \(b\) 小，当且仅当字符串 \(a\) 的字典序小于字符串 \(b\) 的字典序，字典序的定义详见第一题。但是仅仅是找一种方法无法满 足小 D 的好奇心，小 D 更想知道这个游戏有多少种玩法，也就是说，有多少种填数字 的方法满足游戏的要求？

小 D 能力有限，希望你帮助他解决这个问题，即有多少种填 \(0\)、\(1\) 的方法能满足题目要求。由于答案可能很大，你需要输出答案对\(10^9 + 7\)取模的结果。

对于 \(100\%\) 的数据，满足\(n\leq8,m\leq1000000\)

概率与期望

*P1850 换教室 （线性期望DP）

对于刚上大学的牛牛来说，他面临的第一个问题是如何根据实际情况申请合适的课程。

在可以选择的课程中，有 \(2n\) 节课程安排在 \(n\) 个时间段上。在第 \(i\)（\(1 \leq i \leq n\)）个时间段上，两节内容相同的课程同时在不同的地点进行，其中，牛牛预先被安排在教室 \(c_i\) 上课，而另一节课程在教室 \(d_i\) 进行。

在不提交任何申请的情况下，学生们需要按时间段的顺序依次完成所有的 \(n\) 节安排好的课程。如果学生想更换第 \(i\) 节课程的教室，则需要提出申请。若申请通过，学生就可以在第 \(i\) 个时间段去教室 \(d_i\) 上课，否则仍然在教室 \(c_i\) 上课。

由于更换教室的需求太多，申请不一定能获得通过。通过计算，牛牛发现申请更换第 \(i\) 节课程的教室时，申请被通过的概率是一个已知的实数 \(k_i\)，并且对于不同课程的申请，被通过的概率是互相独立的。

学校规定，所有的申请只能在学期开始前一次性提交，并且每个人只能选择至多 \(m\) 节课程进行申请。这意味着牛牛必须一次性决定是否申请更换每节课的教室，而不能根据某些课程的申请结果来决定其他课程是否申请；牛牛可以申请自己最希望更换教室的 \(m\) 门课程，也可以不用完这 \(m\) 个申请的机会，甚至可以一门课程都不申请。

因为不同的课程可能会被安排在不同的教室进行，所以牛牛需要利用课间时间从一间教室赶到另一间教室。

牛牛所在的大学有 \(v\) 个教室，有 \(e\) 条道路。每条道路连接两间教室，并且是可以双向通行的。由于道路的长度和拥堵程度不同，通过不同的道路耗费的体力可能会有所不同。 当第 \(i\)（\(1 \leq i \leq n-1\)）节课结束后，牛牛就会从这节课的教室出发，选择一条耗费体力最少的路径前往下一节课的教室。

现在牛牛想知道，申请哪几门课程可以使他因在教室间移动耗费的体力值的总和的期望值最小，请你帮他求出这个最小值。

对于 \(100\%\) 的数据保证 \(1 \leq n \leq 2000\)，\(0 \leq m \leq 2000\)，\(1 \leq v \leq 300\)，\(0 \leq e \leq 90000\)。

首先看起来是个很套路的线性模型：有个固定n阶段的流程，还有m次技能可以使用，差不多就是分层图的套路。

f[i][j]表示前i节课换了j次课的最小期望，发现不好转移，因为不知道上个阶段这个人出现在哪儿的概率，多一维f[i][j][0/1]来表示这个阶段有没有换课，0就代表他固定在c[i]教室，1代表有k[i]的概率在d[i]，1-k[i]的概率在c[i]

就能转移辣。

#include<bits/stdc++.h>
#define maxn 2010
#define INF 1000000000
using namespace std;
int c[maxn],d[maxn];
double k[maxn],f[maxn][maxn][2];
int dis[maxn][maxn];
int n,m,v,e;
int main()
{
	scanf("%d%d%d%d",&n,&m,&v,&e);
	for(int i=1;i<=n;i++)scanf("%d",&c[i]);
	for(int i=1;i<=n;i++)scanf("%d",&d[i]);
	for(int i=1;i<=n;i++)scanf("%lf",&k[i]);
	for(int i=1;i<=v;i++)for(int j=1;j<=v;j++)if(i!=j)dis[i][j]=dis[j][i]=INF;
	for(int u,v,w;e--;)
	{
		scanf("%d%d%d",&u,&v,&w);
		dis[u][v]=dis[v][u]=min(dis[u][v],w);
	}
	for(int kk=1;kk<=v;kk++)
		for(int i=1;i<=v;i++)
			for(int j=1;j<=v;j++)
				if(dis[i][j]>dis[i][kk]+dis[kk][j])
					dis[i][j]=dis[i][kk]+dis[kk][j];
	for(int i=0;i<=n;i++)for(int j=0;j<=m;j++)f[i][j][0]=f[i][j][1]=1e23;
	f[1][0][0]=f[1][1][1]=0;
	for(int i=2;i<=n;i++)
		for(int j=0;j<=m;j++)
		{
			f[i][j][0]=min(
				f[i-1][j][0]+dis[c[i-1]][c[i]],
				f[i-1][j][1]+dis[c[i-1]][c[i]]*(1-k[i-1])+dis[d[i-1]][c[i]]*k[i-1]);
			if(j>0)f[i][j][1]=min(
				f[i-1][j-1][0]+dis[c[i-1]][c[i]]*(1-k[i])+dis[c[i-1]][d[i]]*k[i],
				f[i-1][j-1][1]+dis[c[i-1]][c[i]]*(1-k[i-1])*(1-k[i])+dis[c[i-1]][d[i]]*(1-k[i-1])*k[i]+dis[d[i-1]][c[i]]*k[i-1]*(1-k[i])+dis[d[i-1]][d[i]]*k[i-1]*k[i]);
		}
	double ans=f[n][m][0];
	for(int i=0;i<=m;i++)
		ans=min(ans,min(f[n][i][0],f[n][i][1]));
	printf("%.2lf",ans);
	return 0;
}

*P4316 绿豆蛙的归宿 （DAG上期望DP）

给出一个有向无环图，起点为1终点为N，每条边都有一个长度，并且从起点出发能够到达所有的点，所有的点也都能够到达终点。绿豆蛙从起点出发，走向终点。
到达每一个顶点时，如果有K条离开该点的道路，绿豆蛙可以选择任意一条道路离开该点，并且走向每条路的概率为 1/K 。
现在绿豆蛙想知道，从起点走到终点的所经过的路径总长度期望是多少？

对于 \(100\%\) 的数据，满足 \(N<=100000，M<=2*N\)

顺：f[i]从1走到i的期望距离

#include <bits/stdc++.h>
using namespace std;
#define add(u, v, w) W[++Si] = w, to[Si] = v, nxt[Si] = idx[u], idx[u] = Si
#define maxn 200010
int to[maxn], idx[maxn], nxt[maxn], W[maxn], Si, du[maxn], in[maxn];
double f[maxn], p[maxn];
queue<int> Q;
int main()
{
    int n, m;
    scanf("%d%d", &n, &m);
    for (int i = 1; i <= m; i++)
    {
        int x, y, z;
        scanf("%d%d%d", &x, &y, &z);
        add(x, y, z);
        du[x]++;
        in[y]++;
    }
    Q.push(1);
    f[1] = 0, p[1] = 1;
    for (; !Q.empty();)
    {
        int u = Q.front();
        Q.pop();
        f[u] /= p[u];
        for (int i = idx[u]; i; i = nxt[i])
        {
            int v = to[i];
            p[v] += p[u] / du[u];
            f[v] += (f[u] + W[i]) * p[u] / du[u];
            in[v]--;
            if (in[v] == 0)
                Q.push(v);
        }
    }
    printf("%.2lf", f[n]);
}

逆：f[i]表示从i走到n的期望距离。u->v时f[u]=sigma(f[v])/du[u]

#include <bits/stdc++.h>
using namespace std;
#define add(u, v, w) W[++Si] = w, to[Si] = v, nxt[Si] = idx[u], idx[u] = Si
#define maxn 200010
int to[maxn], idx[maxn], nxt[maxn], W[maxn], Si, du[maxn], in[maxn];
double f[maxn];
queue<int> Q;
int main()
{
    int n, m;
    scanf("%d%d", &n, &m);
    for (int i = 1; i <= m; i++)
    {
        int x, y, z;
        scanf("%d%d%d", &x, &y, &z);
        add(y, x, z);
        du[x]++;
        in[x]++;
    }
    Q.push(n);
    f[n] = 0;
    for (; !Q.empty();)
    {
        int u = Q.front();
        Q.pop();
        for (int i = idx[u]; i; i = nxt[i])
        {
            int v = to[i];
            f[v] += (f[u] + W[i])/ du[v];
            in[v]--;
            if (in[v] == 0)
                Q.push(v);
        }
    }
    printf("%.2lf", f[1]);
}

*P1365 WJMZBMR打osu! / Easy （线性期望dp）

某一天WJMZBMR在打osu~~~但是他太弱逼了，有些地方完全靠运气:(

我们来简化一下这个游戏的规则

有\(n\)次点击要做，成功了就是o，失败了就是x，分数是按combo计算的，连续\(a\)个combo就有\(a\times a\)分，combo就是极大的连续o。

比如ooxxxxooooxxx，分数就是\(2 \times 2 + 4 \times 4 = 4 +16=20\)。

Sevenkplus闲的慌就看他打了一盘，有些地方跟运气无关要么是o要么是x，有些地方o或者x各有50%的可能性，用?号来表示。

比如oo?xx就是一个可能的输入。
那么WJMZBMR这场osu的期望得分是多少呢？

比如oo?xx的话，?是o的话就是oooxx => 9，是x的话就是ooxxx => 4

期望自然就是\((4+9)/2 =6.5\)了

对于 \(100\%\) 的数据，满足\(n\leq 10^6\)

期望连续长度\(\bar x\)，考虑贡献（增量）

x：\(x^2-x^2=0,\bar x\rightarrow0\)

o：\((x+1)^2-x^2=2x+1,\bar x\rightarrow\bar x+1\)

?：\(\frac{x^2+(x+1)^2}2-x^2=\bar x+0.5,\bar x\rightarrow\frac{\bar x+1}2\)

转移即可

#include <bits/stdc++.h>
#define maxn 1000010
using namespace std;
double p1[maxn], f[maxn];
char qwq[maxn];
int main()
{
    int n;
    scanf("%d", &n);
    cin >> (qwq + 1);
    for (int i = 1; i <= n; i++)
        if (qwq[i] == '?')
        {
            p1[i] = p1[i - 1] / 2 + 0.5;
            f[i] = f[i - 1] + p1[i - 1] + 0.5;
        }
        else if (qwq[i] == 'o')
        {
            p1[i] = p1[i - 1] + 1;
            f[i] = f[i - 1] + 2*p1[i - 1] + 1;
        }
        else
        {
            p1[i] = 0;
            f[i] = f[i - 1];
        }
    printf("%.4lf", f[n]);
}

*P3802 小魔女帕琪 （概率）

现在帕琪与强大的夜之女王，吸血鬼蕾咪相遇了，夜之女王蕾咪具有非常强大的生命力，普通的魔法难以造成效果，只有终极魔法：帕琪七重奏才能对蕾咪造成伤害。帕琪七重奏的触发条件是：连续释放的7个魔法中，如果魔法的属性各不相同，就能触发一次帕琪七重奏。

现在帕琪有7种属性的能量晶体，分别为\(a1,a2,a3,a4,a5,a6,a7\)(均为自然数)，每次释放魔法时，会随机消耗一个现有的能量晶体，然后释放一个对应属性的魔法。

现在帕琪想知道，她释放出帕琪七重奏的期望次数是多少，可是她并不会算，于是找到了学OI的你

对于100%的测试点，\(a1+a2+a3+a4+a5+a6+a7<=10^9\)

设\(n=sigma(ai)\)

对于每一位\(i\)，从\(i\)开始的\(7\)位中，\(7\)个魔法都不相同的概率为\(7! * a1/n * a2/(n-1) * a3/(n-2) * a4/(n-3) *a5/(n-4) * a6/(n-5) * a7/(n-6)\)

由于一共有\(n-6\)个这样的位置\(i\)，所以最后一项除以\(n-6\)可以消去

答案就是\(7! * a1/n * a2/(n-1) * a3/(n-2) * a4/(n-3) *a5/(n-4) * a6/(n-5) * a7\)

\[13!\frac{4}{52}*\frac{4}{51}*\frac{4}{50}*\frac{4}{49}*...*\frac{4}{40}*40 \]