模板库(二) - 数论算法模板

最新推荐文章于 2022-03-31 14:49:34 发布

原创最新推荐文章于 2022-03-31 14:49:34 发布 · 284 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#代码 #模板 #NOIP

模板专栏收录该内容

7 篇文章

订阅专栏

写在前面

“模板库”这一系列文章用来复习 $O I$ 模板
由于时间原因，作者无法一一亲自调试其中的程序，也因如此，有一部分程序来自于互联网，如果您觉得这侵犯了您的合法权益，请联系 $(Q Q 2068926345)$ 删除。
对于给您造成的不便和困扰，我表示深深的歉意。
本系列文章仅用于学习，禁止任何人或组织用于商业用途。
本系列文章中，标记*的为选学算法，在 $N O I P$ 中较少涉及。

数论算法模板

快速幂

【简介】

顾名思义，快速幂就是快速算底数的次幂。其时间复杂度为 $Θ (l o g n)$ ，与朴素的相比效率有了极大的提高。

【代码实现】

#include<cstdio>
int b,p,mod;
inline int Fast_pow(long long b,int p){
    long long ans=1;
    while(p){
        if(p&1) ans=ans*b%mod;
        b=b*b%mod;p>>=1;
    }
    return ans;
}
int main(){
    scanf("%d%d%d",&b,&p,&mod);
    printf("%d^%d mod %d=%d",b,p,mod,Fast_pow(b,p));
    return 0; 
}

复杂度 $Θ (l o g n)$

欧几里得算法

【简介】

欧几里德算法又称辗转相除法，是指用于计算两个正整数的最大公约数。应用领域有数学和计算机两个方面。计算公式。

【代码实现】

#include<cstdio>
int gcd(int a,int b){return b?gcd(b,a%b):a;}
int main(){
    int a,b;
    scanf("%d%d",&a,&b);
    printf("%d",gcd(a,b));
}

复杂度 $Θ (l o g n)$

扩展欧几里得算法

【简介】

扩展欧几里德算法是用来在已知求解一组，使它们满足贝祖等式：（解一定存在，根据数论中的相关定理）。扩展欧几里德常用在求解模线性方程及方程组中。

【代码实现】

#include<cstdio> 
int exgcd(int a,int b,int &x,int &y){
    if(!b)x=1,y=0;
    else exgcd(b,a%b,y,x),y-=x*(a/b);
}
int main(){
    int a,b,x,y;
    scanf("%d%d",&a,&b);
    exgcd(a,b,x,y);
    printf("%d",(x%b+b)%b);
    return 0;
}

复杂度 $Θ (l o g n)$

逆元

【简介】

乘法逆元，是指数学领域群中任意一个元素，都在中有唯一的逆元，具有性质，其中为该群的单位元。

【代码实现】

方法1：求单个逆元

