fork download
copy 
  1. #include <iostream>
  2. #include <vector>
  3. #include <algorithm>
  4. #include <map>
  5. using namespace std;
  6.  
  7. #define COUT(x) cout << #x << " = " << (x) << " (L" << __LINE__ << ")" << endl
  8. template<class T1, class T2> ostream& operator << (ostream &s, pair<T1,T2> P)
  9. { return s << '<' << P.first << ", " << P.second << '>'; }
  10. template<class T> ostream& operator << (ostream &s, vector<T> P)
  11. { for (int i = 0; i < P.size(); ++i) { if (i > 0) { s << " "; } s << P[i]; } return s; }
  12. template<class T> ostream& operator << (ostream &s, vector<vector<T> > P)
  13. { for (int i = 0; i < P.size(); ++i) { s << endl << P[i]; } return s << endl; }
  14.  
  15.  
  16. vector<int> calc(vector<int> f, long long p) {
  17.     int N = (int)f.size();
  18.     vector<int> res(N);
  19.     for (int i = 0; i < N; ++i) res[i] = i;
  20.     for (int i = 0; i < p; ++i) {
  21.         vector<int> tmp(N);
  22.         for (int j = 0; j < N; ++j) tmp[j] = f[res[j]];
  23.         res = tmp;
  24.     }
  25.     return res;
  26. }
  27.  
  28. int sisuu(vector<int> f, int i) {
  29.     int v = f[i];
  30.     int con = 1;
  31.     while (v != i) v = f[v], ++con;
  32.     return con;
  33. }
  34.  
  35. int max_sisuu(vector<int> f) {
  36.     int res = 1;
  37.     for (int i = 0; i < f.size(); ++i) {
  38.         res = max(res, sisuu(f, i));
  39.     }
  40.     return res;
  41. }
  42.  
  43.  
  44. vector<long long> all(int N, long long t) {
  45.     int MOD = 1000000007;
  46.  
  47.     vector<int> f(N, 0), g(N, 0);
  48.     for (int i = 0; i < N; ++i) f[i] = g[i] = i;
  49.  
  50.     map<vector<int>, long long> maf, mag;
  51.     do {
  52.         maf[calc(f, t)]++;
  53.     } while (next_permutation(f.begin(), f.end()));
  54.  
  55.     do {
  56.         mag[calc(g, 2)]++;
  57.     } while (next_permutation(g.begin(), g.end()));
  58.  
  59.     vector<long long> res(N+1, 0);
  60.     for (auto it : maf) {
  61.         int si = max_sisuu(it.first);
  62.  
  63.         for (int k = si; k <= N; ++k) {
  64.             res[k] += it.second * mag[it.first];
  65.         }
  66.     }
  67.  
  68.     return res;
  69. }
  70.  
  71. int sub_all(int N, long long p, long long e) {
  72.     int res = 0;
  73.     vector<int> f(N), ans(N);
  74.     for (int i = 0; i < N; ++i) f[i] = i;
  75.     do {
  76.         auto g = calc(f, p);
  77.         bool ok = true;
  78.         for (int k = 0; k < N/e; ++k) {
  79.             for (int i = 0; i < e; ++i) {
  80.                 if (g[k*e+i] != (k*e) + (i+1)%e) ok = false;
  81.             }
  82.         }
  83.         if (ok) ++res;
  84.     } while (next_permutation(f.begin(), f.end()));
  85.  
  86.     return res;
  87. }
  88.  
  89.  
  90.  
  91.  
  92.  
  93. long long GCD(long long x, long long y) {
  94.     if (y == 0) return x;
  95.     else return GCD(y, x % y);
  96. }
  97.  
  98. // modint: mod 計算を int を扱うように扱える構造体
  99. template<int MOD> struct Fp {
  100.     long long val;
  101.     constexpr Fp(long long v = 0) noexcept : val(v % MOD) {
  102.         if (val < 0) val += MOD;
  103.     }
  104.     constexpr int getmod() { return MOD; }
  105.     constexpr Fp operator - () const noexcept {
  106.         return val ? MOD - val : 0;
  107.     }
  108.     constexpr Fp operator + (const Fp& r) const noexcept { return Fp(*this) += r; }
  109.     constexpr Fp operator - (const Fp& r) const noexcept { return Fp(*this) -= r; }
  110.     constexpr Fp operator * (const Fp& r) const noexcept { return Fp(*this) *= r; }
  111.     constexpr Fp operator / (const Fp& r) const noexcept { return Fp(*this) /= r; }
  112.     constexpr Fp& operator += (const Fp& r) noexcept {
  113.         val += r.val;
  114.         if (val >= MOD) val -= MOD;
  115.         return *this;
  116.     }
  117.     constexpr Fp& operator -= (const Fp& r) noexcept {
  118.         val -= r.val;
  119.         if (val < 0) val += MOD;
  120.         return *this;
  121.     }
  122.     constexpr Fp& operator *= (const Fp& r) noexcept {
  123.         val = val * r.val % MOD;
  124.         return *this;
  125.     }
  126.     constexpr Fp& operator /= (const Fp& r) noexcept {
  127.         long long a = r.val, b = MOD, u = 1, v = 0;
  128.         while (b) {
  129.             long long t = a / b;
  130.             a -= t * b; swap(a, b);
  131.             u -= t * v; swap(u, v);
  132.         }
  133.         val = val * u % MOD;
  134.         if (val < 0) val += MOD;
  135.         return *this;
  136.     }
  137.     constexpr bool operator == (const Fp& r) const noexcept {
  138.         return this->val == r.val;
  139.     }
  140.     constexpr bool operator != (const Fp& r) const noexcept {
  141.         return this->val != r.val;
  142.     }
  143.     friend constexpr ostream& operator << (ostream &os, const Fp<MOD>& x) noexcept {
  144.         return os << x.val;
  145.     }
  146.     friend constexpr Fp<MOD> modpow(const Fp<MOD> &a, long long n) noexcept {
  147.         if (n == 0) return 1;
  148.         auto t = modpow(a, n / 2);
  149.         t = t * t;
  150.         if (n & 1) t = t * a;
  151.         return t;
  152.     }
  153. };
  154.  
  155. // 二項係数ライブラリ
  156. template<class T> struct BiCoef {
  157.     vector<T> fact_, inv_, finv_;
  158.     constexpr BiCoef() {}
  159.     constexpr BiCoef(int n) noexcept : fact_(n, 1), inv_(n, 1), finv_(n, 1) {
  160.         init(n);
  161.     }
  162.     constexpr void init(int n) noexcept {
  163.         fact_.assign(n, 1), inv_.assign(n, 1), finv_.assign(n, 1);
  164.         int MOD = fact_[0].getmod();
  165.         for(int i = 2; i < n; i++){
  166.             fact_[i] = fact_[i-1] * i;
  167.             inv_[i] = -inv_[MOD%i] * (MOD/i);
  168.             finv_[i] = finv_[i-1] * inv_[i];
  169.         }
  170.     }
  171.     constexpr T com(int n, int k) const noexcept {
  172.         if (n < k || n < 0 || k < 0) return 0;
  173.         return fact_[n] * finv_[k] * finv_[n-k];
  174.     }
  175.     constexpr T fact(int n) const noexcept {
  176.         if (n < 0) return 0;
  177.         return fact_[n];
  178.     }
  179.     constexpr T inv(int n) const noexcept {
  180.         if (n < 0) return 0;
  181.         return inv_[n];
  182.     }
  183.     constexpr T finv(int n) const noexcept {
  184.         if (n < 0) return 0;
  185.         return finv_[n];
  186.     }
  187. };
  188.  
  189. const int MOD = 1000000007;
  190. using mint = Fp<MOD>;
  191. BiCoef<mint> bc;
  192.  
  193. // e*k 人を、e 人の k グループに分ける方法の数
  194. mint grouping(int e, int k) {
  195.     return bc.fact(k*e) * bc.finv(k) * modpow(bc.finv(e), k);
  196. }
  197.  
  198. // res[k] = h を位数 e の巡回群の k 個の直積からなる 1 つの順列としたときに
  199. // f^p = h となる f が何個あるか
  200. vector<mint> sub(int N, long long p, long long e) {
  201.     int K = N / e;
  202.  
  203.     // f を構成する巡回群の位数 d として、e = d / GCD(d, p) となる d のみを考える
  204.     vector<long long> D;
  205.     for (long long d = e; d <= N; d += e) {
  206.         if (e * GCD(d, p) == d) D.push_back(d);
  207.     }
  208.  
  209.     // DP
  210.     vector<vector<mint>> dp(D.size()+1, vector<mint>(K+1, 0));
  211.     dp[0][0] = 1;
  212.     for (int i = 0; i < D.size(); ++i) {
  213.         int d = D[i];
  214.         int g = d / e;
  215.         mint mul = bc.fact(g - 1) * modpow(mint(e), g - 1);
  216.         for (int k = 0; k <= K; ++k) {
  217.             mint fac = 1;
  218.             for (int l = 0; k + l * g <= K; ++l) {
  219.                 dp[i + 1][k + l * g] += dp[i][k] *
  220.                     bc.com(k + l * g, k) * grouping(g, l) * fac;
  221.                 fac *= mul;
  222.             }
  223.         }
  224.     }
  225.  
  226.     /*
  227.     COUT("---------------------");
  228.     COUT(N); COUT(p); COUT(e);
  229.     COUT(D);
  230.     COUT(dp);
  231.     */
  232.  
  233.     return dp[D.size()];
  234. }
  235.  
  236. mint solve(int N, long long X, long long Y, long long Z) {
  237.     if (X + Y + Z == 0) {
  238.         auto gdp = sub(N, 2, 1);
  239.         return gdp[N] * modpow(mint(N), N);
  240.     }
  241.     long long t = X - Y + Z;
  242.     if (t < 0) t = -t;
  243.     if (t == 0) {
  244.         auto gdp = sub(N, 2, 1);
  245.         return gdp[N] * bc.fact(N);
  246.     }
  247.  
  248.     vector<vector<mint>> dp(N+1, vector<mint>(N+1, 0));
  249.     dp[0][0] = 1;
  250.     for (int e = 1; e <= N; ++e) {
  251.         const vector<mint>& fdp = sub(N, t, e);
  252.         const vector<mint>& gdp = sub(N, 2, e);
  253.         for (int n = 0; n <= N; ++n) {
  254.             mint fac = 1;
  255.             for (int k = 0; n + k * e <= N; ++k) {
  256.                 dp[e][n + k * e] += dp[e - 1][n] *
  257.                     bc.com(n + k * e, n) * grouping(e, k) * fac * fdp[k] * gdp[k];
  258.                 fac *= bc.fact(e - 1);
  259.  
  260.                 /*
  261.                 COUT("------------");
  262.                 COUT(e);
  263.                 COUT(k);
  264.                 cout << n << " -> " << n + k*e << endl;
  265.  
  266.                 COUT( bc.com(n + k * e, n));
  267.                 COUT(grouping(e, k));
  268.                 COUT(modpow(bc.fact(e - 1), k));
  269.                 COUT(fdp[k]);
  270.                 COUT(gdp[k]);
  271.                 */
  272.             }
  273.         }
  274.     }
  275.  
  276.     COUT(dp);
  277.  
  278.     return dp[N][N];
  279. }
  280.  
  281. int main() {
  282.     bc.init(5100);
  283.     int N;
  284.     long long X, Y, Z;
  285.     while (cin >> N >> X >> Y >> Z) {
  286.         cout << solve(N, X, Y, Z) << endl;
  287.  
  288.         /*
  289.         long long t = X - Y + Z;
  290.         if (t < 0) t = -t;
  291.         cout << all(N, t) << endl;
  292.         */
  293.     }
  294. }
  295.  
Success #stdin #stdout 0s 4172KB
comments (?)
stdin
copy 
8 3 0 0
stdout
copy 
dp = 
1 0 0 0 0 0 0 0 0
1 1 2 12 90 546 6156 81432 942012
1 1 2 12 96 576 6336 83952 1021392
1 1 2 12 96 576 6336 83952 1021392
1 1 2 12 96 576 6336 83952 1026432
1 1 2 12 96 600 6480 84960 1042560
1 1 2 12 96 600 6480 84960 1042560
1 1 2 12 96 600 6480 85680 1048320
1 1 2 12 96 600 6480 85680 1048320
 (L276)
1048320
https://ideone.com/7ApGye
language:
C++14 (gcc 8.3)
created:
5 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!