#130095 (C++14) No.444 旨味の相乗効果

提出ソース

結果

問題	No.444 旨味の相乗効果
ユーザー	Min_25
提出日時	2016-11-12 02:04:14
言語	C++14 (gcc 12.3.0 + boost 1.83.0)
結果	AC
実行時間	4 ms / 2,500 ms
コード長	11,620 bytes
コンパイル時間	2,338 ms
コンパイル使用メモリ	126,436 KB
実行使用メモリ	5,376 KB
最終ジャッジ日時	2024-06-24 20:54:07
合計ジャッジ時間	3,240 ms
ジャッジサーバーID （参考情報）	judge5 / judge1

このコードへのチャレンジ
（要ログイン）

テストケース

テストケース表示

入力	結果	実行時間実行使用メモリ
testcase_00	AC	2 ms 5,248 KB
testcase_01	AC	2 ms 5,376 KB
testcase_02	AC	2 ms 5,376 KB
testcase_03	AC	2 ms 5,376 KB
testcase_04	AC	4 ms 5,376 KB
testcase_05	AC	2 ms 5,376 KB
testcase_06	AC	2 ms 5,376 KB
testcase_07	AC	2 ms 5,376 KB
testcase_08	AC	2 ms 5,376 KB
testcase_09	AC	2 ms 5,376 KB
testcase_10	AC	2 ms 5,376 KB
testcase_11	AC	2 ms 5,376 KB
testcase_12	AC	3 ms 5,376 KB
testcase_13	AC	3 ms 5,376 KB
testcase_14	AC	2 ms 5,376 KB
testcase_15	AC	2 ms 5,376 KB
testcase_16	AC	4 ms 5,376 KB
testcase_17	AC	2 ms 5,376 KB
testcase_18	AC	2 ms 5,376 KB
testcase_19	AC	3 ms 5,376 KB
testcase_20	AC	3 ms 5,376 KB
testcase_21	AC	2 ms 5,376 KB
testcase_22	AC	2 ms 5,376 KB
testcase_23	AC	2 ms 5,376 KB
testcase_24	AC	2 ms 5,376 KB
testcase_25	AC	2 ms 5,376 KB
testcase_26	AC	2 ms 5,376 KB
testcase_27	AC	2 ms 5,376 KB

権限があれば一括ダウンロードができます

ソースコード

raw source code

#pragma GCC optimize ("O3")
#pragma GCC target ("avx")

#include <cstdio>
#include <cassert>
#include <cmath>
#include <cstring>

#include <iostream>
#include <algorithm>
#include <vector>
#include <map>
#include <set>
#include <functional>
#include <tuple>
#include <stack>

#include <initializer_list>

#define _rep(_1, _2, _3, _4, name, ...) name
#define rep2(i, n) rep3(i, 0, n)
#define rep3(i, a, b) rep4(i, a, b, 1)
#define rep4(i, a, b, c) for (int i = int(a); i < int(b); i += int(c))
#define rep(...) _rep(__VA_ARGS__, rep4, rep3, rep2, _)(__VA_ARGS__)

using namespace std;

using i64 = long long;
using u8 = unsigned char;
using u32 = unsigned;
using u64 = unsigned long long;
using f80 = long double;

// for dense
class matrix {
public:
  using R = int;
  using vec = vector<R>;
  using poly = vec;

  struct xor128 {
    xor128(u32 seed=521288629) : z(seed) {}
    u32 operator() () { return next(); }
    u32 operator() (u32 r) { return next() % r; }
    u32 next() {
      u32 t = x ^ (x << 11);
      x = y; y = z; z = w;
      return (w ^= (w >> 19) ^ (t ^ (t >> 8)));
    }
    u32 x = 123456789, y = 362436069, w = 88675123, z = 521288629;
  };

  static void reset_seed(int S) {
    xrand = xor128(S);
  }

private:
  struct matrix64 {
    matrix64(const matrix& rhs) : m(rhs.m_), n(rhs.n_), a(rhs.a, rhs.a + m * n) {}
    i64* operator [] (int i) { return a.data() + n * i; }
    void swap_row(int i, int j) {
      if (i == j) return;
      rep(k, n) swap((*this)[i][k], (*this)[j][k]);
    }
    void print() {
      printf("matrix64(%d, %d):\n", m, n);
      puts("[");
      rep(i, m) {
        printf("  [%d", fixed((*this)[i][0]));
        rep(j, 1, n) printf(", %d", fixed((*this)[i][j]));
        puts("],");
      }
      puts("]");
    }
    int m, n;
    vector<i64> a;
  };

