数学专题学习笔记

开坑于-2022.1.24

update1-2022.1.26

update2-2022.2.4

update3-2022.2.14

update4-2022.2.16

这是一个从约数开始的数学专题，至于能写多长我也不知道，有时间就学一点。

约数#

约数应该没有什么好有疑问的地方，无非就是求一个数的约数，或者将题目转换成与约数有关的内容。

感觉大多数时是被其他算法利用的一个工具。

一些例题#

P2926 [USACO08DEC]Patting Heads

一道比较简单的题目，考虑将问题的答案全部预处理出。

由于 $\text{n}$ 只有 $1e6$ ，所以直接统计每个数的出现次数，然后直接枚举，预处理所有的答案即可。

复杂度为 $O(n\sqrt{n})$。

点击查看代码

#include <bits/stdc++.h>
using namespace std;

int n , a[100010] , cnt[1000010] , ans[1000010];

inline int read()
{
    int asd = 0 , qwe = 1; char zxc;
    while(!isdigit(zxc = getchar())) if(zxc == '-') qwe = -1;
    while(isdigit(zxc)) asd = asd * 10 + zxc - '0' , zxc = getchar();
    return asd * qwe;
}

int main()
{
    n = read();
    for(int i = 1;i <= n;i++)
        a[i] = read() , cnt[a[i]]++;
    for(int i = 1;i <= 1e6 + 10;i++)
    {
        for(int j = i;j <= 1e6 + 10;j += i)
        {
            ans[j] += cnt[i];
        }
    }
    for(int i = 1;i <= n;i++)
        cout << ans[a[i]] - 1 << endl;
    return 0;
}

P1445 [Violet]樱花

根据题意，我们发现直接求答案十分麻烦。

考虑将式子转换一下。

$$\frac{1}{x}+\frac{1}{y}=\frac{1}{n!}$$

$$\frac{x+y}{xy}=\frac{1}{n!}$$

$$(x+y)n!=xy$$

$$xn!=(x-n!)y$$

$$y=\frac{xn!}{x-n!}$$

由于 $\text{y}$ 为正整数。

所以

$$\frac{xn!}{x-n!}>0$$

$$\frac{xn!}{x-n!}=\frac{xn!-n{!}^2+n{!}^2}{x-n!}=\frac{(x-n!)n!+n{!}^2}{x-n!}=n!+\frac{n{!}^2}{x-n!}$$

又因为 $\text{x}$ ， $\text{y}$ 均为正整数。

所以 $x>n!$。

所以问题就变成了求 $n{!}^2$ 的约数个数。

故直接求就可以了。

点击查看代码

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

const int mod = 1e9 + 7;
const int maxn = 1e6 + 10;

int n , ans , cnt;
int prime[maxn] , isprime[maxn];

inline int read()
{
    int asd = 0 , qwe = 1; char zxc;
    while(!isdigit(zxc = getchar())) if(zxc == '-') qwe = -1;
    while(isdigit(zxc)) asd = asd * 10 + zxc - '0' , zxc = getchar();
    return asd * qwe;
}

inline void init()
{
    for(int i = 2;i <= n;i++)
    {
        if(!isprime[i]) prime[cnt++] = i;
        for(int j = 0;prime[j] * i <= n;j++)
        {
            isprime[prime[j] * i] = 1;
            if(i % prime[j] == 0) break;
        }
        
    }
}

signed main()
{
    n = read() , ans = 1;
    init();
    for(int i = 0;i < cnt;i++)
    {
        int x = prime[i] , y = 0;
        for(int j = n;j;j /= x) y += j / x;
        ans = ans * (2 * y + 1) % mod;
    }
    cout << ans << endl;
    return 0;
}

这里用了一个小技巧。

即求大数阶乘的约数。

我们可以将其质因数分解，这样如何求约数就解决了。

至于如何将其质因数分解，我们很容易就能知道，在 $\text{n}$ 个数里，有 $\lfloor \frac{n}{k}\rfloor$ 个是 $\text{k}$ 的倍数。

那么就可以质因数分解了。

P1463 [POI2001][HAOI2007]反素数

由于 $n\le 2\times {10}^9$ 。

所以需要考虑的质数只有九个，而最大的指数只有三十。

我们发现就可以直接暴搜求出答案。

点击查看代码

#include <bits/stdc++.h>
using namespace std;

int n , ans , num;
int prime[10] = {0 , 2 , 3 , 5 , 7 , 11 , 13 , 17 , 19 , 23};

inline int read()
{
    int asd = 0 , qwe = 1; char zxc;
    while(!isdigit(zxc = getchar())) if(zxc == '-') qwe = -1;
    while(isdigit(zxc)) asd = asd * 10 + zxc - '0' , zxc = getchar();
    return asd * qwe;
}

inline void dfs(int u , int last , int sum , int s)
{
    if(s > num || (sum < ans && num == s))
    {
        num = s;
        ans = sum;
    }
    
    if(u == 10) return;
    
    for(int i = 1;i <= last;i++)
    {
        if((long long)sum * prime[u] > n) break;
        sum *= prime[u];
        dfs(u + 1 , i , sum , s * (i + 1));
    }
}

int main()
{
    n = read();
    dfs(1 , 30 , 1 , 1);
    cout << ans;
    return 0;
}

欧拉函数#

欧拉函数是小于等于 $\text{x}$ 的数中与 $\text{x}$ 互质的数的数目。

$$\phi (x)=\prod _{i=1}^n(1-\frac{1}{a_i})$$

其中 $a_i$ 为 $\text{x}$ 的质因子，$\text{n}$ 为 $\text{x}$ 的因子个数。

几个性质#

