头像

Anoxia_3

长春大学


访客:10719

离线:5小时前



Anoxia_3
6小时前

不定方程解的个数、容斥原理、状压、组合计数


题意:需要从N个盒子中总共取出M枝花,每个盒子中的花是相同的,问有多少种取法。


假设每个盒子中的花有无限,记从第i个盒子中取出x[i]枝花,则有:
$x[1] + x[2] + … + x[N] = M$,令$y[i] = x[i] + 1$,

=>$[1] + y[2] + .. + y[N] = M + N$

=>方案数:$C_{N+M-1}^{N-1}$


回到原问题,直接从正面求不好求,那就用补集的思想:合法方案数=总方案数-非法方案数!

总方案数:$C_{N+M-1}^{N-1}$;

非法方案:至少存在从一个从盒子中取花的数量大于盒子中存在的花的数量的取法
记第i个盒子非法的集合为S[i],则所有非法方案数为$|S[1]∩S[2]∩…∩S[N]|$

合法方案数= $C_{N+M-1}^{N-1} + (-1)^{k-1}\sum_{1<=i<=N}{C_{N+M-2-A[i]}^{N-1}} + (-1)^{k-1} \sum_{1<=i<j<=N}{C_{N+M-3-A[i]-A[j]}^{N-1}} + … +$ ,
(其中k是集合个数)


如何方便地去计算非法方案?
容斥原理中项的个数是$2^N$个,对应这N个条件的选或者不选,因此可以把每一种情况看成一个二进制表达式,
1代表选,0代表不选;1的个数是奇数为-,1的个数为偶数为+。


因为N很小,所以C可以直接用定义计算,因为模的是质数,所以除的数可以用快速幂算出逆元转换成乘


#include <iostream>

using namespace std;

typedef long long LL;

const int N = 20 , mod = 1e9 + 7;

LL A[N];
int down = 1;

int qmi(int a, int k)
{
    int res = 1;
    while (k)
    {
        if (k & 1) res = (LL)res * a % mod;
        a = (LL)a * a % mod;
        k >>= 1;
    }
    return res;
}

int C(LL a, LL b)
{
    if (a < b) return 0;
    int up = 1;
    for (LL i = a; i > a - b; i -- ) up = i % mod * up % mod;

    return (LL)up * down % mod; // 费马小定理
}



int main()
{
    LL n , m;
    cin >> n >> m;
    for(int i = 0 ; i < n ; i++)    cin >> A[i];

    for(int i = 1 ; i <= n - 1 ; i++) down = (LL)down * i % mod;
    down = qmi(down , mod - 2);

    int res = 0;
    for(int i = 0 ; i < 1 << n ; i++)
    {
        LL a = m + n - 1 , b = n - 1;
        int sign = 1;
        for(int j = 0 ; j < n ; j++)
        {
            if(i >> j & 1)
            {
                sign *= -1;
                a -= 1 + A[j];
            }
        }

        res = (res + C(a , b) * sign) % mod;
    }

    cout << (res + mod) % mod << endl;
    return 0;
}


活动打卡代码 AcWing 214. Devu和鲜花

Anoxia_3
6小时前

不定方程解的个数、容斥原理、状压、组合计数


题意:需要从N个盒子中总共取出M枝花,每个盒子中的花是相同的,问有多少种取法。


假设每个盒子中的花有无限,记从第i个盒子中取出x[i]枝花,则有:
$x[1] + x[2] + … + x[N] = M$,令$y[i] = x[i] + 1$,

=>$[1] + y[2] + .. + y[N] = M + N$

=>方案数:$C_{N+M-1}^{N-1}$


回到原问题,直接从正面求不好求,那就用补集的思想:合法方案数=总方案数-非法方案数!

总方案数:$C_{N+M-1}^{N-1}$;

非法方案:至少存在从一个从盒子中取花的数量大于盒子中存在的花的数量的取法
记第i个盒子非法的集合为S[i],则所有非法方案数为$|S[1]∩S[2]∩…∩S[N]|$

合法方案数= $C_{N+M-1}^{N-1} + (-1)^{k-1}\sum_{1<=i<=N}{C_{N+M-2-A[i]}^{N-1}} + (-1)^{k-1} \sum_{1<=i<j<=N}{C_{N+M-3-A[i]-A[j]}^{N-1}} + … +$ ,
(其中k是集合个数)


