python3 - <<'PY'
import math
ln2=math.log(2)
def f(x): return 0.0 if x<=0 else -x*math.log(x)/ln2
def H(x): return f(x)+f(1-x)
def obj(r,b,c):
    pA=20-b-c-3/r;pB=3/r+b;a1=6-3/r;a=5-b-c
    if pA<=0 or pB<=0 or a1<0 or a<0 or c<3 or c<4-b*r: return 1e9
    HX=(pA/20)*(f(9/pA)+f(a1/pA)+f(a/pA))+(pB/20)*(f((3/r)/pB)+f(b/pB))
    HY=(pB/20)*H(r)+(c/20)*(f(3/c)+f((1-b*r)/c)+f((c-4+b*r)/c))
    return 2.1822287-HX-HY
def golden(fun,a,b,tol=1e-13):
    g=(math.sqrt(5)-1)/2;c=b-g*(b-a);d=a+g*(b-a);fc=fun(c);fd=fun(d)
    while b-a>tol:
        if fc<fd: b,d,fd=d,c,fc;c=b-g*(b-a);fc=fun(c)
        else: a,c,fc=c,d,fd;d=a+g*(b-a);fd=fun(d)
    return (a+b)/2
v=[0.77006,0.76571,3.60395]
for it in range(20):
    for j,st in enumerate([0.001,0.001,0.002]):
        def g(x,j=j): vv=v.copy();vv[j]=x;return obj(*vv)
        v[j]=golden(g,v[j]-st,v[j]+st)
r,b,c=v;pA=20-b-c-3/r;pB=3/r+b;a1=6-3/r;a=5-b-c
lam=[[9/pA,a1/pA,a/pA],[0,(3/r)/pB,b/pB],[0,0,1]]
mu=[[0,0,1],[0,r,1-r],[3/c,(1-b*r)/c,(c-4+b*r)/c]]
print('v',v,obj(*v))
print('lam',lam,'mu',mu)
def soft(z):
    m=max(z); e=[2**(x-m) for x in z]; s=sum(e); return [x/s for x in e]
def res(xw):
    W=[[xw[0],xw[1],xw[2]],[xw[3],xw[4],xw[5]],[xw[6],xw[7],0.0]]
    out=[]
    for i in range(3):
        wm=[sum(W[j][k]*mu[i][k] for k in range(3)) for j in range(3)]
        lh=soft([-wm[j] for j in range(3)])
        wtl=[sum(W[j][k]*lam[i][j] for j in range(3)) for k in range(3)]
        mh=soft([-wtl[k] for k in range(3)])
        out+=[lh[j]-lam[i][j] for j in range(3)]+[mh[k]-mu[i][k] for k in range(3)]
    return out
def norm(x): return sum(t*t for t in x)
xw=[0.0]*8
for it in range(100):
    r0=res(xw); e0=norm(r0)
    if e0<1e-24: break
    h=1e-5; J=[]
    for j in range(8):
        xp=xw.copy();xp[j]+=h;xm=xw.copy();xm[j]-=h
        rp=res(xp);rm=res(xm);J.append([(rp[i]-rm[i])/(2*h) for i in range(18)])
    JTJ=[[sum(J[a][i]*J[b][i] for i in range(18)) for b in range(8)] for a in range(8)]
    JTr=[sum(J[a][i]*r0[i] for i in range(18)) for a in range(8)]
    A=[row[:]+[-JTr[i]] for i,row in enumerate(JTJ)]
    for col in range(8):
        piv=max(range(col,8),key=lambda rr:abs(A[rr][col]));A[col],A[piv]=A[piv],A[col]
        pv=A[col][col]
        if abs(pv)<1e-14: break
        for jj in range(col,9): A[col][jj]/=pv
        for rr in range(8):
            if rr!=col:
                fac=A[rr][col]
                for jj in range(col,9): A[rr][jj]-=fac*A[col][jj]
    dx=[A[i][8] for i in range(8)]
    t=1.0
    while t>1e-10:
        xn=[xw[i]+t*dx[i] for i in range(8)]
        if norm(res(xn))<e0: xw=xn;break
        t*=0.5
    else: break
print('res',norm(res(xw)),xw)
PY