From owner-freebsd-numerics@FreeBSD.ORG Thu Dec 5 18:23:25 2013 Return-Path: Delivered-To: freebsd-numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [8.8.178.115]) (using TLSv1 with cipher ADH-AES256-SHA (256/256 bits)) (No client certificate requested) by hub.freebsd.org (Postfix) with ESMTPS id B39195B0 for ; Thu, 5 Dec 2013 18:23:25 +0000 (UTC) Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu [128.95.76.21]) (using TLSv1.2 with cipher ECDHE-RSA-AES256-GCM-SHA384 (256/256 bits)) (No client certificate requested) by mx1.freebsd.org (Postfix) with ESMTPS id 777C3119A for ; Thu, 5 Dec 2013 18:23:25 +0000 (UTC) Received: from troutmask.apl.washington.edu (localhost.apl.washington.edu [127.0.0.1]) by troutmask.apl.washington.edu (8.14.7/8.14.7) with ESMTP id rB5INOwT065173 for ; Thu, 5 Dec 2013 10:23:24 -0800 (PST) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.14.7/8.14.7/Submit) id rB5INOmH065172 for freebsd-numerics@freebsd.org; Thu, 5 Dec 2013 10:23:24 -0800 (PST) (envelope-from sgk) Date: Thu, 5 Dec 2013 10:23:24 -0800 From: Steve Kargl To: freebsd-numerics@freebsd.org Subject: Re: dead code in lgamma_r[f]? Message-ID: <20131205182324.GA65135@troutmask.apl.washington.edu> References: <20131205173700.GA64575@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20131205173700.GA64575@troutmask.apl.washington.edu> User-Agent: Mutt/1.5.21 (2010-09-15) X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.17 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Thu, 05 Dec 2013 18:23:25 -0000 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;