如何方便地去计算非法方案?
容斥原理中项的个数是$2^N$个,对应这N个条件的选或者不选,因此可以把每一种情况看成一个二进制表达式,
1代表选,0代表不选;1的个数是奇数为-,1的个数为偶数为+。


因为N很小,所以C可以直接用定义计算,因为模的是质数,所以除的数可以用快速幂算出逆元转换成乘


#include <iostream>

using namespace std;

typedef long long LL;

const int N = 20 , mod = 1e9 + 7;

LL A[N];
int down = 1;

int qmi(int a, int k)
{
    int res = 1;
    while (k)
    {
        if (k & 1) res = (LL)res * a % mod;
        a = (LL)a * a % mod;
        k >>= 1;
    }
    return res;
}

int C(LL a, LL b)
{
    if (a < b) return 0;
    int up = 1;
    for (LL i = a; i > a - b; i -- ) up = i % mod * up % mod;

    return (LL)up * down % mod; // 费马小定理
}



int main()
{
    LL n , m;
    cin >> n >> m;
    for(int i = 0 ; i < n ; i++)    cin >> A[i];

    for(int i = 1 ; i <= n - 1 ; i++) down = (LL)down * i % mod;
    down = qmi(down , mod - 2);

    int res = 0;
    for(int i = 0 ; i < 1 << n ; i++)
    {
        LL a = m + n - 1 , b = n - 1;
        int sign = 1;
        for(int j = 0 ; j < n ; j++)
        {
            if(i >> j & 1)
            {
                sign *= -1;
                a -= 1 + A[j];
            }
        }

        res = (res + C(a , b) * sign) % mod;
    }

    cout << (res + mod) % mod << endl;
    return 0;
}


活动打卡代码 AcWing 208. 开关问题

Anoxia_3
23小时前
#include <iostream>
#include <cstring>

using namespace std;

const int N = 35;

int a[N][N];
int n;

int gauss()
{
    int r = 1 , c = 1;
    for(; c <= n ; c++)
    {
        int t = r;
        while(!a[t][c] && t <= n) t++;

        if(!a[t][c]) continue;

        for(int i = c ; i <= n + 1 ; i++) swap(a[t][i] , a[r][i]);

        for(int i = r + 1 ; i <= n ; i++)
            for(int j = n + 1 ; j >= c ; j--)
                a[i][j] ^= a[i][c] & a[r][j];

        r++;
    }

    int res = 1;
    if(r < n + 1)//如果r == n+1,说明每一个行式子的第c列都存在系数为1的变量,变量只有唯一解
    {
        for(int i = r ; i <= n ; i++)
        {
            if(a[i][n + 1]) return -1;  //出现 0 == !0
            res *= 2;
        }
    }

    return res;
}

int main()
{
    int T;
    cin >> T;

    while(T--)
    {
        memset(a , 0 , sizeof a);
        cin >> n;
        for(int i = 1 ; i <= n ; i++)   cin >> a[i][n + 1];//读入初始状态
        for(int i = 1 ; i <= n ; i++)   //读入最终状态
        {
            int t;
            cin >> t;
            a[i][n + 1] ^= t;   //关注的并不是开关状态,而是是否需要变化
            a[i][i] = 1;    //i开关当然会控制它自己
        }

        int x , y;
        while(cin >> x >> y , x || y) a[y][x] = 1;//y受x的控制

        int t = gauss();

        if(t == -1) puts("Oh,it's impossible~!!");
        else cout << t << endl;
    }
    return 0;
}



$n$维空间求圆心:
记圆心为$(a[0][1],a[0][2],a[0][3],…,a[0][n])$,给定点为$(a[i][1],a[i][2],a[i][3],…a[i][n])$,则可以列出n+1条等式:

形如:$(a[i][1]-a[0][1])^2 + (a[i][2]-a[0][2])^2 + (a[i][3]-a[0][3])^2 + … + (a[i][n]-a[0][n])^2 = R^2- - i$

将第式i-式1,可以消去未知数的平方项,可以得到n条线性方程:
形如:$2(a[i][1] - a[1][1])x[1] + 2(a[i][2] - a[1][2])x[2] +…+ 2(a[i][j] - a[1][j])x[i] +…+2(a[i][n] - a[1][n])x[n]
$