• 若 $\text{x}$ 为质数，显然 $\phi (x)=x-1$

• 若 $m,n$ 互质，则 $\phi (mn)=\phi (m)\times \phi (n)$

• 若 $x=2n,n\mathrm{mod}2=1$，则 $\phi (x)=\phi (n)$

• 若 $\text{x}$ 为质数，则 $\phi (x^k)=x^k-x^{k-1}$

• 若 $x>2$，则 $\phi (x)\mathrm{mod}2=0$

• 小于 $n$ 且与 $n$ 互质的数的总和为 $\phi (n)\times n÷2(n>1)$

• $n=\sum _{d|n}\phi (d)$，即 $\text{n}$ 包括 $\text{n}$ 在内的所有因数的欧拉函数之和等于 $\text{n}$，这一条也被称为欧拉反演。

求法#

有两种求法。

1. $o(\sqrt{n})$ 的时间复杂度求出单个欧拉函数。

点击查看代码

inline int eular(int x)
{
    int res = x;
    for(int i = 2;i * i <= x;i++)
        if(x % i == 0)
        {
            res = res / i * (i - 1);
            while(x % i == 0) x /= i;
        }
    if(x > 1) res = res / x * (x - 1);
    return res;
}

2. $O(n\log (n))$ 的时间复杂度求出$1-n$欧拉函数。

点击查看代码

inline int eular(int x)
{
    for(int i = 1;i <= x;i++) e[i] = i;
    for(int i = 2;i <= x;i++)
        if(e[i] == i)
            for(int j = i;j < x;j++)
                e[j] = e[j] / i * (i - 1);
}

一些例题#

有空了再写题。

高斯消元#

高斯消元主要用于求解线性方程组，方程组主要形式为：

