数论~~~

质数

质数
- Miller-Rabin算法
- 质因子分解
- 质数筛
- - 埃氏筛
  - 欧拉筛
  - 如果只是计数，埃氏筛改进
快速幂
- 乘法快速幂
- 矩阵快速幂
- 1维k阶实战(提醒：最好在mul函数中作乘法时加上（long long）的强制类型转换，或者全部数组换成long long（当没有除余的时候），不然容易越界，而且，即便常规做法不会越界，该做法也可能会越界)
- - 爬楼梯
  - 第n个泰波那契数
  - 多米诺和托米诺平铺
- k维1阶
- - 统计元音字母序列的数目
  - 学生出勤记录II
逆元和除法同余
- 连续数字逆元的线性递推
- 连续阶乘逆元的线性递推
容斥原理
- Coprime Subsequences
- 洛谷1450-硬币购物
- - 题目描述
  - 输入格式
  - 输出格式
  - 样例 #1
  - - 样例输入 #1
    - 样例输出 #1
  - 提示
  - - 数据规模与约定
- 播放列表的数量
高斯消元
- 加法方程组
- - 题目练手
  - - 洛谷4035
    - 洛谷-5027-三角形称重
- 异或方程组
- - 杭电5833-完全平方数
  - 洛谷-2962-令点全为1
  - 外星千足虫-洛谷-2447
- 同余方程组
- - hd-5755
  - poj-2947

质数

Miller-Rabin算法

在这里插入图片描述

#include <bits/stdc++.h>
using namespace std;
typedef __int128 ll;//128位整数
//__128类型要自己写输入输出
//输入
template <typename T> inline T read() {
    T x = 0, f = 1; char ch = 0;
    for (; !isdigit(ch); ch = getchar()) if (ch == '-')f = -1;
    for (; isdigit(ch); ch = getchar())x = (x << 3) + (x << 1) + (ch - '0');
    return x * f;
}
//输出,在本题用不上
template <typename T> inline void write(T x) {
    if (x < 0)putchar('-'), x = -x;
    if (x > 9)write(x / 10);
    putchar(x % 10 + '0');
}
template <typename T> inline void print(T x, char ed = '\n') {
    write(x), putchar(ed);
}
ll t, n;
//快速幂
ll qpow(ll a, ll b, ll mod) {
    ll ret = 1;
    while (b) {
        if (b & 1)ret = (ret * a) % mod;
        a = (a * a) % mod;
        b >>= 1;
    }
    return ret % mod;
}
//miller_rabin
vector<ll> p = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
bool miller_rabin(ll n) {
    if (n < 3 || n % 2 == 0)return n == 2;
    ll u = n - 1, t = 0;
    while (u % 2 == 0)u /= 2, ++t;
    for (auto a : p) {
        if (n == a)return 1;
        if (n % a == 0)return 0;
        ll v = qpow(a, u, n);
        if (v == 1)continue;
        ll s = 1;
        for (; s <= t; ++s) {
            if (v == n - 1)break;
            v = v * v % n;
        }
        if (s > t)return 0;
    }
    return 1;
}
int main() {
    t = read<ll>();
    while (t--) {
        n = read<ll>();
        if (miller_rabin(n))puts("Yes");
        else puts("No");
    }
    return 0;
}

质因子分解

https://leetcode.cn/problems/largest-component-size-by-common-factor/description/
题解
用质因子分解+并查集

class Solution {
public:
    int mem[100001];//记录质因子是否出现，出现在哪
    int fa[20002];//父节点
    int se[20002];//大小
    int ans=1;//记录最大连通块大小
    int find(int x){
        if(fa[x]!=x)fa[x]=find(fa[x]);
        return fa[x];
    }
    void uni(int x,int y){
        int fx=find(x),fy=find(y);
        if(fx!=fy){
            fa[fy]=fx;
            se[fx]+=se[fy];
            ans=max(ans,se[fx]);
        }
    }
    //分解质因子
    void f(int x,int j){
        for(int i=2;i*i<=x;i++){
            if(x%i==0){
                if(mem[i]==-1)mem[i]=j;
                else uni(mem[i],j);
                while(x%i==0)x/=i;
            }
        }
        if(x>1){
            if(mem[x]==-1)mem[x]=j;
            else uni(mem[x],j);
        }
    }
    int largestComponentSize(vector<int>& nums) {
        memset(mem,-1,sizeof mem);
        int n=nums.size();
        for(int i=0;i<n;i++)fa[i]=i,se[i]=1;
        for(int i=0;i<n;i++){
            f(nums[i],i);
        }
        return ans;
    }
};

质数筛

定义：给定整数n,返回1~n范围的所有质数

埃氏筛

从2开始遍历，将2的倍数全记录为合数；再到3,将3的倍数全记录为合数（从33开始，因为32设置过了），如此往下，若是质数，效仿上述过程。直到i*i<=n;

https://leetcode.cn/problems/count-primes/

class Solution {
public:
    bool visited[5000001];
    int countPrimes(int n) {
        n-=1;
        int ans=0;
        for(int i=2;i*i<=n;i++){
            if(!visited[i]){
                for(int j=i*i;j<=n;j+=i){
                    visited[j]=1;
                }
            }
        }
        for(int i=2;i<=n;i++){
            if(!visited[i])ans++;
        }
        return ans;
    }
};

欧拉筛

埃氏筛会重复筛，比如26,34都筛了一次12
欧拉筛争取不重复，即每个合数只被自己的最小质因子筛掉。

class Solution {
public:
    bool visited[5000001];
    int countPrimes(int n) {
        int prime[n/2+1];
        int cnt=0;
        for(int i=2;i<n;i++){
            if(!visited[i])prime[cnt++]=i;
            for(int j=0;j<cnt;j++){
                if(i*prime[j]>=n)break;
                visited[i*prime[j]]=1;
                if(i%prime[j]==0)break;
            }
        }
        return cnt;
    }
};

如果只是计数，埃氏筛改进