  static R gcd(R a, R b) {
    return b ? gcd(b, a % b) : a;
  }
  static R mul_mod(R a, R b) { return i64(a) * b % mod; }
  static i64 add_mod64(i64 a, i64 b) { return (a += b - lmod) < 0 ? a + lmod : a; }
  static i64 sub_mod64(i64 a, i64 b) { return (a -= b) < 0 ? a + lmod : a; }
  static R fixed(R a) { return (a %= mod) < 0 ? a + mod : a; }
  static R mod_inv(R a, R m=0) {
    R b = (m == 0) ? mod : m;
    R s = 1, t = 0;
    while (b) {
      R q = a / b;
      swap(a %= b, b);
      swap(s -= t * q, t);
    }
    if (a != 1) {
      fprintf(stderr, "Error: Modular multiplicative inverse does not exist.\n");
      assert(0);
    }
    return s < 0 ? s + mod : s;
  }

public:
  static void set_mod(R m) {
    mod = m;
    lmod = i64(mod) << 32;
    modq = min(u64(1) << 16, (u64(-1) / mod - 1) / mod);
  }

public:
  matrix(int n) : matrix(n, n) {}
  matrix(int m, int n, bool zfill=true) : m_(m), n_(n) {
    alloc();
    if (zfill) clear();
  }
  matrix(const matrix& m) : matrix(m.m_, m.n_, false) {
    copy(m.a, m.a + m_ * n_, a);
  }
  matrix(initializer_list< initializer_list<R> > l) : 
      matrix(l.size(), l.begin()->size()) {
    int i = 0;
    for (auto& v : l) {
      int j = 0;
      for (auto c : v) (*this)[i][j++] = fixed(c);
      ++i;
    }
  }
  matrix& operator = (matrix&& rhs) {
    if (this == &rhs) return *this;
    this->m_ = rhs.m_;
    this->n_ = rhs.n_;
    free();
    this->a = rhs.a; rhs.a = nullptr;
    return *this;
  }
  ~matrix() { free(); }

  R* operator [] (int i) { return a + i * n_; }
  const R* operator [] (int i) const { return a + i * n_; }

  void clear() { fill(a, a + m_ * n_, 0); }

  void swap_row(int i, int j) {
    if (i == j) return;
    rep(k, n_) swap((*this)[i][k], (*this)[j][k]);
  }

  static matrix random_matrix(int m, int n) {
    auto ret = matrix(m, n, false);
    rep(i, m) rep(j, n) ret[i][j] = xrand(mod);
    return ret;
  }

  static vec random_vector(int n) {
    auto ret = vec(n);
    rep(i, n) ret[i] = xrand(mod);
    return ret;
  }

  matrix transpose() const {
    auto ret = matrix(n_, m_, false);
    rep(i, m_) rep(j, n_) ret[j][i] = (*this)[i][j];
    return ret;
  }

  matrix operator * (const matrix& rhs) const {
    assert(n_ == rhs.m_);
    auto ret = matrix(m_, rhs.n_);
    auto rhsT = rhs.transpose();
    rep(i, m_) mul_mat_vec(rhsT, (*this)[i], ret[i]);
    return ret;
  }
  
  matrix& operator *= (const matrix& rhs) {
    return *this = *this * rhs;
  }

  matrix krylov(int k, vec v) const {
    auto ret = matrix(k + 1, n_);
    rep(i, k) {
      copy(v.begin(), v.end(), ret[i]);
      mul_mat_vec(*this, v.data(), v.data());
    }
    copy(v.begin(), v.end(), ret[k]);
    return ret;
  }

  template <typename T>
  static void vec_sub(i64* a, T* b, int n, R c) {
    rep(i, n) a[i] = sub_mod64(a[i], i64(R(b[i])) * c);
  }

  static vec solve_linear_equations_in_Z_pZ(const matrix& in) {
    matrix64 mat = matrix64(in);
    assert(mat.m + 1 == mat.n);

    int rank = 0;
    rep(i, mat.m) {
      int piv = -1;
      for (int j = mat.m - 1; j >= rank; --j) if (mat[j][i] %= mod) piv = j;
      if (piv < 0) continue;
      if (rank != piv) mat.swap_row(rank, piv);
      rep(k, i + 1, mat.n) mat[rank][k] %= mod;
      auto inv = mod_inv(mat[rank][i]);
      rep(j, rank + 1, mat.m) if (mat[j][i]) {
        vec_sub(mat[j] + i, mat[rank] + i, mat.n - i, mul_mod(mat[j][i], inv));
      }
      ++rank;
    }
    for (int i = rank - 1; i >= 0; --i) if (mat[i][i]) {
      R c = mat[i][mat.m] = mul_mod(mat[i][mat.m] % mod, mod_inv(mat[i][i]));
      rep(j, i) mat[j][mat.m] = sub_mod64(mat[j][mat.m], i64(c) * int(mat[j][i]));
    }
    vec ret(mat.n);
    ret[0] = fixed(1);
    rep(i, rank) ret[ret.size() - 1 - i] = mat[i][mat.m];
    return ret;
  }

