数论练习（博客园）

Yolo@～

已于 2024-09-21 01:01:10 修改

阅读量718

点赞数 17

文章标签： c++ 算法数据结构

于 2024-09-21 00:56:31 首次发布

本文链接：https://blog.csdn.net/2301_80225424/article/details/140923783

版权

1——快速幂

Luogu P1226 【模板】快速幂||取余运算

#include<iostream>
using namespace std;
typedef long long ll;
ll binpow(ll a, ll b, ll p)
{
	ll ans = 1;
	a = a % p;
	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;
	cout << a << '^' << b << " mod " << p << '=' << binpow(a, b, p);
	return 0;
}

2——高精度快速幂

Luogu P1045 [NOIP2003 普及组] 麦森数

#include <iostream>
#include <cstring>
#include <algorithm>
#include <vector>
#include <cmath>
using namespace std;

const int N=500;
int a[N*2],res[N*2],t[N*2];
int p;

void mul(int*a, int*b, int*c){//高精度
  memset(t,0,sizeof(t));
  for(int i=0; i<N; i++)
    for(int j=0; j<N; j++){
      t[i+j] += a[i]*b[j];
      t[i+j+1] += t[i+j]/10;
      t[i+j] %= 10;    
    }
  memcpy(c,t,sizeof(t));
}
void quickpow(int p){//快速幂
  res[0]=1, a[0]=2;
  while(p){
    if(p & 1) mul(res,a,res);
    mul(a,a,a);
    p >>= 1;
  }
  res[0]--; //个位修正
}
int main(){
  cin >> p;
  printf("%d\n",int(p*log10(2))+1);
  quickpow(p);
  for(int i=0, k=499; i<10; i++){
    for(int j=0; j<50; j++, k--)
      printf("%d",res[k]);
    puts("");
  }
  return 0;
}

计算位数

2^p的位数为p*log10(2)+1;

3——矩阵快速幂

Luogu P3390 【模板】矩阵快速幂

#include <iostream>
#include <cstring>
#include <algorithm>
#include <vector>
#include <cmath>
typedef long long ll;
using namespace std;

const ll p = 1e9 + 7;
struct mat {
	ll c[101][101];
}A, res;
ll n, k;
mat mul(mat& a, mat &b)
{
	mat t;//临时矩阵
	memset(t.c, 0, sizeof(t.c));
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n; j++)
		{
			for (int k = 1; k <= n; k++)
			{
				t.c[i][j] = (t.c[i][j] + a.c[i][k] * b.c[k][j]) % p;
			}
		}
	}
	return t;
}
void quickpow(ll k)
{
	for (int i = 1; i <= n; i++)res.c[i][i] = 1;
	while (k)
	{
		if (k & 1)res = mul(res, A);
		A = mul(A, A);
		k >>= 1;
	}
}

int main()
{
	cin >> n >> k;
	for (int i = 1; i <= n; i++)
		for (int j = 1; j <= n; j++)
			cin >> A.c[i][j];
	quickpow(k);
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n; j++)
			cout << res.c[i][j] << ' ';
		cout << '\n';
	}
	return 0;
}

4——矩阵快速幂（斐波那契数列）

Luogu P1962 斐波那契数列

#include<iostream>  
#include<cstring> // 引入cstring库以使用memset  
using namespace std;
typedef long long ll;

struct mat {
    ll c[3][3];
    mat() { memset(c, 0, sizeof(c)); } // 初始化矩阵所有元素为0  
} F, A;

const ll p = 1e9 + 7;

mat operator *(const mat& a, const mat& b) {
    mat t;
    for (int i = 0; i < 2; i++) { // 矩阵大小为2x2，所以索引从0到1  
        for (int j = 0; j < 2; j++) {
            for (int k = 0; k < 2; k++) { // 注意这里应该是k，且k的范围也是0到1  
                t.c[i][j] = (t.c[i][j] + a.c[i][k] * b.c[k][j]) % p;
            }
        }
    }
    return t;
}

void quickpow(ll n) {
    F.c[0][0] = F.c[0][1] = 1; // 初始化F为单位矩阵  
    A.c[0][0] = A.c[0][1] = A.c[1][0] = 1; // 初始化A为给定的转移矩阵  
    A.c[1][1] = 0; // 根据题目要求设置A矩阵的最后一个元素  

    while (n) {
        if (n & 1) F = F * A;
        A = A * A;
        n >>= 1;
    }
}

int main() {
    ll n;
    cin >> n;
    if (n <= 2) {
        cout << "1";
        return 0;
    }
    quickpow(n - 2); // 题目要求的是第n项，但矩阵是从第0项开始计算的，所以需要n-2  
    cout << F.c[0][0]; // 输出矩阵的左上角元素，即第n项的值  
    return 0;
}

Luogu P1939 【模板】矩阵加速（数列）

#include <iostream>
#include <cstring>
#include <algorithm>
using namespace std;

typedef long long LL;
const int mod=1000000007;
struct matrix{
  LL c[4][4];
  matrix(){memset(c, 0, sizeof c);}
} F, A; //F为斐波那契矩阵，A为构造矩阵

matrix operator*(matrix &a, matrix &b){ //矩阵乘法
  matrix t; //临时矩阵
  for(int i=1; i<=3; ++i)
    for(int j=1; j<=3; ++j)
      for(int k=1; k<=3; ++k){
        t.c[i][j] += a.c[i][k]*b.c[k][j]%mod;
        t.c[i][j] %= mod;
      }
  return t;
}
void quickpow(LL n){ //快速幂
  F.c[1][1]=F.c[1][2]=F.c[1][3]=1;
  memset(A.c, 0, sizeof(A.c));  
  A.c[1][1]=A.c[1][2]=A.c[2][3]=A.c[3][1]=1;
  while(n){
    if(n & 1) F = F*A;
    A = A*A;
    n >>= 1;
  }
}
int main(){
  LL T,n; cin>>T;
  while(T--){
    cin>>n;
    if(n<=3){puts("1"); continue;}
    quickpow(n-3);
    cout << F.c[1][1] << endl;    
  }
  return 0;
}

5——最大公约数问题

Luogu P1029 [NOIP2001 普及组] 最大公约数和最小公倍数问题

#include<iostream>
using namespace std;
int gcd(int a, int b)
{
	return (b == 0) ? a : gcd(b, a % b);
}
int main()
{
	int x, y;
	cin >> x >> y;
	int num = 0;
	int t = x * y;
	for (int i = 1; i * i <= t; i++)
		if (t % i == 0 && gcd(i, t / i) == x)
			num += 2;
	if (x == y)num--;
	cout << num;
	return 0;
}

6——判定质数（试除法）

Luogu P5736 【深基7.例2】质数筛

#include<iostream>
using namespace std;
bool prime(int n)
{
	if (n < 2)return 0;
	for (int i = 2; i * i <= n; i++)
		if (n % i == 0)return 0;
	return 1;
}
int main()
{
	int n, m;
	cin >> n;
	while (n--)
	{
		cin >> m;
		if (prime(m))cout << m << ' ';
	}
	return 0;
}

7——质因子分解

Luogu P2043 质因子分解

#include <iostream>
#include <cstring>
#include <algorithm>
using namespace std;

int n;
int a[10001]; //质因子的个数

void decompose(int x){ //分解质因数
  for(int i=2; i*i<=x; i++)
    while(x%i==0) a[i]++, x/=i;
  if(x>1) a[x]++;
}
int main(){
  cin >> n;
  for(int i=2; i<=n; i++) decompose(i);
    for(int i=1;i<=n;i++)
    if(a[i]) cout<<i<<" "<<a[i]<<endl;
  return 0;
}

8——线性筛质数

//埃氏筛法筛质数
#include <iostream>
#include <cstring>
#include <algorithm>
using namespace std;

typedef long long LL;
const int N = 100000010;
int vis[N];  //划掉合数
int prim[N]; //记录质数
int cnt; //质数个数

void Eratosthenes(int n){ //埃氏筛法
  for(LL i=2; i<=n; ++i){
    if(!vis[i]){
      prim[++cnt] = i;
      for(LL j=i*i; j<=n; j+=i)
        vis[j] = 1;
    }
  }
}
int main(){
  int n, q, k;
  scanf("%d %d", &n, &q);
  Eratosthenes(n);
  while(q--){
    scanf("%d", &k);
    printf("%d\n", prim[k]);
  }
  return 0;
}

