python3 <<'PY' import math def F(M,N=200): m11,m12=M[0]; m21,m22=M[1] alpha=m11+m12+m21+m22 if alpha<1e-15: return 0.0 mu1=m11+m12; nu1=m11+m21 best=-1.0 for iq in range(1,N): q=alpha*iq/N; q2=alpha-q lo=max(0.0,q-(alpha-mu1)); hi=min(q,mu1) for ir in range(N+1): r=lo+(hi-lo)*ir/N den=r*(alpha-q)/q-(mu1-r) if abs(den)<1e-14: continue c=(m11*(alpha-q)-(mu1-r)*nu1)/den if c<-1e-8 or c>nu1+1e-8 or c>q+1e-8 or nu1-c>q2+1e-8: continue c=min(max(c,0.0),nu1) r1=r/q; s1=c/q; r2=(mu1-r)/q2; s2=(nu1-c)/q2 val=q*(math.sqrt(max(0,m11*r1*s1))+math.sqrt(max(0,m12*r1*(1-s1)))+math.sqrt(max(0,m21*(1-r1)*s1))+math.sqrt(max(0,m22*(1-r1)*(1-s1))))+q2*(math.sqrt(max(0,m11*r2*s2))+math.sqrt(max(0,m12*r2*(1-s2)))+math.sqrt(max(0,m21*(1-r2)*s2))+math.sqrt(max(0,m22*(1-r2)*(1-s2)))) if val>best: best=val return best M=[[0.25,0.0],[1/12,0.05]] print(F(M,200)) PY