#include<cstdio>
int n,mod;
inline int Fast_pow(long long b,int p){
    long long ans=1;
    while(p){
        if(p&1) ans=ans*b%mod;
        b=b*b%mod;p>>=1;
    }
    return ans;
}
int main(){
    scanf("%d%d",&n,&mod);
    printf("%d",Fast_pow(n,mod-2);
    return 0; 
}

复杂度 $Θ (l o g n)$

方法2: 线性求逆元

#include<cstdio>
long long a[100007]={0,1};
int main(){
    int n,p;
    scanf("%d%d",&n,&p);
    for(int i=2;i<=n;i++) a[i]=(p-(p/i))*a[p%i]%p;
    for(int i=1;i<=n;i++) printf("%lld\n",a[i]);
    return 0;
}

复杂度 $Θ (n)$

筛法求积性函数

【简介】

筛法是一种简单检定素数的算法。据说是古希腊的埃拉托斯特尼 $(E r a t o s t h e n e s$ ，约公元前 $274 ～ 194$ 年)发明的，又称埃拉托斯特尼筛法 $(s i e v e$ $o f$ $E r a t o s t h e n e s)$ ，后来扩展到利用筛法求积性函数。

【一些性质】

积性函数：对于函数 $f (n)$ ，若满足对任意互质的数字 $a$ ， $b$ ， $a * b = n$ ，且 $f (n) = f (a) * f (b)$ ，那么称函数 $f$ 为积性函数。
狄利克雷卷积：对于函数 $f$ ， $g$ ，定义它们的卷积为 $(f∗g)(n)=∑d∣nf(d)g(nd)(f*g)(n)=\sum_{d|n}f(d)g(nd)$ 。

两个积性函数的狄利克雷卷积仍为积性函数。
积性函数都可以用线性筛筛出来

【代码实现】

方法1

#include<cstdio>
#include<cmath>
#include<cstring>
using namespace std;
int is_not[10000007]={1,1};
void prime(int n){
    int m=sqrt(n);
    for(int i=2;i<=m;i++)
        if(!shu[i])for(int j=i*i;j<=n;j+=i)shu[j]=1;	
}
int main(){
    int n,m;
    scanf("%d%d",&n,&m);
    prime(N);
    for(int i=1;i<=m;i++)scanf("%d",&a),printf(shu[a]?"No\n":"Yes\n");
    return 0;
}

复杂度 $Θ ($ $n$ $l o g$ $l o g$ $n)$

方法2

#include<cstdio>
bool is_not[10000007]={1,1};
int prime[664579+7],n,m,a,cnt;
void Prime(int n){
    for(int i=2;i<=n;++i){
        if(!is_not[i])prime[++cnt]=i;
        for(int j=1;j<=cnt && prime[j]*i<=n;++j){
            is_not[i*prime[j]]=1;
            if(i%prime[j]==0)break;
        }
    }
}
int main(){
    scanf("%d%d",&n,&m);
    Prime(n);
    for(int i=1;i<=m;++i)scanf("%d",&a),printf(is_not[a]?"No\n":"Yes\n");
    return 0;
}

复杂度 $Θ (n)$

【利用筛法求一些常见的积性函数】

复杂度 $Θ (n)$

莫比乌斯函数

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

欧拉函数

void getphi(){
    phi[1]=1;
    for(int i=2;i<=N;++i){
        if(!is_not[i])prime[++cnt]=i,phi[i]=i-1;
        for(int j=1;j<=cnt,i*prime[j]<=N;++j){
            is_not[i*prime[j]]=1;
            if(i%prime[j]==0){
                phi[i*prime[j]]=phi[i]*prime[j];break;
            }
            else phi[i*prime[j]]=phi[i]*(prime[j]-1);
        }
    }
}

约数个数

is_not[1]=1,d[1]=1;
for(int i=2;i<=N;++i){
    if(!is_not[i]){
        prime[++cnt]=i;
        d[i]=2;pred[i]=1;
    }
    for(int j=1;j<=num&&i*prime[j]<=N;++j){
        is_not[i*prime[j]]=1;
        if(i%prime[j])d[i*prime[j]]=d[i]*d[prime[j]],pred[i*prime[j]]=1;
        else{
            pred[i*prime[j]]=pred[i]+1;
            d[i*prime[j]]=d[i]/(pred[i]+1)*(pred[i]+2);
            break;
        }
    }
}

约数的和

void Prepare(){
    is_not[1]=1;f[1]=mu[1]=1;
    for(int i=2;i<=N;++i){
        if(!isprime[i]){
            prime[++cnt]=i;f[i]=i+1;mu[i]=-1;
            sumd[i]=1+i;powd[i]=i;
        }
        for(int j=1;j<=num&&i*prime[j]<N;++j){
            is_not[i*prime[j]]=1;
            if(i%prime[j]){
                sumd[i*prime[j]]=prime[j]+1;
                powd[i*prime[j]]=prime[j];
                f[i*prime[j]]=f[i]*f[prime[j]];
            }
            else{
                powd[i*prime[j]]=powd[i]*prime[j];
                sumd[i*prime[j]]=sumd[i]+powd[i*prime[j]];
                f[i*prime[j]]=f[i]/sumd[i]*sumd[i*prime[j]];
                break;
            }
        }
    }
}

*中国剩余定理

【简介】

中国剩余定理是中国古代求解一次同余式组的方法。是数论中一个重要定理。又称中国余数定理。一元线性同余方程组问题最早可见于中国南北朝时期 $($ 公元5世纪 $)$ 的数学著作《孙子算经》卷下第二十六题，叫做“物不知数”问题。
在 $O I$ 中，中国剩余定理主要解决以下问题：

$mn)\begin{cases} x≡a_1(mod\ m_1)\\ x≡a_2(mod\ m_2)\\ x≡a_3(mod\ m_3)\\ \ \ \cdots\\ x≡a_n(mod\ m_n) \end{cases}$

【代码实现】

#include<cstdio>
typedef long long ll;
using namespace std;
ll a[17],b[17],n;
ll exgcd(ll a,ll b,ll &x,ll &y) {
    if(!b){x=1;y=0;return a;}
    long long q=exgcd(b,a%b,y,x);
    y-=a/b*x;
    return q;
}
ll china() {
    ll M=1,res=0;
    for(int i=1;i<=n;i++) M*=a[i];
    for(int i=1;i<=n;i++){
        ll m=M/a[i],p,y;
        ll d=exgcd(m,a[i],p,y);
        res=(res+m*p*b[i])%M;
    }
    if(res<0) res+=M;
    return res;
} 
int main(){
    scanf("%d",&n);
    for(int i=1;i<=n;i++) scanf("%lld%lld",&a[i],&b[i]);
    printf("%lld",china());
}

复杂度 $Θ (n$ $l o g n)$

*BSGS

【简介】

大步小步法 $(B a b y - S t e p - G i a n t - S t e p ，简称 B S G S)$ ，可以较高效的求解形如 $A^x≡B(mod$ $C) (C$ 是素数 $)$ 的同余方程。

【代码实现】