$=(a[i][1]^2 - a[1][1]^2) + (a[i][2]^2 - a[1][2]^2) +…+ (a[i][j]^2 - a[1][j]^2) +..+ (a[i][n]^2 - a[1][n]^2)$

令$b[i][j] = a[i][j] - a[1][j]$ (j∈[1,n]),

$b[i][n+1] = (a[i][1]^2 - a[1][1]^2) + (a[i][2]^2 - a[1][2]^2) +…+ (a[i][j]^2 - a[1][j]^2) +..+ (a[i][n]^2 - a[1][n]^2)$

显然b都是已知数,原始转换成形如:$b[i][1] * x[1] + b[i][2] * x[2] +…+b[i][j] * x[j] +…+ b[i][n] * x[j] = b[i][n+1]$

此时有n条线性方程,n个未知数,可以用高斯消元求得解。

#include <iostream>
#include <cstdio>
#include <cmath>

using namespace std;

const int N = 15;

int n;
double a[N][N] , b[N][N];

void gauss()
{
    for(int r = 1 , c = 1 ; c <= n ; c++)
    {
        int t = r;
        for(int i = r + 1 ; i <= n ; i++)  
            if(fabs(b[t][c]) < fabs(b[i][c])) t = i;

        for(int i = c ; i <= n + 1 ; i++) swap(b[t][i] , b[r][i]);
        for(int i = n + 1 ; i >= c ; i--) b[r][i] /= b[r][c];

        for(int i = r + 1 ; i <= n ; i++)
            for(int j = c + 1 ; j <= n + 1 ; j++)
                b[i][j] -= b[i][c] * b[r][j];

        r++;
    }

    for(int i = n ; i ; i--)//求第i个未知数的解b[i][n + 1]
        for(int j = i + 1  ; j <= n ;j++)//第j个未知数的解b[j][n + 1]
            b[i][n + 1] -= b[i][j] * b[j][n + 1];
}

int main()
{
    cin >> n;

    for(int i = 1 ; i <= n + 1 ; i++)
        for(int j = 1 ; j <= n ; j++)   
            cin >> a[i][j];

    for(int i = 1 ; i <= n ; i++)
        for(int j = 1 ; j <= n ; j++)
        {
            b[i][j] = 2 * (a[i + 1][j] - a[1][j]);
            b[i][n + 1] += a[i + 1][j] * a[i + 1][j] - a[1][j] * a[1][j];
        }   

    gauss();

    for(int i = 1 ; i <= n ; i++)   printf("%.3lf " , b[i][n + 1]);
    return 0;
}



/*
n维空间求圆心:
    记圆心为(a[0][1],a[0][2],a[0][3],...,a[0][n]),给定点为(a[i][1],a[i][2],a[i][3],...a[i][n]),则可以列出n+1条等式:
    (a[1][1]-a[0][1])^2 + (a[1][2]-a[0][2])^2 + (a[1][3]-a[0][3])^2 + ... + (a[1][n]-a[0][n])^2 = R^2---1
    .
    (a[i][1]-a[0][1])^2 + (a[i][2]-a[0][2])^2 + (a[i][3]-a[0][3])^2 + ... + (a[i][n]-a[0][n])^2 = R^2---i
    .
    (a[n+1][1]-a[0][1])^2 + (a[n+1][2]-a[0][2])^2 + (a[n+1][3]-a[0][3])^2 + ... + (a[n+1][n]-a[0][n])^2 = R^2---n+1

将第式i-式1,可以消去未知数的平方项,可以得到n条线性方程:
形如:2(a[i][1] - a[1][1])x[1] + 2(a[i][2] - a[1][2])x[2] +...+ 2(a[i][j] - a[1][j])x[i] +...+ 2(a[i][n] - a[1][n])x[n] 
     = (a[i][1]^2 - a[1][1]^2) + (a[i][2]^2 - a[1][2]^2) +...+ (a[i][j]^2 - a[1][j]^2) +..+ (a[i][n]^2 - a[1][n]^2)

令b[i][j] = a[i][j] - a[1][j] ,
b[i][n+1] = (a[i][1]^2 - a[1][1]^2) + (a[i][2]^2 - a[1][2]^2) +...+ (a[i][j]^2 - a[1][j]^2) +..+ (a[i][n]^2 - a[1][n]^2)
显然b都是已知数,原始转换成形如:b[i][1]*x[1] + b[i][2]*x[2] +...+b[i][j]*x[j] +...+ b[i][n]*x[j] = b[i][n+1]

此时有n条线性方程,n个未知数,可以用高斯消元求得解。

*/
#include <iostream>
#include <cstdio>
#include <cmath>

