check_upper_half.py

"""
Verify Theorem R189: R_d(a) <= 1 for d >= a/2.
Also check the restricted {1,2} path counts.
"""
from collections import defaultdict

def compute_Nd(max_a):
    N = [defaultdict(int) for _ in range(max_a + 1)]
    N[1][0] = 1
    for c in range(1, max_a + 1):
        if not N[c]:
            continue
        n = c*c + 1
        divs = []
        d = 1
        while d*d <= n:
            if n % d == 0:
                divs.append(d)
                if d*d != n:
                    divs.append(n // d)
            d += 1
        children_seen = set()
        for delta in divs:
            child = c + delta
            if child <= max_a and child not in children_seen:
                children_seen.add(child)
                for depth, count in N[c].items():
                    N[child][depth + 1] += count
    return N

def check_upper_half(max_a=60):
    N = compute_Nd(max_a)
    violations = []
    
    for a in range(2, max_a + 1):
        depths = sorted(N[a].keys())
        for i in range(len(depths)-1):
            d = depths[i]
            if d >= a/2:
                if N[a][depths[i+1]] > N[a][d]:
                    violations.append((a, d, N[a][d], N[a][depths[i+1]]))
    
    if violations:
        print("VIOLATIONS of R_d(a) <= 1 for d >= a/2:")
        for a, d, nd, nd1 in violations[:10]:
            print(f"  a={a}, d={d}: N_d={nd}, N_{{d+1}}={nd1}, ratio={nd1/nd:.4f}")
    else:
        print(f"ALL PASS: N_d(a) <= N_{{d+1}}(a) for all d >= a/2, a <= {max_a}")
        print("(i.e., sequence is non-increasing for d >= a/2)")

# Also check restricted {1,2} paths
def compute_restricted(max_a):
    """Count paths using only steps 1 and 2 (2 only when c is odd)."""
    h = [defaultdict(int) for _ in range(max_a + 1)]
    h[1][0] = 1
    for c in range(1, max_a + 1):
        if not h[c]:
            continue
        # Always can go to c+1
        if c+1 <= max_a:
            for depth, count in h[c].items():
                h[c+1][depth+1] += count
        # Can go to c+2 if c is odd (2|c^2+1 iff c odd)
        if c % 2 == 1 and c+2 <= max_a:
            for depth, count in h[c].items():
                h[c+2][depth+1] += count
    return h

def compare_distributions(max_a=30):
    N = compute_Nd(max_a)
    h = compute_restricted(max_a)
    
    print(f"\n{'='*70}")
    print("COMPARISON: Full paths vs {1,2}-only paths")
    print(f"{'='*70}")
    
    for a in [10, 20, 30]:
        ftotal = sum(N[a].values())
        htotal = sum(h[a].values())
        
        mu_full = sum(d*N[a][d] for d in N[a]) / ftotal
        mu_restr = sum(d*h[a][d] for d in h[a]) / htotal if htotal > 0 else 0
        
        print(f"\na={a}: f_full={ftotal}, f_restricted={htotal}, "
              f"mu_full={mu_full:.2f}, mu_restricted={mu_restr:.2f}")

if __name__ == "__main__":
    check_upper_half(60)
    compare_distributions(30)