Skip site navigation (1)Skip section navigation (2)
Date:      Thu, 5 Dec 2013 09:37:00 -0800
From:      Steve Kargl <sgk@troutmask.apl.washington.edu>
To:        freebsd-numerics@freebsd.org
Subject:   dead code in lgamma_r[f]?
Message-ID:  <20131205173700.GA64575@troutmask.apl.washington.edu>

Next in thread | Raw E-Mail | Index | Archive | Help
In reviewing the implementation details for lgamma_r[f],
I noticed that in src/e_lgamma.c (and similar code in
e_lgammaf.c):

	static double sin_pi(double x)
{
...
	} else {
            if(ix>=0x43400000) {
                y = zero; n = 0;                 /* y must be even */
            } else {
                if(ix<0x43300000) z = y+two52;	/* exact */
...
}

double
__ieee754_lgamma_r(double x, int *signgamp)
{
...
	if(hx<0) {
	    if(ix>=0x43300000) 	/* |x|>=2**52, must be -integer */
		return one/zero;
	    t = sin_pi(x);
...
}

sin_pi is only called as shown above.  It is appears that
the 'if(ix>=0x43400000)' block is dead code.  It also appears
that (if I'm reading the code correctly) the 'if(ix<0x43300000)'
is always true.  Thus, if seems that following patch (note
cut-n-paste whitespace munging) should be appropriate:

Index: src/e_lgamma_r.c
===================================================================
--- src/e_lgamma_r.c    (revision 1427)
+++ src/e_lgamma_r.c    (working copy)
@@ -177,15 +177,11 @@
            y   = 2.0*(y - floor(y));           /* y = |x| mod 2.0 */
            n   = (int) (y*4.0);
        } 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;

Am I mssing something here?

-- 
Steve



Want to link to this message? Use this URL: <http://docs.FreeBSD.org/cgi/mid.cgi?20131205173700.GA64575>