using namespace std;

const int N = 15;

int n;
double a[N][N] , b[N][N];

void gauss()
{
    for(int r = 1 , c = 1 ; c <= n ; c++)
    {
        int t = r;
        for(int i = r + 1 ; i <= n ; i++)  
            if(fabs(b[t][c]) < fabs(b[i][c])) t = i;

        for(int i = c ; i <= n + 1 ; i++) swap(b[t][i] , b[r][i]);
        for(int i = n + 1 ; i >= c ; i--) b[r][i] /= b[r][c];

        for(int i = r + 1 ; i <= n ; i++)
            for(int j = c + 1 ; j <= n + 1 ; j++)
                b[i][j] -= b[i][c] * b[r][j];

        r++;
    }

    for(int i = n ; i ; i--)//求第i个未知数的解b[i][n + 1]
        for(int j = i + 1  ; j <= n ;j++)//第j个未知数的解b[j][n + 1]
            b[i][n + 1] -= b[i][j] * b[j][n + 1];
}

int main()
{
    cin >> n;

    for(int i = 1 ; i <= n + 1 ; i++)
        for(int j = 1 ; j <= n ; j++)   
            cin >> a[i][j];

    for(int i = 1 ; i <= n ; i++)
        for(int j = 1 ; j <= n ; j++)
        {
            b[i][j] = 2 * (a[i + 1][j] - a[1][j]);
            b[i][n + 1] += a[i + 1][j] * a[i + 1][j] - a[1][j] * a[1][j];
        }   

    gauss();

    for(int i = 1 ; i <= n ; i++)   printf("%.3lf " , b[i][n + 1]);
    return 0;
}



看到样例中n=3,答案=5,就有可能是卡特兰数。

做法:依次将1~2n内的每个数归入奇数项和偶数项,显然1要排在奇数项的第一位(整个序列的第一位)。

本题结论:奇数项数>=偶数项数

证明:
$\ \ \ $如果奇数项数 < 偶数项数,那么当将一个数归入奇数项时,因为是顺次枚举每个数,那么这个数将比上一个选出归入偶数项的数要大,即如果把这个数放入奇数项,肯定会比之后所有的偶数项都要大,显然会比相邻的偶数项要大,不满足条件,假设不成立。

#include <iostream>

using namespace std;

const int N = 2000010;

typedef long long LL;

int n , mod;
int primes[N] , cnt;
bool st[N];

int qmi(int a, int b)
{
    int res = 1;
    while(b)
    {
        if(b & 1) res = (LL) res * a % mod;
        b >>= 1;
        a = (LL) a * a % mod;
    }
    return res;
}

void init(int n)
{
    for(int i = 2 ; i <= n ; i++)
    {
        if(!st[i]) primes[cnt++] = i;

        for(int j = 0 ; primes[j] <= n / i ; j++)
        {
            st[primes[j] * i] = true;

            if(i % primes[j] == 0) break;
        }
    }
}

int get(int x , int p)
{
    int s = 0;
    while(x) s += x / p , x /= p;
    return s;
}

int C(int a , int b)
{
    int res = 1;
    for(int i = 0 ; i < cnt ; i++)
    {
        int p = primes[i];
        int s = get(a , p) - get(b , p) - get(a - b , p);
        res = (LL) res * qmi(p , s) % mod;
    }
    return res;
}

int main()
{

    cin >> n >> mod;

    init(2 * n);//因为式子中会用到2*n

    cout << (C(2 * n , n) - C(2 * n , n - 1) + mod) % mod << endl;
    return 0;
}


活动打卡代码 AcWing 1316. 有趣的数列

