python3 <<'PY' import math,random a=0.25; b=1.0/12 def bestK(y1,y2,K=4,N=200000): alpha=2*a+2*b-y1-y2; mu1=a+b-y1; nu1=a+b best=-1 for _ in range(N): w=[random.random() for i in range(K)]; sw=sum(w); q=[alpha*v/sw for v in w] u=[random.random() for i in range(K)]; su=sum(u); aa=[mu1*v/su for v in u] if any(aa[i]>q[i]+1e-12 for i in range(K)): continue cc=[0.0]*K rem=nu1 for i in range(K-1): maxc=min(q[i],rem) ci=random.random()*maxc cc[i]=ci; rem-=ci a4=aa[K-1]; q4=q[K-1] if a4<1e-14: continue target=a for i in range(K-1): target-=aa[i]*cc[i]/q[i] c4=target*q4/a4 if c4<-1e-10 or c4>q4+1e-10: continue cc[K-1]=c4 if abs(sum(cc)-nu1)>1e-9: continue val=0.0 for i in range(K): qi=q[i]; ai=aa[i]; ci=cc[i] bi=qi-ai; di=qi-ci val+=math.sqrt(max(0,a*ai*ci))+math.sqrt(max(0,b*ai*di))+math.sqrt(max(0,b*bi*ci))+math.sqrt(max(0,a*bi*di)) if val>best: best=val return best for y1,y2 in [(0.035,0.095),(0.02,0.22),(0.0,0.0)]: print(y1,y2,bestK(y1,y2,4,200000)) PY