Submission #17173142

---

Source Code Expand

#include <algorithm>
#include <array>
#include <cassert>
#include <iostream>
#include <vector>
using namespace std;

template<typename T> ostream& operator<<(ostream &os, const vector<T> &v) { os << '{'; string sep; for (const auto &x : v) os << sep << x, sep = ", "; return os << '}'; }
template<typename T, size_t size> ostream& operator<<(ostream &os, const array<T, size> &arr) { os << '{'; string sep; for (const auto &x : arr) os << sep << x, sep = ", "; return os << '}'; }
template<typename A, typename B> ostream& operator<<(ostream &os, const pair<A, B> &p) { return os << '(' << p.first << ", " << p.second << ')'; }

void dbg_out() { cerr << endl; }
template<typename Head, typename... Tail> void dbg_out(Head H, Tail... T) { cerr << ' ' << H; dbg_out(T...); }

#ifdef NEAL_DEBUG
#define dbg(...) cerr << "(" << #__VA_ARGS__ << "):", dbg_out(__VA_ARGS__)
#else
#define dbg(...)
#endif

template<const int &MOD>
struct _m_int {
    int val;

    _m_int(int64_t v = 0) {
        if (v < 0) v = v % MOD + MOD;
        if (v >= MOD) v %= MOD;
        val = int(v);
    }

    _m_int(uint64_t v) {
        if (v >= MOD) v %= MOD;
        val = int(v);
    }

    _m_int(int v) : _m_int(int64_t(v)) {}
    _m_int(unsigned v) : _m_int(uint64_t(v)) {}

    static int inv_mod(int a, int m = MOD) {
        // https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm#Example
        int g = m, r = a, x = 0, y = 1;

        while (r != 0) {
            int q = g / r;
            g %= r; swap(g, r);
            x -= q * y; swap(x, y);
        }

        return x < 0 ? x + m : x;
    }

    explicit operator int() const { return val; }
    explicit operator unsigned() const { return val; }
    explicit operator int64_t() const { return val; }
    explicit operator uint64_t() const { return val; }
    explicit operator double() const { return val; }
    explicit operator long double() const { return val; }

    _m_int& operator+=(const _m_int &other) {
        val -= MOD - other.val;
        if (val < 0) val += MOD;
        return *this;
    }

    _m_int& operator-=(const _m_int &other) {
        val -= other.val;
        if (val < 0) val += MOD;
        return *this;
    }

    static unsigned fast_mod(uint64_t x, unsigned m = MOD) {
#if !defined(_WIN32) || defined(_WIN64)
        return unsigned(x % m);
#endif
        // Optimized mod for Codeforces 32-bit machines.
        // x must be less than 2^32 * m for this to work, so that x / m fits in an unsigned 32-bit int.
        unsigned x_high = unsigned(x >> 32), x_low = unsigned(x);
        unsigned quot, rem;
        asm("divl %4\n"
            : "=a" (quot), "=d" (rem)
            : "d" (x_high), "a" (x_low), "r" (m));
        return rem;
    }

    _m_int& operator*=(const _m_int &other) {
        val = fast_mod(uint64_t(val) * other.val);
        return *this;
    }

    _m_int& operator/=(const _m_int &other) {
        return *this *= other.inv();
    }

    friend _m_int operator+(const _m_int &a, const _m_int &b) { return _m_int(a) += b; }
    friend _m_int operator-(const _m_int &a, const _m_int &b) { return _m_int(a) -= b; }
    friend _m_int operator*(const _m_int &a, const _m_int &b) { return _m_int(a) *= b; }
    friend _m_int operator/(const _m_int &a, const _m_int &b) { return _m_int(a) /= b; }

    _m_int& operator++() {
        val = val == MOD - 1 ? 0 : val + 1;
        return *this;
    }

    _m_int& operator--() {
        val = val == 0 ? MOD - 1 : val - 1;
        return *this;
    }

    _m_int operator++(int) { _m_int before = *this; ++*this; return before; }
    _m_int operator--(int) { _m_int before = *this; --*this; return before; }

    _m_int operator-() const {
        return val == 0 ? 0 : MOD - val;
    }

    friend bool operator==(const _m_int &a, const _m_int &b) { return a.val == b.val; }
    friend bool operator!=(const _m_int &a, const _m_int &b) { return a.val != b.val; }
    friend bool operator<(const _m_int &a, const _m_int &b) { return a.val < b.val; }
    friend bool operator>(const _m_int &a, const _m_int &b) { return a.val > b.val; }
    friend bool operator<=(const _m_int &a, const _m_int &b) { return a.val <= b.val; }
    friend bool operator>=(const _m_int &a, const _m_int &b) { return a.val >= b.val; }

