double FiveSubtree( double Wab, double Mab, double Wac, double Mac, double Wad, double Mad, double Wae, double Mae, double Wbc, double Mbc, double Wbd, double Mbd, double Wbe, double Mbe, double Wcd, double Mcd, double Wce, double Mce, double Wde, double Mde, double *L) { double Z1[49], Z2[49], Z3[49], B1[7], B2[7], B3[7], t; int i, it, j; Z1[0] = Z2[0] = Wab+Wac+Wad+Wae; Z1[1] = Z2[1] = Wab; Z1[2] = Z2[2] = Wac+Wad+Wae; Z1[3] = Z2[3] = Wac; Z1[4] = Z2[4] = Wad+Wae; Z1[5] = Z2[5] = Wad; Z1[6] = Z2[6] = Wae; B1[0] = B2[0] = Mab+Mac+Mad+Mae; Z1[7] = Z2[7] = Wab; Z1[8] = Z2[8] = Wab+Wbc+Wbd+Wbe; Z1[9] = Z2[9] = Wbc+Wbd+Wbe; Z1[10] = Z2[10] = Wbc; Z1[11] = Z2[11] = Wbd+Wbe; Z1[12] = Z2[12] = Wbd; Z1[13] = Z2[13] = Wbe; B1[1] = B2[1] = Mab+Mbc+Mbd+Mbe; Z1[14] = Z2[14] = Wac+Wad+Wae; Z1[15] = Z2[15] = Wbc+Wbd+Wbe; Z1[16] = Z2[16] = Wac+Wad+Wbc+Wbd+Wae+Wbe; Z1[17] = Z2[17] = Wac+Wbc; Z1[18] = Z2[18] = Wad+Wbd+Wae+Wbe; Z1[19] = Z2[19] = Wad+Wbd; Z1[20] = Z2[20] = Wae+Wbe; B1[2] = B2[2] = Mac+Mad+Mbc+Mbd+Mae+Mbe; Z1[21] = Z2[21] = Wac; Z1[22] = Z2[22] = Wbc; Z1[23] = Z2[23] = Wac+Wbc; Z1[24] = Z2[24] = Wac+Wbc+Wcd+Wce; Z1[25] = Z2[25] = Wcd+Wce; Z1[26] = Z2[26] = Wcd; Z1[27] = Z2[27] = Wce; B1[3] = B2[3] = Mac+Mbc+Mcd+Mce; Z1[28] = Z2[28] = Wad+Wae; Z1[29] = Z2[29] = Wbd+Wbe; Z1[30] = Z2[30] = Wad+Wbd+Wae+Wbe; Z1[31] = Z2[31] = Wcd+Wce; Z1[32] = Z2[32] = Wad+Wbd+Wcd+Wae+Wbe+Wce; Z1[33] = Z2[33] = Wad+Wbd+Wcd; Z1[34] = Z2[34] = Wae+Wbe+Wce; B1[4] = B2[4] = Mad+Mbd+Mcd+Mae+Mbe+Mce; Z1[35] = Z2[35] = Wad; Z1[36] = Z2[36] = Wbd; Z1[37] = Z2[37] = Wad+Wbd; Z1[38] = Z2[38] = Wcd; Z1[39] = Z2[39] = Wad+Wbd+Wcd; Z1[40] = Z2[40] = Wad+Wbd+Wcd+Wde; Z1[41] = Z2[41] = Wde; B1[5] = B2[5] = Mad+Mbd+Mcd+Mde; Z1[42] = Z2[42] = Wae; Z1[43] = Z2[43] = Wbe; Z1[44] = Z2[44] = Wae+Wbe; Z1[45] = Z2[45] = Wce; Z1[46] = Z2[46] = Wae+Wbe+Wce; Z1[47] = Z2[47] = Wde; Z1[48] = Z2[48] = Wae+Wbe+Wce+Wde; B1[6] = B2[6] = Mae+Mbe+Mce+Mde; for( it=0; it<14; it++ ) { for( i=0; i<49; i++ ) Z3[i] = Z2[i]; for( i=0; i<7; i++ ) B3[i] = B2[i]; gausselimi(Z3,B3,L,7); for( i=0; i<7; i++ ) { if( L[i] < MinLen-1e-10 ) { for( j=0; j<7; j++ ) Z2[7*i+j] = 0; Z2[8*i] = 1; B2[i] = MinLen; break; } else if( fabs(L[i]-MinLen)<1e-10 ) { t = -B1[i]; for( j=0; j<7; j++ ) t += Z1[7*i+j]*L[j]; if( t<-1e-10 ) { for( j=0; j<7; j++ ) Z2[7*i+j] = Z1[7*i+j]; B2[i] = B1[i]; break; } } } if( i>=7 ) break; } { double t1, t4, t7, t10, t13, t16, t19, t22, t29, t36; t1 = L[2]+L[4]+L[0]+L[5]; t4 = L[4]+L[2]+L[1]+L[5]; t7 = L[2]+L[4]+L[6]+L[0]; t10 = L[2]+L[4]+L[6]+L[1]; t13 = L[2]+L[3]+L[0]; t16 = L[3]+L[2]+L[1]; t19 = L[4]+L[3]+L[5]; t22 = L[6]+L[4]+L[3]; t29 = L[0]+L[1]; t36 = L[5]+L[6]; return( t1*t1*Wad+t4*t4*Wbd+t7*t7*Wae+t10*t10*Wbe+t13*t13*Wac+ t16*t16*Wbc+t19*t19*Wcd+t22*t22*Wce+t29*t29*Wab+t36*t36*Wde - 2*(t1*Mad+t4*Mbd+t7*Mae+t10*Mbe+t13*Mac+t16*Mbc+t19*Mcd+t22*Mce+ t29*Mab+Mde*t36) ); } }