python3 <<'PY' import math,random a=0.25; b=1.0/12 y1=0.025697868008998203; y2=0.22745484489188828 def Ared(y1,y2,A): return A def bestK(K,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 r=[aa[i]/q[i] for i in range(K)] c=[0.0]*K rem=nu1 for i in range(K-2): ci=random.random()*min(q[i],rem) c[i]=ci; rem-=ci # solve last two for sum rem and m11 a r1=r[K-2]; r2=r[K-1] if abs(r1-r2)<1e-14: continue T=a-sum(r[i]*c[i] for i in range(K-2)) c1=(T-r2*rem)/(r1-r2); c2=rem-c1 c[K-2]=c1; c[K-1]=c2 if any(c[i]<-1e-10 or c[i]>q[i]+1e-10 for i in range(K)): continue val=0.0 for i in range(K): qi=q[i]; ai=aa[i]; ci=c[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 print("K2 approx",0.304319) for K in (3,4,5): print(K,bestK(K,100000)) PY