/*-*-C++-*------------------------------------------------------------------
 * gaol -- Just Another Interval Library
 *--------------------------------------------------------------------------
 * This file is part of the gaol distribution. Gaol was primarily 
 * developed at the Swiss Federal Institute of Technology, Lausanne, 
 * Switzerland, and is now developed at the Institut de Recherche 
 * en Informatique de Nantes, France.
 *
 * Copyright (c) 2001 Swiss Federal Institute of Technology, Switzerland
 * Copyright (c) 2002 Institut de Recherche en Informatique de Nantes, France
 *--------------------------------------------------------------------------
 * gaol is a software distributed WITHOUT ANY WARRANTY. Read the associated 
 * COPYING file for information.
 *--------------------------------------------------------------------------
 * CVS: $Id: gaol_double_op_precise.h,v 1.4 2006/03/22 08:35:12 goualard Exp $
 * Last modified: Tue Feb 21 16:17:08 2006 on caleb
 * By: Frederic Goualard <Frederic.Goualard@lina.univ-nantes.fr>
 *--------------------------------------------------------------------------*/

/*!
  \file   gaol_double_op_precise.h
  \brief  Operations correctly rounded down and up.

  Accurate version using APmathlib

  \author Frédéric Goualard
  \date   2006-02-21
*/


#ifndef __gaol_double_op_precise_h__
#define __gaol_double_op_precise_h__

#include <cmath>
#include "gaol/gaol_mathlib.h"
#include "gaol/gaol_port.h"
#include "gaol/gaol_config.h"
#include "gaol/gaol_fpu.h"

namespace gaol {

  INLINE double previous_float(double d)
  {
    return nextafter(d,-GAOL_INFINITY);
  }

  INLINE double next_float(double d)
  {
    return nextafter(d,GAOL_INFINITY);
  }

  /*!
    \brief pow() correctly rounded down
    
    Assumes a correctly rounded pow(double) when rounding direction
    is to nearest (such as ensured by the IBM math library).
    \caution The current rounding direction must be to nearest.
  */
  INLINE double pow_dn(double d, unsigned int e)
    {
      round_nearest();
      double f=previous_float(upow(d,e));
      round_upward();
      return f;
    }

  /*!
    \brief pow() correctly rounded up

    Assumes a correctly rounded pow(double) when rounding direction
    is to nearest (such as ensured by the IBM math library).
    \caution The current rounding direction must be to nearest.
  */
   INLINE double pow_up(double d, unsigned int e)
    {
      round_nearest();
      double f=next_float(upow(d,e));
      round_upward();
      return f;
    }

  /*!
    \brief nthroot() correctly rounded down
    
    Assumes a correctly rounded pow(double) when rounding direction
    is to nearest (such as ensured by the IBM math library).
    \caution The current rounding direction must be to nearest.
  */
  INLINE double nthroot_dn(double d, double e)
    {
      round_nearest();
      double f=previous_float(upow(d,e));
      round_upward();
      return f;
    }

  /*!
    \brief nthroot() correctly rounded up

    Assumes a correctly rounded pow(double) when rounding direction
    is to nearest (such as ensured by the IBM math library).
    \caution The current rounding direction must be to nearest.
  */
   INLINE double nthroot_up(double d, double e)
    {
      round_nearest();
      double f=next_float(upow(d,e));
      round_upward();
      return f;
    }


  /*!
    \brief Exponential correctly rounded downward

    Assumes a correctly rounded function exp(double) when rounding
    direction is to nearest (such as ensured by the IBM math library).
    \caution The current rounding direction must be to nearest.
  */
  INLINE double exp_dn(double d)
  {
    round_nearest();
    double f=previous_float(uexp(d));
    round_upward();
    return f;
  }

  /*!
    \brief Exponential correctly rounded upward

    Assumes a correctly rounded function exp(double) when rounding
    direction is to nearest (such as ensured by the IBM math library).
    \caution The current rounding direction must be to nearest.
  */
  INLINE double exp_up(double d)
  {
    round_nearest();
    double f=next_float(uexp(d));
    round_upward();
    return f;
  }

  INLINE double cos_dn(double d)
  {
    round_nearest();
    double f=previous_float(ucos(d));
    round_upward();
    return f;
  }

  INLINE double cos_up(double d)
  {
    round_nearest();
    double f=next_float(ucos(d));
    round_upward();
    return f;
  }

  INLINE double tan_dn(double d)
  {
    round_nearest();
    double f=previous_float(utan(d));
    round_upward();
    return f;
  }

  INLINE double tan_up(double d)
  {
    round_nearest();
    double f=next_float(utan(d));
    round_upward();
    return f;
  }

   INLINE double log_dn(double d)
  {
    round_nearest();
    double f=previous_float(ulog(d));
    round_upward();
    return f;
  }

  INLINE double log_up(double d)
  {
    round_nearest();
    double f=next_float(ulog(d));
    round_upward();
    return f;
  }


  INLINE double acos_dn(double d)
  {
    round_nearest();
    double f=previous_float(uacos(d));
    round_upward();
    return f;
  }

  INLINE double acos_up(double d)
  {
    round_nearest();
    double f=next_float(uacos(d));
    round_upward();
    return f;
  }

  INLINE double asin_dn(double d)
  {
    round_nearest();
    double f=previous_float(uasin(d));
    round_upward();
    return f;
  }

  INLINE double asin_up(double d)
  {
    round_nearest();
    double f=next_float(uasin(d));
    round_upward();
    return f;
  }

  INLINE double atan_dn(double d)
  {
    round_nearest();
    double f=previous_float(uatan(d));
    round_upward();
    return f;
  }

  INLINE double atan_up(double d)
  {
    round_nearest();
    double f=next_float(uatan(d));
    round_upward();
    return f;
  }

} // namespace gaol

#endif /* __gaol_double_op_precise_h__ */
