#218155 (C++14) No.303 割れません

提出ソース
結果

問題	No.303 割れません
ユーザー	Min_25
提出日時	2017-11-24 18:31:23
言語	C++14 (gcc 13.2.0 + boost 1.83.0)
結果	AC
実行時間	35 ms / 10,000 ms
コード長	11,956 bytes
コンパイル時間	1,668 ms
コンパイル使用メモリ	106,208 KB
実行使用メモリ	6,728 KB
最終ジャッジ日時	2023-08-18 04:26:29
合計ジャッジ時間	3,775 ms
ジャッジサーバーID （参考情報）	judge11 / judge14
このコードへのチャレンジ
（要ログイン）
テストケース

テストケース表示
入力	結果	実行時間実行使用メモリ
testcase_00	AC	2 ms 4,376 KB
testcase_01	AC	1 ms 4,376 KB
testcase_02	AC	1 ms 4,380 KB
testcase_03	AC	15 ms 4,420 KB
testcase_04	AC	28 ms 5,660 KB
testcase_05	AC	28 ms 5,288 KB
testcase_06	AC	28 ms 5,404 KB
testcase_07	AC	35 ms 6,548 KB
testcase_08	AC	35 ms 6,724 KB
testcase_09	AC	34 ms 6,728 KB
testcase_10	AC	2 ms 4,380 KB
testcase_11	AC	29 ms 5,316 KB
testcase_12	AC	35 ms 6,712 KB
testcase_13	AC	16 ms 4,376 KB
testcase_14	AC	17 ms 4,700 KB
testcase_15	AC	9 ms 4,376 KB
testcase_16	AC	1 ms 4,376 KB
権限があれば一括ダウンロードができます
ソースコード

raw source code
#include <cstdio>
#include <cassert>
#include <cmath>
#include <cstring>

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

#include <tuple>

#define getchar getchar_unlocked
#define putchar putchar_unlocked

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

int get_int() {
  int n, c, sign = 0;
  while ((c = getchar()) < '-');
  if (c == '-') sign = 1, n = 0;
  else n = c - '0';
  while ((c = getchar()) >= '0') n = n * 10 + c - '0';
  return sign ? -n : n;
}

namespace ntt {

using word_t = u64;
using dword_t = __uint128_t;

static constexpr word_t mul_inv(word_t n, int e=6, word_t x=1) {
  return e == 0 ? x : mul_inv(n, e-1, x*(2-x*n));
}

template <word_t mod, word_t prim_root>
class UnsafeMod {
private:
  static const int word_bits = 8 * sizeof(word_t);

public:
  static constexpr word_t inv = mul_inv(mod);
  static constexpr word_t r2 = -dword_t(mod) % mod;
  static constexpr int level = __builtin_ctzll(mod - 1);
  static_assert(inv * mod == 1, "invalid 1/M modulo 2^@.");