/*
看到样例中n=3,答案=5,就有可能是卡特兰数。
在1~2n中依次选出每个数是偶数项还是奇数项,显然1要是奇数项的第一位。
结论:奇数项数>=偶数项数
    证明:
        如果奇数项数 < 偶数项数,那么当将一个数归入奇数项时,因为是顺次枚举每个数,那么这个数将比上一个选出归入偶数项的数要大,
    即如果把这个数放入奇数项,肯定会比之后所有的偶数项都要大,显然会比相邻的偶数项要大,不满足条件,假设不成立。
*/
#include <iostream>

using namespace std;

const int N = 2000010;

typedef long long LL;

int n , mod;
int primes[N] , cnt;
bool st[N];

int qmi(int a, int b)
{
    int res = 1;
    while(b)
    {
        if(b & 1) res = (LL) res * a % mod;
        b >>= 1;
        a = (LL) a * a % mod;
    }
    return res;
}

void init(int n)
{
    for(int i = 2 ; i <= n ; i++)
    {
        if(!st[i]) primes[cnt++] = i;

        for(int j = 0 ; primes[j] <= n / i ; j++)
        {
            st[primes[j] * i] = true;

            if(i % primes[j] == 0) break;
        }
    }
}

int get(int x , int p)
{
    int s = 0;
    while(x) s += x / p , x /= p;
    return s;
}

int C(int a , int b)
{
    int res = 1;
    for(int i = 0 ; i < cnt ; i++)
    {
        int p = primes[i];
        int s = get(a , p) - get(b , p) - get(a - b , p);
        res = (LL) res * qmi(p , s) % mod;
    }
    return res;
}

int main()
{

    cin >> n >> mod;

    init(2 * n);//因为式子中会用到2*n

    cout << (C(2 * n , n) - C(2 * n , n - 1) + mod) % mod << endl;
    return 0;
}



卡特兰数的应用

依据卡塔兰数的推导,找出$(n,m)$关于$y=x+1$的对称点$(m-1,n+1)$,每一条非法路径都对应一条$(0,0)$到$(m-1,n+1)$的路径,因此合法路径=总路径数-非法路径,即$C_{n+m}^{n} - C_{n+m}^{m-1}$

#include <iostream>

using namespace std;

const int N = 10010;

int primes[N] , cnt;
bool st[N];
int a[N] , b[N];

void init(int n)
{
    for(int i = 2 ; i <= n ; i++)
    {
        if(!st[i]) primes[cnt++] = i;

        for(int j = 0 ; primes[j] <= n / i ; j++)
        {
            st[primes[j] * i] = true;

            if(i % primes[j] == 0) break;
        }
    }
}

int get(int n , int p)//n!中p的次数
{
    int res = 0;
    for(int i = n ; i ; i /= p) res += i / p;
    return res;
}

void mul(int a[] , int &len , int b)//a = a * b 高精度*低精度
{
    int t = 0;
    for(int i = 0  ; i < len ; i++)
    {
        t += a[i] * b;
        a[i] = t % 10;
        t /= 10;
    }

    while(t)
    {
        a[len++] = t % 10;
        t /= 10;
    }

}

int C(int a , int b , int c[])//把$C_{a}^{b}$存到c数组中
{
    int len = 1;
    c[0] = 1;

    for(int i = 0 ; i < cnt ; i++)
    {
        int p = primes[i];
        int s = get(a , p) - get(b , p) - get(a - b , p);
        while(s--) mul(c , len , p);
    }

    return len;
}

void sub(int a[] , int al , int b[] , int bl)
{
    for(int i = 0 , t = 0; i < al ; i++)
    {
        a[i] -= t + b[i];
        if(a[i] < 0) t = 1 , a[i] += 10;
        else t = 0;//t代表前一位是否要减一
    }
}

int main()
{
    init(N - 1);

    int n , m;
    cin >> n >> m;
    int al = C(n + m , m , a);
    int bl = C(n + m , m - 1 , b);

    sub(a , al , b , bl);

    int i = al - 1;
    while(!a[i] && i) i--;
    while(i >= 0) cout << a[i--];
    return 0;
}


活动打卡代码 AcWing 1315. 网格

/*
依据卡塔兰数的推导,找出(n,m)关于y=x+1的对称点(m-1,n+1),每一条非法路径都对应一条(0,0)到(m-1,n+1)的路径,
因此合法路径=总路径数-非法路径,即$C_{n+m}^{n} - C_{n+m}^{m-1}$
*/

