[thirdparty/glibc.git] / sysdeps / ieee754 / dbl-64 / e_exp2.c

/* Double-precision floating point 2^x.
   Copyright (C) 1997,1998,2000,2001,2005,2006,2011
   Free Software Foundation, Inc.
   This file is part of the GNU C Library.
   Contributed by Geoffrey Keating <geoffk@ozemail.com.au>

   The GNU C Library is free software; you can redistribute it and/or
   modify it under the terms of the GNU Lesser General Public
   License as published by the Free Software Foundation; either
   version 2.1 of the License, or (at your option) any later version.

   The GNU C Library 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
   Lesser General Public License for more details.

   You should have received a copy of the GNU Lesser General Public
   License along with the GNU C Library; if not, write to the Free
   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
   02111-1307 USA.  */

/* The basic design here is from
   Shmuel Gal and Boris Bachelis, "An Accurate Elementary Mathematical
   Library for the IEEE Floating Point Standard", ACM Trans. Math. Soft.,
   17 (1), March 1991, pp. 26-45.
   It has been slightly modified to compute 2^x instead of e^x.
   */
#include <stdlib.h>
#include <float.h>
#include <ieee754.h>
#include <math.h>
#include <fenv.h>
#include <inttypes.h>
#include <math_private.h>

#include "t_exp2.h"

static const double TWO1023 = 8.988465674311579539e+307;
static const double TWOM1000 = 9.3326361850321887899e-302;

double
__ieee754_exp2 (double x)
{
  static const double himark = (double) DBL_MAX_EXP;
  static const double lomark = (double) (DBL_MIN_EXP - DBL_MANT_DIG - 1);

  /* Check for usual case.  */
  if (__builtin_expect (isless (x, himark), 1))
    {
      /* Exceptional cases:  */
      if (__builtin_expect (! isgreaterequal (x, lomark), 0))
	{
	  if (__isinf (x))
	    /* e^-inf == 0, with no error.  */
	    return 0;
	  else
	    /* Underflow */
	    return TWOM1000 * TWOM1000;
	}

      static const double THREEp42 = 13194139533312.0;
      int tval, unsafe;
      double rx, x22, result;
      union ieee754_double ex2_u, scale_u;
      fenv_t oldenv;

      libc_feholdexcept (&oldenv);
#ifdef FE_TONEAREST
      /* If we don't have this, it's too bad.  */
      libc_fesetround (FE_TONEAREST);
#endif

      /* 1. Argument reduction.
	 Choose integers ex, -256 <= t < 256, and some real
	 -1/1024 <= x1 <= 1024 so that
	 x = ex + t/512 + x1.

	 First, calculate rx = ex + t/512.  */
      rx = x + THREEp42;
      rx -= THREEp42;
      x -= rx;  /* Compute x=x1. */
      /* Compute tval = (ex*512 + t)+256.
	 Now, t = (tval mod 512)-256 and ex=tval/512  [that's mod, NOT %; and
	 /-round-to-nearest not the usual c integer /].  */
      tval = (int) (rx * 512.0 + 256.0);

      /* 2. Adjust for accurate table entry.
	 Find e so that
	 x = ex + t/512 + e + x2
	 where -1e6 < e < 1e6, and
	 (double)(2^(t/512+e))
	 is accurate to one part in 2^-64.  */

      /* 'tval & 511' is the same as 'tval%512' except that it's always
	 positive.
	 Compute x = x2.  */
      x -= exp2_deltatable[tval & 511];

      /* 3. Compute ex2 = 2^(t/512+e+ex).  */
      ex2_u.d = exp2_accuratetable[tval & 511];
      tval >>= 9;
      unsafe = abs(tval) >= -DBL_MIN_EXP - 1;
      ex2_u.ieee.exponent += tval >> unsafe;
      scale_u.d = 1.0;
      scale_u.ieee.exponent += tval - (tval >> unsafe);

      /* 4. Approximate 2^x2 - 1, using a fourth-degree polynomial,
	 with maximum error in [-2^-10-2^-30,2^-10+2^-30]
	 less than 10^-19.  */

      x22 = (((.0096181293647031180
	       * x + .055504110254308625)
	      * x + .240226506959100583)
	     * x + .69314718055994495) * ex2_u.d;
      math_opt_barrier (x22);

      /* 5. Return (2^x2-1) * 2^(t/512+e+ex) + 2^(t/512+e+ex).  */
      libc_fesetenv (&oldenv);

      result = x22 * x + ex2_u.d;

      if (!unsafe)
	return result;
      else
	return result * scale_u.d;
    }
  else
    /* Return x, if x is a NaN or Inf; or overflow, otherwise.  */
    return TWO1023*x;
}
strong_alias (__ieee754_exp2, __exp2_finite)
Commit	Line	Data
63640cb7	1	/* Double-precision floating point 2^x.
0ac5ae23 UD	2	Copyright (C) 1997,1998,2000,2001,2005,2006,2011
0ac5ae23 UD	3	Free Software Foundation, Inc.
63640cb7 UD	4	This file is part of the GNU C Library.
	5	Contributed by Geoffrey Keating <geoffk@ozemail.com.au>
	6
	7	The GNU C Library is free software; you can redistribute it and/or
41bdb6e2 AJ	8	modify it under the terms of the GNU Lesser General Public
	9	License as published by the Free Software Foundation; either
	10	version 2.1 of the License, or (at your option) any later version.
63640cb7 UD	11
	12	The GNU C Library is distributed in the hope that it will be useful,
	13	but WITHOUT ANY WARRANTY; without even the implied warranty of
	14	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
41bdb6e2	15	Lesser General Public License for more details.
63640cb7	16
41bdb6e2 AJ	17	You should have received a copy of the GNU Lesser General Public
	18	License along with the GNU C Library; if not, write to the Free
	19	Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
	20	02111-1307 USA. */
