/*
 * $Id: fvar_fn1.cpp 494 2012-06-13 20:41:16Z johnoel $
 *
 * Author: David Fournier
 * Copyright (c) 2008-2012 Regents of the University of California 
 */
/**
 * \file
 * Description not yet available.
 */
// file fvar_fn.cpp
// math.h functions involving prevariables
#include "fvar.hpp"

#include <stdio.h>
#include <math.h>


void gradfree(dlink *);

//extern prevariable * FRETURN;
//extern int RETURN_INDEX;
//extern dlist * GRAD_LIST;          //js
//extern grad_stack  * GRAD_STACK1;

/**
 * Description not yet available.
 * \param
 */
prevariable& exp(const prevariable& v1)
    {
      double tmp;
      if (++gradient_structure::RETURN_PTR > gradient_structure::MAX_RETURN) gradient_structure::RETURN_PTR = gradient_structure::MIN_RETURN;
      tmp=::exp(v1.v->x);
      gradient_structure::RETURN_PTR->v->x=tmp;
      gradient_structure::GRAD_STACK1->set_gradient_stack(default_evaluation,&(gradient_structure::RETURN_PTR->v->x),
                                    &(v1.v->x),tmp );
      return(*gradient_structure::RETURN_PTR);
    }

/**
 * Description not yet available.
 * \param
 */
prevariable& atan(const prevariable& v1)
    {
      if (++gradient_structure::RETURN_PTR > gradient_structure::MAX_RETURN) gradient_structure::RETURN_PTR = gradient_structure::MIN_RETURN;
      gradient_structure::RETURN_PTR->v->x= ::atan(v1.v->x);
      gradient_structure::GRAD_STACK1->set_gradient_stack(default_evaluation,&(gradient_structure::RETURN_PTR->v->x),
                                    &(v1.v->x),1./(1.+v1.v->x * v1.v->x) );
      return(*gradient_structure::RETURN_PTR);
    }

/**
 * Description not yet available.
 * \param
 */
prevariable& ldexp(const prevariable& v1, const int& exponent)
    {
      if (++gradient_structure::RETURN_PTR > gradient_structure::MAX_RETURN) gradient_structure::RETURN_PTR = gradient_structure::MIN_RETURN;
      gradient_structure::RETURN_PTR->v->x=::ldexp(v1.v->x, exponent);
      gradient_structure::GRAD_STACK1->set_gradient_stack(default_evaluation,&(gradient_structure::RETURN_PTR->v->x),
                                    &(v1.v->x),pow(2.0,exponent) );
      return(*gradient_structure::RETURN_PTR);
    }

/**
 * Description not yet available.
 * \param
 */
prevariable& sqrt(const prevariable& v1)
    {
      double tmp=v1.v->x;
      if (tmp==0.0) 
      {
        cerr << "Attempting to take the derivative of sqrt(prevariable x)"
         " at x=0\n";
        ad_exit(1);
      }
      tmp=::sqrt(tmp);
      if (++gradient_structure::RETURN_PTR > gradient_structure::MAX_RETURN) gradient_structure::RETURN_PTR = gradient_structure::MIN_RETURN;
      gradient_structure::RETURN_PTR->v->x=tmp;
      gradient_structure::GRAD_STACK1->set_gradient_stack(default_evaluation,&(gradient_structure::RETURN_PTR->v->x),
                                    &(v1.v->x),1./(2.*tmp) );
      return(*gradient_structure::RETURN_PTR);
    }

/**
 * Description not yet available.
 * \param
 */
prevariable& sqr(const prevariable& v1)
    {
      double tmp=v1.v->x;
      if (tmp==0.0) 
      {
        cerr << "Attempting to take the derivative of sqrt(prevariable x)"
         " at x=0\n";
        ad_exit(1);
      }
      tmp=::sqrt(tmp);
      if (++gradient_structure::RETURN_PTR > gradient_structure::MAX_RETURN) gradient_structure::RETURN_PTR = gradient_structure::MIN_RETURN;
      gradient_structure::RETURN_PTR->v->x=tmp;
      gradient_structure::GRAD_STACK1->set_gradient_stack(default_evaluation,&(gradient_structure::RETURN_PTR->v->x),
                                    &(v1.v->x),1./(2.*tmp) );
      return(*gradient_structure::RETURN_PTR);
    }

/**
 * Description not yet available.
 * \param
 */
prevariable& tan(const prevariable& v1)
    {
      double t = ::tan(v1.v->x);
      if (++gradient_structure::RETURN_PTR > gradient_structure::MAX_RETURN) gradient_structure::RETURN_PTR = gradient_structure::MIN_RETURN;
      gradient_structure::RETURN_PTR->v->x= t;
      gradient_structure::GRAD_STACK1->set_gradient_stack(default_evaluation,&(gradient_structure::RETURN_PTR->v->x),
                          &(v1.v->x), 1+t*t);
      return(*gradient_structure::RETURN_PTR);
    }

/**
 * Description not yet available.
 * \param
 */
prevariable& tanh(const prevariable& v1)
    {
      double t = ::tanh(v1.v->x);
      if (++gradient_structure::RETURN_PTR > gradient_structure::MAX_RETURN) gradient_structure::RETURN_PTR = gradient_structure::MIN_RETURN;
      gradient_structure::RETURN_PTR->v->x= t;
      gradient_structure::GRAD_STACK1->set_gradient_stack(default_evaluation,&(gradient_structure::RETURN_PTR->v->x),
                          &(v1.v->x), 1-t*t);
      return(*gradient_structure::RETURN_PTR);
    }

/**
 * Description not yet available.
 * \param
 */