首先假设所有的奇数为素数，然后删去非素数的奇数

class Solution {
public:
    bool visited[5000001];
    int countPrimes(int n) {
        n-=1;
        if(n<=1)return 0;
        int ans=(n+1)/2;
        for(int i=3;i*i<=n;i+=2){
            if(!visited[i]){
                for(int j=i*i;j<=n;j+=2*i){
                    if(!visited[j]){
                        visited[j]=1;
                        ans--;
                    }
                }
            }
        }
        return ans;
    }
};

快速幂

在这里插入图片描述

1维2阶：斐波那契数列
2维1阶：例如：dp[i][1]=dp[i-1][1]+dp[i-1][2]

乘法快速幂

洛谷1226

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
int qpow(ll a,ll b,ll p){
    ll ans=1;
    while(b){
        if(b&1)ans=(ans*a)%p;
        a=a*a%p;
        b>>=1;
    }
    return ans;
}
int main(){
    ll a,b,p;
    cin>>a>>b>>p;
    printf("%lld^%lld mod %lld=%d",a,b,p,qpow(a,b,p));
    return 0;
}

矩阵快速幂

斐波那契问题举例：


class Solution {
public:
    vector<vector<int>> mul(vector<vector<int>> a, vector<vector<int>> b) {
        int n = a.size(), m = b[0].size(), k = a[0].size();
        vector<vector<int>> ans(n, vector<int>(m));
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                ans[i][j] = 0;
                for (int c = 0; c < k; c++) {
                    ans[i][j] += a[i][c] * b[c][j];
                }
            }
        }
        return ans;
    }
    vector<vector<int>> qpow(vector<vector<int>> a, int b) {
        int n = a.size();
        vector<vector<int>> ans(n, vector<int>(n));
        for(int i=0;i<n;i++)ans[i][i]=1;
        while (b) {
            if (b & 1)ans = mul(ans, a);
            a = mul(a, a);
            b >>= 1;
        }
        return ans;
    }
    int fib(int n) {
        if(n==0){
            return 0;
        }
        vector<vector<int>>start = {{1, 0}};
        vector<vector<int>> base = {{1, 1}, {1, 0}};
        vector<vector<int>> ans = mul(start, qpow(base, n - 1));
        return ans[0][0];
    }
};

1维k阶实战(提醒：最好在mul函数中作乘法时加上（long long）的强制类型转换，或者全部数组换成long long（当没有除余的时候），不然容易越界，而且，即便常规做法不会越界，该做法也可能会越界)

设k=3;

求start：{F(0),F(1),F(2)};
求base: 使得{F(i-1),F(i-2),F(i-3)} * base == {F(i),F(i-1),F(i-2)}

爬楼梯

https://leetcode.cn/problems/climbing-stairs/description/
核心：
vector<vector>start = {{1, 1}};
vector<vector> base = {{1, 1}, {1, 0}};

class Solution {
public:
    vector<vector<long long int>> mul(vector<vector<long long int>> a, vector<vector<long long int>> b) {
        int n = a.size(), m = b[0].size(), k = a[0].size();
        vector<vector<long long>> ans(n, vector<long long>(m));
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                ans[i][j] = 0;
                for (int c = 0; c < k; c++) {
                    ans[i][j] += a[i][c] * b[c][j];
                }
            }
        }
        return ans;
    }
    vector<vector<long long int>> qpow(vector<vector<long long int>> a, int b) {
        int n = a.size();
        vector<vector<long long int>> ans(n, vector<long long int>(n));
        for(int i=0;i<n;i++)ans[i][i]=1;
        while (b) {
            if (b & 1)ans = mul(ans, a);
            a = mul(a, a);
            b >>= 1;
        }
        return ans;
    }
    int climbStairs(int n) {
        vector<vector<long long int>>start = {{1, 1}};
        vector<vector<long long int>> base = {{1, 1}, {1, 0}};
        vector<vector<long long int>> ans = mul(start, qpow(base, n - 1));
        return ans[0][0];
    }
};

第n个泰波那契数

https://leetcode.cn/problems/n-th-tribonacci-number/description/

class Solution {
public:
    vector<vector<long long int>> mul(vector<vector<long long int>> a, vector<vector<long long int>> b) {
        int n = a.size(), m = b[0].size(), k = a[0].size();
        vector<vector<long long>> ans(n, vector<long long>(m));
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                ans[i][j] = 0;
                for (int c = 0; c < k; c++) {
                    ans[i][j] += a[i][c] * b[c][j];
                }
            }
        }
        return ans;
    }
    vector<vector<long long int>> qpow(vector<vector<long long int>> a, int b) {
        int n = a.size();
        vector<vector<long long int>> ans(n, vector<long long int>(n));
        for(int i=0;i<n;i++)ans[i][i]=1;
        while (b) {
            if (b & 1)ans = mul(ans, a);
            a = mul(a, a);
            b >>= 1;
        }
        return ans;
    }
    int tribonacci(int n) {
        if(n<=1)return n;
        vector<vector<long long int>>start = {{1, 1,0}};
        vector<vector<long long int>> base = {{1, 1,0}, {1, 0,1},{1,0,0}};
        vector<vector<long long int>> ans = mul(start, qpow(base, n - 2));
        return ans[0][0];
    }
};

多米诺和托米诺平铺

https://leetcode.cn/problems/domino-and-tromino-tiling/description/
常规思路

class Solution {
public:
    int N=1e9+7;
    long long a[1001];//有凸出来的
    long long b[1001];//没凸出来的
    int numTilings(int n) {
        if(n<=1)return 1;
        if(n==2)return 2;
        if(n==3)return 5;
        a[1]=1,a[2]=2,b[1]=1,b[2]=2;
        for(int i=3;i<=n;i++){
            a[i]=(b[i-1]+a[i-1])%N;
            b[i]=(b[i-1]+b[i-2]+2*a[i-2])%N;
        }
        return b[n];
    }
};

