SCL--CRT（中国剩余定理）

2015-09-18 22:44:20

总结：复习了一下CRT，整理出这个板子。

核心方程（1）：，几个同余方程模数非互素，可以计算出所有模数的所有素因子的最高幂形 p^k，然后CRT。

核心方程（2）：    

　　有整数解。并且在模下的解是唯一的，解为

        

　 　其中，而为模的逆元。（这部分转自ACdreamer‘s blog）

1：m1，m2 ... mk 互质的情况下

inline ll mulmod(ll x,ll y,ll mod){
    ll res = 0;
    __asm__("movq %1,%%rax\n imulq %2\n idivq %3\n":"=d"(res):"m"(x),"m"(y),
            "m"(mod):"%rax");
    return res;
}

ll Ex_gcd(ll a,ll b,ll &x,ll &y){
    if(b == 0){
        x = 1;
        y = 0;
        return a;
    }
    ll r = Ex_gcd(b,a % b,y,x);
    y -= x * (a / b);
    return r;
}

ll CRT(ll a[],ll m[],ll n){
    ll M = 1;
    for(int i = 1; i <= n; ++i) M *= m[i];
    ll res = 0;
    for(int i = 1; i <= n; ++i){
        ll x,y;
        ll Mi = M / m[i];
        Ex_gcd(Mi,m[i],x,y);
        ll tmp = mulmod(mulmod(Mi,x,M),a[i],M);
        res = (res + tmp) % M;
    }
    return (res + M) % M;
}

 

2：m1，m2 ... mk 不互质的情况下，转载自http://yzmduncan.iteye.com/blog/1323599/

 1 ll Mod;
 2 
 3 ll gcd(ll a,ll b)  {
 4     return b == 0 ? a : gcd(b,a % b);  
 5 }  
 6 
 7 ll Extend_Euclid(ll a,ll b,ll &x,ll &y){  
 8     if(b == 0){
 9         x = 1,y = 0;
10         return a;  
11     }  
12     ll d = Extend_Euclid(b,a % b,x,y);  
13     ll t = x;  
14     x = y;  
15     y = t - a / b * y;  
16     return d;  
17 }  
18 
19 //a在模n乘法下的逆元，没有则返回-1  
20 ll inv(ll a,ll n)  {
21     ll x,y;  
22     ll t = Extend_Euclid(a,n,x,y);  
23     if(t != 1) return -1;
24     return (x % n + n) % n;
25 }  
26 
27 //将两个方程合并为一个  
28 bool merge(ll a1,ll n1,ll a2,ll n2,ll &a3,ll &n3)  {
29     ll d = gcd(n1,n2);  
30     ll c = a2 - a1;  
31     if(c % d) return false;  
32     c = (c % n2 + n2) % n2;  
33     c /= d; 
34     n1 /= d;  
35     n2 /= d;  
36     c *= inv(n1,n2);  
37     c %= n2;  
38     c *= n1 * d;  
39     c += a1;  
40     n3 = n1 * n2 * d;  
41     a3 = (c % n3 + n3) % n3;  
42     return true;  
43 }  
44 
45 //求模线性方程组x=ai(mod ni),ni可以不互质  
46 ll China_Reminder2(int len,ll *a,ll *n)  {
47     ll a1 = a[1],n1 = n[1];
48     ll a2,n2;  
49     for(int i = 2; i <= len; i++){
50         ll aa,nn;  
51         a2 = a[i],n2=n[i];  
52         if(!merge(a1,n1,a2,n2,aa,nn)) return -1;  
53         a1 = aa;  
54         n1 = nn;  
55     }
56     Mod = n1;
57     return (a1 % n1 + n1) % n1;
58 }
59 
60 ll a[1000],b[1000];
61 
62 int main(){
63     .....
64     ll ans = China_Reminder2(N,a,b);
65     return 0;
66 }

 

posted @ 2015-09-18 22:55  Naturain  阅读(229)  评论(0)  编辑  收藏  举报