// assume that F2CAD_USE_Adolc is defined
# include <f2cad/f2cad.hpp>

namespace {
	static bool Recording  = false; // currently recording a function
	static int Domain      = 0;     // domain dimension
	static int Range       = 0;     // range dimension
	static double *Xmemory = 0;     // domain value and differential
	static double **X      = 0;
	static double *Ymemory = 0;     // range value and differential
	static double **Y      = 0;
}

namespace f2cad {

void Independent(integer n, doublereal *x)
{
	// check not two Independents in a row
	assert( ! Recording );

	// make sure the set of independent variables not empty
	assert( n > 0 );

	// check for previously defined function (non-zero pointers)
	if( Xmemory != 0 )
	{	assert( Xmemory != 0 );
		delete [] Xmemory;

		assert( Ymemory != 0 );
		delete [] Ymemory;

		assert( X != 0 );
		delete [] X;

		assert( Y != 0 );
		delete [] Y;
	}

	// allocate memory for independent variable value and differential
	// and store domain value in first column of X
	Xmemory = new double[2 * n];
	X       = new double *[n];
	int j;
	for(j = 0; j < n; j++)
	{	X[j]    = Xmemory + 2 * j;
		X[j][0] = x[j].value();
	}

	// start recording on tape number 1
	Recording = true;
	Domain    = n;
	Range     = 0;
	int tag   = 1;
	int keep  = 1;
	trace_on(tag, keep);

	// declare independent variables
	for(j = 0; j < n; j++)
		x[j] <<= X[j][0];

	return;
}

void Dependent(integer m, doublereal *y)
{
	// check that directly follows call to Independent
	assert( Recording );
	assert( Range == 0 );

	// make sure the set of dependent variables is not empty
	assert( m > 0 );

	// declare the independent variables
	int i;
	double yi;
	for(i = 0; i < m; i++)
		y[i] >>= yi;

	// allocate memory for dependent variable value and differential
	// and store domain value in first column of X
	Ymemory = new double[2 * m];
	Y       = new double *[m];
	for(i = 0; i < m; i++)
		Y[i] = Ymemory + 2 * i;

	// stop the recording
	Range     = m;
	Recording = false;
	trace_off();

}

// linked here from Partial<doublereal>(integer i, integer j)
double Partial(const doublereal &zero, integer i, integer j)
{
	// check that function is defined
	assert( ! Recording );
	assert( Domain > 0 );

	// check dimension of arguments
	int m = Range;
	int n = Domain;
	assert( 0 <= i );
	assert( 0 <= j );
	assert( i < m  );
	assert( j < n );
	assert( zero == doublereal(0) );

	// move data into structure expected by Adolc
	int k;
	for(k = 0; k < Domain; k++)
		X[k][1] = 0.;
	X[j][1] = 1.;

	// compute dy
	int tag  = 1; // tape identifier
	int keep = 1; // keep data in tape
	int d    = 1; // highest derivative degree
	forward(tag, m, n, d, keep, X, Y);

	// return the result
	return Y[i][1];
}

double Value(const doublereal &z)
{	doublereal tmp = z;
	return tmp.value(); 
}

doublereal Sign(const doublereal* x, const doublereal* y)
{	if( *y >= 0. )
		return abs(*x);
	return -abs(*x);
}

} // end f2cad namespace