// #define _GLIBCXX_DEBUG #include // clang-format off std::ostream&operator<<(std::ostream&os,std::int8_t x){return os<<(int)x;} std::ostream&operator<<(std::ostream&os,std::uint8_t x){return os<<(int)x;} std::ostream&operator<<(std::ostream&os,const __int128_t &u){if(!u)os<<"0";__int128_t tmp=u<0?(os<<"-",-u):u;std::string s;while(tmp)s+='0'+(tmp%10),tmp/=10;return std::reverse(s.begin(),s.end()),os< // clang-format off templatestruct make_long{using type= T;}; template<>struct make_long{using type= int16_t;}; template<>struct make_long{using type= uint16_t;}; template<>struct make_long{using type= int32_t;}; template<>struct make_long{using type= uint32_t;}; template<>struct make_long{using type= int64_t;}; template<>struct make_long{using type= uint64_t;}; template<>struct make_long{using type= __int128_t;}; template<>struct make_long{using type= __uint128_t;}; template<>struct make_long{using type= double;}; template<>struct make_long{using type= long double;}; template using make_long_t= typename make_long::type; // clang-format on template class PiecewiseLinearConvex { using D= make_long_t; struct Node { int ch[2]= {0, 0}; T z= 0, x= 0, d= 0, a= 0; D s= 0; size_t sz= 0; friend std::ostream &operator<<(std::ostream &os, const Node &t) { return os << "{z:" << t.z << ",x:" << t.x << ",d:" << t.d << ",a:" << t.a << ",s:" << t.s << ",sz:" << t.sz << ",ch:(" << t.ch[0] << "," << t.ch[1] << ")}"; } }; static inline size_t ni= 1; static inline Node *n= new Node[NODE_SIZE]{Node{}}; static inline void info(int t, int d, std::stringstream &ss) { if (!t) return; info(n[t].ch[0], d + 1, ss); for (int i= 0; i < d; ++i) ss << " "; ss << " ■ " << n[t] << '\n', info(n[t].ch[1], d + 1, ss); } static inline void dump_xs(int t, std::vector &xs) { if (t) push(t), dump_xs(n[t].ch[0], xs), xs.push_back(n[t].x), dump_xs(n[t].ch[1], xs); } static inline void dump_slopes_l(int t, T ofs, std::vector &as) { if (t) push(t), dump_slopes_l(n[t].ch[1], ofs, as), ofs+= n[n[t].ch[1]].a + n[t].d, as.push_back(-ofs), dump_slopes_l(n[t].ch[0], ofs, as); } static inline void dump_slopes_r(int t, T ofs, std::vector &as) { if (t) push(t), dump_slopes_r(n[t].ch[0], ofs, as), ofs+= n[n[t].ch[0]].a + n[t].d, as.push_back(ofs), dump_slopes_r(n[t].ch[1], ofs, as); } static inline int create(T d, T x) { return n[ni].d= d, n[ni].x= x, n[ni].z= 0, ni++; } static inline bool lt(T a, T b) { if constexpr (std::is_floating_point_v) return 1e-15 < b - a; else return a < b; } template static inline int build(Iter bg, Iter ed) { if (bg == ed) return 0; auto md= bg + (ed - bg) / 2; int t= create(md->first, md->second); return n[t].ch[0]= build(bg, md), n[t].ch[1]= build(md + 1, ed), update(t), t; } template static inline void dump(Iter itr, int t) { if (!t) return; push(t); size_t sz= n[n[t].ch[0]].sz; dump(itr, n[t].ch[0]), *(itr + sz)= {n[t].d, n[t].x}, dump(itr + sz + 1, n[t].ch[1]); } static inline void update(int t) { int l= n[t].ch[0], r= n[t].ch[1]; n[t].sz= 1 + n[l].sz + n[r].sz, n[t].a= n[t].d + n[l].a + n[r].a, n[t].s= D(n[t].x) * n[t].d + n[l].s + n[r].s; } template static inline void prop(int &t, T v) { if constexpr (persistent && b) { if (!t) return; n[ni]= n[t], t= ni++; } n[t].z+= v, n[t].s+= D(v) * n[t].a, n[t].x+= v; } static inline void push(int t) { if (n[t].z != 0) prop(n[t].ch[0], n[t].z), prop(n[t].ch[1], n[t].z), n[t].z= 0; } template static inline int join_(int t, int a, int b) { push(a); if constexpr (r) b= join<0>(b, t, n[a].ch[0]); else b= join<0>(n[a].ch[1], t, b); if constexpr (persistent) n[ni]= n[a], a= ni++; if (n[n[a].ch[r]].sz * 4 >= n[b].sz) return n[a].ch[!r]= b, update(a), a; return n[a].ch[!r]= n[b].ch[r], update(a), n[b].ch[r]= a, update(b), b; } template static inline int join(int l, int t, int r) { if constexpr (persistent && b) n[ni]= n[t], t= ni++; if (n[l].sz > n[r].sz * 4) return join_<0>(t, l, r); if (n[r].sz > n[l].sz * 4) return join_<1>(t, r, l); return n[t].ch[0]= l, n[t].ch[1]= r, update(t), t; } static inline std::array split(int t, T x) { if (!t) return {0, 0, 0}; push(t); if (lt(n[t].x, x)) { auto [a, b, c]= split(n[t].ch[1], x); return {join(n[t].ch[0], t, a), b, c}; } else if (lt(x, n[t].x)) { auto [a, b, c]= split(n[t].ch[0], x); return {a, b, join(c, t, n[t].ch[1])}; } return {n[t].ch[0], t, n[t].ch[1]}; } static inline int unite(int l, int r) { if (!l) return r; if (!r) return l; push(l); if constexpr (persistent) n[ni]= n[l], l= ni++; auto [a, b, c]= split(r, n[l].x); return n[l].d+= n[b].d, join<0>(unite(a, n[l].ch[0]), l, unite(n[l].ch[1], c)); } static inline int insert(int t, T x, T d) { if (!t) return n[ni]= Node{{0, 0}, 0, x, d, d, D(x) * d, 1}, ni++; push(t); if constexpr (persistent) n[ni]= n[t], t= ni++; if (lt(x, n[t].x)) return join<0>(insert(n[t].ch[0], x, d), t, n[t].ch[1]); if (lt(n[t].x, x)) return join<0>(n[t].ch[0], t, insert(n[t].ch[1], x, d)); return n[t].d+= d, update(t), t; } static inline int erase(int t, T x, T d) { assert(t); push(t); if (lt(x, n[t].x)) return join(erase(n[t].ch[0], x, d), t, n[t].ch[1]); if (lt(n[t].x, x)) return join(n[t].ch[0], t, erase(n[t].ch[1], x, d)); assert(!lt(n[t].d, d)); if (lt(d, n[t].d)) { if constexpr (persistent) n[ni]= n[t], t= ni++; return n[t].d-= d, update(t), t; } if (n[t].ch[1]) { auto [a, s]= pop<0>(n[t].ch[1]); return join(n[t].ch[0], s, a); } return n[t].ch[0]; } template static inline std::pair pop(int t) { if (push(t); !n[t].ch[r]) return {n[t].ch[!r], t}; auto [a, s]= pop(n[t].ch[r]); if constexpr (r) return {join(n[t].ch[!r], t, a), s}; else return {join(a, t, n[t].ch[!r]), s}; } template static inline bool lgt(T a, T b) { if constexpr (g) return lt(b, a); else return lt(a, b); } template static inline int cut(int t, T x) { if (!t) return t; if (push(t); lgt(n[t].x, x)) return cut(n[t].ch[!r], x); if (lgt(x, n[t].x)) { if constexpr (r) return join(n[t].ch[0], t, cut<1>(n[t].ch[1], x)); else return join(cut<0>(n[t].ch[0], x), t, n[t].ch[1]); } return n[t].ch[!r]; } template static inline D calc_y(int t, T x, T ol, D ou) { for (; t;) { if (push(t); lgt(n[t].x, x)) t= n[t].ch[!r]; else { ol+= n[n[t].ch[!r]].a, ou+= n[n[t].ch[!r]].s; if (!lgt(x, n[t].x)) break; ol+= n[t].d, ou+= D(n[t].x) * n[t].d, t= n[t].ch[r]; } } return D(x) * ol - ou; } template static inline std::array split(int t, T p, T &ol, D &ou) { push(t); T s= ol + n[n[t].ch[!r]].a; if (lt(p, s)) { auto [a, b, c]= split(n[t].ch[!r], p, ol, ou); if constexpr (r) return {a, b, join(c, t, n[t].ch[r])}; else return {join(n[t].ch[r], t, a), b, c}; } ol= s + n[t].d; if (lt(ol, p)) { ou+= n[n[t].ch[!r]].s + D(n[t].x) * n[t].d; auto [a, b, c]= split(n[t].ch[r], p, ol, ou); if constexpr (r) return {join(n[t].ch[!r], t, a), b, c}; else return {a, b, join(c, t, n[t].ch[!r])}; } ou+= n[n[t].ch[!r]].s; return {n[t].ch[0], t, n[t].ch[1]}; } template static inline bool lte(T a, T b) { if constexpr (l) return lt(a, b); else return !lt(b, a); } template static inline std::pair split_cum(int t, T p, T &ol, D &ou) { push(t); T s= ol + n[n[t].ch[!r]].a; if (lte(p, s)) { auto [c, b]= split_cum(n[t].ch[!r], p, ol, ou); if constexpr (l) { if constexpr (r) return {join(c, t, n[t].ch[r]), b}; else return {join(n[t].ch[r], t, c), b}; } else return {c, b}; } ol= s + n[t].d; if (lte(ol, p)) { ou+= n[n[t].ch[!r]].s + D(n[t].x) * n[t].d; auto [a, b]= split_cum(n[t].ch[r], p, ol, ou); if constexpr (l) return {a, b}; else { if constexpr (r) return {join(n[t].ch[!r], t, a), b}; else return {join(a, t, n[t].ch[!r]), b}; } } ou+= n[n[t].ch[!r]].s; return {n[t].ch[!r ^ l], t}; } int mn, lr[2]; bool bf[2]; T o[2], rem, bx[2]; D y; inline D calc_y(T x) { if (!mn) return 0; if (lt(x, n[mn].x)) return -calc_y<0>(lr[0], x, o[0], D(n[mn].x) * o[0]); if (lt(n[mn].x, x)) return calc_y<1>(lr[1], x, o[1], D(n[mn].x) * o[1]); return 0; } inline void slope_eval(bool neg) { T p= neg ? -rem : rem, ol= o[neg]; if (p <= ol) o[neg]-= p, o[!neg]+= p, y+= D(n[mn].x) * rem; else { D ou= D(n[mn].x) * ol; auto [a, b, c]= neg ? split<1>(lr[neg], p, ol, ou) : split<0>(lr[neg], p, ol, ou); o[neg]= ol - p, ol-= n[b].d, ou+= D(n[b].x) * (o[!neg]= p - ol); if (neg) y-= ou, lr[!neg]= join(lr[!neg], mn, a), lr[neg]= c; else y+= ou, lr[!neg]= join(c, mn, lr[!neg]), lr[neg]= a; mn= b; } rem= 0; } template inline void slope_eval_cum() { T p= neg ? -rem : rem, ol= o[neg]; if (lte(p, ol)) o[neg]-= p, o[!neg]+= p, y+= D(n[mn].x) * rem; else { D ou= D(n[mn].x) * ol; auto [a, b]= split_cum(lr[neg], p, ol, ou); o[neg]= ol - p, ol-= n[b].d, ou+= D(n[b].x) * (o[!neg]= p - ol); if constexpr (l) lr[neg]= a; else { if constexpr (neg) lr[!neg]= join(lr[!neg], mn, a); else lr[!neg]= join(a, mn, lr[!neg]); } if constexpr (neg) y-= ou; else y+= ou; mn= b; } rem= 0; } template void add_inf(T x0) { if (bf[r] && !lgt(bx[r], x0)) return; if (assert(!bf[!r] || !lgt(bx[!r], x0)), bf[r]= true, bx[r]= x0; !mn) return; if (lgt(x0, n[mn].x)) return lr[r]= cut(lr[r], x0), void(); D q= n[lr[!r]].s + D(n[mn].x) * o[!r]; T v= o[!r] + n[lr[!r]].a; lr[!r]= cut(lr[!r], x0); if (!r) y-= q, rem+= v; else y+= q, rem-= v; if (lr[!r]) std::tie(lr[r], mn)= pop(lr[!r]), lr[!r]= 0; else mn= lr[r]= 0; o[r]= n[mn].d, o[!r]= 0; } inline void prop(T x) { if constexpr (persistent) mn= create(n[mn].d, n[mn].x); n[mn].x+= x; } public: // f(x) := 0 PiecewiseLinearConvex(): mn(0), lr{0, 0}, bf{0, 0}, o{0, 0}, rem(0), bx{0, 0}, y(0) {} // f(x) := sum max(0, a(x-x0)) PiecewiseLinearConvex(const std::vector> &ramps): PiecewiseLinearConvex() { int m= ramps.size(); if (!m) return; std::vector> w(m); int s= 0, t= 0; for (auto [d, x]: ramps) { if (lt(d, 0)) y-= D(d) * x, rem+= d, d= -d; if (!lt(0, d)) continue; w[s++]= {d, x}; } std::sort(w.begin(), w.begin() + s, [](auto a, auto b) { return a.second < b.second; }); for (int i= 0; i < s; ++i) { if (t && !lt(w[t - 1].second, w[i].second) && !lt(w[i].second, w[t - 1].second)) w[t - 1].first+= w[i].first; else w[t++]= w[i]; } mn= create(w[0].first, w[0].second), o[1]= n[mn].d, lr[1]= build(w.begin() + 1, w.begin() + t); } std::string info() { std::stringstream ss; if (ss << "\n rem:" << rem << ", y:" << y << ", mn:" << mn << ", lr:{" << lr[0] << ", " << lr[1] << "}\n bf[0]:" << bf[0] << ", bf[1]:" << bf[1] << ", bx[0]:" << bx[0] << ", bx[1]:" << bx[1] << "\n " << "o[0]:" << o[0] << ", o[1]:" << o[1] << "\n"; mn) { if (lr[0]) info(lr[0], 1, ss); ss << " ■ " << n[mn] << '\n'; if (lr[1]) info(lr[1], 1, ss); } return ss.str(); } template static inline void rebuild(Args &...plc) { static_assert(std::conjunction_v...>); constexpr size_t m= sizeof...(Args); std::array>, m> ls, rs; std::array, m> mns; int i= 0; (void)(int[]){(mns[i]= {n[plc.mn].d, n[plc.mn].x}, ls[i].resize(n[plc.lr[0]].sz), rs[i].resize(n[plc.lr[1]].sz), dump(ls[i].begin(), plc.lr[0]), dump(rs[i].begin(), plc.lr[1]), ++i)...}; ni= 1, i= 0; (void)(int[]){((plc.mn ? (plc.mn= create(mns[i].first, mns[i].second)) : 0), plc.lr[0]= build(ls[i].begin(), ls[i].end()), plc.lr[1]= build(rs[i].begin(), rs[i].end()), ++i)...}; } static inline void rebuild(std::vector &plcs) { size_t m= plcs.size(); std::vector>> ls(m), rs(m); std::vector> mns(m); for (int i= m; i--;) mns[i]= {n[plcs[i].mn].d, n[plcs[i].mn].x}, ls[i].resize(n[plcs[i].lr[0]].sz), rs[i].resize(n[plcs[i].lr[1]].sz), dump(ls[i].begin(), plcs[i].lr[0]), dump(rs[i].begin(), plcs[i].lr[1]); ni= 1; for (int i= m; i--;) (plcs[i].mn ? (plcs[i].mn= create(mns[i].first, mns[i].second)) : 0), plcs[i].lr[0]= build(ls[i].begin(), ls[i].end()), plcs[i].lr[1]= build(rs[i].begin(), rs[i].end()); } static void reset() { ni= 1; } static bool pool_empty() { if constexpr (persistent) return ni >= NODE_SIZE * 0.8; else return ni + 1000 >= NODE_SIZE; } // f(x) += c void add_const(D c) { y+= c; } // f(x) += ax, / void add_linear(T a) { rem+= a; } // f(x) += max(a(x-x0),b(x-x0)), (a < b) void add_max(T a, T b, T x0) { assert(lt(a, b)); if (bf[0] && x0 <= bx[0]) y-= D(b) * x0, rem+= b; else if (bf[1] && bx[1] <= x0) y-= D(a) * x0, rem+= a; else if (T c= b - a; mn) { if (lt(n[mn].x, x0)) lr[1]= insert(lr[1], x0, c), y-= D(a) * x0, rem+= a; else if (lt(x0, n[mn].x)) lr[0]= insert(lr[0], x0, c), y-= D(b) * x0, rem+= b; else { if constexpr (persistent) mn= create(n[mn].d, n[mn].x); n[mn].d+= c, o[1]+= c, y-= D(a) * x0, rem+= a; } } else mn= create(c, x0), y-= D(a) * x0, rem+= a, o[0]= 0, o[1]= c; } // f(x) += max(0, a(x-x0)) void add_ramp(T a, T x0) { if (lt(0, a)) add_max(0, a, x0); else if (lt(a, 0)) add_max(a, 0, x0); } // f(x) += a|x-x0|, \/ void add_abs(T a, T x0) { if (assert(!lt(a, 0)); lt(0, a)) add_max(-a, a, x0); } // right=false : f(x) += inf (x < x_0), right=true: f(x) += inf (x_0 < x) void add_inf(bool right= false, T x0= 0) { return right ? add_inf<1>(x0) : add_inf<0>(x0); } // f(x) -= max(a(x-x0),b(x-x0)), (a < b) void subtract_max(T a, T b, T x0) { assert(lt(a, b)); if (bf[0] && x0 <= bx[0]) y+= D(b) * x0, rem-= b; else if (bf[1] && bx[1] <= x0) y+= D(a) * x0, rem-= a; else { assert(mn); T c= b - a; if (lt(n[mn].x, x0)) lr[1]= erase(lr[1], x0, c), y+= D(a) * x0, rem-= a; else if (lt(x0, n[mn].x)) lr[0]= erase(lr[0], x0, c), y+= D(b) * x0, rem-= b; else { assert(!lt(n[mn].d, c)); if (lt(c, n[mn].d)) { if constexpr (persistent) mn= create(n[mn].d, n[mn].x); if (n[mn].d-= c, o[0]-= c, y+= D(b) * x0, rem-= b; lt(o[0], 0)) o[1]+= o[0], y+= D(o[0]) * x0, rem-= o[0], o[0]= 0; } else if (lr[1]) std::tie(lr[1], mn)= pop<0>(lr[1]); else if (lr[0]) std::tie(lr[0], mn)= pop<1>(lr[0]); else mn= 0, o[0]= o[1]= 0; } } } // f(x) -= max(0, a(x-x0)) void subtract_ramp(T a, T x0) { if (lt(0, a)) subtract_max(0, a, x0); else if (lt(a, 0)) subtract_max(a, 0, x0); } // f(x) -= a|x-x0|, \/ void subtract_abs(T a, T x0) { if (assert(!lt(a, 0)); lt(0, a)) subtract_max(-a, a, x0); } // f(x) <- f(x-x0) void shift(T x0) { if (bx[0]+= x0, bx[1]+= x0, y-= D(rem) * x0; mn) prop(x0), prop(lr[0], x0), prop(lr[1], x0); } // rev=false: f(x) <- min_{y<=x} f(y), rev=true : f(x) <- min_{x<=y} f(y) void chmin_cum(bool rev= false) { if (bf[0] && bf[1] && !lt(bx[0], bx[1])) y+= D(rem) * bx[0], rem= 0; else if (bool r= lt(rem, 0); r || lt(0, rem)) { T u= (r ? -rem : rem) - o[r] - n[lr[r]].a; if (!lt(u, 0)) { if (r ^ rev) { if (lt(0, u) && bf[r]) { D q= n[lr[r]].s + D(n[mn].x) * o[r] + D(u) * bx[r]; if (r ? y-= q : y+= q; mn) lr[!r]= join(lr[0], mn, lr[1]); o[!r]= u, rem= 0, mn= create(u, bx[r]), lr[!rev]= 0, o[!rev]= 0; } } else { assert(bf[r]); D q= n[lr[r]].s + D(n[mn].x) * o[r] + D(u) * bx[r]; (r ? y-= q : y+= q), rem= 0, mn= lr[0]= lr[1]= 0, o[0]= o[1]= 0; } bf[!rev]= false; return; } if ((r ^ rev)) r ? slope_eval_cum<0, 1>() : slope_eval_cum<0, 0>(); else r ? slope_eval_cum<1, 1>() : slope_eval_cum<1, 0>(); if constexpr (persistent) mn= create(o[rev], n[mn].x); else n[mn].d= o[rev]; } else if (mn) { if (!lt(0, o[rev])) { if (lr[rev]) std::tie(lr[rev], mn)= rev ? pop<0>(lr[rev]) : pop<1>(lr[rev]), o[rev]= n[mn].d; else mn= 0; } else { if constexpr (persistent) mn= create(o[rev], n[mn].x); else n[mn].d= o[rev]; } } bf[!rev]= false, lr[!rev]= 0, o[!rev]= 0; } // f(x) <- min_{lb<=y<=ub} f(x-y). (lb <= ub), \_/ -> \__/ void chmin_slide_win(T lb, T ub) { assert(lb <= ub); if (bf[0] && bf[1] && !lt(bx[0], bx[1])) y+= D(rem) * bx[0], rem= 0; else { if (bool r= lt(rem, 0); r || lt(0, rem)) { T u= (r ? -rem : rem) - o[r] - n[lr[r]].a; if (lt(0, u)) { T b[2]= {lb, ub}; if (bf[r]) { D q= n[lr[r]].s + D(n[mn].x) * o[r] + D(u) * bx[r]; if (r ? y-= q : y+= q; mn) lr[!r]= join(lr[0], mn, lr[1]), prop<0>(lr[!r], b[!r]); lr[r]= 0, rem= 0, o[!r]= u, o[r]= 0, mn= create(u, bx[r] + b[!r]); } else { y-= D(rem) * b[!r]; if (mn) prop(b[!r]), prop(lr[0], b[!r]), prop(lr[1], b[!r]); } bx[0]+= lb, bx[1]+= ub; return; } slope_eval(r); } if (mn) { if (!lt(0, o[0])) prop(ub); else if (!lt(0, o[1])) prop(lb); else lr[1]= join<0>(0, create(o[1], n[mn].x), lr[1]), prop(lb), n[mn].d= o[0], o[1]= 0; prop(lr[0], lb), prop(lr[1], ub); } } bx[0]+= lb, bx[1]+= ub; } std::optional operator()(T x) { if (bf[0] && x < bx[0]) return std::nullopt; if (bf[1] && bx[1] < x) return std::nullopt; return calc_y(x) + D(rem) * x + y; } std::optional min() { bool r= lt(rem, 0); if (!r && !lt(0, rem)) return y; T u= (r ? -rem : rem) - o[r] - n[lr[r]].a; if (lt(0, u)) { if (!bf[r]) return std::nullopt; D q= n[lr[r]].s + D(n[mn].x) * o[r] + D(u) * bx[r]; return r ? y - q : y + q; } return slope_eval(r), y; } std::array argmin() { if (bool r= lt(rem, 0); r || lt(0, rem)) { if (lt(o[r] + n[lr[r]].a, (r ? -rem : rem))) { if (!bf[r]) return {1, 0}; // no solution return {bx[r], bx[r]}; } slope_eval(r); } std::array ret= {bx[0], bx[1]}; int t= mn; if (!t) return ret; bool r= lt(0, o[0]); if (r && lt(0, o[1])) ret[0]= ret[1]= n[t].x; else if (ret[!r]= n[t].x, t= lr[r]; t) { for (; push(t), n[t].ch[!r];) t= n[t].ch[!r]; ret[r]= n[t].x; } else if (!bf[r]) return {1, 0}; // no solution return ret; } size_t size() { return n[lr[0]].sz + n[lr[1]].sz + !!mn; } PiecewiseLinearConvex &operator+=(const PiecewiseLinearConvex &g) { return *this= *this + g; } PiecewiseLinearConvex operator+(PiecewiseLinearConvex g) const { PiecewiseLinearConvex ret= *this; if (g.bf[0]) ret.add_inf(false, g.bx[0]); if (g.bf[1]) ret.add_inf(true, g.bx[1]); if (bf[0]) g.add_inf(false, bx[0]); if (bf[1]) g.add_inf(true, bx[1]); ret.y+= g.y, ret.rem+= g.rem; if (!g.mn) return ret; if (!ret.mn) return ret.mn= g.mn, ret.lr[0]= g.lr[0], ret.lr[1]= g.lr[1], ret.o[0]= g.o[0], ret.o[1]= g.o[1], ret; ret.y+= n[ret.lr[0]].s + D(n[ret.mn].x) * ret.o[0] + n[g.lr[0]].s + D(n[g.mn].x) * g.o[0], ret.rem-= ret.o[0] + n[ret.lr[0]].a + g.o[0] + n[g.lr[0]].a; int t= unite(join(ret.lr[0], ret.mn, ret.lr[1]), join(g.lr[0], g.mn, g.lr[1])); return std::tie(ret.lr[1], ret.mn)= pop<0>(t), ret.lr[0]= 0, ret.o[0]= 0, ret.o[1]= n[ret.mn].d, ret; } std::vector dump_xs() { std::vector xs; if (bf[0]) xs.push_back(bx[0]); dump_xs(lr[0], xs); if (mn) xs.push_back(n[mn].x); dump_xs(lr[1], xs); if (bf[1]) xs.push_back(bx[1]); return xs; } std::vector> dump_xys() { auto xs= dump_xs(); std::vector> xys(xs.size()); for (int i= xs.size(); i--;) xys[i]= {xs[i], operator()(xs[i])}; return xys; } std::vector dump_slopes() { std::vector as; if (mn) as.push_back(-o[0]), dump_slopes_l(lr[0], o[0], as), std::reverse(as.begin(), as.end()), as.push_back(o[1]), dump_slopes_r(lr[1], o[1], as); else as.push_back(0); for (auto &a: as) a+= rem; return as; } }; using namespace std; signed main() { cin.tie(0); ios::sync_with_stdio(0); long long N, K, L, U; cin >> N >> K >> L >> U; PiecewiseLinearConvex f; vector vl[K], vu[K]; for (int i= 0; i < N; ++i) { long long A; cin >> A; long long xl= L > A ? (A - L + 1) / K - 1 : (A - L) / K; f.add_ramp(1, xl); vl[((A - L) % K + K + 1) % K].push_back(xl); long long xu= U < A ? (A - U - 1) / K + 1 : (A - U) / K; f.add_ramp(-1, xu); vu[((A - U) % K + K) % K].push_back(xu); } long long ans= f.min().value(); for (int i= 1; i < K; ++i) { for (long long x: vl[i]) f.subtract_ramp(1, x), f.add_ramp(1, x - 1); for (long long x: vu[i]) f.subtract_ramp(-1, x), f.add_ramp(-1, x - 1); ans= min(ans, (long long)f.min().value()); } cout << ans << '\n'; return 0; }