OUPC2022 G「Restore Test Case」解説

# はじめに
[問題リンク](https://onlinejudge.u-aizu.ac.jp/beta/room.html#OUPC2022/problems/G)
「Restore Test Case」writer の vwxyz です。

ICPCのテストケース形式から思いついた問題です。
考察のほかに高度なアルゴリズム知識と実装が必要な問題です。
（コンテスト開催の数日前までAOJがメモリ16KBに対応していなかったりテストケースが弱かったりで想定よりも簡単な解法が通ってしまったという事情はありますが）テスターは誰も解けておらず、本番でもAC数が少ないと予想していました。実際にACできたのは1チームだけでした。

# 解説
まず、$0\ 0$の行は一番最後に並べることになるので、それ以外の行の並べ方について考えます。
各テストケースの$N\ M$に対応する行をNMタイプ、$L_i\ R_i$に対応する行をLRタイプと呼ぶことにします。

## タイプごとに右側の値は大きい行に並べてよい
テストケースが複数ある場合、それらはどんな順番で並べても構いません。
そこで、テストケースを$M$が大きい順に並べるという条件（・・・①）を付け加えることにします。
LRタイプの行が２つ以上存在する場合について考えます。
同じテストケースに属するとは限らないLRタイプの２つの行$L_1\ R_1$、$L_2\ R_2$（ただし、$R_1\leq R_2$（・・・②）とする）をとってきます。
ここで、条件を満たすように並べたときに$L_1\ R_1$は$L_2\ R_2$よりも先に並べられる（・・・③）と仮定します。
$L_1\ R_1$、$L_2\ R_2$のテストケースにおける$M$の値をそれぞれ$M_1$、$M_2$とおくと、テストケースの制約により
$R_1\leq M_1$
$R_2\leq M_2$（・・・④）
であり、①と③より
$M_2\leq M_1$（・・・⑤）
となります。すると、②と④より
$R_1\leq R_2\leq M_2$
④と⑤より
$R_2\leq M_2\leq M_1$
なので、$L_1\ R_1$の行と$L_2\ R_2$の行は入れ替えても条件が満たされることがわかります。
よって、LRタイプは$R$が大きい順に並べるという条件を付け足せることがわかります。
以上により、NMタイプ、LRタイプだけで見れば各行は右側の値によって降順に並んでいるとしてよい（・・・⑥）ことがわかります。

## 右側の値が大きい行から見て各タイプに分類していく
まず、各行の右側の値が相異なる場合を考えます。
そこで、$0\ 0$以外の$T-1$行を右側の値で降順にソートし、順番にNMタイプとLRタイプのどちらかに分類して上から並べていくことにします。
ただし、NMタイプに分類した場合はそれを下に置いてその下にLRタイプを置くための空行を$N$行確保し、LRタイプに分類した場合は一番上の空行をそれで書き換えることにします。
⑥から、条件を満たす並べ替えが存在するなら適切な分類を行うことでそれを構築できることがわかるので、適切に分類するための条件を考えます。
まず、$A>B$となるような行はNMタイプに分類されることがわかります。他に
LRタイプを置くために空けられた行が無い状態でLRタイプには分類できない（・・・⑦）
すべての行を分類し終わったあとにLRタイプを置くために空けた行が残らない（・・・⑧）
の２つが必要であり、逆に、これらが両方満たされていればよいこともわかります。
ここで、NMタイプが置かれている行数、LRタイプが置かれている行数、LRタイプを置くために確保されている空行の行数の合計を「確保行数」と呼ぶことにします。
すると、NMタイプ、LRタイプに分類することはそれぞれ確保行数に$N+1$、$0$を足すことに相当（・・・⑨）します。
また、今まで見た行数を「暫定行数」とすると
⑦→$暫定行数\leq 確保行数$（・・・⑩）
⑧→確保行数は最終的に$T-1$になる（・・・⑪）
と言い換えられます。

## 条件付きの部分和問題に帰着する
次に、各行の右側の値が相異なるとは限らない場合を考えます。
右側が同じ値のいくつかの行のうちLRタイプを先に見ると$確保行数<暫定行数$になることがあるが、それらの順番を入れ替えることで条件を満たせる場合があります。
そこで、今まで見た行のうち右側の値が今見ている行のそれよりも真に大きいものの数を「暫定行数」とすると、⑩と⑪を満たすように分類すれば条件を満たすように並べられることがわかります。
$0\ 0$以外のすべての行を右側の値で降順にソートしたときの$i$番目（$0\leq i<T$）の行を$A_i\ B_i$、その行を見るときの暫定行数を$P_i$とすると、$P_i$（$0\leq i<T$）は$O(T)$で計算できるので、
$P_i\leq \displaystyle\sum_{i行目を見るまででNMタイプに選んだj}A_j+1$（$i=0,1,\cdots,T-2$）　（・・・⑫）
$T-1=\displaystyle\sum_{NMタイプに選んだj}A_j+1$
となるようにNMタイプを選ぶ問題へと帰着できます。

## 部分和問題をdpで解く（復元する）
$dp[i][j]=i行を見て和がjになっていることがあるか$（$0\leq i\leq T-1,0\leq i\leq T-1$）
とおくと、⑫は
$dp[i][j]=0$（$0\leq i<T-1$,$0\leq j<P_i$）
となり、これを計算していくことで部分和問題を解くことができます。

## dpを高速化する
この$dp$は普通に計算すると$O(T^2)$であり、$1\leq T\leq 200000$なので間に合いません。
そこで、ビットセットを用いると
```
dp[i+1]=dp[i]|dp[i]<<A_i
```
と更新することができ、$O(\frac{T^2}{64})$での計算が可能となります。

## dpのメモリを節約する
この問題では復元を要求されており、このままでは$dp$の計算に必要な空間計算量が$O(T^2)$となってしまいます。
そこで、メモリを省略する方法を考える必要があります。
まず、平方分割のように、$dp[i][j]$をin-placeに計算しながら$i$が$450$（$\fallingdotseq\sqrt{200000}$）の倍数になるようなものだけ保存しておき、その間（$i$が$450$の倍数ではない部分）を復元してメモリから削除するのを繰り返すことで空間計算量を$O(T\sqrt{T})$にすることができます。
しかし、これでは$16KB$には抑えられない（可能性が高い）です。
そこで、Hirschberg's Algorithm（分割統治法のように真ん中の値だけ復元してから左右に分けて解く）を用いると、時間計算量は$O(\frac{T^2}{64})$のままで空間計算量を$O(T)$にすることができます。

## まとめ
以上により、この問題を$O(T\log T+\frac{T^2}{64})$で解けることがわかりました。

# 実装
```cpp
#pragma GCC target("avx2")
#pragma GCC optimize("O3")
#pragma GCC optimize("unroll-loops")
#include <iostream>
#include <bits/stdc++.h>
using namespace std;

struct Bitset {
    private:
    constexpr static unsigned long long mask[] = {
        0x0000000000000000ull, 0x0000000000000001ull, 0x0000000000000003ull, 0x0000000000000007ull,
        0x000000000000000Full, 0x000000000000001Full, 0x000000000000003Full, 0x000000000000007Full,
        0x00000000000000FFull, 0x00000000000001FFull, 0x00000000000003FFull, 0x00000000000007FFull,
        0x0000000000000FFFull, 0x0000000000001FFFull, 0x0000000000003FFFull, 0x0000000000007FFFull,
        0x000000000000FFFFull, 0x000000000001FFFFull, 0x000000000003FFFFull, 0x000000000007FFFFull,
        0x00000000000FFFFFull, 0x00000000001FFFFFull, 0x00000000003FFFFFull, 0x00000000007FFFFFull,
        0x0000000000FFFFFFull, 0x0000000001FFFFFFull, 0x0000000003FFFFFFull, 0x0000000007FFFFFFull,
        0x000000000FFFFFFFull, 0x000000001FFFFFFFull, 0x000000003FFFFFFFull, 0x000000007FFFFFFFull,
        0x00000000FFFFFFFFull, 0x00000001FFFFFFFFull, 0x00000003FFFFFFFFull, 0x00000007FFFFFFFFull,
        0x0000000FFFFFFFFFull, 0x0000001FFFFFFFFFull, 0x0000003FFFFFFFFFull, 0x0000007FFFFFFFFFull,
        0x000000FFFFFFFFFFull, 0x000001FFFFFFFFFFull, 0x000003FFFFFFFFFFull, 0x000007FFFFFFFFFFull,
        0x00000FFFFFFFFFFFull, 0x00001FFFFFFFFFFFull, 0x00003FFFFFFFFFFFull, 0x00007FFFFFFFFFFFull,
        0x0000FFFFFFFFFFFFull, 0x0001FFFFFFFFFFFFull, 0x0003FFFFFFFFFFFFull, 0x0007FFFFFFFFFFFFull,
        0x000FFFFFFFFFFFFFull, 0x001FFFFFFFFFFFFFull, 0x003FFFFFFFFFFFFFull, 0x007FFFFFFFFFFFFFull,
        0x00FFFFFFFFFFFFFFull, 0x01FFFFFFFFFFFFFFull, 0x03FFFFFFFFFFFFFFull, 0x07FFFFFFFFFFFFFFull,
        0x0FFFFFFFFFFFFFFFull, 0x1FFFFFFFFFFFFFFFull, 0x3FFFFFFFFFFFFFFFull, 0x7FFFFFFFFFFFFFFFull,
        0xFFFFFFFFFFFFFFFFull
    };
    
    void correct() {
        if(n % 64) v[bnum-1] &= mask[n % 64];
    }
    public:
    vector<unsigned long long> v;
    int n, bnum;
    Bitset(int n_ = 0) : n(n_) {
        bnum = (n+63) / 64;
        v.resize(bnum);
    }
    int operator[](int i) const {
        return (v[i/64] >> (i%64)) & 1;
    }
    int count() {
        int c = 0;
        for (int i = 0; i < v.size(); i++) {
            c += __builtin_popcountll(v[i]);
        }
        return c;
    }
    int count_range(int l, int r) {
        int c = 0;
        int l2 = (l+63) / 64;
        int r2 = r / 64;
        for(int i = l2; i < r2; i++) {
            c += __builtin_popcountll(v[i]);
        }
        if(l < l2 * 64) {
            for(int i = l % 64; i < 64; i++) c += (v[l2-1] >> i) & 1;
        }
        if(r2 * 64 < r) {
            for(int i = 0; i < r % 64; i++) c += (v[r2] >> i) & 1;
        }
        return c;
    }
    bool all() {
        return count() == n;
    }
    bool any() {
        return count() > 0;
    }
    bool none() {
        return count() == 0;
    }
    void set(int i) {
        v[i / 64] |= 1ull << (i % 64);
    }
    void reset(int i) {
        v[i / 64] &= ~(1ull << (i % 64));
    }
    void flip(int i) {
        v[i / 64] ^= 1ull << (i % 64);
    }
    void resize(int n_) {
        n = n_;
        v.resize((n+63) / 64);
        correct();
    }
    void all_set() {
        fill(v.begin(), v.end(),  ~0ULL);
        correct();
    }
    void all_reset() {
        fill(v.begin(), v.end(), 0);
    }
    void all_flip() {
        for (int i = 0; i < v.size(); i++) {
            v[i] = ~v[i];
        }
        correct();
    }
    Bitset& operator&=(const Bitset& b) {
        for(int i = 0; i < min(bnum, b.bnum); i++) {
            v[i] &= b.v[i];
        }
        return *this;
    }
    Bitset operator&(const Bitset& b) const {
        return Bitset(*this) &= b;
    }
    Bitset& operator|=(const Bitset& b) {
        for(int i = 0; i < min(bnum, b.bnum); i++) {
            v[i] |= b.v[i];
        }
        correct();
        return *this;
    }
    Bitset operator|(const Bitset& b) const {
        return Bitset(*this) |= b;
    }
    Bitset& operator^=(const Bitset& b) {
        for(int i = 0; i < min(bnum, b.bnum); i++) {
            v[i] ^= b.v[i];
        }
        correct();
        return *this;
    }
    Bitset operator^(const Bitset& b) const {
        return Bitset(*this) ^= b;
    }
    Bitset operator~() const {
        Bitset b(*this);
        b.all_flip();
        return b;
    }
    bool operator==(const Bitset& b) const {
        return v == b.v;
    }
    bool operator!=(const Bitset& b) const {
        return v != b.v;
    }
    Bitset& operator<<=(int sz) {
        for(int i = bnum-1; i >= 0; i--) {
            if(i-sz/64 < 0) v[i] = 0;
            else if(i-sz/64-1 < 0 || sz%64 == 0) v[i] = v[i-sz/64] << (sz%64);
            else v[i] = (v[i-sz/64] << (sz%64)) | (v[i-sz/64-1] >> (64-sz%64));
        }
        correct();
        return *this;
    }
    Bitset operator<<(int sz) const {
        return Bitset(*this) <<= sz;
    }
    Bitset& operator>>=(int sz) {
        for(int i = 0; i < bnum; i++) {
            if(i+sz/64 >= bnum) v[i] = 0;
            else if(i+sz/64+1 >= bnum || sz%64 == 0) v[i] = v[i+sz/64] >> (sz%64);
            else v[i] = (v[i+sz/64] >> (sz%64)) | (v[i+sz/64+1] << (64-sz%64));
        }
        return *this;
    }
    Bitset operator>>(int sz) const {
        return Bitset(*this) >>= sz;
    }
    int find_first(){
        for(int i=0;i<v.size();i++){
            if(v[i]){
                return i*64+__builtin_ffsll(v[i])-1;
            }
        }
        assert(0);
        return -1;
    }
};

int main(){
    
	int T;cin>>T;
	vector<vector<int>> AB(200001,vector<int>(0));
	for(int t=0;t<T;t++){
		int a,b;cin>>a>>b;
		AB[b].push_back(a);
	}
	assert(AB[0].size()==1);
	assert(AB[0][0]==0);
	T-=1;
	vector<int> lower_limit(0),A(0),B(0);
	int cnt=0;
	for(int b=200000;b>0;b--){
		for(int a:AB[b]){
			A.push_back(a);
			B.push_back(b);
			lower_limit.push_back(cnt);
		}
		cnt+=AB[b].size();
	}
	AB.clear();
	lower_limit.push_back(T);
	vector<int> partial_sum(T+1);
	partial_sum[0]=0;
	partial_sum[T]=T;
	vector<pair<int,int>> stack={{0,T}};
	while(stack.size()){
		pair<int,int> p=stack[stack.size()-1];
		stack.pop_back();
		int l=p.first;
		int r=p.second;
		int mid=(l+r)/2;
		Bitset dpL(partial_sum[r]-partial_sum[l]+1),dpR(partial_sum[r]-partial_sum[l]+1);
		dpL.set(0);
		int prev=partial_sum[l];
		for(int i=l+1;i<mid+1;i++){
			if(A[i-1]<=B[i-1]){
				dpL|=dpL<<A[i-1]+1;
			}
			else{
				dpL<<=A[i-1]+1;
			}
			for(int x=lower_limit[i-1];x<lower_limit[i];x++){
				if(x>=partial_sum[l]){
					dpL.reset(x-partial_sum[l]);
				}
			}
		}
		dpR.set(partial_sum[r]-max(lower_limit[r],partial_sum[l]));
		for(int i=r-1;i>mid-1;i--){
			dpR<<=max(lower_limit[i+1],partial_sum[l])-max(lower_limit[i],partial_sum[l]);
			if(A[i]<=B[i]){
				dpR|=dpR>>A[i]+1;
			}
			else{
				dpR>>=A[i]+1;
			}
		}
		dpR<<=max(lower_limit[mid]-partial_sum[l],0);
		partial_sum[mid]=partial_sum[l]+(dpL&dpR).find_first();
		if(mid-l>=2){
			stack.push_back({l,mid});
		}
		if(r-mid>=2){
			stack.push_back({mid,r});
		}
	}
	lower_limit.clear();
	int lr=0;
	for(int nm=0;nm<T;nm++){
		if(partial_sum[nm]!=partial_sum[nm+1]){
			cout<<A[nm]<<" "<<B[nm]<<endl;
			for(int cnt=0;cnt<A[nm];cnt++){
				while(partial_sum[lr]!=partial_sum[lr+1]){
					lr+=1;
				}
				cout<<A[lr]<<" "<<B[lr]<<endl;
				lr+=1;
			}
		}
	}
	cout<<0<<" "<<0<<endl;
}
```