#include<cstdio>
#include<map>
#include<cmath>
using namespace std;
__int128 FastPow(__int128 a,__int128 k,__int128 p){
    __int128 ans=1;
    while(k)
    {
        if(k&1) (ans*=a)%=p;
        (a*=a)%=p;
        k>>=1;
    }
    return ans;
}
__int128 BSGS(__int128 a,__int128 b,__int128 p)
{
    map<__int128,__int128>hash;
    hash.clear();
    b%=p;
    __int128 t=(__int128)sqrt((double)p)+1;
    for(__int128 j=0;j<t;j++)
    {
        __int128 val=b*FastPow(a,j,p)%p;
        hash[val]=j;
    }
    a=FastPow(a,t,p);
    if(!a)return(b==0)?1:-1;
    for(__int128 i=0;i<=t;i++)
    {
        __int128 val=FastPow(a,i,p);
        __int128 j=hash.find(val)==hash.end()?-1:hash[val];
        if(j>=0&&i*t-j>=0) return i*t-j;
    }
    return -1;
}
__int128 k,m;
int main(){
    scanf("%lld%lld",&k,&m);
    __int128 ans=BSGS(10,(9*k+1)%m,m);
    printf("%lld",ans);
    return 0;
}

复杂度 $Θ(C)Θ(\sqrt C)$

除法分块

【简介】

求解类似于 $∑i=1nf(i)[ni]\sum_{i=1}^{n}f(i)[\frac{n}{i}]$ 的算式。其中， $f (i)$ 应可以快速计算前缀和。

【代码实现】

#include<cstdio>
int main(){
    long long n,k,ans=0;
    scanf("%lld%lld",&n,&k);
    for(long long l=1,r,t;l<=n;l=r+1)
        r=(t=k/l)?std::min(k/t,n):n,
        ans+=t*(r-l+1)*(r+l)>>1;
    printf("%lld",n*k-ans);
    return 0;
}

复杂度 $Θ(n)Θ(\sqrt n)$

*Pollard-Rho算法

【简介】

有一类问题，要求我们将一个正整数，分解为两个非平凡因子的乘积 $x = a b$ 。
显然我们需要先检测是否为素数如果是素数将无解，可以使用 $M i l l e r - R a b i n$ 算法来进行测试。
$P o l l a r d - R h o$ 是一个非常玄学的方式，用于在 $O(n4)O(\sqrt[4]{n})$ 的期望时间复杂度内计算合数 $n$ 的某个非平凡因子。事实上算法导论给出的是 $O(p)O(\sqrt p)$ ， $p$ 是 $n$ 的某个最小因子，满足 $p$ 与 $np\frac n p$ 互质。但是这些都是期望，未必符合实际。但事实上 $P o l l a r d - R h o$ 算法在实际环境中运行的相当不错。

【代码实现】

#include <cstdio>
#include <algorithm>
#define rep(i, s, t) for(int i = s; i <= t; ++i)
typedef long long ll;
ll H;
ll pls(ll a, ll b, ll p){
    ll res=0;
    for(; b; b>>=1, a=(a+a)%p)
        if(b & 1) res = (res+a)%p;
    return res%p;
}
ll pow(ll a, ll b, ll p){
    ll res = 1;
    for(; b; b>>=1, a=pls(a, a, p))
        if(b&1) res=pls(res, a, p)%p;
    return res % p;
}
bool M(ll p){
    ll x[60] = {0}, s = 20;
    ll rest = p-1, t = 0;
    while(!(rest%2)){
        t++;
        rest >>= 1;
    }
    while(s--){
        ll a = rand()%(p-1)+1;
        x[0] = pow(a, rest, p);
        rep(i, 1, t){
            x[i] = pow(x[i-1], 2, p)%p;
            if(x[i] == 1)
                if((x[i-1] != 1) && (x[i-1] != p-1)) return false;
        }
        if(x[t] ^ 1) return false;
    }
    return true;
}
ll gcd(ll a, ll b){
    while(b){
        ll t = a%b;
        a = b;b = t;
    }
    return a;
}
ll P(ll p){
    ll c = rand()%(p-1)+1, x1 = rand()%(p-1)+1, x2 = x1;
    for(ll i = 2, k = 2; true; ++i){
        x1 = (pls(x1, x1, p) + c)%p;
        ll G = gcd(p, (x2-x1+p)%p);
        if(G > 1 && G < p) return G;
        if(x2 == x1) return p;
        if(i == k) x2 = x1, k <<= 1;
    }
}
void solve(long long n){
    if(n == 1) return;
    if(M(n)){
        H = std::min(H , n);
        return;
    }
    long long p = n;
    while(p == n) p = P(p);
    solve(p);
    solve(n/p);
}
int main(){
    int _;
    scanf("%d", &_);
    while(_--){
        ll p;
        H = 1LL << 54;
        scanf("%lld", &p);
        solve(p);
        if(H ^ p)
            printf("%lld\n", H);
        else puts("Prime");
    }
    return 0;
}

复杂度 $Θ(n4)Θ(\sqrt[4]{n})$