結果

問題 No.187 中華風 (Hard)
コンテスト
ユーザー tonegawa
提出日時 2021-02-07 20:07:10
言語 C++17
(gcc 13.3.0 + boost 1.87.0)
結果
AC  
実行時間 27 ms / 3,000 ms
コード長 7,274 bytes
コンパイル時間 2,262 ms
コンパイル使用メモリ 150,560 KB
最終ジャッジ日時 2025-01-18 13:50:31
ジャッジサーバーID
（参考情報） 		judge2 / judge5
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
other AC * 25
権限があれば一括ダウンロードができます

ソースコード

diff # 

#include <iostream>
#include <string>
#include <vector>
#include <array>
#include <stack>
#include <queue>
#include <deque>
#include <algorithm>
#include <set>
#include <map>
#include <unordered_map>
#include <bitset>
#include <cmath>
#include <functional>
#include <cassert>
#include <iomanip>
#include <numeric>
#include <memory>
#include <random>
#define VV(a, b, c, d) vector<vector<d>>(a, vector<d>(b, c))
#define VVV(a, b, c, d, e) vector<vector<vector<e>>>(a, VV(b, c, d, e));
#define all(obj) (obj).begin(), (obj).end()
typedef long long int ll;
typedef long double ld;
typedef std::vector<ll> v1;
typedef std::vector<v1> v2;
typedef std::vector<v2> v3;
using namespace std;

namespace fact{
  struct MontgomeryReduction{
    __uint128_t MOD, NEG_INV, R2;
    MontgomeryReduction(uint64_t MOD_): MOD(MOD_){
      NEG_INV = 0;
      __uint128_t s = 1, t = 0;
      for(int i=0;i<64;i++){
        if (~t & 1) {
          t += MOD;
          NEG_INV += s;
        }
        t >>= 1;
        s <<= 1;
      }
      R2 = ((__uint128_t)1<<64) % MOD;
      R2 = R2 * R2 % MOD;
    }
    // return x * R % MOD;
    inline uint64_t generate(uint64_t x) const{
      return reduce((__uint128_t)x * R2);
    }
    // return x / R % MOD;
    inline uint64_t reduce(__uint128_t x) const{
      x = (x + ((uint64_t)x * (uint64_t)NEG_INV) * MOD) >> 64;
      return x < MOD? x: x-MOD;
    }
  };
  uint64_t gcd(uint64_t a, uint64_t b){
    if(a == 0) return b;
    if(b == 0) return a;
    int shift = __builtin_ctzll(a | b);
    a >>= __builtin_ctzll(a);
    do {
      b >>= __builtin_ctzll(b);
      if(a > b) std::swap(a, b);
      b -= a;
    }while(b);
    return (a << shift);
  }
  bool is_prime(uint64_t n, const MontgomeryReduction &mr){
    if(n<=1) return false;
    if(n==2) return true;
    if(!(n&1)) return false;
    uint64_t d = n - 1, one = mr.generate(1);
    while(!(d&1)) d >>= 1;
    for(uint64_t a : {2, 325, 9375, 28178, 450775, 9780504, 1795265022}){
      if(n <= a) break;
      uint64_t t = d;
      auto modpow = [&](uint64_t x, uint64_t k){
        uint64_t r = one;
        while(k){
          if(k&1) r = mr.reduce((__uint128_t)r * x);
          x = mr.reduce((__uint128_t)x * x);
          k >>= 1;
        }
        return r;
      };
      uint64_t y = modpow(mr.generate(a), t), n_minus_one = mr.generate(n-1);
      while(t != n-1 && y != one && y != n_minus_one){
        y = mr.reduce((__uint128_t)y * y);
        t <<= 1;
      }
      if(y != n_minus_one && t % 2 == 0) return false;
    }
    return true;
  }