通过常规思路，打表找规律发现递推式f(i)=f(i-1)*2+f(i-3),然后就可以用矩阵快速幂

class Solution {
public:
    int N=1e9+7;
    vector<vector<int>> mul(const vector<vector<int>> &a,const vector<vector<int>> &b) {
        int n = a.size(), m = b[0].size(), k = a[0].size();
        vector<vector<int>> ans(n, vector<int>(m));
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                ans[i][j] = 0;
                for (int c = 0; c < k; c++) {
                    ans[i][j] =(ans[i][j]+(long long)a[i][c] * b[c][j]%N)%N;
                }
            }
        }
        return ans;
    }
    vector<vector<int>> qpow(vector<vector<int>> a, int b) {
        int n = a.size();
        vector<vector<int>> ans(n, vector<int>(n));
        for(int i=0;i<n;i++)ans[i][i]=1;
        while (b) {
            if (b & 1)ans = mul(ans, a);
            a = mul(a, a);
            b >>= 1;
        }
        return ans;
    }
    int numTilings(int n) {
        if(n<=1)return 1;
        vector<vector<int>>start = {{2, 1,1}};
        vector<vector<int>> base = {{2, 1,0}, {0, 0,1},{1,0,0}};
        vector<vector<int>> ans = mul(start, qpow(base, n - 2));
        return ans[0][0];
    }
};

k维1阶

统计元音字母序列的数目

常规做法

class Solution {
public:
    int N=1e9+7;
    int countVowelPermutation(int n) {
        long long dp[n+1][5];
        for(int i=0;i<5;i++)dp[1][i]=1;
        for(int i=2;i<=n;i++){
            dp[i][0]=(dp[i-1][4]+dp[i-1][2]+dp[i-1][1])%N;
            dp[i][1]=(dp[i-1][0]+dp[i-1][2])%N;
            dp[i][2]=(dp[i-1][1]+dp[i-1][3])%N;
            dp[i][3]= dp[i-1][2];
            dp[i][4]=(dp[i-1][2]+dp[i-1][3])%N;
        }
        int ans=0;
        for(int i=0;i<5;i++){
            ans=(ans+dp[n][i])%N;
        }
        return ans;
    }
};

矩阵快速幂
start:{dp[1][0],dp[1][1],dp[1][2],dp[1][3],dp[1][4]}；
base:…………很显然吧

class Solution {
public:
    int N=1e9+7;
    vector<vector<long long int>> mul(vector<vector<long long int>> a, vector<vector<long long int>> b) {
        int n = a.size(), m = b[0].size(), k = a[0].size();
        vector<vector<long long>> ans(n, vector<long long>(m));
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                ans[i][j] = 0;
                for (int c = 0; c < k; c++) {
                    ans[i][j] =(ans[i][j]+a[i][c] * b[c][j]%N)%N;
                }
            }
        }
        return ans;
    }
    vector<vector<long long int>> qpow(vector<vector<long long int>> a, int b) {
        int n = a.size();
        vector<vector<long long int>> ans(n, vector<long long int>(n));
        for(int i=0;i<n;i++)ans[i][i]=1;
        while (b) {
            if (b & 1)ans = mul(ans, a);
            a = mul(a, a);
            b >>= 1;
        }
        return ans;
    }
    int countVowelPermutation(int n) {
        vector<vector<long long int>>start = {{1, 1,1,1,1}};
        vector<vector<long long int>> base = {{0, 1,0,0,0}, {1, 0,1,0,0},{1,1,0,1,1},{0,0,1,0,1},{1,0,0,0,0}};
        vector<vector<long long int>> ans = mul(start, qpow(base, n - 1));
        int cnt=0;
        for(int i=0;i<5;i++){
            cnt=(cnt+ans[0][i])%N;
        }
        return cnt;
    }
};

学生出勤记录II

https://leetcode.cn/problems/student-attendance-record-ii/description/
dp数组说明：

dp[i][k] 长度为i
k=0,不以L结尾，且目前0A
k=1,不以L结尾，且目前1A
k=2,以L结尾，且结尾连续L为1，目前0A;
k=3,以L结尾，且结尾连续L为2，目前0A;
k=4,以L结尾，且结尾连续L为1，目前1A;
k=5,以L结尾，且结尾连续L为2, 目前1A;
常规做法

class Solution {
public:
    int N=1e9+7;
    int checkRecord(int n) {
        long long int dp[n+1][6];
        dp[1][0]=1,dp[1][1]=1,dp[1][2]=1,dp[1][3]=0,dp[1][4]=0,dp[1][5]=0;
        for(int i=2;i<=n;i++){
            dp[i][0]=(dp[i-1][2]+dp[i-1][3]+dp[i-1][0])%N;
            dp[i][1]=(dp[i-1][0]+dp[i-1][1]+p[i-1][2]+dp[i-1][3]+dp[i-1][4]+dp[i-1][5])%N;
            dp[i][2]=dp[i-1][0];
            dp[i][3]=dp[i-1][2];
            dp[i][4]=dp[i-1][1];
            dp[i][5]=dp[i-1][4];
        }
        long long ans=0;
        for(int i=0;i<6;i++){
            ans=(ans+dp[n][i])%N;
        }
        return ans;
    }
};

矩阵快速幂

