/* RCS $Id: q_exp2.c,v 1.2 1997/07/15 14:57:04 ae13 Exp $ */
/*
$Log: q_exp2.c,v $
 * Revision 1.2  1997/07/15  14:57:04  ae13
 * *** empty log message ***
 *
 * Revision 1.1  1997/05/16  13:29:06  ae13
 * Initial revision
 *
 */
/****************************************************************/
/*                                                              */
/*      Entries         : a_real q_exp2(a_real x)               */
/*                        a_real x  ;                           */
/*                                                              */
/*      Arguments       : x = real argument of exp2(x)          */
/*                                                              */
/*      Description     : Fast real exponential function        */
/*                        for double IEEE Floating-Point        */
/*                        Arithmetic                            */
/*                                                              */
/*      Date            : 1997-07-07                            */
/*                                                              */
/*      Author          : W. Hofschuster                        */
/*                                                              */
/****************************************************************/

#ifndef ALL_IN_ONE
#ifdef AIX
#include "/u/p88c/runtime/o_defs.h"
#else
#include "o_defs.h"
#endif
#define local
#endif


#include "q_defs.h"
#include "q_fcth.h"

#ifdef LINT_ARGS
local double q_exp2(double x)
#else
local double q_exp2(x)

double x;
#endif
{
 int j;
 long int n,m;
 double r,q,s;
 double res;

#ifdef PXSCTRAP
 E_SPUSH("q_exp2");
#endif
 /* Step 1: Special cases  */
 if NANTEST(x)                                       /* Test: x=NaN */
      res=q_abortnan(INV_ARG,&x,4);
 else {
 if ((-q_ext1<x) && (x<q_ext1))                      /* |x|<2^-54 */
   res=x+1;
 else
  { if (1023<x) 
      res=q_abortr1(OVER_FLOW,&x,4);                 /* Overflow */
    else
      { if (x<-1022)
	  res=0;                                     /* Result: Underflow */ 
	else
	  {
            if (CUTINT(x)==x)                        /* x is integer */
              { res=1.0; POWER2(res,(long int)x); }
            else {                                   /* x is not integer */
    	      /* Step 2 */
	      if (x>0) n=CUTINT((x*32)+0.5);
	      else     n=CUTINT((x*32)-0.5);         /* round (x)      */
	      j=n % 32;                              /* j=n mod 32     */
	      if (j<0) j+=32;                        /* We force j>=0  */
	      m=(n-j)/32;
	      r=x-n*0.03125;

	      /* Step 3 */	  
              q=(((((q_exc[6]*r+q_exc[5])*r+q_exc[4])*r+q_exc[3])*r+q_exc[2])*r 
                            +q_exc[1])*r+q_exc[0];
	      q=r*q;
	   
	      /* Step 4 */
	      s=q_exld[j]+q_extl[j];
	      res=(q_exld[j]+(q_extl[j]+s*q));
	      POWER2(res,m);
            }
	  }
      }
  }
  }
#ifdef PXSCTRAP
  E_SPOPP("q_exp2");
#endif
  return(res);
}

