Submission #23091215

---

Source Code Expand

#define _USE_MATH_DEFINES
#define _SILENCE_CXX17_C_HEADER_DEPRECATION_WARNING
#define _SILENCE_ALL_CXX17_DEPRECATION_WARNINGS
#define _CRT_SECURE_NO_WARNINGS

#pragma GCC target ("popcnt,bmi2,fma,fma4,avx2")
// #pragma GCC target ("skylake-avx512")
#pragma clang attribute push (__attribute__((target("popcnt,bmi2,fma,fma4,avx2"))), apply_to=function)
// #pragma clang attribute push (__attribute__((target("popcnt,bmi2,fma,fma4,avx2,avx512f,avx512dq,avx512cd,avx512bw,avx512vl"))), apply_to=function)


#if __has_include(<bits/stdc++.h>)
#include <bits/stdc++.h>
#endif // bits/stdc++

#include <any>
#include <array>
#include <cassert>
// #include <charconv>
#include <cinttypes>
#include <codecvt>
#include <condition_variable>
#include <ctgmath>
#include <forward_list>
#include <fstream>
#include <future>
#include <iostream>
#include <list>
#include <map>
#include <queue>
#include <random>
#include <regex>
#include <scoped_allocator>
#include <set>
#include <shared_mutex>
#include <typeindex>
#include <unordered_map>
#include <unordered_set>
#include <valarray>
#include <variant>
#include <nmmintrin.h>
#include <x86intrin.h>

#if __has_include(<atcoder/all>)

#include <atcoder/all>
//using namespace atcoder;
using atcoder::fenwick_tree;
using atcoder::modint;
using atcoder::pow_mod;
using atcoder::inv_mod;
using atcoder::floor_sum;
using UF = atcoder::dsu;

//using mint = atcoder::modint1000000007;
//using mint = atcoder::modint998244353;


#endif // atcoder/all

using namespace std;


#define printe(...) fprintf(stderr, __VA_ARGS__);
#define rep(i,a) for(int i=0;i<(a);i++)
#define reP(i,a) for(int i=0;i<=(a);i++)
#define Rep(i,a) for(int i=(a)-1;i>=0;i--)
#define ReP(i,a) for(int i=(a);i>=0;i--)
#define repft(i,a,b) for(int i=(a);i<(b);i++)
#define repfT(i,a,b) for(int i=(a);i<=(b);i++)
#define Repft(i,a,b) for(int i=(a)-1;i>=(b);i--)
#define RepfT(i,a,b) for(int i=(a);i>=(b);i--)
#define all(c)  begin(c),end(c)
#define FILL(a,v) fill(begin(a),end(a), v)
#define FILL0(a) memset(a,0,sizeof(a))
#define FILL1(a) memset(a,-1,sizeof(a))

#define fs  first
#define sc  second

using schar = signed char;
using uchar = unsigned char;
using ushort = unsigned short;
using uint = unsigned int;
using ull = unsigned long long;
using ll  = long long;
using ull = unsigned long long;
using Pi  = pair<int,int>;
using Pll = pair<ll,ll>;

const int INF  = 1'010'000'000; // 0x3C33'6080
const ll INFLL = 2'002'002'002'002'002'002; //0x1BC8 8A36 B8F3 4052

template <class A, class B> ostream& operator<<(ostream& st, const pair<A, B>& P) { return st << "(" << P.first << "," << P.second << ")"; };
template <class A, class B> pair<A, B> operator+(const pair<A, B>& P, const pair<A, B>& Q) { return pair<A, B>(P.first + Q.first, P.second + Q.second); };
template <class A, class B> pair<A, B> operator-(const pair<A, B>& P, const pair<A, B>& Q) { return pair<A, B>(P.first - Q.first, P.second - Q.second); };
template <typename T> T div_ceil(const T a, const T b) { return (a + b - 1) / b;}

#if _MSC_VER >= 1920 || __INTEL_COMPILER
int bitCount64(unsigned long long bits) {return (int)_mm_popcnt_u64(bits);}
int bitCount32(unsigned int bits)       {return (int)_mm_popcnt_u32(bits);}
#else
int bitCount64(unsigned long long bits) {return (int)__builtin_popcountll(bits);}
int bitCount32(unsigned int bits)       {return (int)__builtin_popcount(bits);}
#endif 


template <typename IdxType = uint32_t, int WORD_SIZE = 64>
struct bit_vector {
    static constexpr uint32_t wsize = WORD_SIZE;
    static_assert(WORD_SIZE == 32 || WORD_SIZE == 64);
    using WordT = typename conditional<WORD_SIZE == 64, uint64_t, uint32_t>::type;

#if _MSC_VER >= 1920 || __INTEL_COMPILER
    inline static int bitCount32(uint32_t bits) {
        return (int)_mm_popcnt_u32(bits);
    }
    inline static int bitCount64(uint64_t bits) {
        return (int)_mm_popcnt_u64(bits);
    }
#else
    inline static int bitCount32(uint32_t bits) {
        return __builtin_popcount(bits);
    }
    inline static int bitCount64(uint64_t bits) {
        return __builtin_popcountll(bits);
    }
#endif

    inline static uint32_t rank(const uint32_t x, const uint32_t i) {
        return bitCount32((uint32_t)_bzhi_u32(x, i));
        //return bitCount32(x & ((1u << i) - 1));
        //return 0;
    }
    inline static uint32_t rank(const uint64_t x, const uint32_t i) {
        //return 0;
        return bitCount64((uint64_t)_bzhi_u64(x, i));
        //return bitCount64(x & ((1ull << i) - 1));
    }
#pragma pack(4)
    struct block_t {
        WordT    nbit;
        IdxType  nsum;
    };
#pragma pack()
    IdxType zeros, n;
    vector<block_t> block;
    bit_vector(const IdxType _n = 0) : n(_n), block(n / wsize + 1) {}
    void reset(const IdxType _n = 0) { n = _n; block.clear(); block.resize(n / wsize + 1); }
    inline bool is_zero(const IdxType i) const {
        return (block[i / wsize].nbit >> (i % wsize)) & 1;
    }
    inline void nset(const IdxType i) {
        block[i / wsize].nbit |= (WordT)1 << (i % wsize);
    }
    inline void nset_block_bit(const IdxType i, const WordT val) {
        block[i].nbit = val;
    }
    inline WordT get_block_nbit(const IdxType i) const {
        return block[i].nbit;
    }
    void build() {
        for (IdxType j = 0; j < n / wsize; ++j) {
            block[j + 1].nsum = block[j].nsum + bitCount64(block[j].nbit);
        }
        zeros = rank0(n);
    }
    inline IdxType rank0(const IdxType i) const {
        auto&& e = block[i / wsize];
        return e.nsum + rank(e.nbit, i % wsize);
    }
    inline IdxType rank1(const IdxType i) const {
        return i - rank0(i);
    }
};

template <class T = uint64_t, typename IdxType = uint32_t, int WORD_SIZE = 64>
struct WaveletMatrix {
    static constexpr uint32_t wsize = WORD_SIZE;
    using WordT = typename bit_vector<IdxType, WORD_SIZE>::WordT;

    IdxType size, bit_len;
    vector<T> a;
    vector<bit_vector<IdxType, WORD_SIZE>> bv;
    WaveletMatrix(IdxType _n = 0) {
        reset(_n);
    }
    void reset(IdxType _n = 0) {
        size = (_n + wsize - 1) & (~(wsize - 1));
        a.resize(size);
        for (IdxType i = _n; i < size; ++i) {
            a[i] = 0;
        }
    }
    static constexpr int get_bit_len(uint64_t val) {
        val |= val >> 1;
        val |= val >> 2;
        val |= val >> 4;
        val |= val >> 8;
        val |= val >> 16;
        val |= val >> 32;
        return bit_vector<>::bitCount64(val);
    }
    inline void set_prebuild(const IdxType i, const T& val) {
        a[i] = val;
    }
    void build(const IdxType data_bit) {
        bit_len = data_bit;
        bv.assign(bit_len, size);
        vector<T> nxt(size);

        T mask = ((T)1 << bit_len) - 1;
        const IdxType n_dw = size / wsize;
        for (IdxType h = bit_len; h--; ) {
            mask >>= 1;

            T* it = &a[0];
            for (IdxType i = 0; i < n_dw; ++i, it += wsize) {
                WordT bits = 0;
                for (int j = 0; j < wsize; ++j) {
                    if (it[j] <= mask) { // 0
                        bits |= ((WordT)1 << j);
                    }
                }
                bv[h].nset_block_bit(i, bits);
            }
            bv[h].build();
            auto its = array{
                begin(nxt),
                begin(nxt) + bv[h].zeros
            };
            for (IdxType i = 0; i < size; ++i) {
                const T v = a[i];
                *(its[v > mask]++) = v & mask;
            }
            swap(a, nxt);
        }
        a.clear();
        a.shrink_to_fit();
    }
    T kth(IdxType l, IdxType r, IdxType k) const {
        T res = 0;
        for (int h = bit_len; h--; ) {
            const IdxType l0 = bv[h].rank0(l);
            const IdxType r0 = bv[h].rank0(r);
            const IdxType d0 = r0 - l0;
            if (k < d0) {
                l = l0;
                r = r0;
            } else {
                k -= d0;
                res |= (T)1 << h;
                l += bv[h].zeros - l0;
                r += bv[h].zeros - r0;
            }
        }
        return res;
    }

