// THE THREE-HALVES ECHO — a farther, parity-aware probe of THE SIBLING HORIZON.
//
// Computes Conway's recurrence
//   a(1)=a(2)=1; a(n)=a(a(n-1))+a(n-a(n-1))
// to 2^K (K<=29 in this instrument), verifies 1<=a(n)<=n and the exact
// clock a(2^k)=2^(k-1), then reproduces the fine- and envelope-grain sibling
// correlations of THE SIBLING HORIZON.
//
// The new object is the parity pair.  For odd k and the following even k+1,
// define
//   center_k = (r_k + r_{k+1})/2,
//   amp_k    = (r_k - r_{k+1})/2  (>0 in the observed tail),
//   m_k      = k + 1/2.
// We report amp*m, amp*m^(3/2), and amp*m^2, and a log-log least-squares
// exponent over the last five complete pairs.  This is descriptive finite
// weather, not an asymptotic proof.
//
// Law-9 gates:
//   * range and clock are asserted;
//   * the five published fine-grain values k=22..26 from THE SIBLING HORIZON
//     must reproduce within 5e-9;
//   * envelope counterparts k=22..26 must reproduce within 5e-9;
//   * the parity sign must hold on every complete pair k=17..K-2 used here.
//
// Evidence class: COMPUTED.
// 2026-07-27 — walker at the open door.

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <numeric>
#include <string>
#include <vector>

using u32 = uint32_t;
using i64 = int64_t;
using ld = long double;

static ld fine_r(const std::vector<u32>& A, int k) {
    const size_t half = size_t(1) << (k - 1);
    const size_t x0 = size_t(1) << k;
    const size_t y0 = half;
    ld sx=0, sy=0, sxx=0, syy=0, sxy=0;
    for (size_t j=0;j<half;++j) {
        const ld x = ld(2)*A[x0+j] - ld(x0+j);
        const ld y = ld(2)*A[y0+j] - ld(y0+j);
        sx += x; sy += y; sxx += x*x; syy += y*y; sxy += x*y;
    }
    const ld n = ld(half);
    const ld vx = sxx - sx*sx/n;
    const ld vy = syy - sy*sy/n;
    const ld cov = sxy - sx*sy/n;
    return cov / std::sqrt(vx*vy);
}

static std::vector<ld> env_profile(const std::vector<u32>& A, size_t lo, size_t hi, int bins=64) {
    const size_t len = hi-lo;
    const size_t width = len / bins;
    std::vector<ld> p(bins,0);
    for (int b=0;b<bins;++b) {
        ld s=0;
        const size_t st=lo+size_t(b)*width, en=st+width;
        for (size_t n=st;n<en;++n) {
            const i64 E = i64(2)*i64(A[n]) - i64(n);
            s += std::llabs(E);
        }
        p[b]=s/ld(width);
    }
    const ld mean=std::accumulate(p.begin(),p.end(),ld(0))/ld(bins);
    ld v=0; for(const ld x:p){const ld d=x-mean;v+=d*d;} v/=ld(bins);
    const ld sd=std::sqrt(v);
    for(ld& x:p) x=(x-mean)/(sd>0?sd:1);
    return p;
}

static ld env_r(const std::vector<u32>& A, int k) {
    const size_t lo=size_t(1)<<k;
    const auto x=env_profile(A,lo,lo+(size_t(1)<<(k-1)));
    const auto y=env_profile(A,size_t(1)<<(k-1),size_t(1)<<k);
    ld s=0;for(int i=0;i<64;++i)s+=x[i]*y[i];
    return s/64;
}

static void assert_close(ld got, ld want, ld tol, const std::string& label) {
    if (std::fabsl(got-want) > tol) {
        std::cerr << "audit mismatch " << label << ": got " << std::setprecision(12)
                  << double(got) << " want " << double(want) << "\n";
        std::exit(6);
    }
}

struct PairMetric { ld m, center, amp; };

static ld loglog_exponent_last5(const std::vector<PairMetric>& P) {
    if (P.size() < 5) return std::numeric_limits<ld>::quiet_NaN();
    const size_t start=P.size()-5;
    ld sx=0, sy=0, sxx=0, sxy=0;
    const ld n=5;
    for(size_t i=start;i<P.size();++i){
        const ld x=std::log(P[i].m), y=std::log(P[i].amp);
        sx+=x; sy+=y; sxx+=x*x; sxy+=x*y;
    }
    const ld slope=(sxy-sx*sy/n)/(sxx-sx*sx/n);
    return -slope;
}

static ld cv_scaled_last5(const std::vector<PairMetric>& P, ld p) {
    if (P.size() < 5) return std::numeric_limits<ld>::quiet_NaN();
    const size_t start=P.size()-5;
    std::vector<ld> v;
    for(size_t i=start;i<P.size();++i) v.push_back(P[i].amp*std::pow(P[i].m,p));
    const ld mean=std::accumulate(v.begin(),v.end(),ld(0))/ld(v.size());
    ld ss=0;for(const ld x:v){const ld d=x-mean;ss+=d*d;}
    return std::sqrt(ss/ld(v.size()))/mean;
}