class Solution {
public:
    int N=1e9+7;
    vector<vector<int>> mul(const vector<vector<int>> &a,  const vector<vector<int>> &b) {
        int n = a.size(), m = b[0].size(), k = a[0].size();
        vector<vector<int>> ans(n, vector<int>(m));
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                ans[i][j] = 0;
                for (int c = 0; c < k; c++) {
                    ans[i][j] =(ans[i][j]+(long long)a[i][c] * b[c][j])%N;
                }
            }
        }
        return ans;
    }
    vector<vector<int>> qpow(vector<vector<int>> a, int b) {
        int n = a.size();
        vector<vector<int>> ans(n, vector<int>(n));
        for(int i=0;i<n;i++)ans[i][i]=1;
        while (b) {
            if (b & 1)ans = mul(ans, a);
            a = mul(a, a);
            b >>= 1;
        }
        return ans;
    }
    int checkRecord(int n) {
        vector<vector<int>>start = {{1, 1,1,0,0,0}};
        vector<vector<int>> base = {{1, 1,1,0,0,0}, {0, 1,0,0,1,0},{1,1,0,1,0,0},{1,1,0,0,0,0},{0,1,0,0,0,1},{0,1,0,0,0,0}};
        vector<vector<int>> ans = mul(start, qpow(base, n - 1));
        int cnt=0;
        for(int i=0;i<6;i++){
            cnt=(cnt+ans[0][i])%N;
        }
        return cnt;
    }
};

逆元和除法同余

在这里插入图片描述

连续数字逆元的线性递推

在这里插入图片描述
洛谷P3811

#include<bits/stdc++.h>
using namespace std;
int a[3000001];
int main(){
    ios::sync_with_stdio(false);
    cin.tie(0);
    int n,p;
    cin>>n>>p;
    a[1]=1;
    cout<<"1\n";
    for(int i=2;i<=n;i++){
        a[i]=(int)(p-(long long)a[p%i]*(p/i)%p);
        cout<<a[i]<<'\n';    
    }
    return 0;

}

连续阶乘逆元的线性递推

在这里插入图片描述

容斥原理

Coprime Subsequences

https://codeforces.com/problemset/problem/803/F

首先无关顺序，总序列数为2的n次方-1
找到最大的数m,令i从m到1遍历，得到以i为最大公约数的序列数，为2**(cnt)-1-dp[i2]-dp[i3]……，结果存在dp[i]中
最终输出dp[1]

#include<bits/stdc++.h>
using namespace std;
#define N 1000000007
int a[100002];//记录每个数的个数
int dp[100002];
//快速幂
int qpow(long long a,int b){
    long long ans=1;
    while(b){
        if(b&1)ans=ans*a%N;
        b>>=1;
        a=a*a%N;
    }
    return ans;
}
int main(){
    int n,x,m=0;
    cin>>n;
    for(int i=0;i<n;i++){
        cin>>x;
        m=max(m,x);
        a[x]++;
    }
    for(int i=m;i>=1;i--){
        long long cnt=a[i];//记录序列中i的倍数的总个数
        for(int j=i*2;j<=m;j+=i){
            cnt+=a[j];  
            dp[i]=((long long)dp[i]-dp[j]+N)%N;
        }
        dp[i]=((long long)dp[i]+qpow(2,cnt)-1+N)%N;
    }
    cout<<dp[1]<<'\n';
    return 0;
}

洛谷1450-硬币购物

题目描述

共有 $4$ 种硬币。面值分别为 $c_1,c_2,c_3,c_4$ 。

某人去商店买东西，去了 $n$ 次，对于每次购买，他带了 $d_i$ 枚 $i$ 种硬币，想购买 $s$ 的价值的东西。请问每次有多少种付款方法。

输入格式

输入的第一行是五个整数，分别代表 $c_1,c_2,c_3,c_4, n$ 。

接下来 $n$ 行，每行有五个整数，描述一次购买，分别代表 $d_1, d_2, d_3, d_4,s$ 。

输出格式

对于每次购买，输出一行一个整数代表答案。

样例 #1

样例输入 #1

1 2 5 10 2
3 2 3 1 10
1000 2 2 2 900

样例输出 #1

4
27

提示

数据规模与约定

对于 $100\%$ 的数据，保证 $\leq c_i, d_i, s \leq 10^5$ ， $\leq n \leq 1000$ 。

题解

dp[i][j]：假设钱随便用的情况下，钱的种类范围从1~i，目标价格为j ，成立的方案数
考虑一种钱脱离实际的方案数，假设目标价格为s,钱的价值为c,实际数量为d,若s-c*(d+1)>=0,
方案数为dp[s-c*(d+1)]，妙不可言！！！
这样就可以运用容斥原理，得到实际的方案数

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll s,dp[100003],c[5],d[5];
ll f(ll x){
    if(x>=0)return dp[x];
    return 0;
}
ll g(int x){
    return (d[x]+1)*c[x];
}
int main(){
    ll n;
    cin>>c[1]>>c[2]>>c[3]>>c[4]>>n;
    //空间压缩优化
    dp[0]=1;
    for(int i=1;i<=4;i++){
        for(int j=c[i];j<=100002;j++){
            dp[j]+=dp[j-c[i]];
        }
    }
    // for(int i=1;i<=4;i++){
    //     dp[i][0]=1;
    //     for(int j=1;j<=100000;j++){
    //         dp[i][j]+=dp[i-1][j];
    //         if(j-c[i]>=0)dp[i][j]+=dp[i][j-c[i]];
    //     }
    // }
    while(n--){
        cin>>d[1]>>d[2]>>d[3]>>d[4]>>s;
        ll ans=dp[s]-f(s-g(1))-f(s-g(2))-f(s-g(3))-f(s-g(4))
        +f(s-g(1)-g(2))+f(s-g(1)-g(3))+f(s-g(1)-g(4))+f(s-g(2)-g(3))+f(s-g(2)-g(4))+f(s-g(3)-g(4))
        -f(s-g(1)-g(2)-g(3))-f(s-g(1)-g(3)-g(4))-f(s-g(4)-g(2)-g(3))-f(s-g(1)-g(2)-g(4))
        +f(s-g(1)-g(2)-g(3)-g(4));
        cout<<ans<<'\n';
    }
    return 0;
}

播放列表的数量

https://leetcode.cn/problems/number-of-music-playlists/description/
先不考虑每首歌至少播一次的条件

