double FourSubtree( double Wab, double Mab, double Wac, double Mac,
 double Wad, double Mad, double Wbc, double Mbc,
 double Wbd, double Mbd, double Wcd, double Mcd,
  double *L)
{ double Z1[25], Z2[25], Z3[25], B1[5], B2[5], B3[5], t;
  int i, it, j;
Z1[0] = Z2[0] = Wab+Wac+Wad;
Z1[1] = Z2[1] = Wab;
Z1[2] = Z2[2] = Wac+Wad;
Z1[3] = Z2[3] = Wac;
Z1[4] = Z2[4] = Wad;
B1[0] = B2[0] = Mab+Mac+Mad;
Z1[5] = Z2[5] = Wab;
Z1[6] = Z2[6] = Wab+Wbc+Wbd;
Z1[7] = Z2[7] = Wbc+Wbd;
Z1[8] = Z2[8] = Wbc;
Z1[9] = Z2[9] = Wbd;
B1[1] = B2[1] = Mab+Mbc+Mbd;
Z1[10] = Z2[10] = Wac+Wad;
Z1[11] = Z2[11] = Wbc+Wbd;
Z1[12] = Z2[12] = Wac+Wad+Wbc+Wbd;
Z1[13] = Z2[13] = Wac+Wbc;
Z1[14] = Z2[14] = Wad+Wbd;
B1[2] = B2[2] = Mac+Mad+Mbc+Mbd;
Z1[15] = Z2[15] = Wac;
Z1[16] = Z2[16] = Wbc;
Z1[17] = Z2[17] = Wac+Wbc;
Z1[18] = Z2[18] = Wac+Wbc+Wcd;
Z1[19] = Z2[19] = Wcd;
B1[3] = B2[3] = Mac+Mbc+Mcd;
Z1[20] = Z2[20] = Wad;
Z1[21] = Z2[21] = Wbd;
Z1[22] = Z2[22] = Wad+Wbd;
Z1[23] = Z2[23] = Wcd;
Z1[24] = Z2[24] = Wad+Wbd+Wcd;
B1[4] = B2[4] = Mad+Mbd+Mcd;
for( it=0; it<10; it++ ) {
   for( i=0; i<25; i++ ) Z3[i] = Z2[i];
   for( i=0; i<5; i++ ) B3[i] = B2[i];
   gausselimi(Z3,B3,L,5);
   for( i=0; i<5; i++ ) {
      if( L[i] < MinLen-1e-10 ) {
           for( j=0; j<5; j++ ) Z2[5*i+j] = 0;
           Z2[6*i] = 1; B2[i] = MinLen;
           break; }
      else if( fabs(L[i]-MinLen)<1e-10 ) {
               t = -B1[i];
               for( j=0; j<5; j++ ) t += Z1[5*i+j]*L[j];
               if( t<-1e-10 ) {
                  for( j=0; j<5; j++ ) Z2[5*i+j] = Z1[5*i+j];
                  B2[i] = B1[i];
                  break; }
           }
      }
   if( i>=5 ) break;
   }
{ double t1, t4, t7, t10, t13, t20;
      t1 = L[2]+L[3]+L[0];
      t4 = L[2]+L[4]+L[0];
      t7 = L[3]+L[2]+L[1];
      t10 = L[4]+L[2]+L[1];
      t13 = L[0]+L[1];
      t20 = L[3]+L[4];
      return( t1*t1*Wac+t4*t4*Wad+t7*t7*Wbc+t10*t10*Wbd+t13*t13*Wab+
	t20*t20*Wcd - 2*(t1*Mac+t4*Mad+t7*Mbc+t10*Mbd+t13*Mab+Mcd*t20) );
      }
}