线性同余方程和矩阵乘法

钊气蓬勃.

已于 2022-10-12 16:08:53 修改

阅读量444

点赞数

文章标签：矩阵算法 c++

于 2022-10-11 17:47:48 首次发布

本文链接：https://blog.csdn.net/m0_63729880/article/details/127260815

版权

线性同余方程

用扩展欧几里得算法来求线性同余方程

ax+by=1

则解除的x,y有多组解为

x=x0+k*b

y=y0-k*a 其中x0和y0是解出来的解,k是整数

则最小正整数解为x0%b,y0%a

扩展欧几里得模板

int exgcd(int a,int b,int &x,int &y)//拓展欧几里得求ax+by=d的方程
{
    if(!b)
    {
        x=1,y=0;
        return a;
    }
    int d=exgcd(b,a%b,y,x);
    y-=a/b*x;
    return d;
}

1.同余方程

203. 同余方程 - AcWing题库

ax=1(mod b)转化为ax+by=1，则直接用拓展欧几里得算法即可求x,y，最小x解为x0%b

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=1e7+10;
int prime[N],cnt;
bool st[N];
int phi[N];
ll s[N];//记录欧拉函数的前缀和
int exgcd(int a,int b,int &x,int &y)//拓展欧几里得求ax+by=d的方程
{
    if(!b)
    {
        x=1,y=0;
        return a;
    }
    int d=exgcd(b,a%b,y,x);
    y-=a/b*x;
    return d;
}
int main()
{
   int a,b;
   cin>>a>>b;
   int x,y;
   exgcd(a,b,x,y);//求x,y
   cout<<((x%b)+b)%b<<endl;//输出mo出来是正的数
    return 0;
}

2.青蛙的约会

222. 青蛙的约会 - AcWing题库

两个青蛙相距b-a米，而每次能缩小跳m-n米，问至少跳多少次使得两个青蛙碰面

设跳了t次，则跳过的距离为t*(m-n)米，又因为是个圈，则设绕了y圈，若相遇则总的米数为

t*(m-n)=(b-a)+y*l，则直接用扩展欧几里得可得t0，则最小的t为t0%l/d

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll exgcd(ll a,ll b,ll &x,ll &y)//拓展欧几里得求ax+by=d的方程
{
    if(!b)
    {
        x=1,y=0;
        return a;
    }
    ll d=exgcd(b,a%b,y,x);
    y-=a/b*x;
    return d;
}
int main()
{
   ll a,b,n,m,l;
   cin>>a>>b>>m>>n>>l;
   ll x,y;
   ll d=exgcd(m-n,l,x,y);
   if((b-a)%d) puts("Impossible");//假如b-a不是公约数
   else
   {
       x*=(b-a)/d;//将x扩大相应的倍数
       ll t=abs(l/d);//将d除过来，让右边为1，才能用t0%b为最小这个定里
       cout<<((x%t)+t)%t<<endl;//输出最小结果
   }
    return 0;
}

3.最幸运的数字

202. 最幸运的数字 - AcWing题库

龟速乘解决乘法爆longlong的情况

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll get_euler(ll c)//获取一个数的欧拉数
{
    ll res=c;
    for(ll i=2;i<=c/i;i++)
        if(c%i==0)
        {
            while(c%i==0) c/=i;
            res=res/i*(i-1);//套欧拉函数的公式
        }
    if(c>1) res=res/c*(c-1);
    return res;
}
ll qmul(ll a,ll k,ll p)//龟速乘
{
    ll res=0;
    while(k)
    {
        if(k&1) res=(res+a)%p;
        a=(a+a)%p;
        k>>=1;
    }
    return res;
}
ll qmi(ll a,ll k,ll p)//快速幂
{
    ll res=1;
    while(k)
    {
        //res=(__int128_t)res*a%p也行
        if(k&1) res=qmul(res,a,p);//因为乘法会爆ll,所以用龟速乘
        //a=(__int128_t)a*a%p也行
        a=qmul(a,a,p);//因为乘法会爆ll,所以用龟速乘
        k>>=1;
    }
    return res;
}
int main()
{
   int T=1;
   ll l;
   while(cin>>l,l)
   {
       int d=1;
       while(l%(d*2)==0&&d*2<=8) d*=2;
       ll c=9*l/d;//获取要整除的C
       ll phi=get_euler(c);//获取他的欧拉数
       ll res=1e18;
       if(c%2==0||c%5==0) res=0;//假如跟10不互质，则没答案
       for(ll d=1;d*d<=phi;d++)//枚举欧拉数的所有约数
            if(phi%d==0)//假如这个是约数了
            {
               if(qmi(10,d,c)==1) res=min(res,d);//判断符不符合条件
               if(qmi(10,phi/d,c)==1) res=min(res,phi/d);//判断符不符合条件
            }
     printf("Case %d: %lld\n",T++,res);
   }
    return 0;
}

