From 4591769b1b0567cafdf0e56be8fc01c4ae1a85be Mon Sep 17 00:00:00 2001 From: Peter Avalos Date: Tue, 3 Jul 2007 18:51:45 +0000 Subject: [PATCH] From FreeBSD's log: Log: Fixed about 50 million errors of infinity ulps and about 3 million errors of between 1.0 and 1.8509 ulps for lgammaf(x) with x between -2**-21 and -2**-70. As usual, the cutoff for tiny args was not correctly translated to float precision. It was 2**-70 but 2**-21 works. Not as usual, having a too-small threshold was worse than a pessimization. It was just a pessimization for (positive) args between 2**-70 and 2**-21, but for the first ~50 million (negative) args below -2**-70, the general code overflowed and gave a result of infinity instead of correct (finite) results near 70*log(2). For the remaining ~361 million negative args above -2**21, the general code gave almost-acceptable errors (lgamma[f]() is not very accurate in general) but the pessimization was larger than for misclassified tiny positive args. Now the max error for lgammaf(x) with |x| < 2**-21 is 0.7885 ulps, and speed and accuracy are almost the same for positive and negative args in this range. The maximum error overall is still infinity ulps. A cutoff of 2**-70 is probably wastefully small for the double precision case. Smaller cutoffs can be used to reduce the max error to nearly 0.5 ulps for tiny args, but this is useless since the general algrorithm for nearly-tiny args is not nearly that accurate -- it has a max error of about 1 ulp. Obtained-from: FreeBSD --- lib/libm/src/e_lgammaf_r.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/libm/src/e_lgammaf_r.c b/lib/libm/src/e_lgammaf_r.c index dbe06c50d4..8ea90217ad 100644 --- a/lib/libm/src/e_lgammaf_r.c +++ b/lib/libm/src/e_lgammaf_r.c @@ -13,7 +13,7 @@ * ==================================================== * * $NetBSD: e_lgammaf_r.c,v 1.6 2002/05/26 22:01:51 wiz Exp $ - * $DragonFly: src/lib/libm/src/e_lgammaf_r.c,v 1.1 2005/07/26 21:15:20 joerg Exp $ + * $DragonFly: src/lib/libm/src/e_lgammaf_r.c,v 1.2 2007/07/03 18:51:45 pavalos Exp $ */ #include @@ -149,7 +149,7 @@ lgammaf_r(float x, int *signgamp) ix = hx&0x7fffffff; if(ix>=0x7f800000) return x*x; if(ix==0) return one/zero; - if(ix<0x1c800000) { /* |x|<2**-70, return -log(|x|) */ + if(ix<0x35000000) { /* |x|<2**-21, return -log(|x|) */ if(hx<0) { *signgamp = -1; return -logf(-x); -- 2.41.0