数学专题学习笔记

开坑于-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 \leq 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}\) 互质的数的数目。

\[φ(x)=\prod_{i=1}^n(1-\frac{1}{a_{i}}) \]

其中 \(a_{i}\) 为 \(\text{x}\) 的质因子，\(\text{n}\) 为 \(\text{x}\) 的因子个数。

几个性质

若 \(\text{x}\) 为质数，显然 \(φ(x)=x-1\)
若 \(m,n\) 互质，则 \(φ(mn)=φ(m)\timesφ(n)\)
若 \(x=2n~,~n \mod 2 = 1\)，则 \(φ(x)=φ(n)\)
若 \(\text{x}\) 为质数，则 \(φ(x^k)= x^k-x^{k-1}\)
若 \(x>2\)，则 \(φ(x) \mod 2 = 0\)
小于 \(n\) 且与 \(n\) 互质的数的总和为 \(φ(n) \times n \div 2(n>1)\)
\(n=\sum_{d|n}φ(d)\)，即 \(\text{n}\) 包括 \(\text{n}\) 在内的所有因数的欧拉函数之和等于 \(\text{n}\)，这一条也被称为欧拉反演。

求法

有两种求法。

\(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;
}

\(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{cases} 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{cases}\]

如何实现

P3389 【模板】高斯消元法

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

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

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

首先，列出增广矩阵。

\[\begin{bmatrix} 1&3&4&5\\ 1&4&7&3\\ 9&3&2&2 \end{bmatrix}\]

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

\[\begin{bmatrix} 9&3&2&2\\ 1&3&4&5\\ 1&4&7&3 \end{bmatrix}\]

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

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

\[\begin{bmatrix} 9&0&0&-8.76\\ 0&3.67&0&19.01\\ 0&0&-1.15&2.76 \end{bmatrix}\]

格式化的结果

\[\begin{bmatrix} a_1&0&0&b_1\\ 0&a_2&0&b_2\\ 0&0&a_3&b_3 \end{bmatrix}\]

最终，输出答案，答案为 \(\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}^n2(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) \]

化成矩阵即：

\[\begin{bmatrix} 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{bmatrix}\]

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

点击查看代码

#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,a~mod~b) \]

拓展欧几里得算法

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

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

因为：

\[\gcd(a,b)=\gcd(b,a\bmod b) \]

所以：

\[p \times a + q\times b=\gcd(a,b)=\gcd(b,a\bmod b) \]

\[=p\times b+q\times a\bmod b=p\times b+q\times (a-a\div b\times b) \]

\[=q\times a+(p-a\div 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'+b \times y'=\gcd(a,b) \]

那么

\[\begin{cases} x_0=x' \times c \div \gcd(a,b)\\ y_0=y' \times c \div \gcd(a,b) \end{cases}\]

即为原方程的一组特解。

而整数通解为：

\[\begin{cases} x_t=x_0 - t \times b'\\ y_t=y_0 - t \times a' \end{cases}\]

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

接着就是，实现细节了。

1.无解

直接判断就可以了。

2.无正整数解

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

3.有正整数解

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

而正整数解的数量就为 \(\displaystyle\frac{x_{max}-x_{min}}{b'}+1\) 或 \(\displaystyle\frac{y_{max}-y_{min}}{a'}+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 \pmod p\)

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

显然有：

\[u,v \leq \sqrt{p} \]

那么：

\[a^l = a^{v\sqrt{p}-w} \equiv n \pmod p \]

\[a^{v\sqrt{p}}\equiv n \times a^w \pmod 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 ≡ b \pmod p\) 中除掉一个 \(d\)，变成\(\frac{a}{d}a^{x-1} \equiv \frac{b}{d} \pmod {\frac{p}{d}}\)。

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

于是问题变成了：

\[\frac{a^k}{D}a^{x-k} \equiv \frac{b}{D} \pmod {\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{cases} x \equiv a_1 \pmod {m_1}\\ x \equiv a_2 \pmod {m_2}\\ x \equiv a_3 \pmod {m_3}\\ ...\\ x \equiv a_n \pmod {m_n} \end{cases}\]

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

拓展中国剩余定理 exCRT

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

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

假设就是合并：

\[\begin{cases} x ≡ a_1\pmod m_1)\\ x ≡ a_2\pmod m_2) \end{cases}\]

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