定义f(n,l,k)：n首歌可选，目标l首歌，至少k首间隔，则
$A_n^{k + 1}{(n - k)^{l - k - 1}} = \frac{{n!{{(n - k)}^{l - k}}}}{(n - k)!}$
现在考虑每首歌至少一次的条件，根据容斥原理，答案为 $\sum\limits_{i = 0}^{n - k - 1} {{{( - 1)}^i} * C_n^i * f(n - i,l,k) =\sum\limits_{i = 0}^{n - k - 1} {{( - 1)}^i}\frac{{n! * {{(n - i - k)}^{l - k}}}}{{i! * (n - i - k)!}}}$

class Solution {
public:
    using ll=long long;
    int N=1e9+7;
    int numMusicPlaylists(int n, int l, int k) {
        ll jie[n+1],inv[n+1];
        jie[0]=1;
        for(int i=1;i<=n;i++)jie[i]=(jie[i-1]*i)%N;
        inv[n]=qpow(jie[n],N-2);
        for(int i=n-1;i>=0;i--)inv[i]=(i+1)*inv[i+1]%N;
        long long ans=0;
        int flag=1;
        for(int i=0;i<n-k;i++){
            ans=(ans+flag*(qpow(n-i-k,l-k)*jie[n]%N*(inv[n-i-k]*inv[i]%N)%N)+N)%N;
            flag=-flag;
        }
        return ans;
    }
    ll qpow(ll a,int b){
        ll ans=1;
        while(b){
            if(b&1)ans=ans*a%N;
            b>>=1;
            a=a*a%N;
        }
        return ans;
    }
};

还可以用dp做
在这里插入图片描述

class Solution {
public:
    long long dp[101][101];
    int N=1e9+7;
    int numMusicPlaylists(int n, int l, int k) {
        dp[0][0]=1;
        for(int i=1;i<=l;i++){
            for(int j=1;j<=min(i,n);j++){
                dp[i][j]+=dp[i-1][j-1]*(n-j+1);
                dp[i][j]+=dp[i-1][j]*max(0,j-k);
                dp[i][j]%=N;
            }
        }
        return dp[l][n];
    }
    
};

高斯消元

在这里插入图片描述
最终可解方程性质特点:自由元所在的行的系数都为0，
主元所在列除了自己全为0
主元的确定有时会依赖于自由元
矛盾：自由元的结果列不为0；
唯一解：无自由元
多解：存在自由元，且对应的结果列为0

加法方程组

模版题洛谷-2455

#include<bits/stdc++.h>
using namespace std;
double sml=1e-7;
double mat[51][52];
void swap1(int x,int y,int n){
    for(int i=1;i<=n+1;i++){
        double tmp=mat[x][i];
        mat[x][i]=mat[y][i];
        mat[y][i]=tmp;
    }
}
void gauss(int n){
    for(int i=1;i<=n;i++){
        int max_row=i;
        for(int j=1;j<=n;j++){
            if(j<i && abs(mat[j][j])>=sml)continue;
            if(abs(mat[j][i])>abs(mat[max_row][i]))max_row=j;
        }
        swap1(i,max_row,n);
        if(abs(mat[i][i])>=sml){
            double tmp=mat[i][i];
            for(int j=i;j<=n+1;j++){
                mat[i][j]/=tmp;
            }
            for(int j=1;j<=n;j++){
                if(i!=j){
                    double rate=mat[j][i];
                    for(int k=i;k<=n+1;k++){
                        mat[j][k]-=mat[i][k]*rate;
                    }
                }
            }
        }
    }
}
int main(){
    int n;
    cin>>n;
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n+1;j++){
            cin>>mat[i][j];
        }
    }
    gauss(n);
    int flag=1;
    for(int i=1;i<=n;i++){
        if(abs(mat[i][i])<sml && abs(mat[i][n+1])>=sml){
            flag=-1;
            break;
        }
        if(abs(mat[i][i])<sml)flag=0;
    }
    if(flag==1){
        for(int i=1;i<=n;i++){
            printf("x%d=%.2lf\n",i,mat[i][n+1]);
        }
    }else{
        cout<<flag<<'\n';
    }
    return 0;
}

题目练手

洛谷4035

通过n个“距离球心相等”的方程组得到球心坐标

#include<bits/stdc++.h>
using namespace std;
double sml=1e-7;
double mat[51][52];
double a[12][12];
void swap1(int x,int y,int n){
    for(int i=1;i<=n+1;i++){
        double tmp=mat[x][i];
        mat[x][i]=mat[y][i];
        mat[y][i]=tmp;
    }
}
void gauss(int n){
    for(int i=1;i<=n;i++){
        int max_row=i;
        for(int j=1;j<=n;j++){
            if(j<i && abs(mat[j][j])>=sml)continue;
            if(abs(mat[j][i])>abs(mat[max_row][i]))max_row=j;
        }
        swap1(i,max_row,n);
        if(abs(mat[i][i])>=sml){
            double tmp=mat[i][i];
            for(int j=i;j<=n+1;j++){
                mat[i][j]/=tmp;
            }
            for(int j=1;j<=n;j++){
                if(i!=j){
                    double rate=mat[j][i];
                    for(int k=i;k<=n+1;k++){
                        mat[j][k]-=mat[i][k]*rate;
                    }
                }
            }
        }
    }
}
int main(){
    int n;
    cin>>n;
    for(int i=1;i<=n+1;i++){
        for(int j=1;j<=n;j++){
            cin>>a[i][j];
        }
    }
    //n个方程组
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n;j++){
            mat[i][j]=2*(a[i][j]-a[i+1][j]);
            mat[i][n+1]+=a[i][j]*a[i][j]-a[i+1][j]*a[i+1][j];
        }
    }
    gauss(n);
    for(int i=1;i<=n;i++){
        printf("%.3lf ",mat[i][n+1]);
    }
    return 0;
}

洛谷-5027-三角形称重