  UnsafeMod() {}
  UnsafeMod(word_t n) : x(init(n)) {};
  static word_t modulus() { return mod; }
  static word_t init(word_t w) { return reduce(dword_t(w) * r2); }
  static word_t reduce(const dword_t w) { 
    return word_t(w >> word_bits) 
         + mod - word_t((dword_t(word_t(w) * inv) * mod) >> word_bits); }
  static UnsafeMod omega() { return UnsafeMod(prim_root).pow((mod - 1) >> level); }
  UnsafeMod& operator += (UnsafeMod rhs) { x += rhs.x; return *this; }
  UnsafeMod& operator -= (UnsafeMod rhs) { x += 3 * mod - rhs.x; return *this; }
  UnsafeMod& operator *= (UnsafeMod rhs) { x = reduce(dword_t(x) * rhs.x); return *this; }
  UnsafeMod operator + (UnsafeMod rhs) const { return UnsafeMod(*this) += rhs; }
  UnsafeMod operator - (UnsafeMod rhs) const { return UnsafeMod(*this) -= rhs; }
  UnsafeMod operator * (UnsafeMod rhs) const { return UnsafeMod(*this) *= rhs; }
  word_t get() const { return reduce(x) % mod; }
  void set(word_t n) { x = n; }
  UnsafeMod pow(word_t e) const {
    UnsafeMod ret = UnsafeMod(1);
    for (UnsafeMod base = *this; e; e >>= 1, base *= base) if (e & 1) ret *= base;
    return ret;
  }
  UnsafeMod inverse() const { return pow(mod - 2); }
  friend ostream& operator << (ostream& os, const UnsafeMod& m) { return os << m.get(); }
  static void debug() { printf("%llu %llu %llu %llu\n", mod, inv, r2, omega().get()); }
  word_t x;
};

template <typename mod_t>
void transform(mod_t* A, int n, const mod_t* roots, const mod_t* iroots) {
  const int logn = __builtin_ctz(n), nh = n >> 1, lv = mod_t::level;
  const mod_t one = mod_t(1), imag = roots[lv - 2];

  mod_t dw[lv - 1]; dw[0] = roots[lv - 3];
  for (int i = 1; i < lv - 2; ++i) dw[i] = dw[i - 1] * iroots[lv - 1 - i] * roots[lv - 3 - i];
  dw[lv - 2] = dw[lv - 3] * iroots[1];

  if (logn & 1) for (int i = 0; i < nh; ++i) {
    mod_t a = A[i], b = A[i + nh];
    A[i] = a + b; A[i + nh] = a - b;
  }
  for (int e = logn & ~1; e >= 2; e -= 2) {
    const int m = 1 << e, m4 = m >> 2;
    mod_t w2 = one;
    for (int i = 0; i < n; i += m) {
      const mod_t w1 = w2 * w2, w3 = w1 * w2;
      for (int j = i; j < i + m4; ++j) {
        mod_t a0 = A[j + m4 * 0] * one, a1 = A[j + m4 * 1] * w2;
        mod_t a2 = A[j + m4 * 2] * w1,  a3 = A[j + m4 * 3] * w3;
        mod_t t02p = a0 + a2, t13p = a1 + a3;
        mod_t t02m = a0 - a2, t13m = (a1 - a3) * imag;
        A[j + m4 * 0] = t02p + t13p; A[j + m4 * 1] = t02p - t13p;
        A[j + m4 * 2] = t02m + t13m; A[j + m4 * 3] = t02m - t13m;
      }
      w2 *= dw[__builtin_ctz(~(i >> e))];
    }
  }
}

template <typename mod_t>
void itransform(mod_t* A, int n, const mod_t* roots, const mod_t* iroots) {
  const int logn = __builtin_ctz(n), nh = n >> 1, lv = mod_t::level;
  const mod_t one = mod_t(1), imag = iroots[lv - 2];

  mod_t dw[lv - 1]; dw[0] = iroots[lv - 3];
  for (int i = 1; i < lv - 2; ++i) dw[i] = dw[i - 1] * roots[lv - 1 - i] * iroots[lv - 3 - i];
  dw[lv - 2] = dw[lv - 3] * roots[1];

  for (int e = 2; e <= logn; e += 2) {
    const int m = 1 << e, m4 = m >> 2;
    mod_t w2 = one;
    for (int i = 0; i < n; i += m) {
      const mod_t w1 = w2 * w2, w3 = w1 * w2;
      for (int j = i; j < i + m4; ++j) {
        mod_t a0 = A[j + m4 * 0], a1 = A[j + m4 * 1];
        mod_t a2 = A[j + m4 * 2], a3 = A[j + m4 * 3];
        mod_t t01p = a0 + a1, t23p = a2 + a3;
        mod_t t01m = a0 - a1, t23m = (a2 - a3) * imag;
        A[j + m4 * 0] = (t01p + t23p) * one; A[j + m4 * 2] = (t01p - t23p) * w1;
        A[j + m4 * 1] = (t01m + t23m) * w2;  A[j + m4 * 3] = (t01m - t23m) * w3;
      }
      w2 *= dw[__builtin_ctz(~(i >> e))];
    }
  }
  if (logn & 1) for (int i = 0; i < nh; ++i) {
    mod_t a = A[i], b = A[i + nh];
    A[i] = a + b; A[i + nh] = a - b;
  }
}

template <typename mod_t>
void convolve(mod_t* A, int s1, mod_t* B, int s2, bool cyclic=false) {
  const int s = cyclic ? max(s1, s2) : s1 + s2 - 1;
  const int size = 1 << (31 - __builtin_clz(2 * s - 1));
  assert(size <= (i64(1) << mod_t::level));
  
  mod_t roots[mod_t::level], iroots[mod_t::level];
  roots[0] = mod_t::omega();
  for (int i = 1; i < mod_t::level; ++i) roots[i] = roots[i - 1] * roots[i - 1];
  iroots[0] = roots[0].inverse();
  for (int i = 1; i < mod_t::level; ++i) iroots[i] = iroots[i - 1] * iroots[i - 1];

  fill(A + s1, A + size, 0); transform(A, size, roots, iroots);
  const mod_t inv = mod_t(size).inverse();
  if (A == B && s1 == s2) {
    for (int i = 0; i < size; ++i) A[i] *= A[i] * inv;
  } else {
    fill(B + s2, B + size, 0); transform(B, size, roots, iroots);
    for (int i = 0; i < size; ++i) A[i] *= B[i] * inv;
  }
  itransform(A, size, roots, iroots);
}

} // namespace ntt