  static vec solve_linear_equations_in_Z_nZ(const matrix& in) {
    matrix64 mat = matrix64(in);
    assert(mat.m + 1 == mat.n);

    int rank = 0;
    rep(i, mat.m) {
      int piv = -1;
      R g = mod;
      rep(j, rank, mat.m) if ( (mat[j][i] %= mod) && g >= 0) {
        R k = gcd(R(mat[j][i]), mod);
        g = (k == 1) ? -1 : gcd(g, k);
        piv = j;
      }
      if (g == mod) continue;
      if (g > 0) {
        // To Do: not optimal
        R g2 = gcd(mat[piv][i], mod);
        for (int j = rank; j < mat.n && g2 != g; ++j) if (mat[j][i]) {
          R g3 = gcd(mat[j][i], g2);
          if (g3 == g2) continue;
          g2 = g3;
          while (gcd(mat[piv][i], mod) != g2) {
            R q = mat[j][i] / mat[piv][i];
            rep(k, i + 1, mat.n) mat[piv][k] %= mod;
            if (q) vec_sub(mat[j] + i, mat[piv] + i, mat.n - i, q);
            mat.swap_row(piv, j);
          }
        }
        assert(g2 == g);
      } else {
        g = 1;
      }
      if (rank != piv) mat.swap_row(rank, piv);
      rep(k, i + 1, mat.n) mat[rank][k] %= mod;
      auto inv = mod_inv(mat[rank][i] / g, mod / g);
      rep(j, rank + 1, mat.m) if (mat[j][i]) {
        vec_sub(mat[j] + i, mat[rank] + i, mat.n - i, mul_mod(mat[j][i] / g, inv));
      }
      ++rank;
    }
    for (int i = rank - 1; i >= 0; --i) if (mat[i][i]) {
      R g = gcd(gcd(R(mat[i][i]), R(mat[i][mat.m] %= mod)), mod);
      R inv = mod_inv(mat[i][i] / g, mod / g);
      R c = mat[i][mat.m] = mul_mod(mat[i][mat.m] / g, inv);
      rep(j, i) mat[j][mat.m] = sub_mod64(mat[j][mat.m], i64(c) * int(mat[j][i]));
    }
    vec ret(mat.n);
    ret[0] = fixed(1);
    rep(i, rank) ret[ret.size() - 1 - i] = mat[i][mat.m];
    return ret;
  }

  static vec linear_recurrence(matrix kry, const vec& v, bool in_Z_pZ) {
    int m = kry.m_, n = kry.n_;
    assert(m == n + 1 && n == int(v.size()));
    rep(i, n) kry[n][i] = (kry[n][i] ? mod - kry[n][i] : 0);

    kry = kry.transpose();
    auto ret = in_Z_pZ ? solve_linear_equations_in_Z_pZ(kry)
                       : solve_linear_equations_in_Z_nZ(kry);
    return ret;
  }

  vec pow_vec(i64 e, const vec& v, bool in_Z_pZ=false) const {
    assert(m_ == n_);
    if (mod == 1) return vec(v.size(), 0);
    int k = min(e, i64(m_));
    auto kry = krylov(k, v);
    if (e <= k) return vec(kry[e], kry[e + 1]);
    
    auto poly = linear_recurrence(kry, v, in_Z_pZ);
    auto rem = x_pow_mod(e, poly);

    reverse(rem.begin(), rem.end());
    tmp64.assign(v.size(), 0);
    rep(i, rem.size()) if (rem[i]) rep(j, v.size()) {
      tmp64[j] = add_mod64(tmp64[j], i64(rem[i]) * kry[i][j]);
    }
    auto ret = vec(v.size());
    rep(i, v.size()) ret[i] = tmp64[i] % mod;
    return ret;
  }

