1. 程式人生 > >[筆記] CRT & exCRT

[筆記] CRT & exCRT

ref play etc 觀察 怎麽 def 答案 lock printf

[筆記] CRT & exCRT

構造法

求多組\(x \equiv r_i (\bmod d_i)\)的解,\(d_i\)互質

余數\((r_i = remainder)\),除數\((d_i=divisor)\)

我們想啊,如果我們能找到一個數 \(k1\equiv1(mod\text{ }3)\)\(5\)\(7\) 的倍數

一個數 $k2\equiv1(mod\text{ }5) $是\(3\)\(7\)的倍數

一個數 $k3\equiv1(mod\text{ }7) $是\(3\)\(5\)的倍數

那麽這樣的話我們的答案就是 \(k1*2+k2*3+k3*2\)

這樣的話就是如何求這三個數了

首先我們求出 \(3,5,7\) 的$ lcm=105$

之後令

\(x1=105/3=35,x2=105/5=21,x3=105/7=15\)

我們可以得到三個線性同余方程

\(35a\equiv1(mod\text{ }3)\)

\(21b\equiv1(mod\text{ }5)\)

\(15c\equiv1(mod\text{ }7)\)

直接上擴歐就行了

所以$ k1=35a=35*2=70$

\(k2=21b=21?1=21\)

\(k3=15c=15?1=15\)

之後我們的答案就是

\((k1?2+k2?3+k3?2)\bmod lcm=23\)

曹沖養豬

板子from wzx(我的構造法過不了,不過有了擴歐也不想用這個方法了)

#include<iostream>
#include<cstdio>
#include<cstring>
#define LL long long
#define re register
#define maxn 25
using namespace std;
inline LL exgcd(LL a,LL b,LL &x,LL &y)
{
    if(!b)
    {
        x=1,y=0;
        return a;
    }
    LL r=exgcd(b,a%b,y,x);
    y-=a/b*x;
    return r;
}
inline LL read()
{
    char c=getchar();
    LL x=0;
    while(c<'0'||c>'9') c=getchar();
    while(c>='0'&&c<='9')
      x=(x<<3)+(x<<1)+c-48,c=getchar();
    return x;
}
LL a[maxn],b[maxn];
LL x,y;
int n;
LL lcm=1;
int main()
{
    n=read();
    for(re int i=1;i<=n;i++)
    {
        b[i]=read();
        a[i]=read();
        lcm*=b[i];
    }
    LL ans=0;
    for(re int i=1;i<=n;i++)
    {
        LL xx=lcm/b[i];
        LL r=exgcd(xx,b[i],x,y);
        x=(x%b[i]+b[i])%b[i];
        ans=(ans+x*xx*a[i]%lcm)%lcm;
    }
    cout<<ans<<endl;
    return 0;
}

擴歐法

\[\begin{cases} x &\equiv& r_1 \pmod{d_1}\x & \equiv& r_2 \pmod{d_2}\&\vdots\x & \equiv& r_n \pmod{d_n} \end{cases}\]

此時如果\(d_1,d_2,\dots,d_n\)不互質怎麽辦?

觀察兩個同余式

  • $x\equiv r_1 \pmod{d_1} $
  • \(x\equiv r_2 \pmod{d_2}\)

可以化為

  • \(x=k_1d_1 + r_1\)
  • \(x=k_2d_2+r_2\)

因此

  • \(k_2d_2+r_2=k_1d_1+r_1\)
  • \(k_2d_2-k_1d_1=r_1-r_2\)

是不是有點像\(ax+bx=c\)

於是我們可以解出來

  • \(k_2 ^{'} d_2+k_1 ^{'} d_1=gcd(d_1,d_2)\)

並乘上\((r_1-r_2)/gcd(d_1,d_2)\)得到之前式子的解

  • \(k_1=k_1 ^{'} \cdot (r_1-r_2)/gcd(d_1,d_2)\)

但是得\(gcd(d_1,d_2)|(r_1-r_2)\),不然是得不到通解的

因為是帶著一個負號的,所以得到

  • \(x_0= -k_1d_1+r_1\)

得到新方程\(x\equiv x_0 \pmod{ lcm(d_1,d_2)}\)

好的,開始,摳細節

  • 首先,最後要取模\(lcm(d1,d2)\),快速乘的時候忘記了,掛掉x1
  • 快速乘裏面不能有負數,掛掉x2
  • 快速乘也要取模\(lcm(d_1,d_2)\),掛掉x3
  • \(lcm(d_1,d_2)\)應該這麽寫,\(d_1 /gcd(d_1,d_2)\cdot d_2\),而不是\(d_1 \cdot d_2 / gcd(d_1,d_2)\),掛掉x4
  • 把除號放裏面,不然有逆元,(其實是因為模了一個不該模的東西)掛掉x5

代碼

#include <iostream>
#include <cstdio>
#include <algorithm>

using namespace std;

typedef long long ll;

const int N = 1e5 + 10;

ll n, x, y;
ll d[N], r[N];

ll mul(ll a, ll b, ll p){//龜速乘
    int f = 1;
    if(a < 0) f = -f, a = -a;
    if(b < 0) f = -f, b = -b;
    ll w = 0;
    while(b){
        if(b & 1)
            w = (w + a) % p;
        b >>= 1;
        a = (a + a) % p;
    }
    return w * f;
}

ll exgcd(ll a, ll b){//擴歐
    ll ans, t;
    if(b == 0){
        x = 1;
        y = 0;
        return a;
    } else {
        ans = exgcd(b, a % b);
        t = x;
        x = y;
        y = t - a / b * y;
    }
        
    return ans;
}

ll exCRT(){
    
    for(int i = 2; i <= n; ++i){
        ll C = r[1] - r[i];
        ll D = exgcd(d[i], d[1]);
        if(C % D) return -1;
        ll k1 = mul(y , c / D, d[1] / D * d[i] );
        ll x0 = mul(-k1 , d[1], d[1] / D * d[i] ) + r[1];
        d[1] = d[1] / D * d[i], r[1] = x0;
        r[1] = (r[1] % d[1] + d[1]) % d[1];//先模一遍,讓絕對值小於,然後處理
    }
    return r[1];
}

int main(){
    scanf("%lld", &n);
    for(int i = 1; i <= n; ++i){
        scanf("%lld %lld", &d[i], &r[i]);
    }
    long long ans = exCRT();
    printf("%lld", ans);
    return 0;
} 

我知道寫的超級醜,比不上fym

[筆記] CRT & exCRT