This is the mail archive of the libc-hacker@sources.redhat.com mailing list for the glibc project.
Note that libc-hacker is a closed list. You may look at the archives of this list, but subscription and posting are not open.
Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
---|---|---|
Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |
I've committed the following patch, Andreas 2001-05-15 Andreas Jaeger <aj@suse.de> * sysdeps/ieee754/ldbl-128/s_expm1l.c: New file, contributed by Stephen L Moshier <moshier@mediaone.net>. * sysdeps/i386/fpu/libm-test-ulps: Adjust for change. * math/libm-test.inc: Add comment with ToDo. * sysdeps/i386/fpu/e_expl.c: Rewritten to C and using a more accurate algorithm. Patch by Stephen L Moshier <moshier@mediaone.net>. * sysdeps/i386/fpu/e_expl.S: Removed. ============================================================ Index: sysdeps/i386/fpu/libm-test-ulps --- sysdeps/i386/fpu/libm-test-ulps 2001/05/14 09:24:51 1.28 +++ sysdeps/i386/fpu/libm-test-ulps 2001/05/15 07:54:47 @@ -484,13 +484,13 @@ # ctan Test "Real part of: ctan (-2 - 3 i) == 0.0037640256415042482 - 1.0032386273536098014 i": -ildouble: 437 -ldouble: 437 +ildouble: 439 +ldouble: 439 Test "Imaginary part of: ctan (-2 - 3 i) == 0.0037640256415042482 - 1.0032386273536098014 i": float: 1 ifloat: 1 -ildouble: 1 -ldouble: 1 +ildouble: 2 +ldouble: 2 Test "Real part of: ctan (0.7 + 1.2 i) == 0.1720734197630349001 + 0.9544807059989405538 i": double: 1 float: 1 @@ -506,13 +506,13 @@ # ctanh Test "Real part of: ctanh (-2 - 3 i) == -0.9653858790221331242 + 0.0098843750383224937 i": -ildouble: 2 -ldouble: 2 +ildouble: 5 +ldouble: 5 Test "Imaginary part of: ctanh (-2 - 3 i) == -0.9653858790221331242 + 0.0098843750383224937 i": float: 1 ifloat: 1 -ildouble: 23 -ldouble: 23 +ildouble: 25 +ldouble: 25 Test "Real part of: ctanh (0 + pi/4 i) == 0.0 + 1.0 i": Test "Imaginary part of: ctanh (0 + pi/4 i) == 0.0 + 1.0 i": float: 1 @@ -551,8 +551,8 @@ float: 12 idouble: 24 ifloat: 12 -ldouble: 4 -ildouble: 4 +ldouble: 12 +ildouble: 12 Test "erfc (9) == 0.41370317465138102381e-36": ldouble: 36 ildouble: 36 ============================================================ Index: math/libm-test.inc --- math/libm-test.inc 2001/05/14 08:14:51 1.33 +++ math/libm-test.inc 2001/05/15 07:54:47 @@ -104,7 +104,12 @@ - Compiler has errors With e.g. gcc 2.7.2.2 the test for cexp fails because of a compiler error. - */ + + + To Do: All parameter should be numbers that can be represented as + exact floating point values. Currently some values cannot be represented + exactly and therefore the result is not the expected result. +*/ #ifndef _GNU_SOURCE # define _GNU_SOURCE ============================================================ Index: sysdeps/ieee754/ldbl-128/s_expm1l.c --- sysdeps/ieee754/ldbl-128/s_expm1l.c created +++ sysdeps/ieee754/ldbl-128/s_expm1l.c Tue May 15 09:50:36 2001 1.1 @@ -0,0 +1,145 @@ +/* expm1l.c + * + * Exponential function, minus 1 + * 128-bit long double precision + * + * + * + * SYNOPSIS: + * + * long double x, y, expm1l(); + * + * y = expm1l( x ); + * + * + * + * DESCRIPTION: + * + * Returns e (2.71828...) raised to the x power, minus one. + * + * Range reduction is accomplished by separating the argument + * into an integer k and fraction f such that + * + * x k f + * e = 2 e. + * + * An expansion x + .5 x^2 + x^3 R(x) approximates exp(f) - 1 + * in the basic range [-0.5 ln 2, 0.5 ln 2]. + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -79,+MAXLOG 100,000 1.7e-34 4.5e-35 + * + */ + +/* Copyright 2001 by Stephen L. Moshier */ + + +#include "math.h" +#include "math_private.h" + +/* exp(x) - 1 = x + 0.5 x^2 + x^3 P(x)/Q(x) + -.5 ln 2 < x < .5 ln 2 + Theoretical peak relative error = 8.1e-36 */ + +static long double + P0 = 2.943520915569954073888921213330863757240E8L, + P1 = -5.722847283900608941516165725053359168840E7L, + P2 = 8.944630806357575461578107295909719817253E6L, + P3 = -7.212432713558031519943281748462837065308E5L, + P4 = 4.578962475841642634225390068461943438441E4L, + P5 = -1.716772506388927649032068540558788106762E3L, + P6 = 4.401308817383362136048032038528753151144E1L, + P7 = -4.888737542888633647784737721812546636240E-1L, + Q0 = 1.766112549341972444333352727998584753865E9L, + Q1 = -7.848989743695296475743081255027098295771E8L, + Q2 = 1.615869009634292424463780387327037251069E8L, + Q3 = -2.019684072836541751428967854947019415698E7L, + Q4 = 1.682912729190313538934190635536631941751E6L, + Q5 = -9.615511549171441430850103489315371768998E4L, + Q6 = 3.697714952261803935521187272204485251835E3L, + Q7 = -8.802340681794263968892934703309274564037E1L, + /* Q8 = 1.000000000000000000000000000000000000000E0 */ +/* C1 + C2 = ln 2 */ + + C1 = 6.93145751953125E-1L, + C2 = 1.428606820309417232121458176568075500134E-6L, +/* ln (2^16384 * (1 - 2^-113)) */ + maxlog = 1.1356523406294143949491931077970764891253E4L, +/* ln 2^-114 */ + minarg = -7.9018778583833765273564461846232128760607E1L, big = 2e4932L; + + +long double +__expm1l (long double x) +{ + long double px, qx, xx; + int32_t ix, sign; + ieee854_long_double_shape_type u; + int k; + + /* Overflow. */ + if (x > maxlog) + return (big * big); + + /* Minimum value. */ + if (x < minarg) + return (4.0 / big - 1.0L); + + /* Detect infinity and NaN. */ + u.value = x; + ix = u.parts32.w0; + sign = ix & 0x80000000; + ix &= 0x7fffffff; + if (ix >= 0x7fff0000) + { + /* Infinity. */ + if (((ix & 0xffff) | u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0) + { + if (sign) + return -1.0L; + else + return x; + } + /* NaN. */ + return (x + x); + } + + /* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */ + xx = C1 + C2; /* ln 2. */ + px = __floorl (0.5 + x / xx); + k = px; + /* remainder times ln 2 */ + x -= px * C1; + x -= px * C2; + + /* Approximate exp(remainder ln 2). */ + px = (((((((P7 * x + + P6) * x + + P5) * x + P4) * x + P3) * x + P2) * x + P1) * x + P0) * x; + + qx = (((((((x + + Q7) * x + + Q6) * x + Q5) * x + Q4) * x + Q3) * x + Q2) * x + Q1) * x + Q0; + + xx = x * x; + qx = x + (0.5 * xx + xx * px / qx); + + /* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2). + + We have qx = exp(remainder ln 2) - 1, so + exp(x) - 1 = 2^k (qx + 1) - 1 + = 2^k qx + 2^k - 1. */ + + px = ldexpl (1.0L, k); + x = px * qx + (px - 1.0); + return x; +} + +weak_alias (__expm1l, expm1l) +#ifdef NO_LONG_DOUBLE +strong_alias (__expm1, __expm1l) weak_alias (__expm1, expm1l) +#endif ============================================================ Index: sysdeps/i386/fpu/e_expl.c --- sysdeps/i386/fpu/e_expl.c created +++ sysdeps/i386/fpu/e_expl.c Tue May 15 09:54:41 2001 1.1 @@ -0,0 +1,75 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>. + */ + +/* + * The 8087 method for the exponential function is to calculate + * exp(x) = 2^(x log2(e)) + * after separating integer and fractional parts + * x log2(e) = i + f, |f| <= .5 + * 2^i is immediate but f needs to be precise for long double accuracy. + * Suppress range reduction error in computing f by the following. + * Separate x into integer and fractional parts + * x = xi + xf, |xf| <= .5 + * Separate log2(e) into the sum of an exact number c0 and small part c1. + * c0 + c1 = log2(e) to extra precision + * Then + * f = (c0 xi - i) + c0 xf + c1 x + * where c0 xi is exact and so also is (c0 xi - i). + * -- moshier@na-net.ornl.gov + */ + +static long double c0 = 1.44268798828125L; +static long double c1 = 7.05260771340735992468e-6L; + +long double +__ieee754_expl (long double x) +{ + long double res, t; + +/* I added the following ugly construct because expl(+-Inf) resulted + in NaN. The ugliness results from the bright minds at Intel. + For the i686 the code can be written better. + -- drepper@cygnus.com. */ + asm ("fxam\n\t" /* Is NaN or +-Inf? */ + "fstsw %%ax\n\t" + "movb $0x45, %%dh\n\t" + "andb %%ah, %%dh\n\t" + "cmpb $0x05, %%dh\n\t" + "je 1f\n\t" /* Is +-Inf, jump. */ + "fldl2e\n\t" /* 1 log2(e) */ + "fmul %%st(1),%%st\n\t" /* 1 x log2(e) */ + "frndint\n\t" /* 1 i */ + "fld %%st(1)\n\t" /* 2 x */ + "frndint\n\t" /* 2 xi */ + "fld %%st(1)\n\t" /* 3 i */ + "fldt c0\n\t" /* 4 c0 */ + "fld %%st(2)\n\t" /* 5 xi */ + "fmul %%st(1),%%st\n\t" /* 5 c0 xi */ + "fsubp %%st,%%st(2)\n\t" /* 4 f = c0 xi - i */ + "fld %%st(4)\n\t" /* 5 x */ + "fsub %%st(3),%%st\n\t" /* 5 xf = x - xi */ + "fmulp %%st,%%st(1)\n\t" /* 4 c0 xf */ + "faddp %%st,%%st(1)\n\t" /* 3 f = f + c0 xf */ + "fldt c1\n\t" /* 4 */ + "fmul %%st(4),%%st\n\t" /* 4 c1 * x */ + "faddp %%st,%%st(1)\n\t" /* 3 f = f + c1 * x */ + "f2xm1\n\t" /* 3 2^(fract(x * log2(e))) - 1 */ + "fld1\n\t" /* 4 1.0 */ + "faddp\n\t" /* 3 2^(fract(x * log2(e))) */ + "fstp %%st(1)\n\t" /* 2 */ + "fscale\n\t" /* 2 scale factor is st(1); e^x */ + "fstp %%st(1)\n\t" /* 1 */ + "fstp %%st(1)\n\t" /* 0 */ + "jmp 2f\n\t" + "1:\ttestl $0x200, %%eax\n\t" /* Test sign. */ + "jz 2f\n\t" /* If positive, jump. */ + "fstp %%st\n\t" + "fldz\n\t" /* Set result to 0. */ + "2:\t\n" + : "=t" (res) : "0" (x) : "ax", "dx"); + return res; +} -- Andreas Jaeger SuSE Labs aj@suse.de private aj@arthur.inka.de http://www.suse.de/~aj
Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
---|---|---|
Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |