math/rational-binomial.hpp

• View this file on GitHub
• Last update: 2023-12-18 23:52:12+09:00
• Include: #include "math/rational-binomial.hpp"

Depends on

• internal/internal-type-traits.hpp
• binary GCD (math-fast/gcd.hpp)
• math/rational.hpp

Required by

• math/bigint-all.hpp
• math/rational-fps.hpp

Verified with

• verify/verify-unit-test/bigrational.test.cpp
• verify/verify-unit-test/rational-number.test.cpp

Code

#pragma once

#include "rational.hpp"

template <typename R = Rational>
struct Binomial_rational {
  vector<R> fc;
  Binomial_rational(int = 0) { fc.emplace_back(1); }
  void extend() {
    int n = fc.size();
    R nxt = fc.back() * n;
    fc.push_back(nxt);
  }
  R fac(int n) {
    if (n < 0) return 0;
    while ((int)fc.size() <= n) extend();
    return fc[n];
  }
  R finv(int n) {
    if (n < 0) return 0;
    return fac(n).inverse();
  }
  R inv(int n) {
    if (n < 0) return -inv(-n);
    return R{1, max(n, 1)};
  }
  R C(int n, int r) {
    if (n < 0 or r < 0 or n < r) return R{0};
    return fac(n) * finv(n - r) * finv(r);
  }
  R operator()(int n, int r) { return C(n, r); }
  template <typename I>
  R multinomial(const vector<I>& r) {
    static_assert(is_integral<I>::value == true);
    int n = 0;
    for (auto& x : r) {
      if (x < 0) return R{0};
      n += x;
    }
    R res = fac(n);
    for (auto& x : r) res *= finv(x);
    return res;
  }

  template <typename I>
  R operator()(const vector<I>& r) {
    return multinomial(r);
  }
  
  R P(int n, int r) {
    if (n < 0 || n < r || r < 0) return R(0);
    return fac(n) * finv(n - r);
  }
  // [x^r] 1 / (1-x)^n
  R H(int n, int r) {
    if (n < 0 || r < 0) return R(0);
    return r == 0 ? 1 : C(n + r - 1, r);
  }
};

#line 2 "math/rational-binomial.hpp"

#line 2 "math/rational.hpp"

#include <cassert>
#include <numeric>
#include <vector>
using namespace std;

#line 2 "internal/internal-type-traits.hpp"

#include <type_traits>
using namespace std;

namespace internal {
template <typename T>
using is_broadly_integral =
    typename conditional_t<is_integral_v<T> || is_same_v<T, __int128_t> ||
                               is_same_v<T, __uint128_t>,
                           true_type, false_type>::type;

template <typename T>
using is_broadly_signed =
    typename conditional_t<is_signed_v<T> || is_same_v<T, __int128_t>,
                           true_type, false_type>::type;

template <typename T>
using is_broadly_unsigned =
    typename conditional_t<is_unsigned_v<T> || is_same_v<T, __uint128_t>,
                           true_type, false_type>::type;

#define ENABLE_VALUE(x) \
  template <typename T> \
  constexpr bool x##_v = x<T>::value;

ENABLE_VALUE(is_broadly_integral);
ENABLE_VALUE(is_broadly_signed);
ENABLE_VALUE(is_broadly_unsigned);
#undef ENABLE_VALUE

#define ENABLE_HAS_TYPE(var)                                   \
  template <class, class = void>                               \
  struct has_##var : false_type {};                            \
  template <class T>                                           \
  struct has_##var<T, void_t<typename T::var>> : true_type {}; \
  template <class T>                                           \
  constexpr auto has_##var##_v = has_##var<T>::value;

#define ENABLE_HAS_VAR(var)                                     \
  template <class, class = void>                                \
  struct has_##var : false_type {};                             \
  template <class T>                                            \
  struct has_##var<T, void_t<decltype(T::var)>> : true_type {}; \
  template <class T>                                            \
  constexpr auto has_##var##_v = has_##var<T>::value;

}  // namespace internal
#line 2 "math-fast/gcd.hpp"

#include <algorithm>
using namespace std;

namespace BinaryGCDImpl {
using u64 = unsigned long long;
using i8 = char;

u64 binary_gcd(u64 a, u64 b) {
  if (a == 0 || b == 0) return a + b;
  i8 n = __builtin_ctzll(a);
  i8 m = __builtin_ctzll(b);
  a >>= n;
  b >>= m;
  n = min(n, m);
  while (a != b) {
    u64 d = a - b;
    i8 s = __builtin_ctzll(d);
    bool f = a > b;
    b = f ? b : a;
    a = (f ? d : -d) >> s;
  }
  return a << n;
}

using u128 = __uint128_t;
// a > 0
int ctz128(u128 a) {
  u64 lo = a & u64(-1);
  return lo ? __builtin_ctzll(lo) : 64 + __builtin_ctzll(a >> 64);
}
u128 binary_gcd128(u128 a, u128 b) {
  if (a == 0 || b == 0) return a + b;
  i8 n = ctz128(a);
  i8 m = ctz128(b);
  a >>= n;
  b >>= m;
  n = min(n, m);
  while (a != b) {
    u128 d = a - b;
    i8 s = ctz128(d);
    bool f = a > b;
    b = f ? b : a;
    a = (f ? d : -d) >> s;
  }
  return a << n;
}

}  // namespace BinaryGCDImpl

long long binary_gcd(long long a, long long b) {
  return BinaryGCDImpl::binary_gcd(abs(a), abs(b));
}
__int128_t binary_gcd128(__int128_t a, __int128_t b) {
  if (a < 0) a = -a;
  if (b < 0) b = -b;
  return BinaryGCDImpl::binary_gcd128(a, b);
}