    _m_int inv() const {
        return inv_mod(val);
    }

    _m_int pow(int64_t p) const {
        if (p < 0)
            return inv().pow(-p);

        _m_int a = *this, result = 1;

        while (p > 0) {
            if (p & 1)
                result *= a;

            p >>= 1;

            if (p > 0)
                a *= a;
        }

        return result;
    }

    friend ostream& operator<<(ostream &os, const _m_int &m) {
        return os << m.val;
    }
};

extern const int MOD = 1e9 + 7;
using mod_int = _m_int<MOD>;


vector<mod_int> inv, factorial, inv_factorial;
int prepared_maximum = -1;

void prepare_factorials(int64_t maximum) {
    static int prepare_calls = 0;

    if (prepare_calls++ == 0) {
        // Make sure MOD is prime, which is necessary for the inverse algorithm below.
        for (int p = 2; p * p <= MOD; p += p % 2 + 1)
            assert(MOD % p != 0);

        inv = {0, 1};
        factorial = inv_factorial = {1, 1};
        prepared_maximum = 1;
    }

    if (maximum > prepared_maximum) {
        inv.resize(maximum + 1);
        factorial.resize(maximum + 1);
        inv_factorial.resize(maximum + 1);

        for (int i = prepared_maximum + 1; i <= maximum; i++) {
            inv[i] = inv[MOD % i] * (MOD - MOD / i);
            factorial[i] = i * factorial[i - 1];
            inv_factorial[i] = inv[i] * inv_factorial[i - 1];
        }

        prepared_maximum = int(maximum);
    }
}

mod_int choose(int64_t n, int64_t r) {
    if (r < 0 || r > n) return 0;
    return factorial[n] * inv_factorial[r] * inv_factorial[n - r];
}

mod_int inv_choose(int64_t n, int64_t r) {
    assert(0 <= r && r <= n);
    return inv_factorial[n] * factorial[r] * factorial[n - r];
}

mod_int permute(int64_t n, int64_t r) {
    if (r < 0 || r > n) return 0;
    return factorial[n] * inv_factorial[n - r];
}

mod_int inv_permute(int64_t n, int64_t r) {
    assert(0 <= r && r <= n);
    return inv_factorial[n] * factorial[n - r];
}


// https://en.wikipedia.org/wiki/Bernoulli_number
vector<mod_int> bernoulli;

void prepare_bernoulli(int maximum) {
    bernoulli.resize(maximum + 1);

    for (int m = 0; m <= maximum; m++) {
        bernoulli[m] = 1;

        for (int k = 0; k < m; k++)
            bernoulli[m] -= choose(m, k) * bernoulli[k] * inv[m - k + 1];
    }
}

using polynomial = vector<mod_int>;

mod_int evaluate(const polynomial &poly, mod_int x) {
    mod_int result = 0, power = 1;

    for (int i = 0; i < int(poly.size()); i++) {
        result += poly[i] * power;
        power *= x;
    }

    return result;
}

polynomial& operator+=(polynomial &a, const polynomial &b) {
    a.resize(max(a.size(), b.size()), 0);

    for (int i = 0; i < int(b.size()); i++)
        a[i] += b[i];

    return a;
}

polynomial operator+(polynomial a, const polynomial &b) {
    return a += b;
}

polynomial& operator-=(polynomial &a, const polynomial &b) {
    a.resize(max(a.size(), b.size()), 0);

    for (int i = 0; i < int(b.size()); i++)
        a[i] -= b[i];

    return a;
}

polynomial operator-(polynomial a, const polynomial &b) {
    return a -= b;
}

polynomial& operator*=(polynomial &poly, mod_int mult) {
    for (int i = 0; i < int(poly.size()); i++)
        poly[i] *= mult;

    return poly;
}

polynomial operator*(polynomial poly, mod_int mult) {
    return poly *= mult;
}

polynomial operator*(mod_int mult, const polynomial &poly) {
    return poly * mult;
}

polynomial operator*(const polynomial &a, const polynomial &b) {
    if (a.empty() || b.empty())
        return {};

    // return NTT::mod_multiply(a, b);
    static const uint64_t U64_BOUND = numeric_limits<uint64_t>::max() - uint64_t(MOD) * MOD;
    vector<uint64_t> result(a.size() + b.size() - 1, 0);

    for (int i = 0; i < int(a.size()); i++)
        for (int j = 0; j < int(b.size()); j++) {
            result[i + j] += uint64_t(a[i]) * uint64_t(b[j]);

            if (result[i + j] > U64_BOUND)
                result[i + j] %= MOD;
        }

    for (uint64_t &x : result)
        if (x >= MOD)
            x %= MOD;

    return vector<mod_int>(result.begin(), result.end());
}