__int128_t解决乘法爆longlong的情况

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll get_euler(ll c)
{
    ll res=c;
    for(ll i=2;i<=c/i;i++)
        if(c%i==0)
        {
            while(c%i==0) c/=i;
            res=res/i*(i-1);
        }
    if(c>1) res=res/c*(c-1);
    return res;
}
ll qmi(ll a,ll k,ll p)
{
    ll res=1;
    while(k)
    {
        if(k&1) res=(__int128_t)res*a%p;
        a=(__int128_t)a*a%p;
        k>>=1;
    }
    return res;
}
int main()
{
   int T=1;
   ll l;
   while(cin>>l,l)
   {
       int d=1;
       while(l%(d*2)==0&&d*2<=8) d*=2;
       ll c=9*l/d;
       ll phi=get_euler(c);
       ll res=1e18;
       if(c%2==0||c%5==0) res=0;
       for(ll d=1;d*d<=phi;d++)
            if(phi%d==0)
            {
               if(qmi(10,d,c)==1) res=min(res,d);
               if(qmi(10,phi/d,c)==1) res=min(res,phi/d);
            }
     printf("Case %d: %lld\n",T++,res);
   }
    return 0;
}

4.曹冲养猪

信息学奥赛一本通（C++版）在线评测系统 (ssoier.cn)

中国剩余定理

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=11;
int A[N],B[N];
int n;
ll exgcd(ll a,ll b,ll &x,ll &y)//用扩展欧几里得求逆元
{
    if(!b)
    {
        x=1,y=0;
        return a;
    }
    ll d=exgcd(b,a%b,y,x);
    y-=a/b*x;
    return d;
}
int main()
{
   cin>>n;
   ll M=1,res=0;
   for(int i=0;i<n;i++)
   {
       scanf("%d%d",&A[i],&B[i]);
       M*=A[i];//处理所有的乘积
   }
   for(int i=0;i<n;i++)
   {
       ll Mi=M/A[i];//获取Mi
       ll ti,x;
       exgcd(Mi,A[i],ti,x);//求Mi的逆元ti
       res+=B[i]*Mi*ti;//图片中有讲答案加上即可
   }
   cout<<((res%M)+M)%M<<endl;//输出正的模余数
    return 0;
}

矩阵乘法

1.斐波那契前 n 项和

1303. 斐波那契前 n 项和 - AcWing题库

信息学奥赛一本通（C++版）在线评测系统 (ssoier.cn)

因为数很大不能直接用公式fn=fn-1+fn-2求，得用矩阵乘法来求

定义一个矩阵Fn=[fn,fn+1]

这题构造一个Fn=【fn,fn+1,sn】的矩阵，是前 fn项和，则A矩阵如下：

然后A的n-1次方用快速幂求出来就行了，logn的时间

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=3;
int n,m;
void mul(int c[],int a[],int b[][N])
{
    int temp[N]={0};
    for(int i=0;i<N;i++)//矩阵乘法
        for(int j=0;j<N;j++)
          temp[i]=(temp[i]+(ll)a[j]*b[j][i])%m;
    memcpy(c,temp,sizeof temp);
}
void mul(int c[][N],int a[][N],int b[][N])
{
    int temp[N][N]={0};
    for(int i=0;i<N;i++)//矩阵乘法
        for(int j=0;j<N;j++)
          for(int k=0;k<N;k++)
              temp[i][j]=(temp[i][j]+(ll)a[i][k]*b[k][j])%m;
    memcpy(c,temp,sizeof temp);
}
int main()
{
   cin>>n>>m;
   int f1[N]={1,1,1};//f1为1时的矩阵
   int a[N][N]={
    {0,1,0},
    {1,1,1},
    {0,0,1}
   };//a矩阵
   n--;//求Fn-1
   while(n)
   {
       if(n&1) mul(f1,f1,a);//res=res*a
       mul(a,a,a);//a=a*a;
       n>>=1;
   }
   cout<<f1[2]<<endl;//输出Fn-1的中的f(n-1)+1=fn
    return 0;
}