标签归档:逆元

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;
}

NOIP2011 -计算系数

题目链接:CH3601

通过逆元计算组合数.

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 quickpow(int m,int n,int mod){
  int res=1;
  while(n){
    if(n&1) res=(ll)res*m%mod;
    m=(ll)m*m%mod;
    n>>=1;
  }
  return res;
}

int factorial(int n,int mod){
      if(!n) return 1;
      int res=n%mod;
      while(--n){
        res=(ll)res*n%mod;
      }
      return res;
}

int main() {
  int a,b,n,m,k;
  const int mod=10007;
  cin>>a>>b>>k>>n>>m;
  a=quickpow(a,n,mod);
  b=quickpow(b,m,mod);
  int x=factorial(k,mod);
  int y1=factorial(n,mod);
  int y2=factorial(k-n,mod);
  y1=quickpow(y1,mod-2,mod);
  y2=quickpow(y2,mod-2,mod);
  int ans=(ll)y1*y2%mod*x%mod*a%mod*b%mod;
  cout<<ans<<endl;
  return 0;
}

POJ 1845- Sumdiv


题目链接:POJ1845

题意:求\(A^B\)的所有约数之和\(\mod 9901\), \(1 \le A,B \le 5*10^7 \).

把A分解质因数,可以得到约数和为:
\((1+p_1+p_1^2+...+p_1^{B*c_1})*(1+p_2+p_2^2+...+p_2^{B*c_2})*...*(1+p_n+p_n^2+...+p_n^{B*c_n})\)

其中每一项都是等比数列求和:\( \frac{p_i^{B*c_i+1}-1}{p_i-1}\),如果分母逆元存在,则用乘以逆元代替除法运算,如果不存在,有\(p_i\mod 9901=1\),则:
\( 1+p_i+p_i^2+...+p_i^{B*c_i}\equiv 1+1+1^2+...+1^{B*c_i}\equiv B*c_i+1 \pmod {9901} \).

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 c[20],p[20];
int m;
void divide(int a){
    m=0;
    for(int i=2;i<=sqrt(a);++i){
      if(a%i==0){
        p[++m]=i;
        c[m]=0;
        while(a%i==0){
          a/=i;
          c[m]++;
        }
      }
    }
    if(a>1){
      p[++m]=a;
      c[m]=1;
    }
}

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

int inv(int mod,int x){
  return quickpow(mod,x,mod-2);
}

int main(){
  int a,b;
  const int mod=9901;
  cin>>a>>b;
  divide(a);
  int ans=1;
  for(int i=1;i<=m;++i){
    int x=(quickpow(mod,p[i],(ll)b*c[i]+1)-1+mod)%mod;
    int y=(p[i]-1)%mod;
    if(y){
      ans=(ll)ans*(x*inv(mod,y))%mod;
    }
    else{
      ans=ans*((ll)b*c[i]+1)%mod;  
    }
  }
  cout<<ans<<endl;
  return 0;
}