import numpy as np, collections
# FOLLOW (the pull): does equidistribution survive climbing the 3-tower?
# branch values v==4 mod 6; census mod 2*3^k for k=3,4,5 (54,162,486).
DEPTH=60
level=np.array([1],dtype=np.uint64)
# collect all branch values across depths once
allbranch=[]
for d in range(DEPTH):
    b=level[(level%6==4)&(level!=4)]
    allbranch.append(b.copy())
    level=np.concatenate((level*2, ((b-1)//3).astype(np.uint64)))
B=np.concatenate(allbranch)
print("branch values total:", B.size)
for M in [54,162,486,1458]:
    res=B%M
    # classes present must all be ==4 mod 6
    classes=[r for r in range(M) if r%6==4]
    counts=np.array([np.count_nonzero(res==c) for c in classes])
    sh=counts/counts.sum()
    print("mod %-5d  %d classes  min*K=%.3f  max*K=%.3f  (K=%d)"%(
        M, len(classes), sh.min()*len(classes), sh.max()*len(classes), len(classes)))

# THE PROPOSED ENGINE, checked directly: is 2 a primitive root mod 3^k?
def order(a,m):
    o=1; x=a%m
    while x!=1: x=(x*a)%m; o+=1
    return o
print("\nis 2 a primitive root mod 3^k?  (order should == phi(3^k)=2*3^(k-1))")
for k in range(1,9):
    m=3**k; phi=2*3**(k-1)
    print("  3^%d=%-7d  ord2=%-7d  phi=%-7d  primitive=%s"%(k,m,order(2,m),phi,order(2,m)==phi))