「学习笔记」初等数论

wdgm4 · 2023-03-12 20:29:14 · 个人记录

前言

更好的阅读体验?

前置知识（这个应该很显然）：\operatorname{lcm}(a,b)=\dfrac{ab}{\gcd(a,b)}

线性筛素数

直接上代码。

const int MAXN=100000008;
bool np[MAXN];
vector<int> prm,pre;
void gg(const int N=100000000){
    pre.resize(N+1);
    for(int i=2;i<=N;i++){
        if(np[i]==false){
            prm.push_back(i);
            pre[i]=i;
        }
        for(auto j:prm) if(i*j<=N){
            int k=i*j;
            pre[k]=j;
            np[k]=true;
            if(i%j==0) break;
        }else break;
    }
}

例题：B3716 分解质因子 3

说明一下，prm 存了所有范围内的素数，pre 存了第 i 个数的最小质因子。

这个板子直接背过就好了。（其实是我懒得讲）

不等式方程（exgcd）

形如 ax+by=c 的方程（a,b 为整数），方程有解当且仅当 \gcd(a,b) \mid c

欧几里得算法：\gcd(a,b)=\gcd(b,a \bmod b)。

递归函数

void gg(int &x,int &y,int a,int b){
    if(a%b==0){
        x=0,y=1;
    }else{
        int xx,yy;
        gg(xx,yy,b,a%b);
        x=yy;y=xx-a/b*yy;
    }
}

裴蜀定理（写于 2023.4.14）

这个是后补的。QWQ

定义

设 a,b 是不全为零的整数，则存在整数 x,y , 使得 ax+by=\gcd(a,b)。

然后就发现前面的式子求的就是裴蜀定理的式子。

然后写一个我自己发现的东西：（只不过我看过的文章都没讲过，可能是因为为太简单了）

求 ax+by=c 的解，其中 a,b 为整数且 \gcd(a,b)|c。

首先求出 ax+by=\gcd(a,b) 的解，设解为 \begin{cases}x=x'\\y=y'\end{cases}

再设 k=\dfrac{c}{\gcd(a,b)}，则 ax+by=c 的解为 \begin{cases}x=kx'\\y=ky'\end{cases}

然后就没了。QWQ

进一步的结论可以去下面的链接。

裴蜀定理

同余

概念

对两个整数 a，b，如果它们除以 d 的余数相同，则称它们 **模

3. 根据积性可得。

例题：P2158 [SDOI2008] 仪仗队。

假如 C君 在仪仗队的左前方，坐标为 (0,0)。

首先可以发现，对于一个坐标 (i,j)，C君 能看到这个位置当且仅当 \gcd(i,j)=1，并且在与 C君 位置相交的对角线左右两边能看到的位置对称。

可以推个柿子。

当 n\le 2 时：

用线性筛欧拉函数即可，记得当 n=1 时特判一下。

代码：

#include<bits/stdc++.h>
#define XD 114514
#define MAXN 40010
using namespace std;
int n;
int phi[MAXN];
bool np[MAXN];
vector<int> prm,pre;
long long ans;
void gg(){
    pre.resize(n+10);
    for(int i=2;i<=n;i++){
        if(np[i]==false){
            phi[i]=i-1;
            pre[i]=i;
            prm.push_back(i);
        }
        for(int j=0;j<prm.size();j++){
            if(i*prm[j]<=n){
                int k=i*prm[j];
                pre[k]=prm[j];
                np[k]=true;
                if(i%prm[j]==0){
                    phi[k]=phi[i]*prm[j];
                    break;
                }else{
                    phi[k]=phi[i]*phi[prm[j]];
                }
            }else break;
        }
    }
}
int main(){
    cin>>n;
    if(n==1) cout<<0;
    else{
        gg();
        for(int i=1;i<=n-1;i++) ans+=phi[i];
        cout<<(ans+1)*2+1;
    }
    return 0;
}

2023.9.9

P2398 GCD SUM

（GCDMAT - GCD OF MATRIX 双倍经验）

求：

运用欧拉反演，可得：

转换，得：

线性欧拉筛即可。（可以用数论分块，也可以不用，因为 n 范围太小了）

不用数论分块：

#include<bits/stdc++.h>
#define XD 114514
#define MAXN 100010
#define int long long
using namespace std;
int ans;
int n,phi[MAXN];
vector<int> prm;
bool vis[MAXN];
void prime(){
    phi[1]=1;
    for(int i=2;i<=n;i++){
        if(!vis[i]){
            prm.push_back(i);
            phi[i]=i-1;
        }
        for(int j=0;j<prm.size();j++){
            if(i*prm[j]<=n){
                int nem=i*prm[j];
                vis[nem]=true;
                if(i%prm[j]==0){
                    phi[nem]=phi[i]*prm[j];
                    break;
                }else{
                    phi[nem]=phi[i]*phi[prm[j]];
                }
            }else break;
        }
    }
}
signed main(){
    cin>>n;
    prime();
    for(int i=1;i<=n;i++){
        ans+=phi[i]*(n/i)*(n/i);
    }
    cout<<ans;
    return 0;
}

用数论分块：

#include<bits/stdc++.h>
#define XD 114514
#define MAXN 100010
#define int long long
using namespace std;
int ans;
int n,sum[MAXN],phi[MAXN];
vector<int> prm;
bool vis[MAXN];
void prime(){
    sum[1]=phi[1]=1;
    for(int i=2;i<=n;i++){
        if(!vis[i]){
            prm.push_back(i);
            phi[i]=i-1;
        }
        for(int j=0;j<prm.size();j++){
            if(i*prm[j]<=n){
                int nem=i*prm[j];
                vis[nem]=true;
                if(i%prm[j]==0){
                    phi[nem]=phi[i]*prm[j];
                    break;
                }else{
                    phi[nem]=phi[i]*phi[prm[j]];
                }
            }else break;
        }
        sum[i]=sum[i-1]+phi[i];
    }
}
signed main(){
    cin>>n;
    prime();
    int r;
    for(int l=1;l<=n;l=r+1){
        int nem=n/l;
        r=n/nem;
        ans+=(sum[r]-sum[l-1])*(n/l)*(n/l);
    }
    cout<<ans;
    return 0;
}

P2303 [SDOI2012] Longge 的问题

手推的式子（很抽象）

单次求欧拉函数即可。

#include<bits/stdc++.h>
#define XD 114514
#define int long long
using namespace std;
int n,ans;
int phi(int x){
    int num=x;
    for(int i=2;i*i<=x;i++){
        if(x%i==0) num=num/i*(i-1);
        while(x%i==0) x/=i;
    }
    if(x>1) num=num/x*(x-1);
    return num;
}
signed main(){
    cin>>n;
    for(int i=1;i*i<=n;i++){
        if(i*i==n) ans+=i*phi(n/i);
        else if(n%i==0) ans+=i*phi(n/i)+(n/i)*phi(i);
    }
    cout<<ans;
    return 0;
}

2023.9.12

UVA11417 GCD / UVA11424 GCD - Extreme (I) /P1390 公约数的和/ SP3871 GCDEX - GCD Extreme /UVA11426 拿行李（极限版） GCD - Extreme (II) （5 倍经验，按数据范围从小到大排序）

求：

转化为：

后面的直接套用上一道题的式子，可以提前线性筛欧拉函数，然后提前把答案都求出来，这样就可以把前两倍经验给切了。QWQ

#include<bits/stdc++.h>
#define XD 114514
#define MAXN 200010
using namespace std;
int n;
bool vis[MAXN];
int phi[MAXN];
vector<int> prm;
void sieve(){
    phi[1]=1;
    for(int i=2;i<=2e5+1;i++){
        if(!vis[i]){
            prm.push_back(i);
            phi[i]=i-1;
        }
        for(int j=0;j<prm.size() and i*prm[j]<=2e5+1;j++){
            vis[i*prm[j]]=true;
            if(i%prm[j]==0){
                phi[i*prm[j]]=phi[i]*prm[j];
                break;
            }
            phi[i*prm[j]]=phi[i]*phi[prm[j]];
        }
    }
}
long long ans[MAXN],num;
int main(){
    ios::sync_with_stdio(false);
    cin.tie(0);cout.tie(0);
    sieve();
    for(int i=1;i<=2e5+1;i++){
        for(int j=1;j*j<=i;j++){
            if(j*j==i) num+=j*phi[j];
            else if(i%j==0) num+=1ll*j*phi[i/j]+(i/j)*phi[j];
        }
        num-=i;
        ans[i]=num;
    }
    while(true){
        cin>>n;if(n==0) break;
        cout<<ans[n]<<"\n";
    }
    return 0;
}

时间复杂度预处理 O(n\sqrt{n})，查询 O(1)。

考虑优化。

抽象的手推式子

看左面的式子。

还是预处理线性筛欧拉函数，查询直接数论分块即可。

在代码里，为防止爆 long long，我把要减的 \frac{(n+1)\times n}{2} 直接在数论分块中算了。QWQ