polynomial& operator*=(polynomial &a, const polynomial &b) {
    return a = a * b;
}

// Returns 1^k + 2^k + ... + x^k as a polynomial in x.
polynomial sum_of_powers(int k) {
    polynomial result(k + 2, 0);

    // https://en.wikipedia.org/wiki/Faulhaber%27s_formula
    for (int a = 0; a <= k; a++)
        result[k - a + 1] = bernoulli[a] * inv_factorial[a] * factorial[k] * inv_factorial[k - a + 1];

    return result;
}

// Returns the sum P(0) + P(1) + ... + P(x) as a polynomial in x.
polynomial partial_sum(const polynomial &poly) {
    if (poly.empty())
        return {};

    polynomial result(poly.size() + 1, 0);
    result[1] = poly[0];

    // To exclude P(0) from the sum, remove the following line.
    result[0] = poly[0];

    // for (int k = 1; k < int(poly.size()); k++) {
    //     // result += sum_of_powers(k) * poly[k];

    //     for (int a = 0; a <= k; a++)
    //         result[k - a + 1] += bernoulli[a] * inv_factorial[a] * factorial[k] * inv_factorial[k - a + 1] * poly[k];
    // }

    vector<mod_int> a_half(poly.size(), 0);
    vector<mod_int> k_half(poly.size(), 0);

    for (int a = 0; a < int(poly.size()); a++)
        a_half[a] = bernoulli[a] * inv_factorial[a];

    reverse(a_half.begin(), a_half.end());

    for (int k = 1; k < int(poly.size()); k++)
        k_half[k] = factorial[k] * poly[k];

    vector<mod_int> product = a_half * k_half;

    for (int i = 1; i <= int(poly.size()); i++)
        result[i] += product[i + poly.size() - 2] * inv_factorial[i];

    return result;
}


// Finds the length of the longest subsequence of values such that compare is true for all consecutive pairs.
template<typename T, typename T_compare>
int longest_increasing_subsequence(vector<T> values, T_compare &&compare) {
    vector<T> best_ending;

    for (T value : values) {
        auto it = lower_bound(best_ending.begin(), best_ending.end(), value, compare);

        if (it == best_ending.end())
            best_ending.push_back(value);
        else
            *it = value;
    }

    return int(best_ending.size());
}


const int POLY = 100;
const int INF = 1e9 + 10;

