P4071 排列计数

原题链接：https://www.luogu.com.cn/problem/P4071

题目描述

求有多少种 (1) 到 (n) 的排列 (a)，满足序列恰好有 (m) 个位置 (i)，使得 (a_i = i)

答案对 (10^9 + 7) 取模

输入格式

本题单测试点内有多组数据

输入的第一行是一个整数 (T)，代表测试数据的整数

以下 (T) 行，每行描述一组测试数据

对于每组测试数据，每行输入两个整数，依次代表 (n) 和 (m)

输出格式

共输出 (T) 行，对于每组测试数据，输出一行一个整数代表答案

输入输出样例

输入 #1

5
1 0
1 1
5 2
100 50
10000 5000

输出 #1

0
1
20
578028887
60695423

说明/提示

数据规模与约定

本题共 (20) 个测试点，各测试点等分，其数据规模如下表

测试点编号	(T=)	(n, m leq)	测试点编号	(T=)	(n, m leq)
(1sim 3)	(10^3)	(8)	(10 sim 12)	(10^3)	(10^3)
(4 sim 6)	(10^3)	(12)	(13 sim 14)	(5 imes 10^5)	(10^3)
(7 sim 9)	(10^3)	(100)	(15 sim 20)	(5 imes 10^5)	(10^6)

分析

首先先选出位置与数值相等的 (m) 个数，那么选择的情况数为 (n choose m)，接着要对剩余的 (n-m) 个数进行错排，我们设 (d[i]) 表示 (i) 个数错排的情况数，现在先让每个数都升序排列，接着拿出第一个数字 (a[1]) ，在后面的 (i-1) 个数字中选择一个与之交换，假设为 (a[k])，选择的情况数为 (i-1) 种，对于每一种选择的情况，(a[k]) 已经实现了错排，而交换后的 (a[1]) 有两种选择，第一种是就安放原位，那么它也实现了错排，接着只需对 (i-2) 个数进行错排即可，第二种情况是它不放在这，那么它还要参与接下来的错排操作，所以要对 (i-1) 个数进行错排，根据加法原理，两种情况的方案数为 (d[i-1]+d[i-2])，又根据乘法原理，(d[i]=(n-1)(d[i-1]+d[i-2]))

接下来考虑时间复杂度，题目中 (N, M le 10^6)，这样组合数就不能通过递推得到，而应该用组合数公式计算：

[{n choose m} = frac{n!}{m!(n-m)!} \ frac{n!}{m!(n-m)!} equiv n!(m!)^{-1}(n-m)!^{-1} pmod P ]

所以我们要提前计算 (i) 的阶乘对 (P) 取模的结果以及其乘法逆元

由于模数 (P=10^9 + 7) 是质数，所以可以根据费马小定理来求乘法逆元：

[a^{P-1} equiv a cdot a^{P-2} equiv 1 pmod P \ a^{-1} equiv a^{P-2} pmod P ]

这样预处理阶乘和递推错排情况数的时间复杂度均为 (O(N))，计算每个阶乘的乘法逆元要使用快速幂，时间复杂度为 (O(N log_2 N))

代码部分

代码中使用 fac[i] 表示 (i! mod P)，inv[i] 表示 ((i!)^{-1} mod P)，d[i] 表示 (i) 个数的错排情况数

#include <cstdio>
using namespace std;

const int N=1e6+10,P=1e9+7;
int t,n,m;
long long ans,fac[N],inv[N],d[N];

long long power_mod(int x,int y)
{
    if(y==0) return 1;
    long long res=power_mod(x,y>>1);
    res=(res*res)%P;
    if(y&1) res=(res*x)%P;
    return res;
}

int main()
{
    fac[0]=1;
    for(int i=1;i<=1e6;i++) fac[i]=fac[i-1]*i%P;
    d[1]=0,d[2]=1;
    for(int i=3;i<=1e6;i++) d[i]=(i-1)*(d[i-1]+d[i-2])%P;
    scanf("%d",&t);
    while(t--)
    {
        scanf("%d%d",&n,&m);
        if(n==m) // special judge
        {
            puts("1");
            continue;
        }
        if(inv[m]==0) inv[m]=power_mod(fac[m],P-2);
        if(inv[n-m]==0) inv[n-m]=power_mod(fac[n-m],P-2);
        ans=fac[n]*inv[m]%P*inv[n-m]%P*d[n-m]%P;
        printf("%lld
",ans);
    }
    return 0;
}