using m64_1 = ntt::UnsafeMod<1121333583512862721, 51>;
using m64_2 = ntt::UnsafeMod<1128298388379402241, 23>; // <= 1.14e18 (sub.D = 3)

template <u64 Modulus>
class FastDiv21 {
  using u128 = __uint128_t;
  static constexpr int s = __builtin_clzll(Modulus);
  static constexpr u64 m = Modulus << s;
  static constexpr u64 v = u64(~(u128(m) << 64) / m);
public:
  pair<u64, u64> divmod(u128 a) const {
    a <<= s;
    u64 ah = a >> 64, al = a;
    u128 q = u128(ah) * v + a;
    u64 qh = u64(q >> 64) + 1, ql = q, r = al - qh * m;
    if (r > ql) --qh, r += m;
    if (r >= m) ++qh, r -= m;
    return {qh, r >> s};
  }
  friend u64 operator % (u128 a, const FastDiv21& b) { return b.divmod(a).second; }
  friend u64 operator / (u128 a, const FastDiv21& b) { return b.divmod(a).first; }
};

using word_t = u64;
using sword_t = i64;
using dword_t = __uint128_t;
constexpr int kWordBits = sizeof(word_t) * 8;

template <word_t Base>
class BigInteger : public vector<word_t> {
  static_assert(word_t(Base) < (word_t(1) << (kWordBits - 1)), "Base is too large.");

public:
  BigInteger() : BigInteger(0) {}
  BigInteger(word_t n) : vector<word_t>(1, n) {}
  BigInteger(size_t s, word_t n) : vector<word_t>(s, n) {}

  void normalize() {
    while (size() > 1 && back() == 0) pop_back();
  }

  BigInteger operator + (const BigInteger& rhs) const {
    // inefficient
    size_t s = max(size(), rhs.size()) + 1;
    BigInteger ret(s, 0); copy(begin(), end(), ret.begin());
    word_t carry = 0;
    for (size_t i = 0; i < rhs.size(); ++i) {
      word_t a = ret[i] + rhs[i] + carry;
      carry = 0;
      if (a >= Base) ++carry, a -= Base;
      ret[i] = a;
    }
    for (size_t i = rhs.size(); carry; ++i) {
      word_t a = ret[i] + carry;
      carry = 0;
      if (a >= Base) ++carry, a -= Base;
      ret[i] = a;
    }
    ret.normalize();
    return ret;
  }

  BigInteger operator - (const BigInteger& rhs) const {
    // inefficient
    assert(size() > rhs.size() || (size() == rhs.size() && back() >= rhs.back()));
    BigInteger ret(size(), 0); copy(begin(), end(), ret.begin());
    word_t carry = 0;
    for (size_t i = 0; i < rhs.size(); ++i) {
      word_t a = ret[i] - rhs[i] - carry;
      carry = 0;
      if (sword_t(a) < 0) a += Base, carry = 1;
      ret[i] = a;
    }
    for (size_t i = rhs.size(); carry; ++i) {
      word_t a = ret[i] - carry;
      carry = 0;
      if (sword_t(a) < 0) a += Base, carry = 1;
      ret[i] = a;
    }
    ret.normalize();
    return ret;
  }