int main() {
    prepare_factorials(POLY + 5);
    prepare_bernoulli(POLY);

    auto piecewise_partial_sum = [](vector<pair<array<int, 2>, polynomial>> piecewise, int start, int end) -> vector<pair<array<int, 2>, polynomial>> {
        auto get_sum = [&](int x) {
            polynomial p;

            for (auto &pr : piecewise)
                if (pr.first[0] <= x) {
                    polynomial sum_poly = partial_sum(pr.second);

                    if (x >= pr.first[1])
                        p += polynomial{evaluate(sum_poly, pr.first[1]) - evaluate(sum_poly, pr.first[0] - 1)};
                    else
                        p += sum_poly - polynomial{evaluate(sum_poly, pr.first[0] - 1)};
                }

            return p;
        };

        vector<array<int, 2>> ranges;

        for (auto &pr : piecewise)
            ranges.push_back(pr.first);

        while (!ranges.empty() && start > ranges.front()[1])
            ranges.erase(ranges.begin());

        while (!ranges.empty() && end < ranges.back()[0])
            ranges.pop_back();

        if (ranges.empty()) {
            ranges = {{start, end}};
        } else {
            if (start < ranges.front()[0])
                ranges.insert(ranges.begin(), array<int, 2>{start, ranges.front()[0] - 1});
            else
                ranges.front()[0] = start;

            if (end > ranges.back()[1])
                ranges.push_back({ranges.back()[1] + 1, end});
            else
                ranges.back()[1] = end;
        }

        vector<pair<array<int, 2>, polynomial>> next_piecewise;

        for (auto &range : ranges)
            next_piecewise.emplace_back(range, get_sum(range[0]));

        return next_piecewise;
    };

    auto solve_ranges = [&](vector<array<int, 2>> ranges) {
        ranges.push_back({INF, INF});
        vector<pair<array<int, 2>, polynomial>> piecewise = {{{0, 0}, {1}}};

        for (auto &range : ranges) {
            piecewise = piecewise_partial_sum(piecewise, range[0], range[1]);
            // dbg(range, piecewise);
        }

        assert(piecewise.size() == 1);
        return evaluate(piecewise[0].second, INF);
    };

    auto solve_ranges_slow = [](vector<array<int, 2>> ranges) {
        int big = 0;

        for (auto &range : ranges)
            big = max(big, range[1]);

        vector<mod_int> dp(big + 1, 0);
        dp[0] = 1;

        for (auto &range : ranges) {
            vector<mod_int> next_dp(big + 1, 0);
            mod_int sum = 0;

            for (int i = 0; i < range[0]; i++)
                sum += dp[i];

            for (int i = range[0]; i <= range[1]; i++) {
                sum += dp[i];
                next_dp[i] = sum;
            }

            swap(dp, next_dp);
        }

        mod_int ways = 0;

        for (int i = 0; i <= big; i++)
            ways += dp[i];

        return ways;
    };

    int N;
    cin >> N;
    vector<int> A(N);

    for (auto &a : A)
        cin >> a;

    vector<int> perm(N);

    for (int i = 0; i < N; i++)
        perm[i] = i;

    mod_int answer = 0;

    do {
        vector<array<int, 2>> ranges;
        int extra = 0;

        for (int i = N - 1; i >= 0; i--) {
            ranges.push_back({1 + extra, A[perm[i]] + extra});

            if (i > 0 && perm[i] > perm[i - 1])
                extra++;
        }

        reverse(ranges.begin(), ranges.end());
        answer += solve_ranges(ranges) * longest_increasing_subsequence(perm, less<int>());
    } while (next_permutation(perm.begin(), perm.end()));

    for (int i = 0; i < N; i++)
        answer /= A[i];

    cout << answer << '\n';
}

Submission Info

Compile Error

./Main.cpp: In function ‘int main()’:
./Main.cpp:439:10: warning: variable ‘solve_ranges_slow’ set but not used [-Wunused-but-set-variable]
  439 |     auto solve_ranges_slow = [](vector<array<int, 2>> ranges) {
      |          ^~~~~~~~~~~~~~~~~
./Main.cpp: In instantiation of ‘_m_int<MOD>::_m_int(uint64_t) [with const int& MOD = MOD; uint64_t = long unsigned int]’:
/usr/include/c++/9/bits/stl_construct.h:75:7:   required from ‘void std::_Construct(_T1*, _Args&& ...) [with _T1 = _m_int<MOD>; _Args = {long unsigned int&}]’
/usr/include/c++/9/bits/stl_uninitialized.h:83:18:   required from ‘static _ForwardIterator std::__uninitialized_copy<_TrivialValueTypes>::__uninit_copy(_InputIterator, _InputIterator, _ForwardIterator) [with _InputIterator = __gnu_cxx::__normal_iterator<long unsigned int*, std::vector<long unsigned int> >; _ForwardIterator = _m_int<MOD>*; bool _TrivialValueTypes = false]’
/usr/include/c++/9/bits/stl_uninitialized.h:140:15:   required from ‘_ForwardIterator std::uninitialized_copy(_InputIterator, _InputIterator, _ForwardIterator) [with _InputIterator = __gnu_cxx::__normal_iterator<long unsigned int*, std::vector<long unsigned int> >; _ForwardIterator = _m_int<MOD>*]’
/usr/include/c++/9/bits/stl_uninitialized.h:307:37:   required from ‘_ForwardIterator std::__uninitialized_copy_a(_InputIterator, _InputIterator, _ForwardIterator, std::allocator<_Tp>&) [with _InputIterator = __gnu_cxx::__normal_iterator<long unsigned int*, std::vector<long unsigned int> >; _ForwardIterator = _m_int<MOD>*; _Tp = _m_int<MOD>]’
/usr/include/c++/9/bits/stl_vector.h:1582:33:   required from ‘void std::vector<_Tp, _Alloc>::_M_range_initialize(_ForwardIterator, _ForwardIterator, std::forward_iterator_tag) [with _ForwardIterator = __gnu_cxx::__normal_iterator<long unsigned int*, std::vector<long unsigned int> >; _Tp = _m_int<MOD>; _Alloc = std::allocator<_m_int<MOD> >]’
/usr/include/c++/9/bits/stl_vector.h:654:4:   required from ‘std::vector<_Tp, _Alloc>::vector(_InputIterator, _InputIter...

Judge Result