X-Git-Url: http://erislabs.net/gitweb/?a=blobdiff_plain;f=lib%2Fexpl.c;h=8dcc963685a87f7a62dc6a68c44005a3efac2693;hb=bbfcd2f1a92c9bdbb8d7d7d0a8a8c6665c316747;hp=c70f445aeb2f51af0a20c0cca989f4039c871107;hpb=a4d796fb141dac5d85328872e2fefbd5c44870e1;p=gnulib.git diff --git a/lib/expl.c b/lib/expl.c index c70f445ae..8dcc96368 100644 --- a/lib/expl.c +++ b/lib/expl.c @@ -1,9 +1,5 @@ -/* Emulation for expl. - Contributed by Paolo Bonzini - - Copyright 2002-2003, 2007, 2009-2012 Free Software Foundation, Inc. - - This file is part of gnulib. +/* Exponential function. + Copyright (C) 2011-2013 Free Software Foundation, Inc. This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -35,113 +31,120 @@ expl (long double x) # include -static const long double C[] = { -/* Chebyshev polynomial coefficients for (exp(x)-1)/x */ -# define P1 C[0] -# define P2 C[1] -# define P3 C[2] -# define P4 C[3] -# define P5 C[4] -# define P6 C[5] - 0.5L, - 1.66666666666666666666666666666666683E-01L, - 4.16666666666666666666654902320001674E-02L, - 8.33333333333333333333314659767198461E-03L, - 1.38888888889899438565058018857254025E-03L, - 1.98412698413981650382436541785404286E-04L, - -/* Smallest integer x for which e^x overflows. */ -# define himark C[6] - 11356.523406294143949491931077970765L, - -/* Largest integer x for which e^x underflows. */ -# define lomark C[7] --11433.4627433362978788372438434526231L, - -/* very small number */ -# define TINY C[8] - 1.0e-4900L, - -/* 2^16383 */ -# define TWO16383 C[9] - 5.94865747678615882542879663314003565E+4931L}; +/* gl_expl_table[i] = exp((i - 128) * log(2)/256). */ +extern const long double gl_expl_table[257]; + +/* A value slightly larger than log(2). */ +#define LOG2_PLUS_EPSILON 0.6931471805599454L + +/* Best possible approximation of log(2) as a 'long double'. */ +#define LOG2 0.693147180559945309417232121458176568075L + +/* Best possible approximation of 1/log(2) as a 'long double'. */ +#define LOG2_INVERSE 1.44269504088896340735992468100189213743L + +/* Best possible approximation of log(2)/256 as a 'long double'. */ +#define LOG2_BY_256 0.00270760617406228636491106297444600221904L + +/* Best possible approximation of 256/log(2) as a 'long double'. */ +#define LOG2_BY_256_INVERSE 369.329930467574632284140718336484387181L + +/* The upper 32 bits of log(2)/256. */ +#define LOG2_BY_256_HI_PART 0.0027076061733168899081647396087646484375L +/* log(2)/256 - LOG2_HI_PART. */ +#define LOG2_BY_256_LO_PART \ + 0.000000000000745396456746323365681353781544922399845L long double expl (long double x) { - /* Check for usual case. */ - if (x < himark && x > lomark) - { - int exponent; - long double t, x22; - int k = 1; - long double result = 1.0; - - /* Compute an integer power of e with a granularity of 0.125. */ - exponent = (int) floorl (x * 8.0L); - x -= exponent / 8.0L; - - if (x > 0.0625) - { - exponent++; - x -= 0.125L; - } - - if (exponent < 0) - { - t = 0.8824969025845954028648921432290507362220L; /* e^-0.25 */ - exponent = -exponent; - } - else - t = 1.1331484530668263168290072278117938725655L; /* e^0.25 */ - - while (exponent) - { - if (exponent & k) - { - result *= t; - exponent ^= k; - } - t *= t; - k <<= 1; - } - - /* Approximate (e^x - 1)/x, using a seventh-degree polynomial, - with maximum error in [-2^-16-2^-53,2^-16+2^-53] - less than 4.8e-39. */ - x22 = x + x*x*(P1+x*(P2+x*(P3+x*(P4+x*(P5+x*P6))))); - - return result + result * x22; - } - /* Exceptional cases: */ - else if (x < himark) - { - if (x + x == x) - /* e^-inf == 0, with no error. */ - return 0; - else - /* Underflow */ - return TINY * TINY; - } - else - /* Return x, if x is a NaN or Inf; or overflow, otherwise. */ - return TWO16383*x; + if (isnanl (x)) + return x; + + if (x >= (long double) LDBL_MAX_EXP * LOG2_PLUS_EPSILON) + /* x > LDBL_MAX_EXP * log(2) + hence exp(x) > 2^LDBL_MAX_EXP, overflows to Infinity. */ + return HUGE_VALL; + + if (x <= (long double) (LDBL_MIN_EXP - 1 - LDBL_MANT_DIG) * LOG2_PLUS_EPSILON) + /* x < (LDBL_MIN_EXP - 1 - LDBL_MANT_DIG) * log(2) + hence exp(x) < 2^(LDBL_MIN_EXP-1-LDBL_MANT_DIG), + underflows to zero. */ + return 0.0L; + + /* Decompose x into + x = n * log(2) + m * log(2)/256 + y + where + n is an integer, + m is an integer, -128 <= m <= 128, + y is a number, |y| <= log(2)/512 + epsilon = 0.00135... + Then + exp(x) = 2^n * exp(m * log(2)/256) * exp(y) + The first factor is an ldexpl() call. + The second factor is a table lookup. + The third factor is computed + - either as sinh(y) + cosh(y) + where sinh(y) is computed through the power series: + sinh(y) = y + y^3/3! + y^5/5! + ... + and cosh(y) is computed as hypot(1, sinh(y)), + - or as exp(2*z) = (1 + tanh(z)) / (1 - tanh(z)) + where z = y/2 + and tanh(z) is computed through its power series: + tanh(z) = z + - 1/3 * z^3 + + 2/15 * z^5 + - 17/315 * z^7 + + 62/2835 * z^9 + - 1382/155925 * z^11 + + 21844/6081075 * z^13 + - 929569/638512875 * z^15 + + ... + Since |z| <= log(2)/1024 < 0.0007, the relative contribution of the + z^13 term is < 0.0007^12 < 2^-120 <= 2^-LDBL_MANT_DIG, therefore we + can truncate the series after the z^11 term. + + Given the usual bounds LDBL_MAX_EXP <= 16384, LDBL_MIN_EXP >= -16381, + LDBL_MANT_DIG <= 120, we can estimate x: -11440 <= x <= 11357. + This means, when dividing x by log(2), where we want x mod log(2) + to be precise to LDBL_MANT_DIG bits, we have to use an approximation + to log(2) that has 14+LDBL_MANT_DIG bits. */ + + { + long double nm = roundl (x * LOG2_BY_256_INVERSE); /* = 256 * n + m */ + /* n has at most 15 bits, nm therefore has at most 23 bits, therefore + n * LOG2_HI_PART is computed exactly, and n * LOG2_LO_PART is computed + with an absolute error < 2^15 * 2e-10 * 2^-LDBL_MANT_DIG. */ + long double y_tmp = x - nm * LOG2_BY_256_HI_PART; + long double y = y_tmp - nm * LOG2_BY_256_LO_PART; + long double z = 0.5L * y; + +/* Coefficients of the power series for tanh(z). */ +#define TANH_COEFF_1 1.0L +#define TANH_COEFF_3 -0.333333333333333333333333333333333333334L +#define TANH_COEFF_5 0.133333333333333333333333333333333333334L +#define TANH_COEFF_7 -0.053968253968253968253968253968253968254L +#define TANH_COEFF_9 0.0218694885361552028218694885361552028218L +#define TANH_COEFF_11 -0.00886323552990219656886323552990219656886L +#define TANH_COEFF_13 0.00359212803657248101692546136990581435026L +#define TANH_COEFF_15 -0.00145583438705131826824948518070211191904L + + long double z2 = z * z; + long double tanh_z = + (((((TANH_COEFF_11 + * z2 + TANH_COEFF_9) + * z2 + TANH_COEFF_7) + * z2 + TANH_COEFF_5) + * z2 + TANH_COEFF_3) + * z2 + TANH_COEFF_1) + * z; + + long double exp_y = (1.0L + tanh_z) / (1.0L - tanh_z); + + int n = (int) roundl (nm * (1.0L / 256.0L)); + int m = (int) nm - 256 * n; + + return ldexpl (gl_expl_table[128 + m] * exp_y, n); + } } #endif - -#if 0 -int -main (void) -{ - printf ("%.16Lg\n", expl (1.0L)); - printf ("%.16Lg\n", expl (-1.0L)); - printf ("%.16Lg\n", expl (2.0L)); - printf ("%.16Lg\n", expl (4.0L)); - printf ("%.16Lg\n", expl (-2.0L)); - printf ("%.16Lg\n", expl (-4.0L)); - printf ("%.16Lg\n", expl (0.0625L)); - printf ("%.16Lg\n", expl (0.3L)); - printf ("%.16Lg\n", expl (0.6L)); -} -#endif