    // [l, r) に含まれる x 未満の値の個数 O(log(max))
    IdxType less_freq(IdxType l, IdxType r, const T x) const {
        // if (x <= 0) return 0;
        // if (l >= r) return 0;
        // if (x >= 1LL << (bit_len)) return r - l;
        IdxType res = 0;

        for (IdxType d = bit_len; d--; ) {
            if (((x >> d) & 1) == 0) {
                l = bv[d].rank0(l);
                r = bv[d].rank0(r);
            } else {
                res += bv[d].rank0(r) - bv[d].rank0(l);
                l += bv[d].zeros - bv[d].rank0(l);
                r += bv[d].zeros - bv[d].rank0(r);
            }
        }
        return res;
    }

    // [l, r) に含まれる [a, b) の範囲の値の個数 O(log(max))
    int range_freq(const IdxType l, const IdxType r, const T a, const T b) const {
        //return less_freq(l, r, b) - less_freq(l, r, a);
        
        if (l >= r) return 0;
        IdxType la = l, lb = l;
        IdxType ra = r, rb = r;
        IdxType res = 0;
        for (IdxType d = bit_len; d--; ) {
            const IdxType la0 = bv[d].rank0(la);
            const IdxType lb0 = bv[d].rank0(lb);
            const IdxType ra0 = bv[d].rank0(ra);
            const IdxType rb0 = bv[d].rank0(rb);
            if (((a >> d) & 1) == 0) {
                la = la0;
                ra = ra0;
            } else {
                res -= ra0 - la0;
                la += bv[d].zeros - la0;
                ra += bv[d].zeros - ra0;
            }
            if (((b >> d) & 1) == 0) {
                lb = lb0;
                rb = rb0;
            } else {
                res += rb0 - lb0;
                lb += bv[d].zeros - lb0;
                rb += bv[d].zeros - rb0;
            }
        }
        return res;
    }
};


uint src[801][801];

template<typename ValT = uint32_t>
struct WM_2DMedian {
    WM_2DMedian() {}

    static_assert(is_same<ValT, uint64_t>() || is_same<ValT, uint32_t>() || is_same<ValT, uint16_t>() || is_same<ValT, uint8_t>());
    static constexpr int val_bit_len = 20;
    using XIdxT = uint16_t;
    using YIdxT = uint16_t;
    using XYIdxT = int32_t;
    bit_vector<XYIdxT> BV[val_bit_len];
    WaveletMatrix<XIdxT> WM[val_bit_len];
private:
    int H, W;
public:
    WM_2DMedian* build(int NN) {
        // H = src.size();
        // W = src[0].size();
        H = W = NN;

        const int HW = H * W;
        const int w_bit_len = WaveletMatrix<>::get_bit_len(W+1); // w=7 [0,6] bl=3; w=8 [0,7] bl=4
        auto res = vector(H, vector<ValT>(W));

        assert(W < 65535); // That is, less than 65534.

        vector<pair<ValT, XIdxT>> buf[2], nxt[2]; // val, Xidx
        buf[0].reserve(HW);
        buf[1].reserve(HW);
        nxt[0].reserve(HW);
        nxt[1].reserve(HW);

        const XIdxT inf = ((XIdxT)1 << w_bit_len) - 1;
        assert(W < inf);

        rep(y, H) rep(x, W) {
            buf[0].emplace_back(src[y][x], x);
        }

        ValT mask = ((ValT)1 << val_bit_len) - 1;
        for (int d = val_bit_len; d--; ) {
            nxt[0].clear();
            nxt[1].clear();
            mask >>= 1;
            WM[d].reset(HW);// = WaveletMatrix<XIdxT>(HW);
            BV[d].reset(HW);// = bit_vector<XYIdxT>(HW);
            int i = -1;
            rep(b, 2) {
                for (const auto p : buf[b]) {
                    ++i;
                    if (p.first <= mask) { // 0
                        WM[d].set_prebuild(i, p.second);
                        nxt[0].emplace_back(p.first, p.second);
                        BV[d].nset(i);
                    } else {               // 1
                        WM[d].set_prebuild(i, inf);
                        nxt[1].emplace_back(p.first & mask, p.second);
                    }
                }
            }
            WM[d].build(w_bit_len);
            BV[d].build();
            buf[0].swap(nxt[0]);
            buf[1].swap(nxt[1]);
        }
        buf[0].clear(); buf[0].shrink_to_fit();
        buf[1].clear(); buf[1].shrink_to_fit();
        nxt[0].clear(); nxt[0].shrink_to_fit();
        nxt[1].clear(); nxt[1].shrink_to_fit();
        return this;
    }
    inline ValT kth(const XIdxT xa, const XIdxT xb, int ya, int yb, XYIdxT k) const { // [xa, xb), [ya, yb)
        ValT res = 0;
        ya *= W;
        yb *= W;
        for (int h = val_bit_len; h--; ) {
            const int d = WM[h].range_freq(ya, yb, xa, xb);
            const XYIdxT top0 = BV[h].rank0(ya);
            const XYIdxT bot0 = BV[h].rank0(yb);
            if (k < d) {
                ya = top0;
                yb = bot0;
            } else {
                k -= d;
                res |= (ValT)1 << h;
                ya += BV[h].zeros - top0;
                yb += BV[h].zeros - bot0;
            }
        }
        return res;

    }
    int run(const int tyo) const {
        int res = INF;
        rep(y, H-tyo+1) {
            const YIdxT ya = y;
            const YIdxT yb = y+tyo;
            const int ky = yb - ya;
            rep(x, W-tyo+1) {
                const XIdxT xa = x;
                const XIdxT xb = x + tyo;
                const int kx = xb - xa;
                const int k = (ky * kx - 1) / 2;
                auto v = kth(xa, xb, ya, yb, k);
                res = min<int>(res, v);
            }
        }
        return res;
    }
};