prevariable& acos(const prevariable& v1)
    {
      if (++gradient_structure::RETURN_PTR > gradient_structure::MAX_RETURN) gradient_structure::RETURN_PTR = gradient_structure::MIN_RETURN;
      gradient_structure::RETURN_PTR->v->x=::acos(v1.v->x);
      gradient_structure::GRAD_STACK1->set_gradient_stack(default_evaluation,&(gradient_structure::RETURN_PTR->v->x),
                          &(v1.v->x),-1./::sqrt(1.- v1.v->x * v1.v->x));
      return(*gradient_structure::RETURN_PTR);
    }

/**
 * Description not yet available.
 * \param
 */
prevariable& asin(const prevariable& v1)
    {
      if (++gradient_structure::RETURN_PTR > gradient_structure::MAX_RETURN) gradient_structure::RETURN_PTR = gradient_structure::MIN_RETURN;
      gradient_structure::RETURN_PTR->v->x=::asin(v1.v->x);
      gradient_structure::GRAD_STACK1->set_gradient_stack(default_evaluation,&(gradient_structure::RETURN_PTR->v->x),
                          &(v1.v->x),1./::sqrt(1.- v1.v->x * v1.v->x));
      return(*gradient_structure::RETURN_PTR);
    }

/**
 * Description not yet available.
 * \param
 */
prevariable& pow(const prevariable& v1, const prevariable& v2)
    {
      if (++gradient_structure::RETURN_PTR > gradient_structure::MAX_RETURN) gradient_structure::RETURN_PTR = gradient_structure::MIN_RETURN;
      double x,y;
      x=::pow(v1.v->x,(v2.v->x)-1);
      y=x* v1.v->x;
      gradient_structure::RETURN_PTR->v->x=y;
      gradient_structure::GRAD_STACK1->set_gradient_stack(default_evaluation,&(gradient_structure::RETURN_PTR->v->x),
                            &(v1.v->x), v2.v->x * x  ,&(v2.v->x), 
                                     y * ::log(v1.v->x) );
      return(*gradient_structure::RETURN_PTR);
    }

/**
 * Description not yet available.
 * \param
 */
prevariable& pow( CGNU_DOUBLE u, const prevariable& v1)
    {
      if (++gradient_structure::RETURN_PTR > gradient_structure::MAX_RETURN) gradient_structure::RETURN_PTR = gradient_structure::MIN_RETURN;
      double y;
      y=::pow(u,(v1.v->x));

      gradient_structure::RETURN_PTR->v->x=y;
      gradient_structure::GRAD_STACK1->set_gradient_stack(default_evaluation,&(gradient_structure::RETURN_PTR->v->x),
                       &(v1.v->x), y * ::log(u) );
 
      return(*gradient_structure::RETURN_PTR);
    }

/**
 * Description not yet available.
 * \param
 */
prevariable& sinh(const prevariable& v1)
    {
      if (++gradient_structure::RETURN_PTR > gradient_structure::MAX_RETURN) gradient_structure::RETURN_PTR = gradient_structure::MIN_RETURN;
      gradient_structure::RETURN_PTR->v->x=::sinh(v1.v->x);
      gradient_structure::GRAD_STACK1->set_gradient_stack(default_evaluation,&(gradient_structure::RETURN_PTR->v->x),
                                    &(v1.v->x),::cosh(v1.v->x) );
      return(*gradient_structure::RETURN_PTR);
    }

/**
 * Description not yet available.
 * \param
 */
prevariable& cosh(const prevariable& v1)
    {
      if (++gradient_structure::RETURN_PTR > gradient_structure::MAX_RETURN) gradient_structure::RETURN_PTR = gradient_structure::MIN_RETURN;
      gradient_structure::RETURN_PTR->v->x=::cosh(v1.v->x);
      gradient_structure::GRAD_STACK1->set_gradient_stack(default_evaluation,&(gradient_structure::RETURN_PTR->v->x),
                                    &(v1.v->x),::sinh(v1.v->x) );
      return(*gradient_structure::RETURN_PTR);
    }

/**
 * Description not yet available.
 * \param
 */
prevariable& atan2(const prevariable& v1, const prevariable& v2)
{
  if (value(v1) == 0 && value(v2) == 0)
  {
    cerr << "Error: The ADMB function \"atan2(y, x)\" is undefined for y and x equal zero.\n"; 
    ad_exit(1);
  }
  if (value(v1) == 0 && value(v2) > 0)
  {
    return atan(v1/v2);
  }
  dvariable x = (sqrt(v2 * v2 + v1 * v1) - v2)/v1;
  return atan(x) * 2.0;
}

/**
 * Description not yet available.
 * \param
 */
prevariable& atan2(const prevariable& v1, CGNU_DOUBLE v2)
{
  if (value(v1) == 0 && v2 == 0)
  {
    cerr << "Error: The ADMB function \"atan2(y, x)\" is undefined for y and x equal zero.\n"; 
    ad_exit(1);
  }
  if (value(v1) == 0 && v2 > 0)
  {
    return atan(v1/v2);
  }
  dvariable x = (sqrt(v2 * v2 + v1 * v1) - v2)/v1;
  return atan(x) * 2.0;
}

/**
 * Description not yet available.
 * \param
 */
prevariable& atan2( CGNU_DOUBLE v1, const prevariable& v2)
{
  if (v1 == 0 && value(v2) == 0)
  {
    cerr << "Error: The ADMB function \"atan2(y, x)\" is undefined for y and x equal zero.\n"; 
    ad_exit(1);
  }
  if (v1 == 0 && value(v2) > 0)
  {
    return atan(v1/v2);
  }
  dvariable x = (sqrt(v2 * v2 + v1 * v1) - v2)/v1;
  return atan(x) * 2.0;
}
