标签归档：EXGCD

BZOJ 1951-古代猪文 (好题,Lucas+CRT+逆元+EXGCD+快速幂)

发表评论

题目链接:BZOJ 1951

通过欧拉定理的推论得到组合数部分的模数后,用Lucas定理计算,但是模数不是质数,对模数分解质因数后发现所有质因子的指数都是1,于是用中国剩余定理合并方程.具体思路可见《算法竞赛进阶指南》172页.

调了好久,才发现是快速幂的参数写错了.不多说了,直接上代码吧...心态崩了

code:

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <map>
#include <string>
#include <vector>

using namespace std;

typedef long long ll;
typedef unsigned long long ull;

int p[4] = {2, 3, 4679, 35617};

int exgcd(int a, int b, int &x, int &y) {
  if (!b) {
    x = 1;
    y = 0;
    return a;
  }
  int d = exgcd(b, a % b, x, y);
  int z = x;
  x = y;
  y = z - (a / b) * y;
  return d;
}

int quickpow(int n, int m, int mod) {
  int res = 1;
  while (m) {
    if (m & 1) res = (ll)res * n % mod;
    n = (ll)n * n % mod;
    m >>= 1;
  }
  return res;
}

int jc[4][36000];
int factorial() {
  for (int i = 0; i < 4; ++i) {
    jc[i][0] = 1;
    for (int j = 1; j < p[i]; ++j) {
      jc[i][j] = (ll)jc[i][j - 1] * j % p[i];
    }
  }
}

int inv(int a, int mod) {
  int x, y;
  exgcd(a, mod, x, y);
  return (x % mod + mod) % mod;
}

int lucas(int n, int m, int i) {
  if (n < m) return 0;
  if (n < p[i] && m < p[i])
    return (ll)jc[i][n] * inv(jc[i][n - m] * jc[i][m] % p[i], p[i]) % p[i];
  return (ll)lucas(n / p[i], m / p[i], i) * lucas(n % p[i], m % p[i], i) % p[i];
}

int CRT(int n, int a[], int mod) {
  int res = 0;
  for (int i = 0; i < n; ++i) {
    int M = mod / p[i];
    int t = inv(M, p[i]);
    res = (res + (ll)a[i] * M % mod * t % mod) % mod;
  }
  return res;
}

int A(int n, int i, int a[]) {
  for (int j = 0; j < 4; ++j) {
    a[j] += lucas(n, i, j);
  }
}

int main() {
  int n, g;
  cin >> n >> g;
  const int mod = 999911658;
  if (g == mod + 1) {
    cout << 0 << endl;
  } else {
    factorial();
    int a[4];
    memset(a, 0, sizeof(a));
    for (int i = 1; i * i <= n; ++i) {
      if (n % i == 0) {
        A(n, i, a);
        if (n / i != i) A(n, n / i, a);
      }
    }
    int ans = CRT(4, a, mod);
    ans = quickpow(g, ans, mod + 1);
    cout << ans << endl;
  }
  return 0;
}

POJ 2891 -Strange Way to Express Integers

发表评论

题目链接: POJ 2891

题意:给出\(n\)对\(a_i,r_i\), 求一个最小正整数\(m\)满足\(m \equiv r_i\pmod a_i\), 或者给出无解.

solution:

\(a_i\)并不保证两两互质, 所以不能用中国剩余定理. 但可以通过数学归纳法, 用拓展欧几里得逐步得出整个线性同余方程组的解.

假设求出前\(k-1\)个方程的一个解\(m\), 记 \(A=lcm( a_1,a_2,...,a_{k-1}) \) ,则 \(m+i*A\) 是通解. 求出满足第\(k\)个方程的\(i\), 就又求出了满足前\(k\)个方程的一个特解.

使用\(n\)次拓展欧几里得算法就能求出整个方程组的解.

code:

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <map>
#include <string>
#include <vector>

using namespace std;

typedef long long ll;
typedef unsigned long long ull;

ll exgcd(ll a, ll b, ll &x, ll &y) {
  if (!b) {
    x = 1;
    y = 0;
    return a;
  }
  ll d = exgcd(b, a % b, x, y);
  ll z = x;
  x = y;
  y = z - y * (a / b);
  return d;
}

ll gcd(ll a, ll b) {
  if (!b) return a;
  return gcd(b, a % b);
}

ll lcm(ll a, ll b) { return a / gcd(a, b) * b; }

int main() {
  int k;
  while (cin >> k) {
    ll a,r;
    cin>>a>>r;
    ll x, y;
    ll A=a;
    exgcd(1, a, x, y);
    x = x%A*r%A;
    int flag = 1;
    for (int i = 1; i < k; ++i) {
      cin>>a>>r;
      if(!flag) continue;
      ll t;
      ll d = exgcd(A, a, t, y);
      if ((r - x) % d != 0) {
        printf("-1\n");
        flag = 0;
        continue;
      }
      t = ( (t* ((r - x) / d)) % a+a )%a ;
      x =x+t*A;
      A = lcm(A, a);
      x=x%A;
    }
    if (flag) cout << x << endl;
  }
  return 0;
}

NOIP 2012-同余方程

发表评论

题目链接 :CH3301

用拓展欧几里得算法求出一组特解\(x_0,y_0\),则\(x_0\)就是原方程的一个解,通解为所有模\(b\)与\(x_0\)同余的整数.题目要求是最小正整数解,可通过取模操作把解的范围移动到\([1,b)\)之间,\([1,b)\)之间的就是最小正整数解.

code:

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <map>
#include <string>
#include <vector>

using namespace std;

typedef long long ll;
typedef unsigned long long ull;

int exgcd(int a,int b,int &x,int &y){
  if(b==0) {x=1;y=0;return a;}
  int d=exgcd(b,a%b,x,y);
  int z=x;x=y;y=z-y*(a/b);
  return d;
}

int main(){
  int a,b;
  cin>>a>>b;
  int x,y;
  exgcd(a,b,x,y);
  cout<<(x%b+b)%b<<endl; //注意x可能是负数
  return 0;
}