显然有：

\(x = k_1m_1 + a_1 = k_2m_2 + a_2\)。

化简后就是 \(k_1m_1 ? k_2m_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^m_{i=1}a_i^{p_i}\) 。

若 \(\forall i, p_i > 1\)，那么 \(μ(n) = 0\)

否则 \(μ(n) = (-1)^m\)。

莫比乌斯反演:

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

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

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

一些例题

P3455 [POI2007]ZAP-Queries

可以推一下式子：

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

\[g(n)=\sum\limits_{n|d}f(d) \]

根据函数意义：

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

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

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

\[f(d)=\sum\limits_{i|d}\mu(\frac{d}{i})g(i) \]

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

\[f(d)=\sum\limits_{i=1}^{\min(a,b)\div 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\limits_{i=1}^n\sum\limits_{j=1}^m[isprime(\gcd(x,y))] \]

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

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

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

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

所以：

\[=\sum\limits_{isprime(d)}^n\sum\limits_{i=1}^{n\div d}\sum\limits_{j=1}^{m\div d}\sum\limits_{k|\gcd(i,j)}\mu(k) \]

枚举 \(k\) :

\[=\sum\limits_{isprime(d)}^n\sum\limits_{k=1}^{n\div k}\mu(k)\sum\limits_{i=1}^{n\div d\div k}\sum\limits_{j=1}^{m\div d\div k} \]

\[=\sum\limits_{isprime(d)}^n\sum\limits_{k=1}^{n\div k}\mu(k)[\frac{n}{d\times k}][\frac{m}{d\times k}] \]

设 \(t = k\times d\)：

\[=\sum\limits_{isprime(d)}^n\sum\limits_{k=1}^{n\div k}\mu(k)[\frac{n}{t}][\frac{m}{t}] \]

\[=\sum\limits_{t=1}^n[\frac{n}{t}][\frac{m}{t}]\sum\limits_{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\limits_{i=1}^{x}\sum\limits_{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\limits_{i=1}^n\sum\limits_{j=1}^mlcm(i,j) \]

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

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

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

令：

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

\[=\sum\limits_{i=1}^{x}\sum\limits_{j=1}^{y}i\times j\sum\limits_{k|\gcd(i,j)}\mu(k) \]

\[=\sum\limits_{k=1}^x\mu(k)k^2\sum\limits_{i=1}^{\frac{x}{k}}\sum\limits_{j=1}^{\frac{y}{k}}i\times j \]

\[=\sum\limits_{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) \div 4 \]

设：

\[sum(x,y)=x\times (x+1)\times y\times (y+1) \div 4 \]

则：

\[=\sum\limits_{k=1}^x\mu(k)k^2\times sum(\frac{x}{k},\frac{y}{k}) \]

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

带回原式：

\[Ans==\sum\limits_{d=1}^nd\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\limits_{x|i}\sum\limits_{y|j}[\gcd(x,y)=1] \]

所以：

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^md(ij) \]

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

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

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

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

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

预处理：

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

所以原式：

\[=\sum\limits_{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;
}

posted @ 2022-01-24 20:18 JiaY19 阅读(15) 评论(0) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部

Loading

JiaY19

数学专题学习笔记

约数

一些例题

欧拉函数

几个性质

求法

一些例题

高斯消元

如何实现

一些例题

欧几里得算法

拓展欧几里得算法

一些例题

BSGS算法

exBSGS

中国剩余定理 CRT

拓展中国剩余定理 exCRT

一些例题

莫比乌斯反演

一些例题

公告