source: CIVL/examples/compare/provesa/head.c

#define epsilon 1.00000e-14

/*****************************************************************************/
/* This set of functions reference triangular elements to an equilateral     */
/* triangle.  The input are the coordinates in the following order:          */
/*      [x1 x2 x3 y1 y2 y3]                                                  */
/* A zero return value indicates success, while a nonzero value indicates    */
/* failure.                                                                  */
/*****************************************************************************/
/* Not all compilers substitute out constants (especially the square root).  */
/* Therefore, they are substituted out manually.  The values below were      */
/* calculated on a solaris machine using long doubles. I believe they are    */
/* accurate.                                                                 */
/*****************************************************************************/

 #define sqrt3    .577350269189625797959429519858         /*1.0/sqrt(3.0) */ 
/*#define tsqrt3  1.15470053837925159591885903972          2.0/sqrt(3.0) */
#define a        .500000000000000000000000000000        /*  1.0/2.0       */
#define b      -1.00000000000000000000000000000        /* -1.0/1.0       */
#define bm1    -2.00000000000000000000000000000        /* -2.0/1.0       */

 double tsqrt3;// = 2.0/sqrt3; 
 //double tsqrt3 = 2.0/sqrt3; 

double o_fcn(double *obj, /* const */ double x[6])
{
 /* static */ double matr[4], f;
 /* static */ double g;
 double retval;

 retval = 0.0;

 /* Calculate M = A*inv(W). */
 matr[0] = x[1] - x[0];
 matr[1] = (2.0*x[2] - x[1] - x[0])*sqrt3;

 matr[2] = x[4] - x[3];
 matr[3] = (2.0*x[5] - x[4] - x[3])*sqrt3;

 /* Calculate det(M). */
 g = matr[0]*matr[3] - matr[1]*matr[2];
 if (g <= epsilon) { *obj = g; retval = 1.0; }
 else {
 /* Calculate norm(M). */
 f = matr[0]*matr[0] + matr[1]*matr[1] + 
     matr[2]*matr[2] + matr[3]*matr[3];

 /* Calculate objective function. */
 (*obj) = a * f / g;
 }
 return retval;
}