#include<bits/stdc++.h>
using namespace std;
double sml=1e-7;
vector<vector<double>> swap1(int x,int y,int n,vector<vector<double>> mat){
    if(x==y)return mat;
    for(int i=1;i<=n+1;i++){
        double tmp=mat[x][i];
        mat[x][i]=mat[y][i];
        mat[y][i]=tmp;
    }
    return mat;
}
vector<vector<double>> gauss(int n,vector<vector<double>> mat){
    for(int i=1;i<=n;i++){
        int max_row=i;
        for(int j=1;j<=n;j++){
            if(j<i && abs(mat[j][j])>=sml)continue;
            if(abs(mat[j][i])>abs(mat[max_row][i]))max_row=j;
        }
        mat=swap1(i,max_row,n,mat);
        if(abs(mat[i][i])>=sml){
            double tmp=mat[i][i];
            for(int j=i;j<=n+1;j++){
                mat[i][j]/=tmp;
            }
            for(int j=1;j<=n;j++){
                if(i!=j){
                    double rate=mat[j][i];//mat[i][i]可以去掉貌似
                    for(int k=i;k<=n+1;k++){
                        mat[j][k]-=mat[i][k]*rate;
                    }
                }
            }
        }
    }
    return mat;
}
int judge(int n,vector<vector<double>> mat){
    int maxt=0;
    double maxv=0;
    int ans=0;
    for(int i=1;i<=n;i++){
        if(abs(mat[i][i])<sml)return 0;
        if(mat[i][n+1]<=0 || mat[i][n+1]!=(int)mat[i][n+1])return 0;
        if(maxv<mat[i][n+1]){
            maxv=mat[i][n+1];
            maxt=1;
            ans=i;
        }else if(maxv==mat[i][n+1]){
            maxt++;
        }
    }
    if(maxt>1)return 0;
    return ans;
}
int main(){
    int n,sum=0,ans=0;
    cin>>n;
    vector<vector<double>> mat(n+2,vector<double>(n+2,0));
    for(int i=1;i<=n+1;i++){
        int cnt,a;
        cin>>cnt;
        for(int j=0;j<cnt;j++){
            cin>>a;
            mat[i][a]=1;
        }
        cin>>a;
        mat[i][n+1]=a;
    }
    vector<vector<double>> tmp;
    for(int i=1;i<=n+1;i++){
        tmp=gauss(n,swap1(i,n+1,n,mat));
        int cur=judge(n,tmp);
        if(cur){
            sum++;
            if(sum==1){
                ans=cur;
            }
        }
        if(sum>1)break;
    }
    if(sum==1)cout<<ans<<'\n';
    else cout<<"illegal\n";
    return 0;
}

异或方程组

杭电5833-完全平方数

难点:1.发现异或方程组 2.总方案数：由于自由元的确定后，主元便随之确定，且自由元之间互不影响，所以总方案数为2**（自由元数）-1（减去全不选的情况）

#include <bits/stdc++.h>
using namespace std;
#define N 1000000007
int mat[306][306];
int primes[] = {0, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61,67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511, 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657, 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889, 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997, 1999};
long long power2[301];
void swap1(int x, int y, int n)
{
    for (int i = 1; i <= n + 1; i++)
    {
        int tmp = mat[x][i];
        mat[x][i] = mat[y][i];
        mat[y][i] = tmp;
    }
}
void gauss(int n)
{
    for (int i = 1; i <= n; i++)
    {
        for (int j = 1; j <= n; j++)
        {
            if (j < i && mat[j][j] == 1)
                continue;
            if (mat[j][i] == 1)
            {
                swap1(j, i, n);
                break;
            }
        }
        if (mat[i][i] == 1)
        {
            for (int j = 1; j <= n; j++)
            {
                if (i != j && mat[j][i] == 1)
                {
                    for (int k = i; k <= n + 1; k++)
                    {
                        mat[j][k] ^= mat[i][k];
                    }
                }
            }
        }
    }
}
int solve(int n)
{
    for (int i = 1; i <= 303; i++)
    {
        for (int j = 1; j <= 304; j++)
        {
            mat[i][j] = 0;
        }
    }
    long long x;
    for (int i = 1; i <= n; i++)
    {
        cin >> x;
        for (int j = 1; j <= 303 && x != 0; j++)
        {
            while (x % primes[j] == 0)
            {
                mat[j][i] ^= 1;
                x /= primes[j];
            }
        }
    }
    gauss(303);
    int maincnt = 0;
    for (int i = 1; i <= 303; i++)
    {
        if (mat[i][i])maincnt++;
    }
    return power2[n - maincnt] - 1;
}
int main()
{
    power2[0] = 1;
    for (int i = 1; i < 301; i++)
    {
        power2[i] = power2[i - 1] * 2 % N;
    }
    int t, n,qcnt=0;
    cin>>t;
    while (t--)
    {
        cin >> n;
        qcnt++;
        printf("Case #%d:\n%d\n",qcnt,solve(n));
    }
    return 0;
}

洛谷-2962-令点全为1

难点：操作最少的方案数
用dfs

