/*========================================================================
//
//      mathfp/s_exp.c
//
//      Part of the standard mathematical function library
//
//========================================================================
// ####ECOSGPLCOPYRIGHTBEGIN####                                            
// -------------------------------------------                              
// This file is part of eCos, the Embedded Configurable Operating System.   
// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2004 Free Software Foundation, Inc.
// Copyright (C) 2004 eCosCentric Limited                                   
//
// eCos is free software; you can redistribute it and/or modify it under    
// the terms of the GNU General Public License as published by the Free     
// Software Foundation; either version 2 or (at your option) any later      
// version.                                                                 
//
// eCos is distributed in the hope that it will be useful, but WITHOUT      
// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or    
// FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License    
// for more details.                                                        
//
// You should have received a copy of the GNU General Public License        
// along with eCos; if not, write to the Free Software Foundation, Inc.,    
// 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.            
//
// As a special exception, if other files instantiate templates or use      
// macros or inline functions from this file, or you compile this file      
// and link it with other works to produce a work based on this file,       
// this file does not by itself cause the resulting work to be covered by   
// the GNU General Public License. However the source code for this file    
// must still be made available in accordance with section (3) of the GNU   
// General Public License v2.                                               
//
// This exception does not invalidate any other reasons why a work based    
// on this file might be covered by the GNU General Public License.         
// -------------------------------------------                              
// ####ECOSGPLCOPYRIGHTEND####                                              
//========================================================================
//#####DESCRIPTIONBEGIN####
//
// Author(s):     jlarmour
// Contributors:  
// Date:          2004-02-30
// Purpose:       
// Description:   Derived from newlib libm
//                http://sources.redhat.com/newlib
//                in turn based on SunSoft fdlibm.
//                See copyright/license information below.
// Usage:         #include <math.h>
//
//####DESCRIPTIONEND####
//======================================================================*/


/* @(#)z_exp.c 1.0 98/08/13 */
/******************************************************************
 * The following routines are coded directly from the algorithms
 * and coefficients given in "Software Manual for the Elementary
 * Functions" by William J. Cody, Jr. and William Waite, Prentice
 * Hall, 1980.
 ******************************************************************/

/*
FUNCTION
        <<exp>>, <<expf>>---exponential
INDEX
        exp
INDEX
        expf

ANSI_SYNOPSIS
        #include <math.h>
        double exp(double <[x]>);
        float expf(float <[x]>);

TRAD_SYNOPSIS
        #include <math.h>
        double exp(<[x]>);
        double <[x]>;

        float expf(<[x]>);
        float <[x]>;

DESCRIPTION
        <<exp>> and <<expf>> calculate the exponential of <[x]>, that is,
        @ifnottex
        e raised to the power <[x]> (where e
        @end ifnottex
        @tex
        $e^x$ (where $e$
        @end tex
        is the base of the natural system of logarithms, approximately 2.71828).

RETURNS
        On success, <<exp>> and <<expf>> return the calculated value.
        If the result underflows, the returned value is <<0>>.  If the
        result overflows, the returned value is <<HUGE_VAL>>.  In
        either case, <<errno>> is set to <<ERANGE>>.

PORTABILITY
        <<exp>> is ANSI C.  <<expf>> is an extension.

*/

/*****************************************************************
 * Exponential Function
 *
 * Input:
 *   x - floating point value
 *
 * Output:
 *   e raised to x.
 *
 * Description:
 *   This routine returns e raised to the xth power.
 *
 *****************************************************************/

#include <float.h>
#include "mathincl/fdlibm.h"
#include "zmath.h"

#ifndef _DOUBLE_IS_32BITS

static const double INV_LN2 = 1.4426950408889634074;
static const double LN2 = 0.6931471805599453094172321;
static const double p[] = { 0.25, 0.75753180159422776666e-2,
                     0.31555192765684646356e-4 };
static const double q[] = { 0.5, 0.56817302698551221787e-1,
                     0.63121894374398504557e-3,
                     0.75104028399870046114e-6 };

double
_DEFUN (exp, (double),
        double x)
{
  int N;
  double g, z, R, P, Q;

  switch (numtest (x))
    {
      case NAN:
        errno = EDOM;
        return (x);
      case INF:
        errno = ERANGE;
        if (ispos (x))
          return (z_infinity.d);
        else
          return (0.0);
      case 0:
        return (1.0);
    }

  /* Check for out of bounds. */
  if (x > BIGX || x < SMALLX)
    {
      errno = ERANGE;
      return (x);
    }

  /* Check for a value too small to calculate. */
  if (-z_rooteps < x && x < z_rooteps)
    {
      return (1.0);
    }

  /* Calculate the exponent. */
  if (x < 0.0)
    N = (int) (x * INV_LN2 - 0.5);
  else
    N = (int) (x * INV_LN2 + 0.5);

  /* Construct the mantissa. */
  g = x - N * LN2;
  z = g * g;
  P = g * ((p[2] * z + p[1]) * z + p[0]);
  Q = ((q[3] * z + q[2]) * z + q[1]) * z + q[0];
  R = 0.5 + P / (Q - P);

  /* Return the floating point value. */
  N++;
  return (ldexp (R, N));
}

#endif /* _DOUBLE_IS_32BITS */

// EOF s_exp.c