#include<iostream>
#include<cstring>
#include<algorithm>
using namespace std;
const int N = 1e8 + 10;
int vis[N];
int prime[N];
int cnt;

void get_prim(int n) {
	for (int i = 2; i <= n; i++)
	{
		if (!vis[i])prime[++cnt] = i;
		for (int j = 1; i * prime[j] <= n; j++)
		{
			vis[i * prime[j]] = 1;
			if (i % prime[j] == 0)break;
		}
	}
}
int main()
{
	int n, q, k;
	scanf("%d %d", &n, &q);
	get_prim(n);
	while (q--)
	{
		scanf("%d", &k);
		printf("%d\n", prime[k]);
	}
	return 0;
}

9——筛法求欧拉函数

#include<iostream>
using namespace std;
int phi(int n)
{
	int res = n;
	for (int i = 2; i * i <= n; i++)
	{
		if (n % i == 0)
		{
			res = res / i * (i - 1);
			while (n % i == 0)n /= i;
		}
	}
	if (n > 1)res = res / n * (n - 1);
	return res;
}
int main()
{
	int n;
	cin >> n;
	cout << phi(n) << '\n';
	return 0;
}

#include <iostream>
using namespace std;

const int N = 1000010;
int p[N], vis[N], cnt;
int phi[N];

void get_phi(int n){//筛法求欧拉函数
  phi[1] = 1;
  for(int i=2; i<=n; i++){
    if(!vis[i]){
      p[cnt++] = i;
      phi[i] = i-1;
    }
    for(int j=0; i*p[j]<=n; j++){
      int m = i*p[j];
      vis[m] = 1;
      if(i%p[j] == 0){
        phi[m] = p[j]*phi[i];
        break;
      }
      else
        phi[m]=(p[j]-1)*phi[i];
    }
  }
}
int main(){
  int n;
  cin >> n;
  get_phi(n);
  for(int i=1; i<=n; i++)
    printf("%d\n", phi[i]);
  return 0;
}

10——筛法求约数个数

#include<iostream>
using namespace std;
const int N = 1e6 + 10;

int p[N], notprime[N], cnt;
int a[N];//记录i的最小质因子的次数
int d[N];//记录约数个数
void get_d(int n)
{
	d[1] = 1;
	for (int i = 2; i <= n; i++)
	{
		if (!notprime[i]) {
			p[++cnt] = i;
			a[i] = 1;
			d[i] = 2;
		}
		for (int j = 1; i * p[j] <= n; j++) {
			int m = i * p[j];
			notprime[m] = 1;
			if (i % p[j] == 0)
			{
				a[m] = a[i] + 1;
				d[m] = d[i] / a[m] * (a[m] + 1);
				break;
			}
			else
			{
				a[m] = 1;
				d[m] = d[i] * 2;
			}
		}
	}
}
int main()
{
	int n;
	cin >> n;
	get_d(n);
	for (int i = 1; i <= n; i++)
		cout << d[i] << '\n';
	return 0;
}

11——筛法求约数和

#include <iostream>
using namespace std;

const int N = 1000010;
int p[N], vis[N], cnt;
//g[i]表示i的最小质因子的1+p^1+...+p^k
int g[N], f[N];//f[i]表示i的约数和

void get_f(int n){ //筛法求约数和
  g[1] = f[1] = 1;
  for(int i=2; i<=n; i++){
    if(!vis[i]){
      p[++cnt] = i;
      g[i] = f[i] = i+1;
    }
    for(int j=1; i*p[j]<=n; j++){
      int m = i*p[j]; 
      vis[m] = 1;
      if(i%p[j] == 0){
        g[m] = g[i]*p[j]+1;
        f[m] = f[i]/g[i]*g[m];
        break;
      } 
      else{
        g[m] = p[j]+1;
        f[m] = f[i]*g[m];
      }
    }
  }
}
int main(){
  int n;
  cin >> n;
  get_f(n);
  for(int i=1; i<=n; i++)
    printf("%d\n",f[i]);
  return 0;
}

11——筛法求莫比乌斯函数

#include <iostream>
using namespace std;

const int N = 1000010;
int p[N], vis[N], cnt;
int mu[N];

void get_mu(int n){//筛法求莫比乌斯函数
  mu[1] = 1;
  for(int i=2; i<=n; i++){
    if(!vis[i]){
      p[++cnt] = i;
      mu[i] = -1;
    }
    for(int j=1; i*p[j]<=n; j++){
      int m = i*p[j]; 
      vis[m] = 1;
      if(i%p[j] == 0){
        mu[m] = 0;
        break;
      } 
      else
        mu[m] = -mu[i];
    }
  }
}
int main(){
  int n;
  cin >> n;
  get_mu(n);
  for(int i=1; i<=n; i++)
    printf("%d\n",mu[i]);
  return 0;
}

12——扩展欧拉定理

核心逻辑概述

这段代码的主要目的是计算并输出一个模幂运算的结果，即 abmodm，其中 b 是一个非常大的数（以字符串形式给出）。为了实现这一点，代码使用了几个关键的数学概念和编程技巧：

欧拉函数（Euler's Totient Function）：用于计算与模数 m 互质的数的个数，这在处理模幂运算中的降幂问题时非常有用。
降幂：由于 b 非常大，直接计算 ab 是不现实的。因此，代码通过 depow 函数将 b 转换为一个更小的等效值（在模 ϕ(m) 的意义下），这里 ϕ(m) 是 m 的欧拉函数。如果 a 和 m 互质，根据欧拉定理，这可以大大减少计算量。
快速幂：qpow 函数使用快速幂算法来高效地计算 abmodm。快速幂算法通过将指数 b 分解为二进制形式，并迭代地应用平方和乘法来减少乘法操作的次数。

代码逻辑细节

输入处理：代码通过 scanf 读取整数 a 和 m，以及一个表示大整数 b 的字符串 s。
计算欧拉函数：get_phi 函数计算并返回 m 的欧拉函数 ϕ(m)。
降幂处理：depow 函数遍历字符串 s，将 b 转换为整数（如果可能的话），并在这个过程中检查是否需要应用降幂。如果 b≥ϕ(m)，则通过取模和可能的加法来调整 b。
模幂运算：qpow 函数使用快速幂算法计算 abmodm，其中 b 可能是通过 depow 函数处理后的值。
输出：最后，代码使用 printf 输出模幂运算的结果。

13——裴蜀定理

Luogu P4549 【模板】裴蜀定理

#include<iostream>
#include<cmath>
using namespace std;

int gcd(int a, int b)
{
	return (b == 0) ? a : gcd(b, a % b);
}
int main()
{
	int n, a, s;
	cin >> n;
	cin >> a;
	s = a;
	for (int i = 2; i <= n; i++)
	{
		cin >> a;
		s = gcd(s, abs(a));
	}
	cout << s;
	return 0;
}

不定方程扩展欧几里得算法

Luogu P1516 青蛙的约会

#include<iostream>
#include<cmath>
using namespace std;
typedef long long ll;
ll ex_gcd(ll a, ll b, ll& x, ll& y)
{
	if (b == 0)
	{
		x = 1; y = 0; return a;
	}
	ll d = ex_gcd(b, a % b, x, y);
	ll t = x;
	x = y;
	y = t - a / b * y;
	return d;
}
int main()
{
	ll m, n, x, y, l, x0, y0;
	cin >> x >> y >> m >> n >> l;
	ll a = n - m;
	ll b = l;
	ll c = x - y;
	if (a < 0)
	{
		a = -a;
		c = -c;
	}
	ll g = ex_gcd(a, b, x0, y0);
	if (c % g != 0) { cout << "Impossible"; return 0; }
		x0 = x0 * (c / g);
		cout << ((x0 %(b / g) + (b / g)) % (b / g));
	return 0;
}

同余方程乘法逆元扩欧算法

//乘法逆元
#include <iostream>
#include <cstring>
#include <algorithm>
using namespace std;
    
int exgcd(int a,int b,int &x,int &y){
  if(b == 0){
    x = 1, y = 0; return a;
  }
  int x1, y1, d;
  d = exgcd(b, a%b, x1, y1);
  x = y1, y = x1-a/b*y1;
  return d;
}
int main(){
  int a, m, x, y;
  scanf("%d%d", &a, &m);
  exgcd(a, m, x, y);
  printf("%d", (x%m+m)%m);
  return 0;
}

//同余方程
#include <iostream>
#include <cstring>
#include <algorithm>
using namespace std;
    
int exgcd(int a,int b,int &x,int &y){
  if(b == 0){x = 1, y = 0; return a;}
  int x1, y1, d;
  d = exgcd(b, a%b, x1, y1);
  x = y1, y = x1-a/b*y1;
  return d;
}
int main(){
  int a, b, m, x, y;
  scanf("%d%d%d", &a, &b, &m);
  int d = exgcd(a, m, x, y);
  if(b%d == 0) 
    printf("%d", 1ll*x*b/d%m);
  else puts("none");
  return 0;
}

同余式乘法逆元费马小定理

#include<iostream>
using namespace std;

typedef long long LL;
int a, p;

int quickpow(LL a, int b, int p){
  int res = 1;
  while(b){
    if(b & 1) res = res*a%p;
    a = a*a%p;
    b >>= 1;
  }
  return res;
}
int main(){
  cin >> a >> p;
  if(a % p)
    printf("%d\n",quickpow(a,p-2,p));
  return 0;
}

中国剩余定理

#include<iostream>
using namespace std;
typedef long long ll;

ll n, a[11], b[11];
ll ex_gcd(ll a, ll b, ll& x, ll& y)
{
	if (b == 0)
	{
		x = 1; y = 0; return a;
	}
	ll d, x1, y1;
	d=ex_gcd(b, a % b, x1, y1);
	x = y1;
	y = x1 - a / b * y1;
	return d;
}
ll CRT(ll m[], ll r[]) {
	ll M = 1, ans = 0;
	for (int i = 1; i <= n; i++)M *= m[i];
	for (int i = 1; i <= n; i++)
	{
		ll mi = M / m[i], x, y;
		ex_gcd(mi, m[i], x, y);
		ans = (ans + r[i] * mi * x %M) % M;
	}
	return (ans % M + M) % M;
}
int main()
{
	cin >> n;
	for (int i = 1; i <= n; i++)
	{
		cin >> a[i] >> b[i];
	}
	cout << CRT(a, b) << '\n';
	return 0;
}

扩展中国剩余定理

Luogu P4777 【模板】扩展中国剩余定理（EXCRT）

#include <iostream>
#include <cstring>
#include <algorithm>
using namespace std;

typedef __int128 LL;
const int N = 100005;
LL n, m[N], r[N];

LL exgcd(LL a,LL b,LL &x,LL &y){
  if(b==0){x=1, y=0; return a;}
  LL d, x1, y1;
  d = exgcd(b, a%b, x1, y1);
  x = y1, y = x1-a/b*y1;
  return d;
}
LL EXCRT(LL m[], LL r[]){
  LL m1, m2, r1, r2, p, q;
  m1 = m[1], r1 = r[1];
  for(int i=2; i<=n; i++){
    m2 = m[i], r2 = r[i];
    LL d = exgcd(m1,m2,p,q);
    if((r2-r1)%d){return -1;}
    p=p*(r2-r1)/d; //特解
    p=(p%(m2/d)+m2/d)%(m2/d);
    r1 = m1*p+r1;
    m1 = m1*m2/d;
  }
  return (r1%m1+m1)%m1;
}
int main(){
  scanf("%lld", &n);
  for(int i = 1; i <= n; ++i)
    scanf("%lld%lld", m+i, r+i);
  printf("%lld\n",EXCRT(m,r));
  return 0;
}

威尔逊定理

#include <iostream>
using namespace std;

typedef long long LL;
const int N = 1000001;
const int mx = 3000008;
int s[N],p[N],vis[mx],t,n;

void get_prim(){
  for(LL i = 2; i < mx; ++i)
    if(!vis[i]){
      if((i-7)%3 == 0) 
        p[(i-7)/3] = 1;
      for(LL j=i*i; j<mx; j+=i) 
        vis[j] = 1;
    }
}
int main(){
  get_prim();
  for(int i=2; i<N; ++i)
    s[i] = s[i-1]+p[i];
  scanf("%d", &t);
  while(t--){
    scanf("%d", &n);
    printf("%d\n", s[n]);
  }
}

BSGC算法

Luogu P3846 [TJOI2007] 可爱的质数/【模板】BSGS

#include<iostream>  
#include<cmath>
#include<algorithm>
#include<unordered_map>
#include<string>
using namespace std;
typedef long long ll;

ll bsgs(ll a, ll b, ll p) {
	a %= p;
	b %= p;
	if (b == 1)return 0;//x=0;
	ll m = ceil(sqrt(p));//sqrt(p)+1;
	ll t = b;
	unordered_map<int, int>hash;
	hash[b] = 0;
	for (int j = 1; j < m; j++) {
		t = t * a % p;//求b*a^j;
		hash[t] = j;
	}
	ll mi = 1;
	for (int i = 1; i <= m; i++)
	{
		mi = mi * a % p;//求(a^m)
	}
	t = 1;
	for (int i = 1; i <= m; i++) {
		t = t * mi % p;//求(a^m)^i
		if (hash.count(t))
			return i * m - hash[t];
	}
	return -1;//无解；
}
int main()
{
	int a, b, p;
	cin >> p >> a >> b;
	int res = bsgs(a, b, p);
	if (res == -1)cout << "no solution";
	else cout << res << '\n';
	return 0;
}

扩BSGS算法

#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <unordered_map>
using namespace std;

typedef long long LL;

LL gcd(LL a, LL b){
  return b==0?a:gcd(b,a%b);
}
LL exbsgs(LL a, LL b, LL p){
  a %= p; b %= p;
  if(b==1||p==1)return 0;//x=0
  
  LL d, k=0, A=1;
  while(true){
    d = gcd(a,p);
    if(d==1) break;
    if(b%d) return -1; //无解
    k++; b/=d; p/=d;
    A = A*(a/d)%p; //求a^k/D
    if(A==b) return k;
  }
  
  LL m=ceil(sqrt(p));
  LL t = b;
  unordered_map<int,int> hash;
  hash[b] = 0;
  for(int j = 1; j < m; j++){
    t = t*a%p; //求b*a^j
    hash[t] = j;
  }
  LL mi = 1;
  for(int i = 1; i <= m; i++)
    mi = mi*a%p; //求a^m
  t = A;
  for(int i=1; i <= m; i++){
    t = t*mi%p; //求(a^m)^i
    if(hash.count(t)) 
      return i*m-hash[t]+k;
  }
  return -1; //无解
}
int main(){
  LL a, p, b;
  while((scanf("%lld%lld%lld",&a,&p,&b)!=EOF)&&a){
    LL res = exbsgs(a, b, p);
    if(res == -1) puts("No Solution");
    else printf("%lld\n",res);
  }
  return 0;
}

模板线性方程组高斯消元法

Luogu P3389 【模板】高斯消元法

//高斯——约旦消元法

#include <iostream>
#include <algorithm>
#include <cmath>
using namespace std;

const int N=110;
const double eps=1e-6;
int n;
double a[N][N]; //增广矩阵

bool Gauss_Jordan(){
  for(int i=1; i<=n; ++i){    //第i主元
    for(int k=i; k<=n; ++k)   //换非0行
      if(fabs(a[k][i])>eps)
        {swap(a[k],a[i]); break;}
    if(fabs(a[i][i])<eps) return 0;
    
    for(int k=1; k<=n; ++k){  //对角化
      if(k==i) continue;
      for(int j=n+1; j>=i; --j) 
        a[k][j]-=a[k][i]/a[i][i]*a[i][j];
    }
  }
  for(int i=1; i<=n; ++i) 
    a[i][n+1]/=a[i][i];     //除以主元
  return 1;                 //存在唯一解
}
int main(){
  scanf("%d",&n);
  for(int i=1; i<=n; i++)
    for(int j=1; j<=n+1; j++) 
      scanf("%lf",&a[i][j]);
  if(Gauss_Jordan()) 
    for(int i=1; i<=n; i++) 
      printf("%.2lf\n",a[i][n+1]);
  else puts("No Solution");
}

#include<iostream>  
#include<cmath>
#include<algorithm>
#include<unordered_map>
#include<cstdio>
#include<string>
using namespace std;
typedef long long ll;

