python3 <<'PY' import math a=0.25; b=1.0/12 def Aval(y1,y2,q,r): alpha=2*a+2*b-y1-y2; mu1=a+b-y1; nu1=a+b; q2=alpha-q if q<=0 or q2<=0: return -1 lo=max(0.0,q-(alpha-mu1)); hi=min(q,mu1) if rhi+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-10 or c>q+1e-10 or c>nu1+1e-10 or nu1-c>q2+1e-10: return -1 r2=mu1-r; c2=nu1-c return (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)))+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)))) def best2(y1,y2,N=400): alpha=2*a+2*b-y1-y2; mu1=a+b-y1; mu2=alpha-mu1 best=-1; arg=None for iq in range(1,N): q=alpha*iq/N lo=max(0.0,q-mu2); hi=min(q,mu1) for ir in range(N+1): r=lo+(hi-lo)*ir/N v=Aval(y1,y2,q,r) if v>best: best=v; arg=(q,r) return best,arg def opty3(y1,y2): best=-1; ba=None for j in range(2001): y3=a*j/2000 beta=a+b-y3; B=math.sqrt(beta)*(b+math.sqrt(max(0,a*(a-y3)))) 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))) if B+C>best: best=B+C; ba=y3 return best,ba for y1,y2 in [(0.035,0.095),(0.035,0.1),(0.03,0.09),(0.04,0.1)]: A,arg=best2(y1,y2,300) bc,y3=opty3(y1,y2) print(y1,y2,"A",A,"bc",bc,"y3",y3,"tot",A+bc,"arg",arg) PY