#include <stdio.h>
#include <math.h>
#ifdef _CIVL
#include <civlc.cvh>
#endif
#define epsilon 1.0e-14

/*****************************************************************************/
/* 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   5.77350e-01        /*  1.0/sqrt(3.0) */
#define tsqrt3  1.15470053837925159591885903972e+00        /*  2.0/sqrt(3.0) */
#define a       5.00000000000000000000000000000e-01        /*  1.0/2.0       */
#define b      -1.00000000000000000000000000000e-00        /* -1.0/1.0       */
#define bm1    -2.00000000000000000000000000000e-00        /* -2.0/1.0       */

/*****************************************************************************/
/* 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.                                                                  */

/*****************************************************************************/
/* Optimal derivative calculation courtesy of Paul Hovland (at least we      */
/* think it is optimal).  The original code provided was modified to         */
/* reduce the number of flops and intermediate variables, and improve the    */
/* locality of reference.                                                    */
/*                                                                           */
/* This requires 46 flops.  The function only requires 23 flops.             */
/*****************************************************************************/

int g_fcn(double *obj, double g_obj[6], const double x[6])
{
 static double matr[4], f;
 static double adj_m[4], g;
 static double loc1;

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

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

 /* Calculate det(M). */
 g = matr[0]*matr[3] - matr[1]*matr[2];

 #ifdef _gAssumption
 $assume(g > epsilon);
 #endif
    if (g <= epsilon) {
#ifndef _FIX_G
 *obj = g;
#endif
        return 1;
    }
    
    
 /* Calculate norm(M). */
 f = matr[0]*matr[0] + matr[1]*matr[1] + 
     matr[2]*matr[2] + matr[3]*matr[3];

 /* Calculate objective function. */
 /* loc1 = a * pow(g, b); */
 #ifdef _CIVL
 $assume(g != 0);
 #endif
 loc1 = a / g;
 *obj = f * loc1;

 /* Calculate the derivative of the objective function. */
 f = 2.0*loc1;
 g = b*(*obj)/g; 

 adj_m[0] = matr[0]*f + matr[3]*g;
 adj_m[1] = matr[1]*f - matr[2]*g;
 adj_m[2] = matr[2]*f - matr[1]*g;
 adj_m[3] = matr[3]*f + matr[0]*g;

 loc1 = 0.577350F*adj_m[1];
 
 g_obj[0] = -adj_m[0] - loc1;
 g_obj[1] = adj_m[0] - loc1;
 g_obj[2] = 2.0*loc1;

 loc1 = 0.577350F*adj_m[3];
 g_obj[3] = -adj_m[2] - loc1;
 g_obj[4] = adj_m[2] - loc1;
 g_obj[5] = 2.0*loc1;
 return 0;
}