  uint64_t rho(uint64_t n){
    if((n&1)==0) return 2;
    for(auto p:{3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41}){
      if(n%p==0) return p;
      if(n<p) break;
    }
    MontgomeryReduction mr(n);
    if(is_prime(n, mr)) return n;
    uint64_t m = 1LL << ((64 - __builtin_clzll(n)) / 8);
    uint32_t c = 0;
    auto f = [&](uint64_t num){
      uint64_t tmp = mr.generate(num);
      tmp = mr.reduce(mr.reduce((__uint128_t)tmp * tmp)) + c;
      return tmp >= n ? tmp - n : tmp;};

    while(true){
      c = (3 * c + 1)%99;
      uint64_t y = 2, r = 1, q = 1, g = 1, ys, x;
      while(g==1){
        x = y;
        for(int i=0;i<r;i++) y = f(y);
        uint64_t k = 0;
        while(k < r && g == 1){
          q = mr.generate(q);
          ys = y;
          for(int j=0;j<min(m, r-k);j++){
            y = f(y);
            q = mr.reduce((__uint128_t)q * mr.generate(x>y?x-y:y-x));
          }
          g = gcd(mr.reduce(q), n);
          k += m;
        }
        r <<= 1;
      }
      if(g==n){
        while(g==1){
          ys = f(ys);
          g = gcd(n, (x>y?x-y:y-x));
        }
      }
      if(g<n) return g;
    }
  }
  vector<ll> factorize(ll n) {
    if(n == 1) return {};
    ll x = rho(n);
    if(x == n) return {x};
    vector<ll> le = factorize(x);
    vector<ll> ri = factorize(n / x);
    le.insert(le.end(), ri.begin(), ri.end());
    return le;
  }
  vector<ll> PrimeFactor(ll N){
    vector<ll> p = factorize(N);
    sort(all(p));
    decltype(p)::iterator result = std::unique(p.begin(), p.end());
    p.erase(result, p.end());
    return p;
  }
  vector<ll> DivisorRestore(ll N){
    vector<ll> p = factorize(N);
    sort(all(p));
    vector<vector<ll>> x(1, {1});
    ll num = 1, idx = 0;
    for(int i=0;i<p.size();i++){
      num *= p[i];
      x[idx].push_back(num);
      if(i!=p.size()-1&&p[i+1]!=p[i]){
        x.push_back(vector<ll>{1});
        num = 1;
        idx++;
      }
    }
    ll l = 0, r = 1;
    vector<ll> ret{1};
    for(int i=0;i<x.size();i++){
      for(auto e:x[i]){
        for(int j=l;j<r;j++){
          ret.push_back(ret[j] * e);
        }
      }
      l = r;
      r = ret.size();
    }
    return vector<ll>(ret.begin()+l, ret.end());
  }
}

template<typename T>
T mod_inv(T a, T MOD){
  T u[] = {a, 1, 0}, v[] = {MOD, 0, 1}, t;
    while(*v){
      t = *u / *v;
      swap(u[0] -= t * v[0], v[0]);
      swap(u[1] -= t * v[1], v[1]);
      swap(u[2] -= t * v[2], v[2]);
    }
    u[1] %= MOD;
    return (u[1] < 0) ? (u[1] + MOD) : u[1];
}
template<int MOD>
int mod_inv(int a){
  int u[] = {a, 1, 0}, v[] = {MOD, 0, 1}, t;
  while(*v){
    t = *u / *v;
    swap(u[0] -= t * v[0], v[0]);
    swap(u[1] -= t * v[1], v[1]);
    swap(u[2] -= t * v[2], v[2]);
  }
  u[1] %= MOD;
  return (u[1] < 0) ? (u[1] + MOD) : u[1];
}

//A_i = x_i mod m_i
//-> lcm(m)以下のmodで全ての条件を満たす数を返す
//xが全て0の場合は0
//unsafe: 互いに素かつ解が必ず存在する O(N^2)
ll garner_unsafe(const vector<ll>&x, vector<ll> m, ll MOD){
  int n = x.size();
  assert(n == m.size());
  vector<ll> kp(n+1, 0), rmult(n+1, 1);
  m.push_back(MOD);
  for(int i=0;i<n;i++){
    ll y = ((x[i] - kp[i]) * mod_inv<ll>(rmult[i], m[i])) % m[i];
    if(y < 0) y += m[i];
    for(int j=i+1;j<=n;j++){
      kp[j] = (kp[j] + rmult[j] * y) % m[j];
      rmult[j] = (rmult[j] * m[i]) % m[j];
    }
  }
  return kp[n];
}

//A_i = x_i mod m_i
//-> lcm(m)以下のmodで全ての条件を満たす数を返す
//xが全て0の場合は0
//互いに素でない場合、 解が存在しない場合も使える O((Nlog(M_MAX))^2 + N * M_MAX^(1/4))
ll garner(const vector<ll> x, vector<ll> m, ll MOD){
  int n = x.size();
  assert(n == m.size());
  unordered_map<ll, vector<pair<ll, ll>>> mp;
  bool f = false;
  for(int i=0;i<n;i++){
    if(x[i]) f = true;
    auto prime = fact::PrimeFactor(m[i]);
    for(auto p:prime){
      ll num = 1, mi = m[i];
      while(mi % p == 0){
        mi /= p;
        num *= p;
      }
      mp[p].emplace_back(x[i] % num, num);
    }
  }
  if(!f){
    ll ans = 1;
    for(auto &k:mp){
      ll max_p = -1, mo = -1;
      for(auto pa:k.second) if(pa.second > max_p) max_p = pa.second, mo = pa.first;
      ans = (ans * max_p) % MOD;
    }
    return ans;
  }
  vector<ll> x2, m2;
  for(auto &k:mp){
    ll max_p = -1, mo = -1;
    for(auto pa:k.second) if(pa.second > max_p) max_p = pa.second, mo = pa.first;
    for(auto pa:k.second) if(mo % pa.second != pa.first) return -1;
    x2.push_back(mo);
    m2.push_back(max_p);
  }
  return garner_unsafe(x2, m2, MOD);
}

int main(){
  int n;std::cin >> n;
  vector<ll> x(n), m(n);
  for(int i=0;i<n;i++) std::cin >> x[i] >> m[i];
  std::cout << garner(x, m, 1000000007) << '\n';
}
0