Skip site navigation (1)Skip section navigation (2)
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>