#include <bits/stdc++.h>
using namespace std;
#define N 1000000007
int op[40];//记录每个自由元是否操作
int n;
int mat[40][40];
int ans=0;//记录最小操作数
void swap1(int x, int y, int n)
{
    for (int i = 1; i <= n + 1; i++)
    {
        int tmp = mat[x][i];
        mat[x][i] = mat[y][i];
        mat[y][i] = tmp;
    }
}
void gauss(int n)
{
    for (int i = 1; i <= n; i++)
    {
        for (int j = 1; j <= n; j++)
        {
            if (j < i && mat[j][j] == 1)
                continue;
            if (mat[j][i] == 1)
            {
                swap1(j, i, n);
                break;
            }
        }
        if (mat[i][i] == 1)
        {
            for (int j = 1; j <= n; j++)
            {
                if (i != j && mat[j][i] == 1)
                {
                    for (int k = i; k <= n + 1; k++)
                    {
                        mat[j][k] ^= mat[i][k];
                    }
                }
            }
        }
    }
}
void dfs(int x,int num){//x代表遍历到了编号几，num代表以及操作了几个点
    if(num>=ans)return;//剪枝
    if(x==0){
        ans=num;
        return;
    }
    if(mat[x][x]==0){
        op[x]=0;
        dfs(x-1,num);
        op[x]=1;
        dfs(x-1,num+1);
    }else{
        int cur=mat[x][n+1];
        for(int i=x+1;i<=n;i++){
            if(op[i])cur^=mat[x][i];
        }
        dfs(x-1,num+cur);
    }
}
int main()
{
    ios::sync_with_stdio(false);
    cin.tie(0);
    int m,a,b;
    cin>>n>>m;
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n;j++){
            mat[i][j]=0;
        }
        mat[i][i]=1;
        mat[i][n+1]=1;
    }
    for(int i=1;i<=m;i++){
        cin>>a>>b;
        mat[a][b]=1;
        mat[b][a]=1;
    }
    gauss(n);
    int flag=0;
    for(int i=1;i<=n;i++){
        if(mat[i][i]==0){
            flag=1;
            break;
        }
        ans+=mat[i][n+1];
    }
    if(flag){
        ans=n;
        dfs(n,0);
        cout<<ans<<'\n';
    }else{
        cout<<ans<<'\n';
    }
    return 0;
}

外星千足虫-洛谷-2447

难点：1.规模庞大，有2000,需采用位图存储 2.如何得到需要的最少记录数（“即到第k次记录时就能够确定所有虫子种类”）：用need变量，在gauss中每次找寻最大值时，更新need为max(j,need)

#include <bits/stdc++.h>
using namespace std;
const int MAXN=2002;
const int BIT=64;
const int MAXM=MAXN/BIT+1;
long long mat[MAXN][MAXM];
int need=0;
int n,m,max1;
void swap1(int x, int y, int bits){
    for(int i=0;i<=bits/BIT;i++){
        long long tmp=mat[x][i];
        mat[x][i]=mat[y][i];
        mat[y][i]=tmp;
    }
}
int get1(int row,int col){
    return ((mat[row][col/BIT]>>(col%BIT))&1);
}
void set1(int row,int col,int v){
    if(v==0)mat[row][col/BIT]&=(~(1LL<<(col%BIT)));
    else mat[row][col/BIT]|=(1LL<<(col%BIT));
}

void exo(int x,int y,int bits){
    for(int i=0;i<=bits/BIT;i++){
        mat[y][i]^=mat[x][i];
    }
}
void gauss(int nn)
{
    for (int i = 1; i <= nn; i++)
    {
        for (int j = i; j <= nn; j++)
        {
            if (get1(j,i))
            {
                need=max(need,j);
                swap1(j, i, nn+2);
                break;
            }
        }
        if (get1(i,i))
        {
            for (int j = 1; j <= nn; j++)
            {
                if (i != j && get1(j,i))exo(i,j,nn+2);
            }
        }
        else{
            return;
        }
    }
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(0);
    int m;
    string s;
    int b;
    cin>>n>>m;
    max1=max(n,m);
    for(int i=1;i<=m;i++){
        cin>>s>>b;
        for(int j=0;j<n;j++){
            set1(i,j+1,s[j]-'0');
        }
        set1(i,max1+1,b);
    }
    gauss(max1);
    int flag=0;
    for(int i=n;i>=1;i--){
        if(get1(i,i)==0){
            flag=1;
            break;
        }
    }
    if(flag){
        cout<<"Cannot Determine\n";
    }else{
        cout<<need<<'\n';
        for(int i=1;i<=n;i++){
            cout<<(get1(i,max1+1) ? "?y7M#" : "Earth")<<'\n';
        }
    }
    return 0;
}

同余方程组

模版

void gauss1(int n){
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n;j++){
            if(j<i && mat[j][j]!=0)continue;
            if(mat[j][i]){
                swap1(j,i,n);
                break;
            }
        }
        if(mat[i][i]){
            for(int j=1;j<=n;j++){
                if(j!=i && mat[j][i]){
                    int gcd1=gcd(mat[i][i],mat[j][i]);
                    int a=mat[i][i]/gcd1;
                    int b=mat[j][i]/gcd1;
                    if(j<i && mat[j][j]){
                        for(int k=j;k<i;k++){
                            mat[j][k]=(mat[j][k]*a)%MOD;
                        }
                    }
                    for(int k=i;k<=n+1;k++){
                        mat[j][k]=((mat[j][k]*a-mat[i][k]*b)%MOD+MOD)%MOD;
                    }
                }
            }
        }
    }
    //把能求解的主元系数统一变成1
    for(int i=1;i<=n;i++){
        if(mat[i][i]){
            bool flag=0;
            for(int j=i+1;j<=n;j++){
                if(mat[i][j]){
                    flag=1;
                    break;
                }
            }
            if(!flag){
                mat[i][n+1]=mat[i][n+1]*inv[mat[i][i]]%MOD;
                mat[i][i]=1;
            }
        }
    }
}

hd-5755

在这里插入图片描述

#include <bits/stdc++.h>
using namespace std;
const int MOD=3;
int inv[MOD];
int ans;
int mat[904][904];
int g[32][32];
int dx[]={-1,0,1,0};
int dy[]={0,-1,0,1};
void swap1(int x, int y, int n){
    for (int i = 1; i <= n + 1; i++){
        int tmp = mat[x][i];
        mat[x][i] = mat[y][i];
        mat[y][i] = tmp;
    }
}
int gcd(int a,int b){
    return b==0 ? a : gcd(b,a%b);
}

