python3 <<'PY'
import math
a=0.25; b=1.0/12
y1=0.023739012910611847; y2=0.2300466238859805; q=0.38975347654236847; r=0.2910741941684634
mu1=a+b-y1; nu1=a+b; alpha=2*a+2*b-y1-y2; q2=alpha-q
den=r*(alpha-q)/q-(mu1-r); c=(a*(alpha-q)-(mu1-r)*nu1)/den
r2=mu1-r; c2=nu1-c; y3=a; gam=y1+y2+y3
comps=[]
comps.append((q,[r/q,1-r/q,0],[c/q,0,1-c/q]))
comps.append((q2,[r2/q2,1-r2/q2,0],[c2/q2,0,1-c2/q2]))
comps.append((b,[0,0,1],[0,1,0]))
comps.append((gam,[y1/gam,y2/gam,y3/gam],[0,0,1]))
Pgen=[[0.0]*3 for _ in range(3)]
for w,rr,ss in comps:
  for i in range(3):
    for j in range(3):
      Pgen[i][j]+=w*rr[i]*ss[j]
P=[[a,0,b],[b,0,a],[0,b,a]]
print("c",c,"q2",q2,"bounds",0<=c<=q,0<=c2<=q2)
for i in range(3): print([round(Pgen[i][j]-P[i][j],12) for j in range(3)])
aff=0
for w,rr,ss in comps:
  t=0
  for i in range(3):
    for j in range(3):
      t+=math.sqrt(max(0,P[i][j]*rr[i]*ss[j]))
  aff+=w*t
  print("comp",w,"aff contrib",w*t)
print("aff",aff,"cost",1-aff)
PY