$$\{\begin{array}{l}{a}_{1,1}{x}_1+a_{2,1}{x}_2+...a_{n,1}{x}_n=b_1\\ a_{1,2}{x}_1+a_{2,2}{x}_2+...a_{n,2}{x}_n=b_2\\ ...\\ ...\\ a_{1,n}{x}_1+a_{2,n}{x}_2+...a_{n,n}{x}_n=b_n\end{array}$$

如何实现#

P3389 【模板】高斯消元法

这里不介绍朴素的高斯消元了。

介绍一下高斯——约旦消元。

我们以样例为例，模拟一下消元步骤。

首先，列出增广矩阵。

$$\left[\begin{array}{cccc}1& 3& 4& 5\\ 1& 4& 7& 3\\ 9& 3& 2& 2\end{array}\right]$$

将此时遍历到的一列中最大的数设为主元。

$$\left[\begin{array}{cccc}9& 3& 2& 2\\ 1& 3& 4& 5\\ 1& 4& 7& 3\end{array}\right]$$

然后对其他所有的行，进行消元。

不断循环上步骤，最终将矩阵变成下矩阵。

$$\left[\begin{array}{cccc}9& 0& 0& -8.76\\ 0& 3.67& 0& 19.01\\ 0& 0& -1.15& 2.76\end{array}\right]$$

格式化的结果

$$\left[\begin{array}{cccc}{a}_1& 0& 0& b_1\\ 0& a_2& 0& b_2\\ 0& 0& a_3& b_3\end{array}\right]$$

最终，输出答案，答案为 $\frac{b_i}{a_i}$

点击查看代码

#include <bits/stdc++.h>
using namespace std;

int n;
double a[110][110];

inline int read()
{
    int s = 0 , f = 1; char ch;
    while(!isdigit(ch = getchar())) if(ch == '-') f = -1;
    while(isdigit(ch)) s = s * 10 + ch - '0' , ch = getchar();
    return s * f;
}

int main()
{
    n = read();
    for(int i = 1;i <= n;i++)
    {
        for(int j = 1;j <= n + 1;j++) 
            a[i][j] = read();
    }
    for(int i = 1;i <= n;i++)
    {
        int m = i;
        for(int j = i + 1;j <= n;j++)
        {
            if(fabs(a[j][i]) > fabs(a[m][i]))
                m = j;
        }
        for(int j = 1;j <= n + 1;j++) swap(a[i][j] , a[m][j]);
        if(int(a[i][i]) == 0)
        {
            puts("No Solution");
            return 0;
        }
        for(int j = 1;j <= n;j++)
        {
            if(j != i)
            {
                double d = a[j][i] / a[i][i];
                for(int k = i + 1;k <= n + 1;k++)
                {
                    a[j][k] -= a[i][k] * d;
                }
            }
        }
    }
    for(int i = 1;i <= n;i++)
        printf("%.2lf\n" , a[i][i]);
    return 0;
}

一些例题#

P2455 [SDOI2006]线性方程组

一道十分套路的题。

特判一下无解与无穷解。

特判方法可以根据方程的特性，若是无穷解，则出现零是，最右边必然全是零，否则无解。

其余的就是板子了。

点击查看代码

#include <bits/stdc++.h>
using namespace std;

int n , cnt;
double a[100][100];

inline int read()
{
    int s = 0 , f = 1; char ch;
    while(!isdigit(ch = getchar())) if(ch == '-') f = -1;
    while(isdigit(ch)) s = s * 10 + ch - '0' , ch = getchar();
    return s * f;
}

int main()
{
    n = read() , cnt = 1;
    for(int i = 1;i <= n;i++)
        for(int j = 1;j <= n + 1;j++)
            a[i][j] = read();
    for(int i = 1;i <= n;i++)
    {
        int m = cnt;
        for(int j = cnt + 1;j <= n;j++)
            if(fabs(a[j][i]) > fabs(a[m][i])) m = j;
        if(a[m][i] == 0) continue;   
        for(int j = 1;j <= n + 1;j++)
            swap(a[cnt][j] , a[m][j]);
        for(int j = 1;j <= n;j++)
            if(cnt != j)
            {
                double w = a[j][i] / a[cnt][i];
                for(int k = i + 1;k <= n + 1;k++)
                    a[j][k] -= a[cnt][k] * w;
            }
        ++cnt;
    }
    if(cnt <= n)
    {
        while(cnt <= n)
        {
            // cout << a[cnt][n] << endl;
            if(a[cnt++][n + 1] != 0) 
            {
                puts("-1");
                return 0;
            }
        }
            
        puts("0");
        return 0;
    }
    for(int i = 1;i <= n;i++)
        printf("x%d=%.2lf\n" , i , (int(a[i][n + 1] / a[i][i] * 100) == 0 ? 0 : a[i][n + 1] / a[i][i]));
    return 0;
}

P4035 [JSOI2008]球形空间产生器

这道题主要就是列出线性方程，其余的就可以直接套板子。

由于一个球体的球心到球面的任意一个点距离相等，因此，我们需要求的 $(x_1,x_2,x_3,x_4...x_n)$ 需要满足：

$$\sum _{j=1}^n(a_{i,j}-x_j{)}^2=C$$

我们可以相邻两个作差。

$$\sum _{j=1}^{n}2(a_{i,j}-a_{i+1,j})x_j=\sum _{j=1}^n(a_{i,j}^2-a_{i+1,j}^2)(i=1,2,3...n)$$

化成矩阵即：

$$\left[\begin{array}{ccccc}2(a_{1,1}-a_{2,1})& 2(a_{1,2}-a_{2,2})& ...& 2(a_{1,n}-a_{2,n})& \sum _{j=1}^n(a_{1,j}^2-a_{2,j}^2)\\ 2(a_{2,1}-a_{3,1})& 2(a_{2,2}-a_{3,2})& ...& 2(a_{2,n}-a_{3,n})& \sum _{j=1}^n(a_{2,j}^2-a_{3,j}^2)\\ ...\\ 2(a_{n,1}-a_{n+1,1})& 2(a_{n,2}-a_{n+1,2})& ...& 2(a_{n,n}-a_{n+1,n})& \sum _{j=1}^n(a_{n,j}^2-a_{n+1,j}^2)\end{array}\right]$$

最后，就可以直接套板子了。

点击查看代码

#include <bits/stdc++.h>
using namespace std;

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

inline int read()
{
    int s = 0 , f = 1; char ch;
    while(!isdigit(ch = getchar())) if(ch == '-') f = -1;
    while(isdigit(ch)) s = s * 10 + ch - '0' , ch = getchar();
    return s * f;
}

int main()
{
    n = read();
    for(int i = 1;i <= n + 1;i++)
        for(int j = 1;j <= n;j++)
            cin >> a[i][j];
    for(int i = 1;i <= n;i++)
        for(int j = 1;j <= n;j++)
            b[i][j] = 2 * (a[i][j] - a[i + 1][j]),
            b[i][n + 1] += (a[i][j] * a[i][j] - a[i + 1][j] * a[i + 1][j]);
    for(int i = 1;i <= n;i++)
    {
        int m = i;
        for(int j = i + 1;j <= n;j++)
            if(fabs(b[j][i]) > fabs(b[m][i])) m = j;
        for(int j = 1;j <= n + 1;j++)
            swap(b[i][j] , b[m][j]);
        for(int j = 1;j <= n;j++)
            if(i != j)
            {
                double w = b[j][i] / b[i][i];
                for(int k = i + 1;k <= n + 1;k++)
                    b[j][k] -= w * b[i][k];
            }
    }
    for(int i = 1;i <= n;i++)
        printf("%.3lf " , b[i][n + 1] / b[i][i]);
    return 0;
}

欧几里得算法#

非常简单的一个算法。

$$gcd(a,b)=gcd(b,amodb)$$

拓展欧几里得算法#

扩展欧几里得算法是用来在已知 $(a,b)$ 时，求解一组 $(p,q)$ ，使得 $p\times a+q\times b=gcd(a,b)$。

可以简单的阐述一下这个算法：

因为：

$$gcd(a,b)=gcd(b,amodb)$$

所以：

$$p\times a+q\times b=gcd(a,b)=gcd(b,amodb)$$

$$=p\times b+q\times amodb=p\times b+q\times (a-a÷b\times b)$$

$$=q\times a+(p-a÷b\times q)\times b$$

程序只需要不断递归操作，最后当 $\text{b}$ 减少到 $\text{0}$ 时，就可以得出 $p=1,q=0$ ，然后递归回去就可以求出 $\text{p,q}$ 了。

点击查看代码

inline int exgcd(int a , int b , int &x , int &y)
{
    if(!b)
    {
        x = 1, y = 0;
        return a;
    }
    int res = exgcd(b , a % b , x , y) , tmp = x;
    x = y , y = tmp - a / b * y;
    return res;
}

一些例题#

【模板】二元一次不定方程 (exgcd)

模板题。

首先，必然可知，当且仅当 $gcd(a,b)|c$ 时，方程有整数解

然后，我们可以让方程变为拓展欧几里德算法可以求解的方程。

我们先求出：

$$a\times x^{\prime }+b\times y^{\prime }=gcd(a,b)$$

那么

$$\{\begin{array}{l}{x}_0=x^{\prime }\times c÷gcd(a,b)\\ y_0=y^{\prime }\times c÷gcd(a,b)\end{array}$$

即为原方程的一组特解。

而整数通解为：

$$\{\begin{array}{l}{x}_t=x_0-t\times b^{\prime }\\ y_t=y_0-t\times a^{\prime }\end{array}$$

其中 $x,y$ 的最小增减幅度分别为 ${\displaystyle \frac{b}{gcd(a,b)},{\displaystyle \frac{a}{gcd(a,b)}}}$ 即 $b^{\prime },a^{\prime }$。

接着就是，实现细节了。

1.无解

直接判断就可以了。

2.无正整数解

根据二元一次方程的通解可得，$x,y$ 的增减性相反，所以若存在正整数解，则当 $x$ 取得最小正整数解时， $y$ 取得最大正整数解，若此时 $y\le 0$,则说明该方程无正整数解。

3.有正整数解

同样，根据二元一次方程通解，求出 $x,y$ 的最大值和最小值。

而正整数解的数量就为 ${\displaystyle \frac{x_{max}-x_{min}}{b^{\prime }}+1}$ 或 ${\displaystyle \frac{y_{max}-y_{min}}{a^{\prime }}+1}$

点击查看代码

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

int t;

inline int read()
{
    int asd = 0 , qwe = 1; char zxc;
    while(!isdigit(zxc = getchar())) if(zxc == '-') qwe = -1;
    while(isdigit(zxc)) asd = asd * 10 + zxc - '0' , zxc = getchar();
    return asd * qwe;
}

inline int exgcd(int a , int b , int &x , int &y)
{
    if(!b)
    {
        x = 1, y = 0;
        return a;
    }
    int res = exgcd(b , a % b , x , y) , tmp = x;
    x = y , y = tmp - a / b * y;
    return res;
}

signed main()
{
    t = read();
    while(t--)
    {
        int a = read() , b = read() , c = read();
        int x = 0 , y = 0 , z = __gcd(a , b);
        if(c % z != 0)
        {
            puts("-1");
            continue;
        }
        exgcd(a , b , x , y);
        x = x * c / z , y = y * c / z;
        a /= z , b /= z;
        int t , xmax , xmin , ymax , ymin;
        t = ceil((double)(1 - x) / b);
        xmin = x + t * b , ymax = y - t * a;
        t = ceil((double)(1 - y) / a);
        xmax = x - t * b , ymin = y + t * a;
        if(ymax <= 0) cout << xmin << " " << ymin << endl;
        else cout << (xmax - xmin) / b + 1 << " " << xmin << " " << ymin << " " << xmax << " " << ymax << endl;
    }
    return 0;
}

P1516 青蛙的约会

几乎与上一题一模一样，只需要将式子改一下就可以了。

$$(n-m)\times t+l\times p=x-y$$

求出最小的 $t$ 即可。

点击查看代码

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

int t;

inline int read()
{
    int asd = 0 , qwe = 1; char zxc;
    while(!isdigit(zxc = getchar())) if(zxc == '-') qwe = -1;
    while(isdigit(zxc)) asd = asd * 10 + zxc - '0' , zxc = getchar();
    return asd * qwe;
}

inline int exgcd(int a , int b , int &x , int &y)
{
    if(!b)
    {
        x = 1, y = 0;
        return a;
    }
    int res = exgcd(b , a % b , x , y) , tmp = x;
    x = y , y = tmp - a / b * y;
    return res;
}

signed main()
{
    int a = read() , b = read() , m = read() , n = read() , l = read();
    int x = 0 , y = 0 , z;
    if((a - b) % (z = exgcd(n - m , l , x , y)) != 0)
        puts("Impossible");
    else cout << ((a - b) / z % (l / z) * x % (l / z) + l) % (l / z);
    return 0;
}

BSGS算法#

BSGS算法，全称Baby-Step-Giant-Step算法，又称大步小步算法。

它用于计算离散对数。

即一个最小的非负整数 $l$，满足 $a^l\equiv n(\mathrm{mod}p)$

我们可以运用整除分块的思想，将 $l$ 设为 $v\sqrt{p}-w$。

显然有：

$$u,v\le \sqrt{p}$$

那么：

$$a^l=a^{v\sqrt{p}-w}\equiv n(\mathrm{mod}p)$$

$$a^{v\sqrt{p}}\equiv n\times a^w(\mathrm{mod}p)$$

那么我们可以把右式的值存到一个哈希表里，每次遍历 $v$ 时直接查表就可以了。

复杂度 $O(\sqrt{p})$

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

模板题

点击查看代码

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#define int long long

using namespace std;
using namespace __gnu_pbds;

int n , b , p , s;
gp_hash_table<int , int> q;

inline int read()
{
    int s = 0 , f = 1; char ch;
    while(!isdigit(ch = getchar())) if(ch == '-') f = -1;
    while(isdigit(ch)) s = s * 10 + ch - '0' , ch = getchar();
    return s * f;
}

inline int power(int x , int y , int mod)
{
    int res = 1;
    while(y)
    {
        if(y & 1) res = res * x % mod;
        x = x * x % mod , y /= 2;
    }
    return res % mod;
}

signed main()
{
    p = read() , b = read() , n = read() , s = ceil(sqrt(p));
    for(int i = 1;i <= s;i++)
        q[(n * power(b , i , p) % p)] = i;
    for(int i = 1;i <= s;i++)
    {
        int res = power(b , i * s , p);
        if(q.find(res) != q.end())
        {
            // cout << i << " " << s << " " << res << " " << q[res] << endl;
            cout << i * s - q[res];
            return 0;
        }
    }
    cout << "no solution";
    return 0;
}

exBSGS#

拓展 $BSGS$ 算法，即 $ExBSGS$。

在 $BSGS$ 中，存在限制 $gcd(p,a)=1$，如果两者不互质的话需要使用拓展 $BSGS$ 算法。

继续考虑 $gcd(p,a)$，设其等于 $d$。

我们在原方程 $a^x\equiv b(\mathrm{mod}p)$ 中除掉一个 $d$，变成$\frac{a}{d}{a}^{x-1}\equiv \frac{b}{d}(\mathrm{mod}\frac{p}{d})$。

重复这个过程，直到 $gcd$ 为 1。或者再这个过程中 $b$ 不能整除 $d$，则无解。

于是问题变成了：

$$\frac{a^k}{D}{a}^{x-k}\equiv \frac{b}{D}(\mathrm{mod}\frac{p}{D})$$

此时 $gcd$ 为 1，可以使用 $BSGS$ 求解。

【模板】扩展 BSGS/exBSGS

这道题卡常属实有点狠。

直接套板子就可以了。

但要注意，因为这道题太过毒瘤，快速幂要写成递推式。

点击查看代码

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#define int long long

using namespace std;
using namespace __gnu_pbds;

int n , b , p , s , k , l;
gp_hash_table<int , int> q;

inline int read()
{
    int s = 0 , f = 1; char ch;
    while(!isdigit(ch = getchar())) if(ch == '-') f = -1;
    while(isdigit(ch)) s = s * 10 + ch - '0' , ch = getchar();
    return s * f;
}

inline int power(int x , int y , int mod)
{
    int res = 1;
    while(y)
    {
        if(y & 1) res = res * x % mod;
        x = x * x % mod , y /= 2;
    }
    return res % mod;
}

inline void exbsgs()
{
    q.clear() , b %= p , n %= p;
    if(n == 1 || p == 1)
    {
        puts("0");
        return;
    }
    l = 1 , k = 0;
    while(1)
    {
        int res = __gcd(b , p);
        if(res == 1) break;
        if(n % res != 0)
        {
            puts("No Solution");
            return;
        }
        p /= res , n /= res , l = (l * (b / res)) % p , k++;
        if(l == n)
        {
            printf("%lld\n" , k);
            return;
        }
    }
    s = ceil(sqrt(p));
    int ret = 1;
    for(int i = 1;i <= s;i++)
    {
        ret = ret * b % p;
        q[(n * ret % p)] = i;
    }
    ret = power(b , s , p);
    int sumres = 1;
    for(int i = 1;i <= s;i++)
    {
        sumres = sumres * ret % p;
        int res = (sumres * l) % p;
        if(q.find(res) != q.end())
        {
            printf("%lld\n" , i * s - q[res] + k);
            return;
        }
    }
    puts("No Solution");
}

signed main()
{
    b = read() , p = read() , n = read();
    while(b || p || n) exbsgs() , b = read() , p = read() , n = read();
    return 0;
}

至于例题大多都直接套板子就可以做出来，这里就不多放题了。

中国剩余定理 CRT#

求解方程组

$$\{\begin{array}{l}x\equiv a_1(\mathrm{mod}{m}_1)\\ x\equiv a_2(\mathrm{mod}{m}_2)\\ x\equiv a_3(\mathrm{mod}{m}_3)\\ ...\\ x\equiv a_n(\mathrm{mod}{m}_n)\end{array}$$

但原版的 $CRT$ 限制巨大，需要 $m$ 全都互质，所以这里就不再介绍了，直接就介绍 $exCRT$

拓展中国剩余定理 exCRT#

此时不保证所有的 $m_i$ 互质。那么显然最终的答案就是在模 $lcm(m_1,m_2,...,m_n)$ 意义下进行的。

考虑如何合并两个方程的答案。

假设就是合并：

$$\{\begin{array}{l}x\equiv a_1(\mathrm{mod}m{)}_1)\\ x\equiv a_2(\mathrm{mod}m{)}_2)\end{array}$$

首先肯定要合并模数为 $lcm(m_1,m_2)$ 。

显然有：

$x=k_1{m}_1+a_1=k_2{m}_2+a_2$。

化简后就是 $k_1{m}_1?k_2{m}_2=a_2?a_1$

$m_1,m_2,a_1,a_2$ 都是常数，这玩意就是一个 $exgcd$ 求解了。

如果 $exgcd$ 无解显然整个方程组都无解，否则求解后得到 $x$，合并模
数之后就得到了一个新的方程，然后接着把所有方程都合并就好了。

一些例题#

【模板】中国剩余定理(CRT)/曹冲养猪

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

点击查看代码

#include <bits/stdc++.h>
#define int __int128
using namespace std;

const int maxn = 100010;

int n , ans1 , ans2 , a[maxn] , b[maxn];

inline int read()
{
    int asd = 0 , qwe = 1; char zxc;
    while(!isdigit(zxc = getchar())) if(zxc == '-') qwe = -1;
    while(isdigit(zxc)) asd = asd * 10 + zxc - '0' , zxc = getchar();
    return asd * qwe;
}

inline int lcm(int a , int b)
{
    return a / __gcd(a , b) * b;
}

inline void exgcd(int a , int b , int &x , int &y)
{
    if(!b)
    {
        x = 1 , y = 0;
        return;
    }
    exgcd(b , a % b , x , y);
    int tmp = x;
    x = y , y = tmp - a / b * y;
}

signed main()
{
    n = read();
    for(int i = 1;i <= n;i++) a[i] = read() , b[i] = read() % a[i];
    ans1 = a[1] , ans2 = b[1];
    for(int i = 2;i <= n;i++)
    {
        int x = 0 , y = 0 , gcd = __gcd(ans1 , a[i]);
        exgcd(ans1 , a[i] , x , y);
        int c = b[i] - ans2;
        x = x * c / gcd % (a[i] / gcd);
        if(x < 0) x += a[i] / gcd;
        int mod = lcm(ans1 , a[i]); 
        ans2 = (ans1 * x + ans2) % mod;
        if(ans2 < 0) ans2 += mod;
        ans1 = mod;
    }
    cout << (long long)(ans2 % ans1);
    return 0;
}

莫比乌斯反演#

这个东西包含的内容很多，会时不时更新，并不能一次性学完，可能开始的内容及其基础，后面可能才会更加深入。

莫比乌斯函数:

设 $n=\prod _{i=1}^m{a}_i^{p_i}$ 。

若 $\mathrm{\forall }i,p_i>1$，那么 $\mu (n)=0$

否则 $\mu (n)=(-1{)}^m$。

莫比乌斯反演:

$g(n)=\sum _{d|n}f(d)$

$f(n)=\sum _{d|n}g(d)\mu (\frac{n}{d})$

这个东西最主要的是多写题，多运用，熟悉的就更快。

一些例题#

P3455 [POI2007]ZAP-Queries

可以推一下式子：

$$f(d)=\sum _{i=1}^a\sum _{j=1}^b[gcd(i,j)=d]$$

$$g(n)=\sum _{n|d}f(d)$$

根据函数意义：

$$g(d)=\sum _{i=1}^a\sum _{j=1}^b[d|gcd(i,j)]$$

$$g(d)=\sum _{i=1}^{a÷d}\sum _{j=1}^{b÷d}[1|gcd(i,j)]$$

$$g(d)=[\frac{a}{d}]\times [\frac{b}{d}]$$

$$f(d)=\sum _{i|d}\mu (\frac{d}{i})g(i)$$

$$f(d)=\sum _{i|d}\mu (\frac{d}{i})[\frac{a}{i}][\frac{b}{i}]$$

$$f(d)=\sum _{i=1}^{min(a,b)÷d}\mu (i)[\frac{a}{i\times d}][\frac{b}{i\times d}]$$

这个式子就可以用整除分块来做了。

点击查看代码

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

const int maxn = 50010;

int n , ans , tot , mu[maxn] , sum[maxn] , prime[maxn];

bool vis[maxn] , isprime[maxn];

inline int read()
{
    int asd = 0 , qwe = 1; char zxc;
    while(!isdigit(zxc = getchar())) if(zxc == '-') qwe = -1;
    while(isdigit(zxc)) asd = asd * 10 + zxc - '0' , zxc = getchar();
    return asd * qwe;
}

inline void get_mu()
{
    int s = 50000;
    mu[1] = 1;
    for(int i = 2;i <= s;i++)
    {
        if(!vis[i]) mu[i] = -1 , vis[i] = 1 , prime[++tot] = i , isprime[i] = 1;
        for(int j = 1;j <= tot && i * prime[j] <= s;j++)
        {
            vis[prime[j] * i] = 1;
            if(i % prime[j] == 0) break;
            else mu[prime[j] * i] = -mu[i];
        }
    }
    for(int i = 1;i <= s;i++) sum[i] = sum[i - 1] + mu[i];
}

signed main()
{
    get_mu() , n = read();
    for(int i = 1;i <= n;i++)
    {
        int a = read() , b = read() , d = read() , x = min(a , b);
        ans = 0;
        for(int l = 1 , r = 0;l <= x;l = r + 1)
            r = min(a / (a / l) , b / (b / l)),
            ans += (a / (l * d)) * (b / (l * d)) * (sum[r] - sum[l - 1]);
        cout << ans << endl;
    }
    return 0;
}

YY的GCD

虽然又是一道莫比乌斯反演模板题，但这个蒟蒻还是没能独立做出来。

但至少还是做出来了。

同样的，推一下式子:

我们设：$n<m$

$$\sum _{i=1}^n\sum _{j=1}^m[isprime(gcd(x,y))]$$

$$=\sum _{isprime(d)}^n\sum _{i=1}^n\sum _{j=1}^m[gcd(x,y)=d]$$

$$=\sum _{isprime(d)}^n\sum _{i=1}^{n÷d}\sum _{j=1}^{m÷d}[gcd(x,y)=1]$$

我们知道莫比乌斯函数的性质：

$$\sum _{d|n}\mu (d)=[n=1]$$

所以：

$$=\sum _{isprime(d)}^n\sum _{i=1}^{n÷d}\sum _{j=1}^{m÷d}\sum _{k|gcd(i,j)}\mu (k)$$

$$=\sum _{isprime(d)}^n\sum _{i=1}^{n÷d}\sum _{j=1}^{m÷d}\sum _{k|gcd(i,j)}\mu (k)$$

枚举 $k$ :

$$=\sum _{isprime(d)}^n\sum _{k=1}^{n÷k}\mu (k)\sum _{i=1}^{n÷d÷k}\sum _{j=1}^{m÷d÷k}$$

$$=\sum _{isprime(d)}^n\sum _{k=1}^{n÷k}\mu (k)[\frac{n}{d\times k}][\frac{m}{d\times k}]$$

设 $t=k\times d$：

$$=\sum _{isprime(d)}^n\sum _{k=1}^{n÷k}\mu (k)[\frac{n}{t}][\frac{m}{t}]$$

$$=\sum _{t=1}^n[\frac{n}{t}][\frac{m}{t}]\sum _{d|t,isprime(d)}\mu (\frac{t}{d})$$

前面一堆直接整除分块，后面的预处理前缀和即可。

复杂度$O(t\sqrt{n})$

点击查看代码

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

const int maxn = 10000010;

int t , n , m , cnt , ans , h[maxn] , mu[maxn] , sum[maxn] , prime[maxn];

bool vis[maxn];

inline int read()
{
    int asd = 0 , qwe = 1; char zxc;
    while(!isdigit(zxc = getchar())) if(zxc == '-') qwe = -1;
    while(isdigit(zxc)) asd = asd * 10 + zxc - '0' , zxc = getchar();
    return asd * qwe;
}

inline void get_mu()
{
	int s = 1e7; mu[1] = 1;
	for(int i = 2;i <= s;i++)
	{
		if(!vis[i]) mu[i] = -1 , prime[++cnt] = i , vis[i] = 1;
		for(int j = 1;j <= cnt && i * prime[j] <= s;j++)
		{
			vis[i * prime[j]] = 1;
			if(i % prime[j] == 0) break;
			mu[i * prime[j]] = -mu[i];
		}
	}
	for(int i = 1;i <= cnt;i++)
		for(int j = 1;j * prime[i] <= s;j++) h[j * prime[i]] += mu[j];
	for(int i = 1;i <= s;i++) sum[i] = sum[i - 1] + h[i];
}

signed main()
{
    t = read() , get_mu();
    while(t--)
    {
    	n = read() , m = read() , ans = 0;
    	if(n > m) swap(n , m);
    	for(int l = 1 , r = 0;l <= n;l = r + 1)
    	{
    		r = min(n / (n / l) , m / (m / l));
    		ans += (n / l) * (m / l) * (sum[r] - sum[l - 1]);
		}
		printf("%lld\n" , ans);
	}
    return 0;
}

[HAOI2011]Problem b

这一道题可以看做第一题的拓展：

$$f(d,x,y)=\sum _{i=1}^x\sum _{j=1}^y[gcd(i,j)=d]$$

那么根据容斥：

$$Ans=f(k,b,d)-f(k,a-1,d)-f(k,b,c-1)+f(k,a-1,c-1)$$

直接整除分块求解即可。

点击查看代码

#include <bits/stdc++.h>
using namespace std;

const int maxn = 50010;

int n , cnt , mu[maxn] , sum[maxn] , prime[maxn];

bool vis[maxn];

inline int read()
{
    int asd = 0 , qwe = 1; char zxc;
    while(!isdigit(zxc = getchar())) if(zxc == '-') qwe = -1;
    while(isdigit(zxc)) asd = asd * 10 + zxc - '0' , zxc = getchar();
    return asd * qwe;
}

inline void get_mu()
{
	int s = 50000; mu[1] = 1;
	for(int i = 2;i <= s;i++)
	{
		if(!vis[i]) mu[i] = -1 , prime[++cnt] = i , vis[i] = 1;
		for(int j = 1;j <= cnt && i * prime[j] <= s;j++)
		{
			vis[i * prime[j]] = 1;
			if(i % prime[j] == 0) break;
			mu[i * prime[j]] = -mu[i];
		}
	}
	for(int i = 1;i <= s;i++) sum[i] = sum[i - 1] + mu[i];
}

inline int solve(int k , int x , int y)
{
	int ans = 0;
	if(x > y) swap(x , y);
	for(int l = 1 , r = 0;l <= x;l = r + 1)
	{
		r = min(x / (x / l) , y / (y / l));
		ans += (sum[r] - sum[l - 1]) * (x / (l * k)) * (y / (l * k));
	}
	return ans;
}

int main()
{
    n = read() , get_mu();
    while(n--)
    {
    	int a = read() , b = read() , c = read() , d = read() , k = read();
    	cout << solve(k , b , d) - solve(k , a - 1 , d) - solve(k , b , c - 1) + solve(k , a - 1 , c - 1) << endl;
	}
    return 0;
}

[国家集训队]Crash的数字表格 / JZPTAB

推式子，推式子(〃'▽'〃)：

$$\sum _{i=1}^n\sum _{j=1}^{m}lcm(i,j)$$

$$=\sum _{i=1}^n\sum _{j=1}^m\frac{i\times j}{gcd(i,j)}$$

$$=\sum _{i=1}^n\sum _{j=1}^m\sum _{d|i,d|j,gcd(\frac{i}{d},\frac{j}{d})=1}\frac{i\times j}{d}$$

$$=\sum _{d=1}^{n}d\sum _{i=1}^{\frac{n}{d}}\sum _{j=1}^{\frac{m}{d}}[gcd(i,j)=1]\times i\times j$$

令：

$$g(x,y)=\sum _{i=1}^x\sum _{j=1}^y[gcd(i,j)=1]\times i\times j$$

$$=\sum _{i=1}^x\sum _{j=1}^{y}i\times j\sum _{k|gcd(i,j)}\mu (k)$$

$$=\sum _{k=1}^x\mu (k)k^2\sum _{i=1}^{\frac{x}{k}}\sum _{j=1}^{\frac{y}{k}}i\times j$$

$$=\sum _{k=1}^x\mu (k)k^2\times \frac{x}{k}\times (\frac{x}{k}+1)\times \frac{y}{k}\times (\frac{y}{k}+1)÷4$$

设：

$$sum(x,y)=x\times (x+1)\times y\times (y+1)÷4$$

则：

$$=\sum _{k=1}^x\mu (k)k^2\times sum(\frac{x}{k},\frac{y}{k})$$

整除分块 $O(\sqrt{n})$ 求解即可。

带回原式：

$$Ans==\sum _{d=1}^{n}d\times g(\frac{n}{d},\frac{m}{d})$$

太棒了，你会发现又是一个整除分块。

时间复杂度： $O(n)$

点击查看代码

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

const int mod = 20101009;
const int maxn = 10000010;

int n , m , cnt , sum , inv4 , mu[maxn] , sum1[maxn] , sum2[maxn] , prime[maxn];

bool vis[maxn];

inline int read()
{
    int asd = 0 , qwe = 1; char zxc;
    while(!isdigit(zxc = getchar())) if(zxc == '-') qwe = -1;
    while(isdigit(zxc)) asd = asd * 10 + zxc - '0' , zxc = getchar();
    return asd * qwe;
}

inline void get_mu()
{
	int s = 1e7; mu[1] = 1;
	for(int i = 2;i <= s;i++)
	{
		if(!vis[i]) mu[i] = -1 , vis[i] = 1 , prime[++cnt] = i;
		for(int j = 1;j <= cnt && prime[j] * i <= s;j++)
		{
			vis[i * prime[j]] = 1;
			if(i % prime[j] == 0) break;
			mu[i * prime[j]] = -mu[i];
		}
	}
	for(int i = 1;i <= s;i++) sum2[i] = (sum2[i - 1] + i) % mod;
	for(int i = 1;i <= s;i++) sum1[i] = ((sum1[i - 1] + mu[i] * i * i % mod) % mod + mod) % mod;
}

inline int power(int x , int y)
{
	int res = 1;
	while(y)
	{
		if(y & 1) res = res * x % mod;
		x = x * x % mod , y /= 2; 
	}
	return res;
}

inline int inv(int x)
{
	return power(x , mod - 2);
}

inline int calc(int x , int y)
{
	return x * (x + 1) % mod * y % mod * (y + 1) % mod * inv4 % mod;
}

inline int solve(int x , int y)
{
	int ans = 0;
	if(x > y) swap(x , y);
	for(int l = 1 , r = 0;l <= x;l = r + 1)
	{
		r = min(x / (x / l) , y / (y / l));
		ans = (ans + ((sum1[r] - sum1[l - 1]) % mod + mod) % mod * calc(x / l , y / l) % mod) % mod; 
	}
	return ans;
} 

signed main()
{
    n = read() , m = read();
    inv4 = inv(4) , get_mu();
    if(n > m) swap(n , m);
    for(int l = 1 , r = 0;l <= n;l = r + 1)
    {
    	r = min(n / (n / l) , m / (m / l));
    	sum = (sum + (sum2[r] - sum2[l - 1]) * ((solve(n / l , m / l) % mod + mod) % mod) % mod) % mod;
	}
	cout << (sum % mod + mod);
    return 0;
}

[SDOI2015]约数个数和

上来就直接推式子：

首先，必然会有：

$$d(ij)=\sum _{x|i}\sum _{y|j}[gcd(x,y)=1]$$

所以：

$$\sum _{i=1}^n\sum _{j=1}^{m}d(ij)$$

$$=\sum _{i=1}^n\sum _{j=1}^m\sum _{x|i}\sum _{y|j}[gcd(x,y)=1]$$

$$=\sum _{i=1}^n\sum _{j=1}^m\sum _{x|i}\sum _{y|j}\sum _{d|gcd(x,y)}\mu (d)$$

$$=\sum _{x=1}^n\sum _{y=1}^m[\frac{n}{x}][\frac{m}{y}]\sum _{d|gcd(x,y)}\mu (d)$$

$$=\sum _{d=1}^n\mu (d)\sum _{x=1}^{\frac{n}{d}}\sum _{y=1}^{\frac{m}{d}}[\frac{n}{xd}][\frac{m}{xd}]$$

$$=\sum _{d=1}^n\mu (d)\sum _{x=1}^{\frac{n}{d}}[\frac{n}{xd}]\sum _{y=1}^{\frac{m}{d}}[\frac{m}{xd}]$$

预处理：

$$sum(x)=\sum _{i=1}^x[\frac{x}{i}]$$

所以原式：

$$=\sum _{d=1}^n\mu (d)sum[\frac{n}{d}]sum[\frac{m}{d}]$$

好像现在这样就可以直接前缀和+整除分块了。

点击查看代码

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

const int maxn = 50010;

int t , n , m , cnt , ans , g[maxn] , mu[maxn] , sum[maxn] , prime[maxn];

bool vis[maxn];

inline int read()
{
    int asd = 0 , qwe = 1; char zxc;
    while(!isdigit(zxc = getchar())) if(zxc == '-') qwe = -1;
    while(isdigit(zxc)) asd = asd * 10 + zxc - '0' , zxc = getchar();
    return asd * qwe;
}

inline void get_mu()
{
	int s = 50000; mu[1] = 1;
	for(int i = 2;i <= s;i++)
	{
		if(!vis[i]) mu[i] = -1 , prime[++cnt] = i , vis[i] = 1;
		for(int j = 1;j <= cnt && prime[j] * i <= s;j++)
		{
			vis[i * prime[j]] = 1;
			if(i % prime[j] == 0) break;
			mu[i * prime[j]] = -mu[i];
		}
	}
	for(int i = 1;i <= s;i++) sum[i] = sum[i - 1] + mu[i];
	for(int i = 1;i <= s;i++) for(int l = 1 , r = 0;l <= i;l = r + 1)
		r = i / (i / l) , g[i] += (i / l) * (r - l + 1);
}

signed main()
{
    t = read() , get_mu();
    while(t--)
    {
    	n = read() , m = read() , ans = 0;
    	if(n > m) swap(n , m);
    	for(int l = 1 , r = 0;l <= n;l = r + 1)
    	{
    		r = min(n / (n / l) , m / (m / l));
    		ans += (sum[r] - sum[l - 1]) * g[n / l] * g[m / l];
		}
		cout << ans << endl;
	}
    return 0;
}