// (c)2009 Exstrom Laboratories LLC

function expectation(pv,rv){var n=pv.length;var s=0.0;for(var k=0;k<n;++k)s+=pv[k]*rv[k];return s;}
function dgsingle(rv,pv,f){var n=pv.length;var s=0.0;for(var i=0;i<n;++i)
s+=pv[i]*rv[i]/(1.0+rv[i]*f);return s;}
function dgidouble(rv1,rv2,pv1,pv2,f1,f2){var n1=pv1.length;var n2=pv2.length;var s1=0.0;var s2=0.0;var d=0.0;for(var i=0;i<n1;++i)
for(var j=0;j<n2;++j){d=pv1[i]*pv2[j]/(1.0+rv1[i]*f1+rv2[j]*f2);s1+=rv1[i]*d;s2+=rv2[j]*d;}
return[s1,s2];}
function single_iter(rv,pv){var n=pv.length;var passback=function(x){var pk,rk,d1,d2,d3;var s1=0.0;var s2=0.0;var s3=0.0;for(var k=0;k<n;++k){pk=pv[k];rk=rv[k];d1=1+rk*x;d2=d1==0?rk*100000:rk/d1;d3=pk*d2;s1+=d2;s2+=d3;s3+=d2*d3;}
return x-s2/(s1*s2-s3);}
return passback;}
function idouble_iter(rv1,rv2,pv1,pv2){var n1=pv1.length;var n2=pv2.length;var passback=function(v0){var pk,rk;var pl,sl;var d1;var Rkl,Skl,pklRkl,pklSkl;var s1=0.0;var s2=0.0;var s3=0.0;var s4=0.0;var s5=0.0;var s6=0.0;var s7=0.0;var x=v0[0];var y=v0[1];for(var k=0;k<n1;++k){var pk=pv1[k];var rk=rv1[k];for(var l=0;l<n2;++l){pl=pv2[l];sl=rv2[l];d1=1+rk*x+sl*y;if(d1==0)d1=0.000001;Rkl=rk/d1;Skl=sl/d1;pklRkl=pk*pl*Rkl;pklSkl=pk*pl*Skl;s1+=Rkl;s2+=Skl;s3+=pklRkl;s4+=pklSkl;s5+=pklRkl*Rkl;s6+=pklSkl*Skl;s7+=pklRkl*Skl;}}
var fx=s1*s3-s5;var fy=s2*s3-s7;var gx=s1*s4-s7;var gy=s2*s4-s6;var det=fx*gy-fy*gx;var dx=(s4*fy-s3*gy)/det;var dy=(s3*gx-s4*fx)/det;return[x+dx,y+dy]};return passback;}
function kfsingle(rv,pv){var fiter=single_iter(rv,pv);var rmax=rv.max();var rmin=rv.min();if(rmax<=0)return-1000000.0;if(rmin>=0)return 1000000.0;var xmax,xmin,x;if(expectation(pv,rv)>=0){xmax=-1.0/rmin;xmin=0.0;x=xmax/2.0;}
else{xmax=0.0;xmin=-1.0/rmax;x=xmin/2.0;}
while((xmax-xmin)>1.0e-6){if(dgsingle(rv,pv,x)>0)
xmin=x;else
xmax=x;x=(xmax+xmin)/2.0;}
x0=x;x1=fiter(x);while(Math.abs(x1-x0)>1.0e-10){x0=x1;x1=fiter(x1);}
return x1;}
function kfidouble(rv1,rv2,pv1,pv2){var fiter=idouble_iter(rv1,rv2,pv1,pv2);var fs1=kfsingle(rv1,pv1);var fs2=kfsingle(rv2,pv2);if(Math.abs(fs1)==1000000.0)
if(Math.abs(fs2)==1000000.0)
return[fs1,fs2];else
return[fs1,0];if(Math.abs(fs2)==1000000.0)
return[0,fs2];var xmax,xmin,x;if(fs1>0){xmax=fs1;xmin=0;}
else{xmax=0;xmin=fs1;}
x=fs1/2.0;var ymax,ymin,y;if(fs2>0){ymax=fs2;ymin=0;}
else{ymax=0;ymin=fs2;}
y=fs2/2.0;while(((xmax-xmin)>1.0e-10)&&((ymax-ymin)>1.0e-10)){var fv=dgidouble(rv1,rv2,pv1,pv2,x,y);if(fv[0]>0)
xmin=x;else
xmax=x;if(fv[1]>0)
ymin=y;else
ymax=y;x=(xmax+xmin)/2.0;y=(ymax+ymin)/2.0;}
var v0=[x,y];var v1=fiter([x,y]);while((Math.abs(v1[0]-v0[0])>1.0e-12)||(Math.abs(v1[1]-v0[1])>1.0e-12)){v0=v1;v1=fiter(v1);}
return v1;}
function multi_bond_iter(prbv,retv,retb){var n=prbv.length;var passback=function(x){var pk,rk,d0,d1,d2;var s1=0;var s2=0;var s3=0;for(var k=0;k<n;++k){pk=prbv[k];rk=retv[k];d0=1+retb+(rk-retb)*x;d1=d0==0?(rk-retb)*100000:(rk-retb)/d0;d2=pk*d1;s1+=d1;s2+=d2;s3+=d1*d2;}
return x-s2/(s1*s2-s3);}
return passback;}
function frac_multi_bond(prbv,retv,retb){var fiter=multi_bond_iter(prbv,retv,retb);var gmax=Number.NEGATIVE_INFINITY;var lmax=Number.POSITIVE_INFINITY;var xmin=0;var xmax=0;var xstart=0;var x0,x1;for(var i=0;i<retv.length;++i){if(retv[i]>gmax)gmax=retv[i];if(retv[i]<lmax)lmax=retv[i];}
if(expectation(prbv,retv)>=retb)
xmax=lmax>=retb?Number.POSITIVE_INFINITY:(1+retb)/(retb-lmax);else
xmin=gmax<=retb?Number.NEGATIVE_INFINITY:(1+retb)/(retb-gmax);x1=(xmin+xmax)/2;if(Math.abs(x1)!=Number.POSITIVE_INFINITY)
do{x0=x1;x1=fiter(x1);if(x1>xmax)x1=xmax;if(x1<xmin)x1=xmin;}while(Math.abs(x1-x0)>0.00001);return x1;}