/**
 * @brief binary GCD
 */
#line 10 "math/rational.hpp"

// T : 値, U : 比較用
template <typename T, typename U>
struct RationalBase {
  using R = RationalBase;
  using Key = T;
  T x, y;
  RationalBase() : x(0), y(1) {}
  template <typename T1>
  RationalBase(const T1& _x) : RationalBase<T, U>(_x, T1{1}) {}
  template <typename T1, typename T2>
  RationalBase(const pair<T1, T2>& _p)
      : RationalBase<T, U>(_p.first, _p.second) {}
  template <typename T1, typename T2>
  RationalBase(const T1& _x, const T2& _y) : x(_x), y(_y) {
    assert(y != 0);
    if (y == -1) x = -x, y = -y;
    if (y != 1) {
      T g;
      if constexpr (internal::is_broadly_integral_v<T>) {
        if constexpr (sizeof(T) == 16) {
          g = binary_gcd128(x, y);
        } else {
          g = binary_gcd(x, y);
        }
      } else {
        g = gcd(x, y);
      }
      if (g != 0) x /= g, y /= g;
      if (y < 0) x = -x, y = -y;
    }
  }
  // y = 0 の代入も認める
  static R raw(T _x, T _y) {
    R r;
    r.x = _x, r.y = _y;
    return r;
  }
  friend R operator+(const R& l, const R& r) {
    if (l.y == r.y) return R{l.x + r.x, l.y};
    return R{l.x * r.y + l.y * r.x, l.y * r.y};
  }
  friend R operator-(const R& l, const R& r) {
    if (l.y == r.y) return R{l.x - r.x, l.y};
    return R{l.x * r.y - l.y * r.x, l.y * r.y};
  }
  friend R operator*(const R& l, const R& r) { return R{l.x * r.x, l.y * r.y}; }
  friend R operator/(const R& l, const R& r) { return R{l.x * r.y, l.y * r.x}; }
  R& operator+=(const R& r) { return (*this) = (*this) + r; }
  R& operator-=(const R& r) { return (*this) = (*this) - r; }
  R& operator*=(const R& r) { return (*this) = (*this) * r; }
  R& operator/=(const R& r) { return (*this) = (*this) / r; }
  R operator-() const { return raw(-x, y); }
  R inverse() const {
    assert(x != 0);
    R r = raw(y, x);
    if (r.y < 0) r.x = -r.x, r.y = -r.y;
    return r;
  }
  R pow(long long p) const {
    R res{1}, base{*this};
    while (p) {
      if (p & 1) res *= base;
      base *= base;
      p >>= 1;
    }
    return res;
  }
  friend bool operator==(const R& l, const R& r) {
    return l.x == r.x && l.y == r.y;
  };
  friend bool operator!=(const R& l, const R& r) {
    return l.x != r.x || l.y != r.y;
  };
  friend bool operator<(const R& l, const R& r) {
    return U{l.x} * r.y < U{l.y} * r.x;
  };
  friend bool operator<=(const R& l, const R& r) { return l < r || l == r; }
  friend bool operator>(const R& l, const R& r) {
    return U{l.x} * r.y > U{l.y} * r.x;
  };
  friend bool operator>=(const R& l, const R& r) { return l > r || l == r; }
  friend ostream& operator<<(ostream& os, const R& r) {
    os << r.x;
    if (r.x != 0 && r.y != 1) os << "/" << r.y;
    return os;
  }

  // T にキャストされるので T が bigint の場合は to_ll も要る
  T to_mint(T mod) const {
    assert(mod != 0);
    T a = y, b = mod, u = 1, v = 0, t;
    while (b > 0) {
      t = a / b;
      swap(a -= t * b, b);
      swap(u -= t * v, v);
    }
    return U((u % mod + mod) % mod) * x % mod;
  }
};

using Rational = RationalBase<long long, __int128_t>;
#line 4 "math/rational-binomial.hpp"

template <typename R = Rational>
struct Binomial_rational {
  vector<R> fc;
  Binomial_rational(int = 0) { fc.emplace_back(1); }
  void extend() {
    int n = fc.size();
    R nxt = fc.back() * n;
    fc.push_back(nxt);
  }
  R fac(int n) {
    if (n < 0) return 0;
    while ((int)fc.size() <= n) extend();
    return fc[n];
  }
  R finv(int n) {
    if (n < 0) return 0;
    return fac(n).inverse();
  }
  R inv(int n) {
    if (n < 0) return -inv(-n);
    return R{1, max(n, 1)};
  }
  R C(int n, int r) {
    if (n < 0 or r < 0 or n < r) return R{0};
    return fac(n) * finv(n - r) * finv(r);
  }
  R operator()(int n, int r) { return C(n, r); }
  template <typename I>
  R multinomial(const vector<I>& r) {
    static_assert(is_integral<I>::value == true);
    int n = 0;
    for (auto& x : r) {
      if (x < 0) return R{0};
      n += x;
    }
    R res = fac(n);
    for (auto& x : r) res *= finv(x);
    return res;
  }

  template <typename I>
  R operator()(const vector<I>& r) {
    return multinomial(r);
  }
  
  R P(int n, int r) {
    if (n < 0 || n < r || r < 0) return R(0);
    return fac(n) * finv(n - r);
  }
  // [x^r] 1 / (1-x)^n
  R H(int n, int r) {
    if (n < 0 || r < 0) return R(0);
    return r == 0 ? 1 : C(n + r - 1, r);
  }
};

Back to top page