Date: Thu, 5 Dec 2013 10:23:24 -0800 From: Steve Kargl <sgk@troutmask.apl.washington.edu> To: freebsd-numerics@freebsd.org Subject: Re: dead code in lgamma_r[f]? Message-ID: <20131205182324.GA65135@troutmask.apl.washington.edu> In-Reply-To: <20131205173700.GA64575@troutmask.apl.washington.edu> References: <20131205173700.GA64575@troutmask.apl.washington.edu>
next in thread | previous in thread | raw e-mail | index | archive | help
On Thu, Dec 05, 2013 at 09:37:00AM -0800, Steve Kargl wrote: > In reviewing the implementation details for lgamma_r[f], > I noticed that in src/e_lgamma.c (and similar code in > e_lgammaf.c): These should be e_lgamma_r.c and e_lgammaf_r.c. I've only done some limited testing, but the patch following my sig seems to work. Note the first part of the diff allows us to remove the cast to float in e_lgamma_r.c, which minimizes the diff between e_lgamma_r.c and e_lgammaf_r.c (and e_lgammal_r.c :). -- Steve Index: src/e_lgamma_r.c =================================================================== --- src/e_lgamma_r.c (revision 1427) +++ src/e_lgamma_r.c (working copy) @@ -173,19 +173,15 @@ */ z = floor(y); if(z!=y) { /* inexact anyway */ - y *= 0.5; - y = 2.0*(y - floor(y)); /* y = |x| mod 2.0 */ - n = (int) (y*4.0); + y /= 2; + y = 2*(y - floor(y)); /* y = |x| mod 2.0 */ + n = (int)(y*4); } else { - if(ix>=0x43400000) { - y = zero; n = 0; /* y must be even */ - } else { - if(ix<0x43300000) z = y+two52; /* exact */ - GET_LOW_WORD(n,z); - n &= 1; - y = n; - n<<= 2; - } + z = y+two52; /* exact */ + GET_LOW_WORD(n,z); + n &= 1; + y = n; + n<<= 2; } switch (n) { case 0: y = __kernel_sin(pi*y,zero,0); break; Index: src/e_lgammaf_r.c =================================================================== --- src/e_lgammaf_r.c (revision 1427) +++ src/e_lgammaf_r.c (working copy) @@ -106,19 +106,15 @@ */ z = floorf(y); if(z!=y) { /* inexact anyway */ - y *= (float)0.5; - y = (float)2.0*(y - floorf(y)); /* y = |x| mod 2.0 */ - n = (int) (y*(float)4.0); + y /= 2; + y = 2*(y - floorf(y)); /* y = |x| mod 2.0 */ + n = (int)(y*4); } else { - if(ix>=0x4b800000) { - y = zero; n = 0; /* y must be even */ - } else { - if(ix<0x4b000000) z = y+two23; /* exact */ - GET_FLOAT_WORD(n,z); - n &= 1; - y = n; - n<<= 2; - } + z = y+two23; /* exact */ + GET_FLOAT_WORD(n,z); + n &= 1; + y = n; + n<<= 2; } switch (n) { case 0: y = __kernel_sindf(pi*y); break;
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20131205182324.GA65135>