const int N = 110;
const double eps = 1e-6;
int n;
double a[N][N];//增广矩阵

int gauss()
{
	int c, r;//当前列，当前行
	for (c = r = 0; c < n; c++)
	{
		//1.找到c列的最大行t
		int t = r;
		for (int i = r; i < n; i++)
			if (fabs(a[i][c]) > fabs(a[t][c]))t = i;
		if (fabs(a[t][c]) < eps)continue;//c列已0化
		//2.把最大行换到最上面
		for (int i = c; i < n + 1; i++)swap(a[t][i], a[r][i]);
		//3.把当前行的第一个数变为一
		for (int i = n; i >= c; i--)a[r][i] /= a[r][c];
		//4.把c当前列c下面的所有数全部消为0
		for (int i = r + 1; i < n; i++)
			if (fabs(a[i][c]) > eps)
				for (int j = n; j >= c; j--)
					a[i][j] -= a[i][c] * a[r][j];
		r++;//这一行做完换下一行
	}
	if (r < n)//提前变成梯形矩阵
	{
		for (int i = r; i < n; i++)
			if (fabs(a[i][n]) > eps)//a[i][n]==bi
				return 2;//左边等于0，右边不等于0，无解
		return 1;//无穷多解
	}
	//5.唯一解，从下往上回代，得到方程的解
	for (int i = n - 1; i >= 0; i--)
		for (int j = i + 1; j < n; j++)
			a[i][n] -= a[i][j] * a[j][n];
	return 0;
}
int main()
{
	cin >> n;
	for (int i = 0; i < n; i++)
		for (int j = 0; j <= n; j++)
			cin >> a[i][j];
	int t = gauss();
	if (t)cout << "No Solution";
	else for (int i = 0; i < n; i++)
		printf("%.2lf\n", a[i][n]);
	return 0;
}

#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;
const int N = 110;
const double eps = 1e-6;
int n;
double a[N][N];//增广矩阵
int gauss()
{
	for (int i = 1; i <= n; i++)//第i主元
	{                             
		for(int k=i;k<=n;k++)//换非零行
			if (fabs(a[k][i]) > eps)
			{
				swap(a[k], a[i]); break;
			}
		if (fabs(a[i][i]) < eps)return 0;
		for (int j = n + 1; j >= i; j--)//变1
			a[i][j] /= a[i][i];
		for (int k = i + 1; k <= n; k++)//变0
			for (int j = n + 1; j >= i; j--)
				a[k][j] -= a[k][i] * a[i][j];
	}
	for (int i = n - 1; i >= 1; i--)
	{
		for (int j = i + 1; j <= n; j++)
			a[i][n + 1] -= a[i][j] * a[j][n + 1];
		return 1;//存在唯一解
	}
}
int main()
{
	cin >> n;
	for (int i = 1; i <= n; i++)
		for (int j = 1; j <= n; j++)
			cin >> a[i][j];
	if (gauss())
		for (int i = 1; i < n + 1; i++)
			printf("%.2lf\n", a[i][n + 1]);
	else cout << "No Solution";
	return 0;
}

模板矩阵求逆

Luogu P4783 【模板】矩阵求逆

#include<iostream>
#include<cstdio>
#include<cmath>
#define LL long long
using namespace std;

const int N=405,P=1e9+7;
int n;
LL a[N][N<<1];

LL quickpow(LL a, LL b){
  LL ans = 1;
  while(b){
    if(b & 1) ans = ans*a%P;
    a = a*a%P;
    b >>= 1;
  }
  return ans;
}
bool Gauss_Jordan(){    
  for(int i=1;i<=n;++i){ //枚举主元的行列
    int r = i;
    for(int k=i; k<=n; ++k) //找非0行
      if(a[k][i]) {r=k; break;}
    if(r!=i) swap(a[r],a[i]); //换行
    if(!a[i][i]) return 0;  
    
    int x=quickpow(a[i][i],P-2); //求逆元
    for(int k=1; k<=n; ++k){ //对角化
      if(k == i) continue;
      int t=a[k][i]*x%P;
      for(int j=i; j<=2*n; ++j) 
        a[k][j]=((a[k][j]-t*a[i][j])%P+P)%P;
    } 
    for(int j=1; j<=2*n; ++j) //除以主元即乘以主元的逆
      a[i][j]=(a[i][j]*x%P);
  }
  return 1;
}
int main(){
  scanf("%d",&n);
  for(int i=1; i<=n; ++i)
    for(int j=1; j<=n; ++j)
      scanf("%lld",&a[i][j]),a[i][i+n]=1;
  if(Gauss_Jordan())
    for(int i=1; i<=n; ++i){
      for(int j=n+1; j<=2*n; ++j) 
        printf("%lld ",a[i][j]);
      puts("");
    }
  else puts("No Solution");
  return 0;
}

求组合数递推法杨辉三角

#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;

const int N = 2010, p = 1e9 + 7;
int c[N][N];
void init()
{
	for (int i = 0; i < N; i++)c[i][0] = 1;
	for (int i = 1; i < N; i++)
		for (int j = 1; j <= i; j++)
			c[i][j] = (c[i - 1][j] + c[i - 1][j - 1]) % p;
}
int main()
{
	init();
	int T, n, m;
	cin >> T;
	while (T--)
	{
		cin >> n >> m;
		cout << c[n][m];
	}
	return 0;
}

求组合数快速幂

#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;

typedef long long ll;
const long long p = 1e9 + 7;
const ll N = 1e5 + 10;
ll f[N], g[N];
int quickpow(int a, int b)
{
	ll res = 1;
	while (b)
	{
		if (b & 1)res *= a % p;
		a *= a % p;
		b >>= 1;
	}
	return res;
}

void init()
{
	f[0] = g[0] = 1;
	for (int i = 1; i < N; i++)
	{
		f[i] = f[i - 1] * i % p;
		g[i] = g[i - 1] * quickpow(i, p - 2) % p;
	}
}

ll getc(ll n, ll m)
{
	return f[n] * g[m] % p * g[n - m] % p;
}
int main()
{
	ll n, m, t;
	init();
	cin >> t;
	while (t--)
	{
		cin >> n >> m;
		cout << getc(n, m) << '\n';
	}
	return 0;
}

求组合数卢卡斯（Lucas)定理

Luogu P3807 【模板】卢卡斯定理/Lucas 定理

#include<bits/stdc++.h>
using namespace std;
long long n,m,p;
long long pow(long long a,long long k,long long p)
{
    long long ret=1;
    for(;k;k>>=1,a=(a*a)%p)
        if(k&1)ret=(ret*a)%p;
    return ret;
}
long long c(long long n,long long m,long long p)
{
    if(m>n)return 0;
    if(n>=p||m>=p)
        return (c(n/p,m/p,p)*c(n%p,m%p,p))%p;
    long long ans1=1,ans2=1;
    for(int i=m+1;i<=n;i++)ans1=(ans1*i)%p;
    for(int i=1;i<=n-m;i++)ans2=(ans2*i)%p;
    return (ans1*pow(ans2,p-2,p))%p;
}
//这个函数首先检查 m>n 的情况，如果是，则组合数为0。
如果 n≥p 或 m≥p，则使用 Lucas 定理递归地计算。
否则，直接计算 C(n,m)modp，其中 ans1 存储了分子部分（即 n×(n−1)×⋯×(m+1)），ans2 存储了分母部分（即 1×2×⋯×(n−m)），然后返回 (ans1 * pow(ans2, p-2, p)) % p，这里 pow(ans2, p-2, p) 计算了 ans2 模 p 的逆元。

int main()
{
    int t;
    scanf("%d",&t);
    for(int k=1;k<=t;k++)
    {
        scanf("%lld%lld%lld",&n,&m,&p);
        printf("%lld\n",c(n+m,m,p));
    }
}

#include<iostream>  
#include<algorithm>  
#include<cmath>  
using namespace std;
typedef long long ll;
const int N = 1e5 + 10;
int f[N], g[N];

ll quickpow(ll a, ll b, ll p) {
    ll res = 1;
    a %= p; // 确保 a 在模 p 的范围内  
    while (b) {
        if (b & 1) res = (res * a) % p;
        a = (a * a) % p;
        b >>= 1;
    }
    return res;
}

