lib/glibc
1
0
Fork 0
mirror of https://sourceware.org/git/glibc.git synced 2025-12-06 12:01:08 +03:00
Code Activity

Fix lgammaf spurious underflows (bug 18220).

Browse Source 
The flt-32 implementation of lgammaf produces spurious underflow
exceptions for some large arguments, because of calculations involving
x^-2 multiplied by small constants.  This patch fixes this by
adjusting the threshold for a simpler computation to 2**26 (the error
in the simpler computation is on the order of 0.5 * log (x), for a
result on the order of x * log (x)).

Tested for x86_64 and x86.

	[BZ #18220]
	* sysdeps/ieee754/flt-32/e_lgammaf_r.c (__ieee754_lgammaf_r): Use
	2**26 not 2**58 as threshold for returning x * (log (x) - 1).
	* math/auto-libm-test-in: Add another test of lgamma.
	* math/auto-libm-test-out: Regenerated.
This commit is contained in:
Joseph Myers
2015-05-15 17:21:08 +00:00
parent b2fb252408
commit ff069f024a
5 changed files with 38 additions and 4 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
ChangeLog
View File

@@ -1,3 +1,11 @@
2015-05-15 Joseph Myers <joseph@codesourcery.com>
[BZ #18220]
* sysdeps/ieee754/flt-32/e_lgammaf_r.c (__ieee754_lgammaf_r): Use
2**26 not 2**58 as threshold for returning x * (log (x) - 1).
* math/auto-libm-test-in: Add another test of lgamma.
* math/auto-libm-test-out: Regenerated.
2015-05-15 Wilco Dijkstra <wdijkstr@arm.com>
* stdio-common/printf_fp.c (___printf_fp): Use abs.

2
NEWS
View File

@@ -17,7 +17,7 @@ Version 2.22
17999, 18007, 18019, 18020, 18029, 18030, 18032, 18036, 18038, 18039,
18042, 18043, 18046, 18047, 18068, 18080, 18093, 18100, 18104, 18110,
18111, 18125, 18128, 18138, 18185, 18196, 18197, 18206, 18210, 18211,
18217, 18247, 18287, 18319, 18333, 18346, 18397, 18409.
18217, 18220, 18247, 18287, 18319, 18333, 18346, 18397, 18409.
* Cache information can be queried via sysconf() function on s390 e.g. with
_SC_LEVEL1_ICACHE_SIZE as argument.

1
math/auto-libm-test-in
View File

@@ -1757,6 +1757,7 @@ lgamma 0.5
lgamma -0.5
lgamma 0.7
lgamma 1.2
lgamma 0x3.8p56
lgamma 0x1p-5
lgamma -0x1p-5
lgamma 0x1p-10

25
math/auto-libm-test-out
View File

@@ -139612,6 +139612,31 @@ lgamma 1.2
= lgamma tonearest ldbl-128ibm 0x1.33333333333333333333333333p+0L : -0x1.5db138c7d70c72cb0b57c9e83d8p-4L 1 : inexact-ok
= lgamma towardzero ldbl-128ibm 0x1.33333333333333333333333333p+0L : -0x1.5db138c7d70c72cb0b57c9e83dp-4L 1 : inexact-ok
= lgamma upward ldbl-128ibm 0x1.33333333333333333333333333p+0L : -0x1.5db138c7d70c72cb0b57c9e83dp-4L 1 : inexact-ok
lgamma 0x3.8p56
= lgamma downward flt-32 0x3.8p+56f : 0x8.8bdd4p+60f 1 : inexact-ok
= lgamma tonearest flt-32 0x3.8p+56f : 0x8.8bdd4p+60f 1 : inexact-ok
= lgamma towardzero flt-32 0x3.8p+56f : 0x8.8bdd4p+60f 1 : inexact-ok
= lgamma upward flt-32 0x3.8p+56f : 0x8.8bdd5p+60f 1 : inexact-ok
= lgamma downward dbl-64 0x3.8p+56 : 0x8.8bdd41bf4484p+60 1 : inexact-ok
= lgamma tonearest dbl-64 0x3.8p+56 : 0x8.8bdd41bf44848p+60 1 : inexact-ok
= lgamma towardzero dbl-64 0x3.8p+56 : 0x8.8bdd41bf4484p+60 1 : inexact-ok
= lgamma upward dbl-64 0x3.8p+56 : 0x8.8bdd41bf44848p+60 1 : inexact-ok
= lgamma downward ldbl-96-intel 0x3.8p+56L : 0x8.8bdd41bf4484605p+60L 1 : inexact-ok
= lgamma tonearest ldbl-96-intel 0x3.8p+56L : 0x8.8bdd41bf4484605p+60L 1 : inexact-ok
= lgamma towardzero ldbl-96-intel 0x3.8p+56L : 0x8.8bdd41bf4484605p+60L 1 : inexact-ok
= lgamma upward ldbl-96-intel 0x3.8p+56L : 0x8.8bdd41bf4484606p+60L 1 : inexact-ok
= lgamma downward ldbl-96-m68k 0x3.8p+56L : 0x8.8bdd41bf4484605p+60L 1 : inexact-ok
= lgamma tonearest ldbl-96-m68k 0x3.8p+56L : 0x8.8bdd41bf4484605p+60L 1 : inexact-ok
= lgamma towardzero ldbl-96-m68k 0x3.8p+56L : 0x8.8bdd41bf4484605p+60L 1 : inexact-ok
= lgamma upward ldbl-96-m68k 0x3.8p+56L : 0x8.8bdd41bf4484606p+60L 1 : inexact-ok
= lgamma downward ldbl-128 0x3.8p+56L : 0x8.8bdd41bf44846050819264e2d568p+60L 1 : inexact-ok
= lgamma tonearest ldbl-128 0x3.8p+56L : 0x8.8bdd41bf44846050819264e2d57p+60L 1 : inexact-ok
= lgamma towardzero ldbl-128 0x3.8p+56L : 0x8.8bdd41bf44846050819264e2d568p+60L 1 : inexact-ok
= lgamma upward ldbl-128 0x3.8p+56L : 0x8.8bdd41bf44846050819264e2d57p+60L 1 : inexact-ok
= lgamma downward ldbl-128ibm 0x3.8p+56L : 0x8.8bdd41bf44846050819264e2d4p+60L 1 : inexact-ok
= lgamma tonearest ldbl-128ibm 0x3.8p+56L : 0x8.8bdd41bf44846050819264e2d4p+60L 1 : inexact-ok
= lgamma towardzero ldbl-128ibm 0x3.8p+56L : 0x8.8bdd41bf44846050819264e2d4p+60L 1 : inexact-ok
= lgamma upward ldbl-128ibm 0x3.8p+56L : 0x8.8bdd41bf44846050819264e2d8p+60L 1 : inexact-ok
lgamma 0x1p-5
= lgamma downward flt-32 0x8p-8f : 0x3.72d02cp+0f 1 : inexact-ok
= lgamma tonearest flt-32 0x8p-8f : 0x3.72d03p+0f 1 : inexact-ok

6
sysdeps/ieee754/flt-32/e_lgammaf_r.c
View File

@@ -219,15 +219,15 @@ __ieee754_lgammaf_r(float x, int *signgamp)
case 3: z *= (y+(float)2.0); /* FALLTHRU */
r += __ieee754_logf(z); break;
}
/* 8.0 <= x < 2**58 */
} else if (ix < 0x5c800000) {
/* 8.0 <= x < 2**26 */
} else if (ix < 0x4c800000) {
t = __ieee754_logf(x);
z = one/x;
y = z*z;
w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
r = (x-half)*(t-one)+w;
} else
/* 2**58 <= x <= inf */
/* 2**26 <= x <= inf */
r = x*(__ieee754_logf(x)-one);
if(hx<0) r = nadj - r;
return r;