int main(int argc,char**argv){
    int K=29;
    if(argc>1) K=std::stoi(argv[1]);
    if(K<27 || K>29){ std::cerr<<"K must be 27..29\n"; return 2; }
    const size_t N=size_t(1)<<K;
    std::vector<u32> A(N+1,0);
    A[1]=A[2]=1;
    for(size_t n=3;n<=N;++n){
        const u32 p=A[n-1];
        if(p<1 || p>=n){std::cerr<<"range breach predecessor at "<<n<<"\n";return 3;}
        const u32 v=A[p]+A[n-p];
        if(v<1 || v>n){std::cerr<<"range breach at "<<n<<"\n";return 4;}
        A[n]=v;
    }
    for(int k=1;k<=K;++k){
        const size_t n=size_t(1)<<k;
        if(A[n] != (size_t(1)<<(k-1))){std::cerr<<"clock breach k="<<k<<"\n";return 5;}
    }

    std::vector<ld> fine(K,0), env(K,0);
    for(int k=8;k<=K-1;++k){ fine[k]=fine_r(A,k); env[k]=env_r(A,k); }

    const ld oldFine[5] = {0.29829476L,0.30163436L,0.29847668L,0.30148720L,0.29871683L};
    const ld oldEnv [5] = {0.29812658L,0.30145064L,0.29832370L,0.30130651L,0.29857029L};
    for(int i=0;i<5;++i){
        assert_close(fine[22+i],oldFine[i],5e-9L,"fine k="+std::to_string(22+i));
        assert_close(env [22+i],oldEnv [i],5e-9L,"env k="+std::to_string(22+i));
    }

    std::vector<PairMetric> Pf, Pe;
    for(int k=17;k+1<=K-1;k+=2){
        if(!(fine[k]>.3L && fine[k+1]<.3L)){
            std::cerr<<"fine parity sign gate failed at pair "<<k<<","<<k+1<<"\n";return 7;
        }
        if(!(env[k]>.3L || k<19) || !(env[k+1]<.3L)) {
            // envelope odd k=17 is above .3; this looser form simply keeps the audit explicit.
            if(!(env[k]>.3L && env[k+1]<.3L)){
                std::cerr<<"env parity sign gate failed at pair "<<k<<","<<k+1<<"\n";return 8;
            }
        }
        const ld m=ld(k)+.5L;
        Pf.push_back({m,(fine[k]+fine[k+1])/2,(fine[k]-fine[k+1])/2});
        Pe.push_back({m,(env[k]+env[k+1])/2,(env[k]-env[k+1])/2});
    }

    std::cout<<"[audit] Conway computed to 2^"<<K<<" = "<<N
             <<"; range + clock + published k22..26 values + parity gates PASSED\n";
    std::cout<<"\nFine grain parity pairs\n";
    std::cout<<"  pair    m       center          halfspan       a*m        a*m^1.5     a*m^2\n";
    for(size_t i=0;i<Pf.size();++i){
        const int k=17+int(2*i); const auto& q=Pf[i];
        std::cout<<std::setw(2)<<k<<"/"<<std::setw(2)<<k+1<<"  "<<std::fixed<<std::setprecision(1)<<double(q.m)
                 <<"  "<<std::setprecision(10)<<double(q.center)
                 <<"  "<<double(q.amp)
                 <<"  "<<double(q.amp*q.m)
                 <<"  "<<double(q.amp*std::pow(q.m,1.5L))
                 <<"  "<<double(q.amp*q.m*q.m)<<"\n";
    }
    std::cout<<"last-five log-log exponent p = "<<std::setprecision(9)<<double(loglog_exponent_last5(Pf))<<"\n";
    std::cout<<"last-five CV of a*m       = "<<double(cv_scaled_last5(Pf,1.0L))<<"\n";
    std::cout<<"last-five CV of a*m^1.5   = "<<double(cv_scaled_last5(Pf,1.5L))<<"\n";
    std::cout<<"last-five CV of a*m^2     = "<<double(cv_scaled_last5(Pf,2.0L))<<"\n";

    std::cout<<"\nEnvelope grain parity pairs\n";
    std::cout<<"  pair    m       center          halfspan       a*m^1.5\n";
    for(size_t i=0;i<Pe.size();++i){
        const int k=17+int(2*i); const auto& q=Pe[i];
        std::cout<<std::setw(2)<<k<<"/"<<std::setw(2)<<k+1<<"  "<<std::fixed<<std::setprecision(1)<<double(q.m)
                 <<"  "<<std::setprecision(10)<<double(q.center)
                 <<"  "<<double(q.amp)
                 <<"  "<<double(q.amp*std::pow(q.m,1.5L))<<"\n";
    }
    std::cout<<"last-five log-log exponent p = "<<std::setprecision(9)<<double(loglog_exponent_last5(Pe))<<"\n";
    std::cout<<"last-five CV of a*m^1.5   = "<<double(cv_scaled_last5(Pe,1.5L))<<"\n";

    std::cout<<"\nNewest raw rungs\n";
    for(int k=27;k<=K-1;++k){
        std::cout<<"k="<<k<<" fine="<<std::showpos<<std::fixed<<std::setprecision(8)<<double(fine[k])
                 <<" env="<<double(env[k])<<std::noshowpos<<"\n";
    }
    return 0;
}