#include <iostream>

using namespace std;

const int N = 10010;

int primes[N] , cnt;
bool st[N];
int a[N] , b[N];

void init(int n)
{
    for(int i = 2 ; i <= n ; i++)
    {
        if(!st[i]) primes[cnt++] = i;

        for(int j = 0 ; primes[j] <= n / i ; j++)
        {
            st[primes[j] * i] = true;

            if(i % primes[j] == 0) break;
        }
    }
}

int get(int n , int p)//n!中p的次数
{
    int res = 0;
    for(int i = n ; i ; i /= p) res += i / p;
    return res;
}

void mul(int a[] , int &len , int b)//a = a * b 高精度*低精度
{
    int t = 0;
    for(int i = 0  ; i < len ; i++)
    {
        t += a[i] * b;
        a[i] = t % 10;
        t /= 10;
    }

    while(t)
    {
        a[len++] = t % 10;
        t /= 10;
    }

}

int C(int a , int b , int c[])//把$C_{a}^{b}$存到c数组中
{
    int len = 1;
    c[0] = 1;

    for(int i = 0 ; i < cnt ; i++)
    {
        int p = primes[i];
        int s = get(a , p) - get(b , p) - get(a - b , p);
        while(s--) mul(c , len , p);
    }

    return len;
}

void sub(int a[] , int al , int b[] , int bl)
{
    for(int i = 0 , t = 0; i < al ; i++)
    {
        a[i] -= t + b[i];
        if(a[i] < 0) t = 1 , a[i] += 10;
        else t = 0;//t代表前一位是否要减一
    }
}

int main()
{
    init(N - 1);

    int n , m;
    cin >> n >> m;
    int al = C(n + m , m , a);
    int bl = C(n + m , m - 1 , b);

    sub(a , al , b , bl);

    int i = al - 1;
    while(!a[i] && i) i--;
    while(i >= 0) cout << a[i--];
    return 0;
}



法一、
在[L,R]中找出k个不降序的数,等价于在[0,R-L]中找出k个不降序的数,

即:$0<=a[1]<=a[2]<=a[3]<=…<=a[k]<=R-L$

=> 令$x[1]=a[1] , x[2]=a[2]-a[1] ,
x[i]=x[i]-x[i-1]$,使得$x[1]+x[2]+…+x[k]=a[k]<=R-L$,的解的组数,其中每个$x$都是非负整数

=> 令$y[i]=x[i]+1$,使得$y[1]+y[2]+…+y[k]=a[k]+k<=R-L+k$的解的组数,
其中每个$y$都是正整数,利用隔板法,因为是不等式,所以解的组数是$C^{k}_{R-L+k}$


法二、
令$b[1] = a[1] , b[2] = a[2] + 1$ , $b[3] = a[3] + 2 , … , b[i] = a[i] + i - 1$,

那么原不等式就转换成严格小于的形式:$L <= b[1] < b[2] < b[3] < … < b[k] <= R+k-1$,

解的组数即$C^{k}_{R-L+k}$ ,

枚举k,最终结果是 $ΣC_{k}^{R-L+k}$ , 化简后: $C_{R-L+n+1}^{R-L+1}+1$

#include <iostream>

using namespace std;

typedef long long LL;

const int p = 1e6 + 3;

int qmi(int a , int b)
{
    int res = 1;
    while(b)
    {
        if(b & 1) res = (LL) res * a % p;
        b >>= 1;
        a = (LL) a * a % p;
    }
    return res;
}

int C(int a, int b)
{
    if (a < b) return 0;

    int down = 1, up = 1;
    for (int i = a, j = 1; j <= b; i --, j ++ )
    {
        up = (LL)up * i % p;
        down = (LL)down * j % p;
    }

    return (LL)up * qmi(down, p - 2) % p;
}

int lucas(int a , int b)
{
    if(a < p && b < p) return C(a , b);
    return (LL)lucas(a / p , b / p) * C(a % p , b % p) % p;
}

int main()
{
    int T;
    cin >> T;


    int l , r , n;
    while(T--)
    {
        cin >> n >> l >> r;

        cout << (lucas(r - l + n + 1 , r - l + 1) + p - 1) % p << endl;
    }
    return 0;
}