void init(ll p) {
    f[0] = g[0] = 1;
    for (int i = 1; i < N; i++) {
        f[i] = (f[i - 1] * i) % p;
        g[i] = (g[i - 1] * quickpow(i, p - 2, p)) % p;
    }
}

ll getc(ll n, ll m, ll p) {
    if (m > n)return 0;
    return (f[n] * g[m] % p * g[n-m] % p) % p;
}
ll  Lucas(ll n, ll m, ll p)
{
    if (m == 0)return 1;
    return Lucas(n / p, m / p, p) * getc(n % p, m % p, p) % p;
}

int main() {
    ll n, m, p, t;
    cin >> t;
    while (t--) {
        cin >> n >> m >> p;
        init(p); // 每次循环前根据新的 p 初始化  
        cout << Lucas(m+n,n,p) << '\n';
    }
    return 0;
}

求组合数高精度线性筛

#include<iostream>  
#include<algorithm>  
#include<cmath>  
using namespace std;
typedef long long ll;

const int N = 10010;
int prim[N], vis[N], cnt;

void get_prim(int n) {//筛素数
	for (int i = 2; i <= n; i++)
	{
		if (!vis[i])prim[++cnt] = i;
		for (int j = 0; i * prim[j] <= n; j++)
		{
			vis[i * prim[j]] = 1;
			if (i % prim[j])break;
		}
	}
}
int get(int n, int p) {//n的阶乘中p的个数
	int s = 0;
	while (n)s += n / p, n /= p;
	return s;
}
int getps(int n, int m, int p)//c中p的个数
{
	return get(n, p) - get(m, p) - get(n - m, p);
}
void mul(int c[], int p, int& len)
{
	int t = 0;
	for (int i = 0; i < len; i++)
	{
		t += c[i] * p;
		c[i] = t % 10;
		t /= 10;
	}
	while (t)
	{
		c[len++] = t % 10;
		t /= 10;
	}
}
int main()
{
	int n, m;
	cin >> n >> m;
	get_prim(n);
	int c[N], len = 1; c[0] = 1;
	for (int i = 0; i < cnt; i++)
	{
		int p = prim[i];
		int s = getps(n, m, p);
		while (s--)
		{
			mul(c, p, len);
		}
	}
	for (int i = len - 1; i > 0; i--)
		cout << c[i];
	return 0;
}

隔板法方程的解

Luogu P1771 方程的解

#include <iostream>
#include <cstring>
#include <algorithm>
using namespace std;

const int N = 150, P = 1000;
int C[1000][100][N];

int qpow(int a, int b){//快速幂
  int res = 1;
  while(b){
    if(b & 1) res = res*a%P;
    a = a*a%P;
    b >>= 1;
  }
  return res;
}
void add(int c[],int a[],int b[]){//高精度
  for(int i = 0; i < N; i++){
    c[i] += a[i]+b[i];
    c[i+1] += c[i]/10;
    c[i] %= 10;
  }
}
void getC(int n, int k){//组合数
  for(int i=0; i<n; i++)
    for(int j=0; j<=i && j<k; j++)
      if(j == 0) C[i][j][0] = 1;
      else add(C[i][j],C[i-1][j],C[i-1][j-1]);  
}
int main(){
  int k, x;
  cin >> k >> x;
  int n = qpow(x%P, x);
  getC(n, k);
  int i = N-1;
  while(C[n-1][k-1][i]==0) i--;
  while(i>=0) printf("%d", C[n-1][k-1][i--]);
  return 0;
}

容斥原理集合的并

#include <iostream>
#include <cstring>
#include <algorithm>
using namespace std;

typedef long long LL;
const int N = 20;
int n, m, prim[N];

int calc(){ //容斥原理
  int res = 0;
  for(int i=1; i<1<<m; i++){//枚举状态
    int t = 1, sign = -1;
    for(int j=0; j<m; j++) //过滤状态
      if(i & 1<<j){
        if((LL)t*prim[j] > n){
          t = 0; break;
        }
        t *= prim[j]; //质数的积
        sign = -sign;
      }
    if(t) res += n/t*sign; //交集的和     
  }
  return res;
}
int main(){
  cin >> n >> m;
  for(int i=0; i<m; i++) cin>>prim[i];
  cout << calc();
  return 0;
}

这段代码实现了一个基于容斥原理的算法，用于计算在给定范围内的整数中，不能被给定的质数集合中任何一个质数整除的数的个数。具体来说，它解决了以下问题：给定一个整数 n 和一个包含 m 个质数的数组 prim，求在 1 到 n（包含 n）之间，有多少个整数不能被 prim 数组中的任何一个质数整除。

代码详细分析

变量定义

LL：定义了长整型 long long 的别名，用于处理可能的大数。
N：定义了数组 prim 的最大可能长度，这里设为 20。
n, m：分别表示整数范围的上限和质数数组的长度。
prim[N]：存储质数的数组。

calc 函数

这个函数是实现容斥原理的核心部分。

初始化结果：res 初始化为 0，用于存储最终的答案。
枚举状态：通过循环 for(int i=1; i<1<<m; i++) 枚举所有可能的质数组合状态。这里使用了位运算技巧，1<<m 是 2^m，表示所有质数可能的组合数量（包括空集）。i 从 1 开始，因为 0 表示空集，这里不需要考虑。
过滤状态：对于每个状态 i，通过循环 for(int j=0; j<m; j++) 检查该状态对应的质数是否应该被包含在当前组合中。如果 i & (1<<j) 非零，表示第 j 个质数被选中。
- 检查积是否超过 n：如果当前质数的乘积 t 乘以第 j 个质数后超过 n，则无需继续计算，将 t 设为 0 并跳出循环。
- 更新乘积和符号：如果第 j 个质数被选中，则将其乘到 t 上，并改变符号 sign。sign 初始为 -1，每次选择质数时取反，以符合容斥原理的加减交替规则。
计算并累加结果：如果 t 不为 0（即当前组合有效），则计算该组合对应的贡献，即 n/t * sign。这里 n/t 表示当前组合（即能被这些质数整除的数）的个数，sign 决定了是加还是减这个贡献。
返回结果：最终返回 res，即所有有效组合贡献的累加结果。

main 函数

输入 n 和 m，然后输入 m 个质数到 prim 数组中。
调用 calc 函数计算答案，并输出结果。

注意事项

代码中使用了 long long 类型来避免整数溢出，特别是在计算质数乘积和 n/t 时。
使用了位运算来高效地枚举所有质数的组合状态。
需要注意的是，当 n 非常大或 m 和质数乘积较大时，这段代码的效率可能会受到影响，因为需要枚举所有可能的质数组合。然而，对于题目中给定的限制（N=20），这段代码是足够高效的。
容斥原理集合的交
Luogu P1450 [HAOI2008] 硬币购物

#include <iostream>
#include <cstring>
#include <algorithm>
using namespace std;

typedef long long LL;
int c[4],d[4],n,s;
LL f[100005];

void pack_pre(){ //完全背包预处理
    f[0] = 1;
    for(int i=0; i<4; i++)
      for(int j=c[i]; j<100005; j++)
        f[j] += f[j-c[i]];
}
LL calc(LL s){ //容斥原理
  LL res = 0;
  for(int i=1; i<1<<4; i++){//枚举状态
    LL t = 0, sign = -1;
    for(int j=0; j<4; j++) //过滤状态
      if(i & 1<<j){
        t += c[j]*(d[j]+1);
        sign = -sign;      
      }
    if(s>=t) res += f[s-t]*sign;
  }
  return f[s]-res;
}
int main(){
    for(int i=0; i<4; i++) scanf("%d",&c[i]);
    pack_pre(); 
    scanf("%d", &n);
    while(n--){
        for(int i=0; i<4; i++) scanf("%d",&d[i]);
        scanf("%d",&s);
        printf("%lld\n", calc(s));
    }
    return 0;
}

这段代码实现了一个基于完全背包和容斥原理的算法，用于解决一个特定的问题，通常这类问题涉及到在给定的容量限制下，使用有限数量的不同物品（每种物品有固定的重量和数量限制），计算能够组成特定总重量的方案数。不过，这里的实现方式略有不同，它首先通过完全背包预处理出所有可能的重量组合的方案数，然后通过容斥原理排除掉不符合特定数量限制的组合。

