python3 <<'PY'
import math,random
a=0.25; b=1.0/12
def G(y1,y2,q,r):
    alpha=2*a+2*b-y1-y2; mu1=a+b-y1; mu2=b+a-y2; nu1=a+b
    q2=alpha-q
    if y1<-1e-12 or y1>b+1e-12 or y2<-1e-12 or y2>a+1e-12 or q<=1e-14 or q2<=1e-14: return -1
    lo=max(0.0,q-mu2); hi=min(q,mu1)
    if r<lo-1e-12 or r>hi+1e-12: return -1
    den=r*q2/q-(mu1-r)
    if abs(den)<1e-14: return -1
    c=(a*q2-(mu1-r)*nu1)/den
    if c<-1e-12 or c>q+1e-12 or c>nu1+1e-12 or nu1-c>q2+1e-12: return -1
    r2=mu1-r; c2=nu1-c
    t1=math.sqrt(max(0,a*r*c))+math.sqrt(max(0,b*r*(q-c)))+math.sqrt(max(0,b*(q-r)*c))+math.sqrt(max(0,a*(q-r)*(q-c)))
    t2=math.sqrt(max(0,a*r2*c2))+math.sqrt(max(0,b*r2*(q2-c2)))+math.sqrt(max(0,b*(q2-r2)*c2))+math.sqrt(max(0,a*(q2-r2)*(q2-c2)))
    A=t1+t2
    y3=a; gam=y1+y2+y3
    C=math.sqrt(gam)*(math.sqrt(max(0,b*y1))+math.sqrt(max(0,a*y2))+math.sqrt(max(0,a*y3)))
    B=b**1.5
    return A+B+C
# boundary opt
y1b=0.024697867630391723; y2b=0.23045484476568612
best=-1; arg=None
for it in range(300000):
  y1=max(0,min(b,y1b+random.uniform(-0.01,0.01)))
  y2=max(0,min(a,y2b+random.uniform(-0.01,0.01)))
  alpha=2*a+2*b-y1-y2; mu1=a+b-y1; mu2=b+a-y2
  # choose q fraction of alpha
  q=alpha*(1-random.random()*0.02)
  q2=alpha-q
  lo=max(0.0,q-mu2); hi=min(q,mu1)
  if hi<lo: continue
  r=lo+(hi-lo)*random.random()
  g=G(y1,y2,q,r)
  if g>best: best=g; arg=(y1,y2,q,r)
print("rand",best,arg)
# pattern search
if arg:
  v=list(arg); step=[0.002,0.002,0.002,0.002]
  for round_ in range(200):
    improved=False
    for j in range(4):
      for sgn in (1,-1):
        w=v[]; w[j]+=sgn*step[j]
        g=G(*w)
        if g>G(*v): v=w; improved=True
    if not improved:
      step=[s*0.6 for s in step]
    if max(step)<1e-12: break
  print("ref",G(*v),v,"q2",2*a+2*b-v[0]-v[1]-v[2])
PY