void gauss(int n){
    for(int i=1;i<=n;i++){
        //寻找阶段
        for(int j=1;j<=n;j++){
            if(j<i && mat[j][j]!=0)continue;
            if(mat[j][i]){
                swap1(j,i,n);
                break;
            }
        }
        //消除阶段
        if(mat[i][i]){
            for(int j=1;j<=n;j++){
                if(j!=i && mat[j][i]){
                    int gcd1=gcd(mat[i][i],mat[j][i]);
                    int a=mat[i][i]/gcd1;
                    int b=mat[j][i]/gcd1;
                    if(j<i && mat[j][j]){
                        // for(int k=j;k<i;k++){
                        //     mat[j][k]=(mat[j][k]*a)%MOD;
                        // }
                        //和一般的gauss不同，因为这题令自由元全为0，故主元不会被自由元影响，所以无需上面的操作
                        mat[j][j]=(mat[j][j]*a)%MOD;
                    }
                    for(int k=i;k<=n+1;k++){
                        mat[j][k]=((mat[j][k]*a-mat[i][k]*b)%MOD+MOD)%MOD;
                    }
                }
            }
        }
    }
    ans=0;
    //和一般的gauss不同把能求解的主元系数统一变成1，因为这题令自由元全为0，故主元不会被自由元影响，于是可以将所有的主元的系数化为1
    for(int i=1;i<=n;i++){
        if(mat[i][i]){
            mat[i][n+1]=mat[i][n+1]*inv[mat[i][i]]%MOD;
            ans+=mat[i][n+1];
        }
    }
}
int main(){
    //逆元数组
    inv[1]=1;
    for(int i=2;i<=MOD;i++){
        inv[i]=(int)(MOD-(long long)inv[MOD%i]*(MOD/i)%MOD);
    }
    int n,m,t;
    cin>>t;
    while(t--){
        cin>>n>>m;
        //初始化
        for(int i=1;i<=n*m;i++){
            for(int j=1;j<=n*m+1;j++){
                mat[i][j]=0;
            }
        }
        //读入数据
        for(int i=1;i<=n;i++){
            for(int j=1;j<=m;j++){
                cin>>g[i][j];
            }
        }
        //设置mat数组
        for(int i=1;i<=n;i++){
            for(int j=1;j<=m;j++){
                mat[(i-1)*m+j][(i-1)*m+j]=2;
                for(int k=0;k<4;k++){
                    int a=i+dx[k],b=j+dy[k];
                    if(a>=1 && a<=n && b>=1 && b<=m){
                        mat[(i-1)*m+j][(a-1)*m+b]=1;
                    }
                }
                mat[(i-1)*m+j][n*m+1]=3-g[i][j]%3;
            }
        }
        gauss(n*m);
        cout<<ans<<'\n';
        for(int i=1;i<=n*m;i++){
            if(mat[i][i]){
                for(int k=0;k<mat[i][n*m+1];k++){
                    cout<<(i-1)/m+1<<' '<<(i-1)%m+1<<'\n';
                }
            }
        }
    }
    return 0;
}

poj-2947

在这里插入图片描述

#include <bits/stdc++.h>
using namespace std;
const int MOD=7;
map<string,int> mapp;
int inv[MOD];
int n,m,s,k,tool;
string l,r;
int mat[904][904];
void swap1(int x, int y, int n){
    for (int i = 1; i <= n + 1; i++){
        int tmp = mat[x][i];
        mat[x][i] = mat[y][i];
        mat[y][i] = tmp;
    }
}
int gcd(int a,int b){
    return b==0 ? a : gcd(b,a%b);
}
void gauss(int n){
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n;j++){
            if(j<i && mat[j][j]!=0)continue;
            if(mat[j][i]){
                swap1(j,i,n);
                break;
            }
        }
        if(mat[i][i]){
            for(int j=1;j<=n;j++){
                if(j!=i && mat[j][i]){
                    int gcd1=gcd(mat[i][i],mat[j][i]);
                    int a=mat[i][i]/gcd1;
                    int b=mat[j][i]/gcd1;
                    if(j<i && mat[j][j]){
                        for(int k=j;k<i;k++){
                            mat[j][k]=(mat[j][k]*a)%MOD;
                        }
                    }
                    for(int k=i;k<=n+1;k++){
                        mat[j][k]=((mat[j][k]*a-mat[i][k]*b)%MOD+MOD)%MOD;
                    }
                }
            }
        }
    }
    //把能求解的主元系数统一变成1
    for(int i=1;i<=n;i++){
        if(mat[i][i]){
            bool flag=0;
            for(int j=i+1;j<=n;j++){
                if(mat[i][j]){
                    flag=1;
                    break;
                }
            }
            if(!flag){
                mat[i][n+1]=mat[i][n+1]*inv[mat[i][i]]%MOD;
                mat[i][i]=1;
            }
        }
    }
}

int main(){
    mapp["MON"]=1,mapp["TUE"]=2,mapp["WED"]=3,mapp["THU"]=4,mapp["FRI"]=5,mapp["SAT"]=6,mapp["SUN"]=0;
    //逆元数组
    inv[1]=1;
    for(int i=2;i<=MOD;i++){
        inv[i]=(int)(MOD-(long long)inv[MOD%i]*(MOD/i)%MOD);
    }
    while(1){
        cin>>n>>m;
        if(!n)break;
        s=max(n,m);
        for(int i=1;i<=s;i++){
            for(int j=1;j<=s+1;j++){
                mat[i][j]=0;
            }
        }
        for(int i=1;i<=m;i++){
            cin>>k;
            cin>>l>>r;
            while(k--){
                cin>>tool;
                mat[i][tool]=(mat[i][tool]+1)%MOD;
            }
            mat[i][s+1]=(mapp[r]-mapp[l]+1+MOD)%MOD;
        }
        gauss(s);
        int flag=1;
        for(int i=1;i<=s;i++){
            if(!mat[i][i] && mat[i][s+1]){
                flag=-1;
                break;
            }
            if(i<=n && mat[i][i]==0)flag=0;
        }
        if(flag==-1){
            cout<<"Inconsistent data.";
        }else if(flag==0){
            cout<<"Multiple solutions.";
        }else{
            for(int i=1;i<=n;i++){
                if(mat[i][s+1]<3)mat[i][s+1]+=7;
                cout<<mat[i][s+1]<<' ';
            }
        }
        cout<<'\n';
    }
    return 0;
}