代码分析

全局变量定义：
- c[4]：存储四种物品的重量。
- d[4]：存储四种物品的最大使用数量。
- n：表示接下来要处理的查询数量。
- s：表示目标总重量。
- f[100005]：用于存储完全背包问题的解，即重量从0到100004时，能够组成的方案数。
pack_pre 函数：
- 这个函数通过完全背包的解法预处理出f数组。对于每种物品，它遍历所有可能的重量（从该物品的重量开始），并更新f数组，表示在加入当前物品后，能够组成的总重量及其对应的方案数。
calc 函数：
- 这个函数使用容斥原理来计算在给定的数量限制下，能够组成目标重量s的方案数。
- 它枚举所有可能的状态（即哪些物品被考虑在内），对于每种状态，计算如果仅使用这些物品能够组成的最大重量t，并根据容斥原理的奇偶性（通过sign变量控制）来累加或减去对应的方案数。
- 最后，从f[s]（即不考虑数量限制时，能够组成重量s的方案数）中减去通过容斥原理计算出的不符合数量限制的方案数，得到最终结果。
main 函数：
- 首先读取四种物品的重量。
- 调用pack_pre函数进行完全背包的预处理。
- 读取查询数量n，并对于每个查询，读取四种物品的数量限制和目标重量s，然后调用calc函数计算并输出结果。

卡特兰数（栈）

Luogu P1044 [NOIP2003 普及组] 栈

//f[n]=f[n-1]*(4*n-2)/(n+1)
#include<iostream>
using namespace std;
typedef long long ll;
int n;
ll f[20];
int main()
{
	cin >> n;
	f[0] = 1;
	for (int i = 1; i <= n; i++)
		f[i] = f[i - 1] * (4 * i - 2) / (i + 1);
	cout << f[n];
	return 0;
}

整除分块（数论分块）

Luogu P2261 [CQOI2007]余数求和

#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;
typedef long long ll;

int main()
{
	ll n, k, l, r, res;
	cin >> n >> k;
	res = n * k;
	for (l = 1; l <= n; l = r + 1)
	{
		if (k / l == 0)break;
		r = min(k / (k / l), n);
		res -= (k / l) * (r - l + 1) * (l + r) / 2;
	}
	cout << res;
	return 0;
}
//

普通生成函数

HDU - 2152 Fruit

#include <iostream>
#include <cstring>
#include <algorithm>
using namespace std;

int n,m;
int a[110],b[110];//存幂次
int C[110],D[210];//存系数

int calc(){
  for(int i=0;i<=m;++i) C[i]=D[i]=0;
  //填充第1项的系数
  for(int i=a[1];i<=b[1];++i) C[i]=1;
  //从第2项开始枚举
  for(int i=2;i<=n;++i){
    //计算x^(j+k)的系数
    for(int j=0;j<=m;++j)
      for(int k=a[i];k<=b[i];++k)
        D[j+k]+=C[j];
    //转存C，清空D
    for(int j=0;j<=m;++j)
      C[j]=D[j], D[j]=0;
  }
  return C[m];
}

int main(){
  while(~scanf("%d%d",&n,&m)){
    for(int i=1; i<=n; ++i)
      scanf("%d%d",&a[i],&b[i]);
    printf("%d\n",calc());
  }
  return 0;
}

HDU - 1085 Holding Bin-Laden Captive!

#include <iostream>
#include <cstring>
#include <algorithm>
using namespace std;

int n,a[4],b[4],C[8005],D[8005];

void calc(int m){
  for(int i=0; i<=m; ++i) C[i]=D[i]=0;
  for(int i=0; i<=a[1]; ++i) C[i]=1;
  for(int i=2; i<=3; ++i){
    for(int j=0; j<=m; ++j)
      for(int k=0; k<=a[i]*b[i]&&j+k<=m; k+=b[i])
        D[j+k]+=C[j];
    for(int j=0; j<=m; ++j)
      C[j]=D[j],D[j]=0;
  }
}
int main(){
  while(scanf("%d%d%d",&a[1],&a[2],&a[3])&&(a[1]||a[2]||a[3])){
    b[2]=2,b[3]=5;
    int m=a[1]*1+a[2]*2+a[3]*5;
    calc(m);
    int x=1;
    while(x<=m && C[x]) x++;
    printf("%d\n",x);
  }
  return 0;
}

这段代码的主要思路是解决一个特定类型的组合问题，具体来说，它涉及到了使用有限数量的几种物品（每种物品有不同的重量和价值）来组成不超过某个总重量的所有可能方式中，能够组成的最小不能达到的重量。这里的物品及其限制条件是以一种特定的方式给出的：

输入：代码首先通过scanf函数从标准输入读取三种物品的数量（a[1], a[2], a[3]），这里a[i]表示第i种物品的数量。注意，代码中a数组的索引是从1开始的，这在C++中并不常见，但在这个特定上下文中是允许的（尽管可能会导致一些混淆）。此外，对于第二和第三种物品，它们的“价值”或“重量”增量是固定的（b[2]=2和b[3]=5），而第一种物品被视为每个物品的重量为1（这在代码中未显式声明，但从上下文可以推断出来）。
计算总重量限制：然后，代码计算了一个总重量限制m，这是通过将每种物品的数量乘以它们各自的重量（或价值）增量得到的。
动态规划解决方案：接下来，calc函数使用动态规划的思想来计算一个数组C（在过程中使用了辅助数组D来避免覆盖C中的旧值），其中C[i]表示能够组成重量i的方案数。注意，这里的“方案数”实际上是通过标记哪些重量是可达的（用非零值表示）来隐含的。
- 初始化时，将C[0]设为1，表示重量0总是可达的（不需要任何物品）。
- 然后，对于每种物品，代码遍历所有可能的重量，并更新C数组，以反映加入当前物品后哪些新的重量变得可达。
- 注意，对于第二和第三种物品，由于它们的重量增量不是1，代码需要以一种特殊的方式遍历可能的重量增加量（k+=b[i]）。
查找最小不可达重量：最后，main函数中的循环查找最小的x（从1开始），使得C[x]为0，即重量x是不可达的。这个x就是答案，表示使用给定数量的物品无法组成的最小重量。
输出：对于每组输入，代码输出这个最小不可达重量。

注意：代码中的C和D数组大小被设置为8005，这是基于总重量限制m的最大可能值（假设每种物品的数量都足够大）来选择的。然而，这个大小可能是一个保守的估计，具体取决于输入数据的范围。如果输入数据保证总重量限制远小于8005，那么可以减小数组的大小以节省内存。

此外，代码中a数组的索引从1开始，这可能会导致一些混淆，特别是对于那些习惯于从0开始索引的C++程序员来说。在更复杂或更长的代码中，最好坚持使用从0开始的索引，以提高代码的可读性和可维护性。

指数生成函数

HDU - 1521 排列组合（指数生成函数）

生成函数的应用

POJ 3734 Blocks

#include <cstdio>
using namespace std;

const int P=10007;
int qpow(int a, int b){
  int res=1;
  while(b){
    if(b & 1) res=res*a%P;
    a=a*a%P;
    b>>=1;
  }
  return res;
}
int main(){
  int t,n,ans;
  scanf("%d", &t);
  while(t--){
    scanf("%d", &n);
    ans=(qpow(4,n-1)+qpow(2,n-1))%P;
    printf("%d\n", ans);
  }
  return 0;
}

狄利克雷卷积

和式的变换

G38 和式的变换_哔哩哔哩_bilibili

Luogu P3455 [POI2007]ZAP-Queries

#include <iostream>  
#include <cstring>  
#include <algorithm>  
using namespace std;

typedef long long ll;
const int N = 50005;
int vis[N], prime[N], mu[N], cnt;

void init() {
    mu[1] = 1;
    for (int i = 2; i < N; i++) {
        if (!vis[i]) {
            prime[++cnt] = i;
            mu[i] = -1;
        }
        for (int j = 1; j <= cnt && i * prime[j] < N; j++) {
            vis[i * prime[j]] = 1;
            if (i % prime[j] == 0) {
                mu[i * prime[j]] = 0;
                break;
            }
            mu[i * prime[j]] = -mu[i];
        }
    }
    for (int i = 1; i < N; i++) mu[i] += mu[i - 1];
}