  BigInteger operator * (const BigInteger& rhs) const {
    int s1 = size(), s2 = rhs.size(), s = s1 + s2 - 1;
    int ntt_size = 1 << (31 - __builtin_clz(2 * s - 1));
    vector<m64_1> f1(ntt_size); copy(begin(), end(), f1.begin());
    if (this != &rhs) {
      vector<m64_1> g1(ntt_size); copy(rhs.begin(), rhs.end(), g1.begin());
      ntt::convolve(f1.data(), s1, g1.data(), s2);
    } else {
      ntt::convolve(f1.data(), s1, f1.data(), s1);
    }
    vector<m64_2> f2(ntt_size); copy(begin(), end(), f2.begin());
    if (this != &rhs) {
      vector<m64_2> g2(ntt_size); copy(rhs.begin(), rhs.end(), g2.begin());
      ntt::convolve(f2.data(), s1, g2.data(), s2);
    } else {
      ntt::convolve(f2.data(), s1, f2.data(), s1);
    }
    BigInteger ret(s1 + s2, 0);
    auto fdiv = FastDiv21<Base>();
    const auto p1 = m64_1::modulus(), p2 = m64_2::modulus();
    const auto inv = m64_2(p1).inverse();
    dword_t carry = 0;
    for (int i = 0; i < s1 + s2 - 1; ++i) {
      auto r1 = f1[i].get(), r2 = f2[i].get(); // not optimal
      auto prod = r1 + dword_t((m64_2(r2 + p2 - r1) * inv).get()) * p1;
      prod += carry;
      word_t ph = prod >> kWordBits, pl = prod;
      word_t qh = ph / Base, r = ph % Base;
      word_t ql; tie(ql, r) = fdiv.divmod(dword_t(r) << kWordBits | pl);
      carry = dword_t(qh) << kWordBits | ql;
      ret[i] = r;
    }
    ret[s1 + s2 - 1] = carry;
    ret.normalize();
    return ret;
  }

  BigInteger pow(word_t e) const {
    if (e == 0) return BigInteger(1);
    BigInteger ret = *this;
    for (int mask = (1 << __lg(e)) >> 1; mask; mask >>= 1) {
      ret = ret * ret;
      if (mask & e) ret = ret * (*this);
    }
    return ret;
  }
};

constexpr u64 ten_pow(int e, u64 x=1) {
  return e <= 0 ? x : ten_pow(e - 1, x * 10);
}

template <int Digits>
class DecimalBigInteger : public BigInteger< ten_pow(Digits) > {
  using BigInt = BigInteger< ten_pow(Digits) >;
public:
  DecimalBigInteger(int a) : BigInt(a) {}
  DecimalBigInteger(const BigInt& b) : BigInt(b) {}
  void print() const {
    printf("%llu", BigInt::back());
    char str[Digits + 1] = {};
    for (int i = BigInt::size() - 2; i >= 0; --i) {
      auto a = (*this)[i];
      for (int j = 0; j < Digits; ++j) {
        str[Digits - 1 - j] = a % 10 + '0';
        a /= 10;
      }
      for (int j = 0; j < Digits; ++j) putchar(str[j]);
    }
    puts("");
  }
};
using BigInt = DecimalBigInteger<15>;

pair<BigInt, BigInt> fib(int n) {
  if (n <= 2) {
    return {BigInt((n + 1) >> 1), BigInt((n + 2) >> 1)};
  }
  auto res = fib((n >> 1) - 1);
  auto &a = res.first, &b = res.second;
  BigInt aa = a * a, bb = b * b, bb2 = bb + bb;
  BigInt c = bb + aa, d = (bb2 + bb2) - aa;
  if (n & 2) d = d - 2;
  else d = d + 2;
  if (n & 1) {
    return {d, d + d - c};
  } else {
    return {d - c, d};
  }
}

void solve() {
  int N;
  while (~scanf("%d", &N)) {
    if (N == 2) {
      puts("3");
      puts("INF");
    } else {
      BigInt ans(0);
      if (N & 1) {
        ans = fib(N).first;
      } else {
        auto res = fib(N / 2 - 1);
        ans = res.first * res.second;
        ans = ans + ans;
      }
      printf("%d\n", N);
      ans.print();
    }
  }
}

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

結果

テストケース

ソースコード