uint32_t a[801][801];

void solve(long long n, long long k){
    WM_2DMedian<> WM;
    

    vector<pair<uint32_t, Pi>> v;
    v.reserve(n*n);
    rep(i, n){
        rep(j, n){
            v.emplace_back(a[i][j], Pi(i,j));
        }
    }
    sort(all(v));
    rep(i, n*n) {
        int y = v[i].second.first;
        int x = v[i].second.second;
        src[y][x] = i;
    }

    WM.build(n);

    /*
    int m = (k * k - 1) / 2;
    int res = INF;
    reP(i, n - k) {
        reP(j, n - k) {
            int v = WM.kth(i, i + k + 1, j, j + k + 1, m);
            cerr << v << " ";
            res = min(res, v);
        }
        cerr << endl;
    }
    cout << res << endl;
    //*/
    auto m = WM.run(k);

    /*
    rep(i, n){
        rep(j,n){
            cerr << b[i][j] << " ";
        }
        cerr << endl;
    }
    cerr << "m" << m << endl;
    //*/

    cout << a[v[m].second.first][v[m].second.second] << endl;
}


// ~/.local/lib/python3.6/site-packages/atcodertools/tools/templates/default_template.cpp
// Generated by 2.3.0 https://github.com/kyuridenamida/atcoder-tools  (tips: You use the default template now. You can remove this line by using your custom template)
int main(){
    long long N;
    scanf("%lld", &N);
    long long K;
    scanf("%lld", &K);
    for(int i = 0 ; i < N ; i++){
        for(int j = 0 ; j < N ; j++){
            scanf("%d", &a[i][j]);
        }
    }
    solve(N, K);
    return 0;
}
#pragma clang attribute pop

Submission Info

Compile Error

./Main.cpp:6:13: warning: unknown pragma ignored [-Wunknown-pragmas]
#pragma GCC target ("popcnt,bmi2,fma,fma4,avx2")
            ^
In file included from ./Main.cpp:46:
In file included from /opt/ac-library/atcoder/all:6:
In file included from /opt/ac-library/atcoder/maxflow:1:
/opt/ac-library/atcoder/maxflow.hpp:38:13: warning: unused variable 'm' [-Wunused-variable]
        int m = int(pos.size());
            ^
/opt/ac-library/atcoder/maxflow.hpp:53:13: warning: unused variable 'm' [-Wunused-variable]
        int m = int(pos.size());
            ^
In file included from ./Main.cpp:46:
In file included from /opt/ac-library/atcoder/all:7:
In file included from /opt/ac-library/atcoder/mincostflow:1:
/opt/ac-library/atcoder/mincostflow.hpp:39:13: warning: unused variable 'm' [-Wunused-variable]
        int m = int(pos.size());
            ^
In file included from ./Main.cpp:46:
In file included from /opt/ac-library/atcoder/all:9:
In file included from /opt/ac-library/atcoder/scc:1:
/opt/ac-library/atcoder/scc.hpp:17:13: warning: unused variable 'n' [-Wunused-variable]
        int n = internal.num_vertices();
            ^
In file included from ./Main.cpp:46:
In file included from /opt/ac-library/atcoder/all:11:
In file included from /opt/ac-library/atcoder/string:1:
/opt/ac-library/atcoder/string.hpp:179:14: warning: unused variable 'd' [-Wunused-variable]
    for (int d : s) {
             ^
./Main.cpp:92:10: warning: unused variable 'INFLL' [-Wunused-const-variable]
const ll INFLL = 2'002'002'002'002'002'002; //0x1BC8 8A36 B8F3 4052
         ^
7 warnings generated.

Judge Result