ll calc(int n, int m, int k) {
    if (n > m) swap(n, m);
    n /= k;
    m /= k;
    ll ans = 0;
    for (int l = 1, r; l <= min(n, m); l = r + 1) {
        r = min(n / (n / l), m / (m / l));
        ans += (ll)(mu[r] - mu[l - 1]) * (n / l) * (m / l);
    }
    return ans;
}

int main() {
    init();
    int t, n, m, k;
    cin >> t;
    while (t--) {
        cin >> n >> m >> k;
        cout << calc(n, m, k) << endl;  // 添加 endl 以确保输出换行  
    }
    return 0;
}

该代码实现了一个高效的算法来计算在给定范围内（由n和m定义），且每个数都是k的倍数的所有数对中，互质的数对的数量。这里使用了莫比乌斯反演（Möbius Inversion）和数论分块（也称为整除分块或除法分块）技术来优化计算。下面是代码的详细思路与分析：

1. 初始化（`init`函数）

莫比乌斯函数（mu）：mu[i]定义为如果i是1个或多个不同质数的乘积，则mu[i]为-1；如果i含有某个质数的平方因子，则mu[i]为0；特别地，mu[1] = 1。这个函数在计算互质数对时非常有用。
质数筛选：通过埃拉托斯特尼筛法（Sieve of Eratosthenes）的变种来筛选质数，并将mu数组初始化为对应的值。
前缀和：计算mu数组的前缀和，以便在后续的分块计算中快速求和。

2. 计算函数（`calc`函数）

参数处理：首先处理n和m，将它们都除以k，因为题目要求的是k的倍数的数对，所以这一步是为了将问题规模缩小到k的倍数范围内。
分块计算：利用数论分块技术，将[1, min(n, m)]范围内的数分成多个块，每个块内的数除以n或m的结果相同，从而减少计算量。
莫比乌斯函数的应用：对于每个块[l, r]，利用mu数组的前缀和快速计算[l, r]范围内所有数的莫比乌斯函数值的和，然后乘以(n/l) * (m/l)，即这个块内所有可能的数对组合数。这样做是因为在数论中，gcd(i, j) = 1的数对(i, j)的数量可以通过莫比乌斯函数来表示和计算。

3. 主函数（`main`）

输入处理：读取测试用例数量t，然后对于每个测试用例，读取n、m和k的值。
调用计算函数：对于每组n、m和k，调用calc函数计算并输出互质的数对的数量。

复杂度分析

初始化：质数筛选和莫比乌斯函数的计算都是O(N log log N)的，其中N是mu数组的大小。
计算函数：对于每个查询，数论分块技术将查询分解为O(sqrt(min(n, m)))个块，每个块内的计算是常数时间的（通过前缀和快速计算），因此总的时间复杂度是O(sqrt(min(n, m)))。

Luogu P2257 YY的GCD

#include <iostream>
#include <cstdio>
#include <cstring>
using namespace std;
#define LL long long
const int N = 10000010;
int vis[N],p[N],mu[N],cnt;
int F[N];

void init(){
  mu[1]=1;
  for(int i=2;i<N;i++){
    if(!vis[i]) p[++cnt]=i,mu[i]=-1;
    for(int j=1;i*p[j]<N;j++){
      vis[i*p[j]]=1;
      if(i%p[j]==0) break;
      mu[i*p[j]]=-mu[i];
    }
  }
  for(int i=1;i<=cnt;i++)
    for(int j=p[i];j<N;j+=p[i])
      F[j]+=mu[j/p[i]];
  for(int i=1;i<N;i++) F[i]+=F[i-1];
}
LL calc(int n,int m){
  if(n>m) swap(n,m);  
  LL ans=0;
  for(int l=1,r;l<=n;l=r+1){
    r=min(n/(n/l),m/(m/l));
    ans+=1ll*(F[r]-F[l-1])*(n/l)*(m/l);
  }
  return ans;
}
int main(){
  init();
  int n,m,T; 
  scanf("%d",&T);
  while(T--){
    scanf("%d%d",&n,&m);
    printf("%lld\n",calc(n,m));
  }
}

这段代码实现了一个高效的算法来计算在给定范围内（由n和m定义）内所有数对中，满足gcd(x, y) = 1（即x和y互质）的数对的数量。这里使用了莫比乌斯反演（Möbius Inversion）和数论分块（也称为整除分块或除法分块）技术来优化计算。下面是代码的详细思路分析：

1. 初始化（`init`函数）

莫比乌斯函数（mu）：mu[i]是莫比乌斯函数，定义为如果i是1个或多个不同质数的乘积，则mu[i] = -1；如果i含有某个质数的平方因子，则mu[i] = 0；特别地，mu[1] = 1。这个函数在计算互质数对时非常有用。
质数筛选：通过埃拉托斯特尼筛法（Sieve of Eratosthenes）的变种来筛选质数，并同时计算莫比乌斯函数mu的值。
计算F数组：F[j]表示所有小于等于j的数中，其质因数分解后只含有一个质因子的数的莫比乌斯函数值之和。这个数组是为了后续计算互质数对数量时能够快速累加。具体来说，对于每个质数p[i]，遍历其所有倍数j，并将mu[j/p[i]]累加到F[j]上。这样，F[j]就包含了所有小于等于j的、只含有一个质因子的数的莫比乌斯函数值之和。
前缀和：最后，计算F数组的前缀和，以便在后续的分块计算中快速求和。

2. 计算函数（`calc`函数）

参数处理：首先处理n和m，确保n <= m，因为互质数对的计算是对称的，且结果只与n和m的相对大小有关。
数论分块：利用数论分块技术，将[1, min(n, m)]范围内的数分成多个块，每个块内的数除以n或m的结果相同，从而减少计算量。具体来说，对于每个块[l, r]，计算(n/l) * (m/l)，即这个块内所有可能的数对组合数，并乘以F[r] - F[l-1]，即这个块内所有满足条件的数对的莫比乌斯函数值之和（注意这里利用了莫比乌斯反演的性质）。
累加结果：将所有块的结果累加，得到最终的答案。

3. 主函数（`main`）

读取测试用例：首先读取测试用例的数量T，然后对于每个测试用例，读取n和m的值。
调用计算函数：对于每组n和m，调用calc函数计算并输出互质的数对的数量。

Luogu P3327 [SDOI2015]约数个数和

#include <iostream>
#include <cstdio>
#include <cstring>
using namespace std;
#define LL long long

const int N = 50010;
int vis[N],p[N],mu[N],cnt;
LL F[N];

void init(){
  mu[1]=1;
  for(int i=2;i<N;i++){
    if(!vis[i]) p[++cnt]=i,mu[i]=-1;
    for(int j=1;i*p[j]<N;j++){
      vis[i*p[j]]=1;
      if(i%p[j]==0) break;
      mu[i*p[j]]=-mu[i];
    }
  }
  for(int i=1;i<N;i++) mu[i]+=mu[i-1];
  for(int i=1;i<N;i++)
    for(int l=1,r;l<=i;l=r+1){
      r=(i/(i/l));
      F[i]+=1ll*(r-l+1)*(i/l);
    }
}
LL calc(int n,int m){
  if(n>m) swap(n,m);  
  LL ans=0;
  for(int l=1,r;l<=n;l=r+1){
    r=min(n/(n/l),m/(m/l));
    ans+=(mu[r]-mu[l-1])*F[n/l]*F[m/l];
  }
  return ans;
}
int main(){
  init();
  int n,m,T; 
  scanf("%d",&T);
  while(T--){
    scanf("%d%d",&n,&m);
    printf("%lld\n",calc(n,m));
  }
}

