当前位置：首页 > news >正文

数论~~~

news 来源：原创 2025/6/6 12:17:16

质数

质数
- 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;}}}}}//把能求解的主元系数统一变成1for(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，故主元不会被自由元影响，于是可以将所有的主元的系数化为1for(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;}}}}}//把能求解的主元系数统一变成1for(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;
}