#include<bits/stdc++.h>
#define XD 114514
#define MAXN 4000010
using namespace std;
int n;
bool vis[MAXN];
int phi[MAXN];
long long num[MAXN];
vector<int> prm;
void sieve(){
    phi[1]=1;num[1]=1;
    for(int i=2;i<=4e6;i++){
        if(!vis[i]){
            prm.push_back(i);
            phi[i]=i-1;
        }
        num[i]=num[i-1]+phi[i];
        for(int j=0;j<prm.size() and i*prm[j]<=4e6;j++){
            vis[i*prm[j]]=true;
            if(i%prm[j]==0){
                phi[i*prm[j]]=phi[i]*prm[j];
                break;
            }
            phi[i*prm[j]]=phi[i]*phi[prm[j]];
        }
    }
}
long long ans;
int main(){
    ios::sync_with_stdio(false);
    cin.tie(0);cout.tie(0);
    sieve();
    while(true){
        ans=0;
        cin>>n;
        if(n==0) break;
        int r;
        for(int l=1;l<=n;l=r+1){
            r=n/(n/l);
            ans+=(1ll*(r+1)*r/2-1ll*(l-1)*l/2)*(1ll*num[n/l]-1);
        }
        cout<<ans<<"\n";
    }

    return 0;
}

说句闲话，这篇题解的式子（是我上面推的式子的方法）有点问题。QWQ

2023.10.3

膜运算

这是一个裸题，因为显然可以用欧几里得定理搞成原题。

而我才不会告诉你我想的有亿点复杂。

推式子。

2023.10.9

发现一个可能有用的等式：

欧拉定理

对任意正整数 a,m，\gcd(a,m)=1，一定有：a^{\varphi(m)} \equiv 1 \pmod m。

逆元

逆元补充 （by 2023.5.11）

定义

对于任意正整数 x 满足 x \times x^{-1} \equiv 1 \pmod p，则称 x^{-1} 为 x 在模 p 同余下的逆元（其实就跟普通的 x \times x^{-1}=1 中的 x^{-1} 差不多，只不过这里的逆元是在模 p 同余下，而且在模 p 同余下，x\times x^{-1}=1，然后就可以得出 x^k \times x^{-1} \equiv x^{k-1} \pmod p）

为了方便理解，我就先把下文中所有的字母下面带 0 的都为这个字母代表的数在模某一个数同余下的逆元。（因为我在听课时老师是用 x^{-1} 表示 x 在模 p 同余下的逆元QWQ）

性质

逆元存在性定理：x 在模 p 同余下存在逆元当且仅当 \gcd(x,p)=1。

推论：当且仅当模数 p 是质数时，[1, p − 1] 内所有整数都存在模 p 下的逆元。 0 没有逆元。

逆元唯一性定理：模 p 同余下，一个整数 x 的逆元若存在，则唯一。

定理：在模质数 p 同余下， [1, p − 1] 内所有整数的逆元互不相同。

定理：一个数的逆元的逆元等于它自身。

计算逆元

先上例题。

P1082 [NOIP2012 提高组] 同余方程

根据逆元定义可得这是一道求逆元板子题（显然）。

因为题目说输入数据保证一定有解，根据逆元存在性定理可得 \gcd(a,b)=1。

开始求逆元喽。

（由于 k 为任意整数，我们不知道它是正数还是负数，就把前面的正负号省略了）

然后可以发现，这不就是一个不定方程吗。（因为 a,b 为已知量，只需求 x）

在前面我们知道 \gcd(a,b)=1，所以可以用递归求解即可。

最后记得把 x 取模到 [1,b-1] 的范围即可。QWQ

代码如下。

#include<bits/stdc++.h>
#define XD 114514
#define int long long
using namespace std;
inline void gg(int &x,int &y,int a,int b){
    if(a%b==0){
        x=0,y=1;
    }else{
        int xx,yy;
        gg(xx,yy,b,a%b);
        x=yy;
        y=xx-a/b*yy;
    }
}
signed main(){
    int a,b,x,y;
    cin>>a>>b;
    gg(x,y,a,b);
    cout<<(x%b+b)%b;
    return 0;
}

当然也可以像下面这样写。

#include<bits/stdc++.h>
#define XD 114514
#define int long long
using namespace std;
inline void gg(int &x,int &y,int a,int b){
    if(a%b==0){
        x=0,y=1;
        return;
    }
    gg(y,x,b,a%b);
    y-=a/b*x;
}
signed main(){
    int a,b,x,y;
    cin>>a>>b;
    gg(x,y,a,b);
    cout<<(x%b+b)%b;
    return 0;
}

补充

写于 2023.5.25

由于扩展中国剩余定理中有一些比较难懂的东西，结果发现是前面学的有问题，所以在这里来个补充。