63640cb7 UD	21
	22	/* The basic design here is from
	23	Shmuel Gal and Boris Bachelis, "An Accurate Elementary Mathematical
	24	Library for the IEEE Floating Point Standard", ACM Trans. Math. Soft.,
	25	17 (1), March 1991, pp. 26-45.
	26	It has been slightly modified to compute 2^x instead of e^x.
	27	*/
63640cb7 UD	28	#include <stdlib.h>
	29	#include <float.h>
	30	#include <ieee754.h>
	31	#include <math.h>
	32	#include <fenv.h>
	33	#include <inttypes.h>
	34	#include <math_private.h>
	35
	36	#include "t_exp2.h"
	37
d38f1dba UD	38	static const double TWO1023 = 8.988465674311579539e+307;
d38f1dba UD	39	static const double TWOM1000 = 9.3326361850321887899e-302;
63640cb7 UD	40
	41	double
	42	__ieee754_exp2 (double x)
	43	{
	44	static const double himark = (double) DBL_MAX_EXP;
601d2942	45	static const double lomark = (double) (DBL_MIN_EXP - DBL_MANT_DIG - 1);
63640cb7 UD	46
63640cb7 UD	47	/* Check for usual case. */
99ce7b04	48	if (__builtin_expect (isless (x, himark), 1))
63640cb7	49	{
99ce7b04 UD	50	/* Exceptional cases: */
	51	if (__builtin_expect (! isgreaterequal (x, lomark), 0))
	52	{
	53	if (__isinf (x))
	54	/* e^-inf == 0, with no error. */
	55	return 0;
	56	else
	57	/* Underflow */
	58	return TWOM1000 * TWOM1000;
	59	}
	60
63640cb7 UD	61	static const double THREEp42 = 13194139533312.0;
	62	int tval, unsafe;
	63	double rx, x22, result;
	64	union ieee754_double ex2_u, scale_u;
	65	fenv_t oldenv;
	66
d38f1dba	67	libc_feholdexcept (&oldenv);
63640cb7 UD	68	#ifdef FE_TONEAREST
63640cb7 UD	69	/* If we don't have this, it's too bad. */
d38f1dba	70	libc_fesetround (FE_TONEAREST);
63640cb7 UD	71	#endif
	72
	73	/* 1. Argument reduction.
	74	Choose integers ex, -256 <= t < 256, and some real
	75	-1/1024 <= x1 <= 1024 so that
	76	x = ex + t/512 + x1.
	77
	78	First, calculate rx = ex + t/512. */
	79	rx = x + THREEp42;
	80	rx -= THREEp42;
	81	x -= rx; /* Compute x=x1. */
	82	/* Compute tval = (ex*512 + t)+256.
	83	Now, t = (tval mod 512)-256 and ex=tval/512 [that's mod, NOT %; and
	84	/-round-to-nearest not the usual c integer /]. */
	85	tval = (int) (rx * 512.0 + 256.0);
	86
	87	/* 2. Adjust for accurate table entry.
	88	Find e so that
	89	x = ex + t/512 + e + x2
	90	where -1e6 < e < 1e6, and
	91	(double)(2^(t/512+e))
	92	is accurate to one part in 2^-64. */
	93
	94	/* 'tval & 511' is the same as 'tval%512' except that it's always
	95	positive.
	96	Compute x = x2. */
	97	x -= exp2_deltatable[tval & 511];
	98
	99	/* 3. Compute ex2 = 2^(t/512+e+ex). */
	100	ex2_u.d = exp2_accuratetable[tval & 511];
	101	tval >>= 9;
	102	unsafe = abs(tval) >= -DBL_MIN_EXP - 1;
	103	ex2_u.ieee.exponent += tval >> unsafe;
	104	scale_u.d = 1.0;
	105	scale_u.ieee.exponent += tval - (tval >> unsafe);
	106
	107	/* 4. Approximate 2^x2 - 1, using a fourth-degree polynomial,
	108	with maximum error in [-2^-10-2^-30,2^-10+2^-30]
	109	less than 10^-19. */
	110
	111	x22 = (((.0096181293647031180
	112	* x + .055504110254308625)
	113	* x + .240226506959100583)
	114	* x + .69314718055994495) * ex2_u.d;
d38f1dba	115	math_opt_barrier (x22);
63640cb7 UD	116
63640cb7 UD	117	/* 5. Return (2^x2-1) * 2^(t/512+e+ex) + 2^(t/512+e+ex). */
d38f1dba	118	libc_fesetenv (&oldenv);
63640cb7 UD	119
	120	result = x22 * x + ex2_u.d;
	121
	122	if (!unsafe)
	123	return result;
	124	else
	125	return result * scale_u.d;
	126	}
63640cb7 UD	127	else
	128	/* Return x, if x is a NaN or Inf; or overflow, otherwise. */
	129	return TWO1023*x;
	130	}
0ac5ae23	131	strong_alias (__ieee754_exp2, __exp2_finite)