QOJ.ac

QOJ

ID题目提交者结果用时内存语言文件大小提交时间测评时间
#455068#3140. DivModulotimicsWA 0ms3620kbC++204.3kb2024-06-25 19:23:252024-06-25 19:23:26

Judging History

This is the latest submission verdict.

  • [2024-06-25 19:23:26]
  • Judged
  • Verdict: WA
  • Time: 0ms
  • Memory: 3620kb
  • [2024-06-25 19:23:25]
  • Submitted

answer

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef vector<int> vi;
typedef pair<int,int> pii;
typedef pair<double, double> pdd;
#define pb push_back
#define mp make_pair
#define fs first
#define sc second
#define rep(i, from, to) for (int i = from; i < (to); ++i)
#define all(x) x.begin(), x.end()
#define sz(x) (int)(x).size()
#define FOR(i, to) for (int i = 0; i < (to); ++i)
typedef vector<vector<int> > vvi;
typedef vector<ll> vll;
typedef vector<vll> vvll;
typedef vector<pair<int, int> > vpi;
typedef pair<ll,ll> pll;
typedef vector<string> vs;
#define Nmax 200200
#define MOD 1000000007
ll M,N,D;
vll a, m;

class Sieve {
public:
  vi fp; // first prime (i == fp[i] means it's prime itself)
  vi primes; // all primes <= N
  vi nump; // num unique primes

  Sieve(int N) {
    fp.resize(N+10);
    nump.resize(N+10);
    for(int i=2;i<=N;++i) {
      if(!fp[i]) {
        primes.pb(i);
        for(int j=i;1LL*i*j<=N;++j) {
          if(!fp[i*j]) fp[i*j] = i;
          ++nump[i*j];
        }
        fp[i] = i;
      }
    }
  }
  
  bool isPrime(int x) {
    return fp[x] == x;
  }

  vpi getPrimes(int x) { // prime decomp of x
    vpi ret;
    while(x > 1) {
      int p = fp[x], cnt = 0;
      while(x % p == 0) {
        x /= p; ++cnt;
      }
      ret.pb(mp(p, cnt));
    }
    return ret;
  }
};

ll fexp(ll N, ll p) {
  ll cur = 1;
  ll ret = 0;
  while(cur <= N/p) {
    cur *= p;
    ret += N/cur;
  }
  return ret;
}

ll pw(ll x, ll y, ll mod) {
  if(!y) return 1;
  if(y%2) return (x * pw(x, y-1, mod)) % mod;
  ll z = pw(x, y/2, mod);
  return (z*z) % mod;
}

ll divy(ll N, ll p, ll pe) {
  ll per = 1; // period
  for(ll i=1;i<pe;++i) {
    if(i%p) per = (per * i) % pe;
  }

  ll ret = pw(per, N/pe, pe);

  //cout<<N<<" "<<per<<" "<<N/pe<<" "<<N%pe<<endl;

  for(ll i=1;i<=N%pe;++i) { // last section
    if(i%p) ret = (ret * i) % pe;
  }
  if(N >= p) ret = (ret * divy(N/p, p, pe)) % pe;
  return ret;
}

// Euclid extins
ll gcd(ll a, ll b, ll &x, ll &y) {
  if(!b) { x=1; y=0; return a; } 
  else {
    ll x0, y0, d = gcd(b,a%b,x0,y0);
    x = y0; y = x0 - a/b *y0; return d;
  }
}
pll euclid(ll a, ll b, ll c) { //ax - by = c;
  ll x, y, sol1, sol2; ll d = gcd(a,b,x,y);
  if(c%d) return mp(0,0); //no sol
  else { sol1 = (c/d)*x; sol2 = -(c/d)*y; }
  //only if want minimal
  while(sol1 < 0 || sol2 < 0) { sol1 += b/d; sol2 += a/d;}
  while(sol1 >= b/d || sol2 >= a/d) { sol1 -= b/d; sol2 -= a/d;}
  return mp(sol1,sol2);
}

// Chinese remainder theorem x = a1 mod m1, x = a2 mod m2
pll crt(ll a1, ll m1, ll a2, ll m2) {
  ll d = gcd(m1,m2);
  pll p = euclid(m1,m2,d);
  ll s = p.fs, t = -p.sc;
  if (a1 % d != a2 % d) return make_pair(0, -1);
  ll x = (s * a2 * m1 + t * a1 * m2) % (m1*m2);
  x = (x + m1*m2) % (m1*m2);
  return mp(x/d, m1*m2/d);
}

// x = a[i] (mod m[i]) also returns period (lcm)
pll chinese(vll &a, vll &m) {
  pll ret = mp(a[0], m[0]);
  for (int i = 1; i < a.size(); ++i) {
    ret = crt(ret.fs, ret.sc, a[i], m[i]); if (ret.sc == -1) break;
  }
  return ret;
}

ll inv(ll N, ll mod) { // gcd(N, mod) = 1
  // want a so that a * N - mod * b = 1
  pii e = euclid(N, mod, 1);

  if((e.fs * N) % mod != 1) cout<<"CALL THE POLICE "<<endl;

  return e.fs;
}

int main() {
  cin.sync_with_stdio(0);
  cin.tie(0);

  cin>>M>>N>>D;
  Sieve sv(D);

  vpi v = sv.getPrimes(D);
  vll w;
  
  ll fit = 1000000000000000000LL;

  for(auto x: v) {
    ll p = x.fs, e = x.sc;
    ll k = fexp(M, p) - fexp(N, p) - fexp(M-N, p);
    w.pb(k);
    fit = min(fit, k/e);
  }

  for(int i=0;i<sz(v);++i) {
    ll p = v[i].fs, e = v[i].sc;
    ll tot = 1; FOR(i, e) tot *= p;

    ll k = w[i] - fit*e;

    if(k >= e) {
      a.pb(0);
      m.pb(tot);
      cout<<0<<" "<<tot<<endl;
      continue;
    }

    ll mdp = 1; // how much it is mod p^e = p^k * divy(comb, p^(e-k))
    FOR(i, k) mdp *= p;
    ll md = 1;
    FOR(i, e-k) md *= p;

    mdp = (mdp * divy(M, p, md)) % tot;
    mdp = mdp * inv(divy(N, p, md), md) % tot;
    mdp = mdp * inv(divy(M-N, p, md), md) % tot;
    // x / p^e = mdp mod tot

    ll fac = inv(D/tot, tot);
    fac = pw(fac, fit, tot);

    mdp = (mdp * fac ) % tot;

    a.pb(mdp);
    m.pb(tot);
    //cout<<mdp<<" "<<tot<<endl;
  }

  pii sol = chinese(a, m);

  FOR(i, sz(a)) {
    if(sol.fs % m[i] != a[i]) cout<<"CALL THE PRIEST"<<endl;
  }

  cout<<sol.fs<<"\n";
  
  return 0;
}

// A = 2*ab + 2*ac + TAB + TAC
// B = 2*ba +  

詳細信息

Test #1:

score: 100
Accepted
time: 0ms
memory: 3620kb

input:

9 2 3

output:

1

result:

ok single line: '1'

Test #2:

score: 0
Accepted
time: 0ms
memory: 3604kb

input:

5 2 5

output:

2

result:

ok single line: '2'

Test #3:

score: -100
Wrong Answer
time: 0ms
memory: 3620kb

input:

6 3 6

output:

0 2
2

result:

wrong answer 1st lines differ - expected: '2', found: '0 2'