初始化（init函数）：
- 莫比乌斯函数（mu）：mu[i]是莫比乌斯函数，定义为如果i是k个不同质数的乘积，则mu[i] = (-1)^k；如果i含有平方质因子，则mu[i] = 0。特别地，mu[1] = 1。这个数组用于后续计算互质数对。
- 质数筛选：通过埃拉托斯特尼筛法（Sieve of Eratosthenes）的变种来筛选质数，并同时计算mu数组的值。
- 前缀和：计算mu数组的前缀和，以便在后续计算中快速获取区间和。
- 计算F数组：F[i]数组用于存储1到i中所有数的约数个数之和（更具体地说，是Σ(i/j)，其中j是i的约数）。这个数组通过数论分块技术预先计算好，以便在后续计算中快速使用。
计算函数（calc函数）：
- 首先处理n和m，确保n <= m，因为结果是对称的。
- 使用数论分块技术来减少计算量。对于每个块[l, r]，其中l是块的起始位置，r是块的结束位置，且在该块内所有数除以n（或m）的结果相同，从而可以统一计算。
- 对于每个块，计算(mu[r] - mu[l-1])，这表示在区间[l, r]内所有数的莫比乌斯函数值之和。
- 然后计算F[n/l] * F[m/l]，这表示在[1, n/l]和[1, m/l]范围内，分别选取一个数作为a和b的约数个数乘积，这个乘积与gcd(a, b) = 1的条件通过莫比乌斯函数相关联。
- 将上述两部分相乘并累加到答案ans中。
主函数（main）：
- 读取测试用例数量T，然后对于每个测试用例，读取n和m的值。
- 调用calc函数计算并输出互质的数对的数量。

莫比乌斯反演

Luogu P1829 [国家集训队]Crash的数字表格

#include <algorithm>
#include <cstdio>
using namespace std;
const int N=10000010;
const int P=20101009;
int vis[N],p[N],mu[N],S[N],cnt;

void init(){
  mu[1]=1;
  for(int i=2; i<N; ++i){
    if(!vis[i])p[++cnt]=i,mu[i]=-1;
    for(int j=1; i*p[j]<N; ++j){
      vis[i*p[j]]=1;
      if(i%p[j] == 0) break;
      mu[i*p[j]]=-mu[i];
    }
  }
  for(int i=1;i<N;++i)
    S[i]=(S[i-1]+1LL*mu[i]*i*i%P+P)%P;
}
int G(int n, int m){
  return (1LL*n*(n+1)/2%P)*(1LL*m*(m+1)/2%P)%P;
}
int F(int n, int m){
  int res=0;
  for(int l=1,r; l<=n; l=r+1){
    r=min(n/(n/l),m/(m/l));
    res=(res+1LL*(S[r]-S[l-1])*G(n/l,m/l)%P+P)%P;
  }
  return res;
}
int calc(int n, int m){
  if(n>m) swap(n,m);
  int res=0;
  for(int l=1,r; l<=n; l=r+1){
    r=min(n/(n/l),m/(m/l));
    res=(res+1LL*(r-l+1)*(l+r)/2%P*F(n/l,m/l)%P)%P;
  }
  return res;
}
int main(){
  init();  
  int n,m;
  scanf("%d%d", &n, &m);
  printf("%d\n", calc(n,m));
}

Luogu P3704 [SDOI2017]数字表格

#include<iostream>
#include<cstdio>
#include<cstdlib>
using namespace std;
#define P 1000000007
#define N 1000010
int mu[N],p[N],vis[N],cnt;
int f[N],g[N],F[N];

int qpow(int a,int b){
  int res=1;
  while(b){
    if(b&1) res=1ll*res*a%P;
    a=1ll*a*a%P;
    b>>=1;
  }
  return res;
}
void init(){
  mu[1]=1;
  for(int i=2;i<N;++i){
    if(!vis[i]) p[++cnt]=i, mu[i]=-1;
    for(int j=1;i*p[j]<N;++j){
      vis[i*p[j]]=1; 
      if(i%p[j]==0) break; mu[i*p[j]]=-mu[i];
    }
  }
  f[1]=g[1]=F[0]=F[1]=1;
  for(int i=2;i<N;++i){
    f[i]=(f[i-1]+f[i-2])%P; g[i]=qpow(f[i],P-2); F[i]=1;    
  }
  for(int i=1;i<N;++i)
    for(int j=i;j<N;j+=i)
      if(mu[j/i]) F[j]=1ll*F[j]*(mu[j/i]==1?f[i]:g[i])%P;
  for(int i=2;i<N;++i) F[i]=1ll*F[i]*F[i-1]%P; //前缀积
}
int calc(int n, int m){
  if(n>m) swap(n,m); int r,s,ans=1;
  for(int l=1;l<=n;l=r+1){
    r=min(n/(n/l),m/(m/l));
    s=1ll*F[r]*qpow(F[l-1],P-2)%P; //区间积
    ans=1ll*ans*qpow(s,1ll*(n/l)*(m/l)%(P-1))%P;
  } 
  return ans;
}
int main(){
  init();
  int T,n,m;
  scanf("%d",&T);
  while(T--){
    scanf("%d%d",&n,&m);
    printf("%d\n",calc(n,m));
  }
}

杜教筛

Luogu P4213 【模板】杜教筛（Sum）

#include <algorithm>
#include <cstdio>
#include <cstring>
#include <map>
using namespace std;
#define LL long long
const int N=2000010;
LL vis[N],pm[N],mu[N],phi[N],cnt;
map<LL,LL> mp_mu, mp_phi;

void init(){
  mu[1]=phi[1]=1;
  for(int i=2; i<N; i++){
    if(!vis[i])pm[++cnt]=i,mu[i]=-1,phi[i]=i-1;
    for(int j=1;i*pm[j]<N;j++){
      int p=pm[j]; vis[i*p]=1;
      if(i%p==0){phi[i*p]=phi[i]*p; break;}
      mu[i*p]=-mu[i]; phi[i*p]=phi[i]*(p-1);
    }
  }
  for(int i=1;i<N;i++)mu[i]+=mu[i-1],phi[i]+=phi[i-1]; 
}
LL Sphi(LL n){
  if(n<N) return phi[n]; //预处理的剪枝
  if(mp_phi[n]) return mp_phi[n]; //记忆化剪枝
  LL ans = n*(n+1)/2;
  for(LL l=2,r; l<=n; l=r+1){
    r=n/(n/l);
    ans-=Sphi(n/l)*(r-l+1);
  }
  return mp_phi[n]=ans;
}
LL Smu(LL n){
  if(n<N) return mu[n];
  if(mp_mu[n]) return mp_mu[n];
  LL ans = 1;
  for(LL l=2,r; l<=n; l=r+1){
    r=n/(n/l);
    ans-=Smu(n/l)*(r-l+1);
  }
  return mp_mu[n]=ans;
}
int main(){
  init();
  LL T, n; scanf("%lld", &T);
  while(T--){
    scanf("%lld", &n);
    printf("%lld %lld\n", Sphi(n), Smu(n));
  }
}

Luogu P3768 简单的数学题

#include<bits/stdc++.h>
#define N 5000000
#define LL long long
using namespace std;

int vis[N],p[N],phi[N],cnt;
LL P,n,sum[N],inv;
map<LL,LL> mp;

LL qpow(LL a,LL b){
  LL res=1;
  while(b){
    if(b&1) res=res*a%P;
    a=a*a%P;
    b>>=1;
  }
  return res;
}
void init(){
  inv=qpow(6,P-2);
  phi[1]=1;
  for(int i=2;i<N;i++){
    if(!vis[i])p[++cnt]=i,phi[i]=i-1;
    for(int j=1; i*p[j]<N; j++){
      vis[i*p[j]] = 1;
      if(i%p[j]==0){phi[i*p[j]] = p[j]*phi[i];break;}
      phi[i*p[j]]=(p[j]-1)*phi[i];
    }
  }
  for(int i=1;i<N;i++)
    sum[i]=(sum[i-1]+1LL*i*i%P*phi[i]%P)%P;
}
LL S2(LL x){x%=P;return x*(x+1)%P*(2*x+1)%P*inv%P;}
LL S3(LL x){x%=P;return(x*(x+1)/2)%P*((x*(x+1)/2)%P)%P;}
LL Sn(LL x){
  if(x<N) return sum[x];
  if(mp[x]) return mp[x];
  LL res=S3(x);
  for(LL l=2,r;l<=x;l=r+1){
    r=x/(x/l);
    res=(res-(S2(r)-S2(l-1))%P*Sn(x/l)%P+P)%P;
  }
  return mp[x]=res;
}
LL calc(){
  LL ans=0;
  for(LL l=1,r;l<=n;l=r+1){
    r=n/(n/l);
    ans=(ans+(Sn(r)-Sn(l-1))%P*S3(n/l)%P+P)%P;
  }
  return ans;
}
int main(){
  scanf("%lld%lld",&P,&n);
  init();
  printf("%lld\n",calc());
  return 0;
}