// THE SIBLING HORIZON — extend THE SIBLING SIGN on Conway alone.
//
// Computes Conway's recurrence
//   a(1)=a(2)=1; a(n)=a(a(n-1))+a(n-a(n-1))
// to 2^K, verifies 1<=a(n)<=n and a(2^k)=2^(k-1), then asks the exact
// same sibling-correlation questions as TOOLS/sibling_sign.py at farther rungs.
// Fine grain uses E(n)=2a(n)-n (a factor-2 version of e(n)=a(n)-n/2,
// which leaves Pearson r unchanged). Envelope grain uses 64-bin means of |E|.
//
// Evidence class: COMPUTED.  The goal is not to prove a limit but to decide
// whether the apparent +0.30 pinning bends toward 1/3, toward 0.30, or elsewhere.
// Iris, 2026-07-26.

#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) {
        ld x = ld(2)*A[x0+j] - ld(x0+j);
        ld y = ld(2)*A[y0+j] - ld(y0+j);
        sx += x; sy += y; sxx += x*x; syy += y*y; sxy += x*y;
    }
    ld n = ld(half);
    ld vx = sxx - sx*sx/n;
    ld vy = syy - sy*sy/n;
    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 m = len / bins;
    std::vector<ld> p(bins,0);
    for (int b=0;b<bins;++b) {
        ld s=0;
        size_t st=lo+size_t(b)*m, en=st+m;
        for (size_t n=st;n<en;++n) {
            i64 E = i64(2)*i64(A[n]) - i64(n);
            s += std::llabs(E);
        }
        p[b]=s/ld(m);
    }
    ld mean=std::accumulate(p.begin(),p.end(),ld(0))/ld(bins);
    ld v=0; for(ld x:p){ld d=x-mean;v+=d*d;} v/=ld(bins);
    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) {
    size_t lo=size_t(1)<<k;
    auto x=env_profile(A,lo,lo+(size_t(1)<<(k-1)));
    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;
}

int main(int argc,char**argv){
    int K=25;
    if(argc>1) K=std::stoi(argv[1]);
    if(K<8 || K>29){ std::cerr<<"K must be 8..29\n"; return 2; }
    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){
        u32 p=A[n-1];
        if(p<1 || p>=n){std::cerr<<"range breach predecessor at "<<n<<"\n";return 3;}
        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){
        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::cout<<"[audit] Conway computed to 2^"<<K<<" = "<<N<<"; range + clock passed\n";
    std::cout<<"  k       fine_r       env_r      fine-0.30     fine-1/3\n";
    std::vector<ld> tail;
    for(int k=8;k<=K-1;++k){
        ld rf=fine_r(A,k); ld re=env_r(A,k); tail.push_back(rf);
        std::cout<<std::setw(3)<<k<<"  "<<std::showpos<<std::fixed<<std::setprecision(8)
                 <<std::setw(12)<<double(rf)<<"  "<<std::setw(12)<<double(re)
                 <<"  "<<std::setw(12)<<double(rf-ld(.30))
                 <<"  "<<std::setw(12)<<double(rf-ld(1.0L/3.0L))<<std::noshowpos<<"\n";
    }
    if(tail.size()>=8){
        ld m=0;for(size_t i=tail.size()-8;i<tail.size();++i)m+=tail[i];m/=8;
        std::cout<<"tail mean last 8 fine r = "<<std::fixed<<std::setprecision(10)<<double(m)<<"\n";
    }
}