先看一道题。（这是我自己出的QWQ）

解 ax \equiv b \pmod p

按上面的思路可得 ax+pk=b

由裴蜀定理可得，如果 \gcd(a,p) \nmid b 则无解。

设 \gcd(a,p)=g

用 exgcd 可求 au+pv=g

即 x=\dfrac{bu}{g}，k=\dfrac{bv}{g}

这说明 x_0\equiv \dfrac{bu}{g}\pmod p 是 ax\equiv b \pmod p 的解。

设 x_1 是 ax\equiv b \pmod p 的其他解。

则 ax_1 \equiv ax_0 \pmod p

所以 p \mid ax_1-ax_0

这蕴含 \dfrac{p}{g} \mid \dfrac{a(x_1-x_0)}{g}

因为 \dfrac{p}{g}，\dfrac{a}{g} 没有公约数。

所以 \dfrac{p}{g} \mid x_1-x_0

即 x_1=x_0+k\dfrac{p}{g}（k \in \mathbb{Z}）

而相差 p 的倍数的解被规定相同，所以恰好有 g 个不同的解。

即取 k=0,1,...,g-1。

其实学完扩展中国剩余定理再回来看可以感觉到其实 ax+by=c 也不一定只有一个解或无解。（只不过在 OI 中我们只需要求差不多最小正整数解即可QWQ）

好像很显然的样子QWQ。

欧拉定理的继续

欧拉定理有一个前提条件就是 \gcd(a,m)=1，根据逆元存在性定理就可得出这个 a 在模 m 同余下一定有逆元。让我们推个式子。QWQ

然后我们就得到了另一个求逆元的方法。（虽然不是很好用）

但是当 m 是质数时，a^{\varphi(m)} \equiv 1 \pmod m 就可以转化成 a^{m-1} \equiv 1 \pmod m 于是费马小定理出现了。

费马小定理

对任意整数 a 和质数 p，a^{p-1} \equiv 1 \pmod p。

根据式子可以推出：

然后我们就又双叒叕得到了一个求逆元的方法（用快速幂求一下 a^{p-2} \bmod p 即可）。

上个题QWQ P2613 【模板】有理数取余

求 \frac{a}{b} \bmod 19260817 的值。

让我们算一算吧。

所以我们只需算出 b 模 19260817 同余下的逆元。由于 19260817 为一个质数，用费马小定理求即可。

而无解的情况就是 b 是 19260817 倍数，那么 b \bmod 19260817 一定等于 0。

直接上代码。

#include<bits/stdc++.h>
#define XD 114514
#define mod 19260817
using namespace std;
inline long long read(){
    char ch=getchar(),h;long long w=0;
    while(ch>'9' or ch<'0') h=ch,ch=getchar();
    while(ch>='0' and ch<='9') w=(w*10+(ch-'0'))%mod,ch=getchar();
    if(h=='-') w=-w;
    return w;
}
long long a,b;
inline long long gg(long long x,long long k){
    //求 x 在模 mod 同余下的逆元，即求 x^(mod-2)%mod 的值
    //这里 k 就表示 mod-2 的值 
    long long ans=1;
    while(k){
        if(k&1){
            ans*=x;ans%=mod;
        }
        x*=x;x%=mod;
        k>>=1;
    }
    return ans;
}//快速幂是啥偶就不讲了QWQ 
int main(){
    a=read(),b=read();
    if(b==0) cout<<"Angry!\n";
    else cout<<a*gg(b,mod-2)%mod;
    return 0;
}

线性求逆元

设 inv_i 表示 i 的逆元，有递推式：

边界为 inv_1 = 1。

又要开始推式子了。

再解释一下，我们刚开始设 p=i \times k +r，其中 r \in [1,n)，显然，k=\lfloor \frac{p}{i} \rfloor，r=p \bmod i，这样应该能看懂了吧。QWQ

P3811 【模板】乘法逆元

下面是代码。

#include<bits/stdc++.h>
#define XD 114514
#define MAXN 3000010
using namespace std;
int n,p,inv[MAXN];
int main(){
    ios::sync_with_stdio(false);
    cin>>n>>p;
    inv[1]=1;
    for(int i=2;i<=n;i++){
//      inv[i]=(1ll*-(p/i)*inv[p%i]%p+p)%p;
        inv[i]=1ll*(p-p/i)*inv[p%i]%p;//这两个算式等价，只不过这个少取模一次 
        //实测，在P3811中第一个用398ms，第二个用338ms
    }
    for(int i=1;i<=n;i++) cout<<inv[i]<<'\n';
    return 0;
}