  static poly poly_mul(const poly& f, const poly& g) {
    int s1 = f.size(), s2 = g.size(), s = s1 + s2 - 1;
    tmp64.assign(s, 0);
    rep(i, s1) if (f[i]) rep(j, s2) {
      tmp64[i + j] = add_mod64(tmp64[i + j], i64(f[i]) * int(g[j]));
    }
    poly ret(s);
    rep(i, s) ret[i] = tmp64[i] % mod;
    return ret;
  }

  static poly poly_rem(const poly& f, const poly& g) {
    int s1 = f.size(), s2 = g.size();
    if (s1 < s2) return f;
    assert(g[0] == 1);
    tmp64.resize(s1);
    copy(f.begin(), f.end(), tmp64.begin());
    rep(i, s1 - s2 + 1) {
      int c = tmp64[i] % mod;
      if (c) rep(j, 1, s2) tmp64[i + j] = sub_mod64(tmp64[i + j], i64(c) * int(g[j]));
      tmp64[i] = c;
    }
    poly ret(s2 - 1);
    rep(i, s2 - 1) ret[i] = tmp64[s1 - s2 + 1 + i] % mod;
    return ret;
  }

  static poly x_pow_mod(i64 e, const poly& f) {
    if (e == 0) return poly(1, 1);
    poly ret = poly({1, 0});
    i64 mask = (i64(1) << __lg(e)) >> 1;
    while (mask) {
      ret = poly_mul(ret, ret);
      if (e & mask) ret.push_back(0);
      ret = poly_rem(ret, f);
      mask >>= 1;
    }
    return ret;
  }

  static void mul_mat_vec(const matrix& m, const R* v, R* res) {
    buff.resize(m.m_);
    rep(i, m.m_) buff[i] = dot(m[i], v, m.n_);
    copy(buff.begin(), buff.begin() + m.m_, res);
  }

  void print() const {
    printf("matrix(%d, %d):\n", m_, n_);
    puts("[");
    rep(i, m_) {
      printf("  [%d", (*this)[i][0]);
      rep(j, 1, n_) printf(", %d", (*this)[i][j]);
      puts("],");
    }
    puts("]");
  }

private:
  void alloc() {
    a = new R[m_ * n_];
  }
  void free() {
    delete [] a; 
  }

  static int dot(const int* a, const int* b, int n) {
    if (n == 0) return 0;
    int k = min(n, modq) & ~3;
    int r = (k ? n % k : n), q = (k ? n / k : 0);
    i64 s0 = 0, s1 = 0;
    if (k >= 12) {
      rep(i, q) {
        u64 s = 0;
        rep(j, i * k, (i + 1) * k) s += i64(a[j]) * b[j];
        s0 += u32(s), s1 += u32(s >> 32);
      }
      if (r) {
        u64 s = 0;
        rep(j, q * k, n) s += i64(a[j]) * b[j];
        s0 += u32(s), s1 += u32(s >> 32);
      }
    } else {
      rep(j, n) {
        i64 s = i64(a[j]) * b[j];
        s0 += u32(s), s1 += u32(s >> 32);
      }
    }
    return (((s1 % mod) << 32) + s0) % mod;
  }

  static int mod, modq;
  static i64 lmod;
  static vector<int> buff;
  static vector<i64> tmp64;
  static xor128 xrand;

  int m_, n_;
  R* a;
};
vector<int> matrix::buff;
vector<i64> matrix::tmp64;
int matrix::mod, matrix::modq;
i64 matrix::lmod;
matrix::xor128 matrix::xrand;

void solve() {
  const int mod = 1e9 + 7;
  matrix::set_mod(mod);
  auto pow_mod = [&] (int b, i64 e) {
    int ret = 1;
    for (; e; e >>= 1, b = i64(b) * b % mod) if (e & 1) ret = i64(ret) * b % mod;
    return ret;
  };
  int A[128];
  int n; i64 c;
  while (~scanf("%d %lld", &n, &c)) {
    matrix mat(n, n, true);
    int ans = 0;
    rep(i, n) {
      scanf("%d", &A[i]);
      rep(j, i + 1) mat[i][j] = A[i];
      ans = (ans + mod - pow_mod(A[i], c)) % mod;
    }
    auto res = mat.pow_vec(c - 1, matrix::vec(A, A + n));
    rep(i, n) ans = (ans + res[i]) % mod;
    printf("%d\n", ans);
  }
}

int main() {
  auto beg = clock();
  solve();
  auto end = clock();
  fprintf(stderr, "%.3f sec\n", double(end - beg) / CLOCKS_PER_SEC);
}

yukicoder

結果

テストケース

ソースコード