说句闲话，我刚开始做这道题时一直 \color{black}{TLE} 了两个点，最后把换行符 endl 改为 '\n' 就 \color{green}{AC} 了，这个故事告诉我们，珍爱生命，远离 endl。

B3717 组合数问题

首先知道：

显然：

之后就把 998244353 当做 p 就行了不然打一个九位数太麻烦了。

然后求 (m!)^{-1}，显然：

所以先线性求一遍逆元，然后预处理出所有阶乘的逆元即可。

上代码。

#include<bits/stdc++.h>
#define XD 114514

using namespace std;
int T,N;
const int p=998244353;
int ans;
int main(){
    ios::sync_with_stdio(false);
    cin>>T>>N;
    vector<int> inv(N+1),a(N+1),b(N+1);
    inv[1]=1;//inv_i 表示 i 的逆元 
    a[1]=1;//a_i 表示 i!的值 
    b[0]=b[1]=1;//b_i 表示 i! 的逆元的值 
    for(int i=2;i<=N;i++){
        inv[i]=1ll*(p-p/i)*inv[p%i]%p;
        a[i]=1ll*a[i-1]*i%p;
        b[i]=1ll*b[i-1]*inv[i]%p;
    }
    while(T--){
        int n,m,num=1;
        cin>>n>>m;
        num=1ll*a[n]*b[m]%p*b[n-m]%p;
        ans^=num;
    }
    cout<<ans;
    return 0;
}

威尔逊定理（on 2023.7.11）

定理

若 p 为素数，那么：

证明

当 p=2 时，(p-1)! \equiv 1 \equiv -1 \pmod p

取一个大于 2 的质数 p，那么对于任意正整数 x （1 \le x \le p-1），它一定有一个逆元，且这 p-1 个数的逆元不同。对于这些逆元，逆元是它本身的是 1 和 p-1。我们可以把剩下的整数 2,3,4,...,p-1 分成 (p-3)/2 份，每一份的乘积都模 p 余 1，这说明：

再在两边同时乘上 1 和 p-1，得：

证毕。

它的逆定理为：

正整数 p 是质数的充要条件为：

(p-1)! \equiv -1 \pmod p

威尔逊定理用的很少，不过扩展卢卡斯定理可能要用。

中国剩余定理（CRT）

它是用来解决一个问题的。

问题

给定 n 个线性同余方程：

其中保证 m_1,m_2,...,m_n 两两互质。求 x 的通解。

定理

上面的方程组有解，且按如下方式构造：

记 M=\prod\limits^n_{i=1}m_i ，M_i=\frac{M}{m_i}，t_i 为 M_i 在模 m_i 同余下的逆元。

则方程组的唯一通解是：

x \equiv \sum\limits^n_{i=1}a_i t_i M_i \pmod M

证明

显然：

且，

（能够保证 M_i 在模 m_i 同余下一定有逆元）

另外，对于一个 j（j \ne i），由于 m_i | M_j，于是 M_j \equiv 0 \pmod{m_i}。

这就验证了 x=\sum\limits^n_{i=1} a_iM_it_i 是方程组的一组特解。

运用小学知识（其实就是我懒得证明）可得每两组相邻特解之间相差为 M 的倍数，于是 x+kM 也是一个特解，由此，通解唯一性得证。

例题

P1495 【模板】中国剩余定理（CRT）/ 曹冲养猪

模板题，直接上代码。

#include<bits/stdc++.h>
#define XD 114514
#define MAXN 20
#define int long long
using namespace std;
int n;
void gg(int &x,int &y,int a,int b){//求逆元
    if(a%b==0){
        x=0,y=1;
        return;
    }
    gg(y,x,b,a%b);
    y-=a/b*x;
}
int a[MAXN],b[MAXN],M[MAXN],m=1,t[MAXN],ans;
signed main(){
    cin>>n;
    for(int i=1;i<=n;i++){
        cin>>a[i]>>b[i];//x=ai*k+bi -> x mod ai=bi mod ai
        m*=a[i];
    }
    for(int i=1;i<=n;i++){
        M[i]=m/a[i];
        int x=0,y=0;
        gg(x,y,M[i],a[i]);
        t[i]=(x%a[i]+a[i])%a[i];
    }
    for(int i=1;i<=n;i++){
        ans+=b[i]*t[i]%m*M[i]%m;
        ans%=m;
    }
    cout<<ans;
    return 0;
}

P3868 [TJOI2009] 猜数字

然后套模板即可。

拓展中国剩余定理

等我学完再写。QWQ

终于学完力。QWQ（on 2023.5.23）

求

就可以用 exgcd 求出 b_1h_1+b_2h_2=g 的解。