From owner-freebsd-numerics@FreeBSD.ORG Fri Dec 6 01:55:28 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 68E124FC for ; Fri, 6 Dec 2013 01:55:28 +0000 (UTC) Received: from mail110.syd.optusnet.com.au (mail110.syd.optusnet.com.au [211.29.132.97]) by mx1.freebsd.org (Postfix) with ESMTP id 28A211E44 for ; Fri, 6 Dec 2013 01:55:27 +0000 (UTC) Received: from c122-106-156-23.carlnfd1.nsw.optusnet.com.au (c122-106-156-23.carlnfd1.nsw.optusnet.com.au [122.106.156.23]) by mail110.syd.optusnet.com.au (Postfix) with ESMTPS id 62215783839; Fri, 6 Dec 2013 12:55:19 +1100 (EST) Date: Fri, 6 Dec 2013 12:55:16 +1100 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: Bruce Evans Subject: Re: dead code in lgamma_r[f]? In-Reply-To: <20131206102724.W1033@besplex.bde.org> Message-ID: <20131206114426.I1329@besplex.bde.org> References: <20131205173700.GA64575@troutmask.apl.washington.edu> <20131205182324.GA65135@troutmask.apl.washington.edu> <20131205195225.GA65732@troutmask.apl.washington.edu> <20131206102724.W1033@besplex.bde.org> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed X-Optus-CM-Score: 0 X-Optus-CM-Analysis: v=2.1 cv=DstvpgP+ c=1 sm=1 tr=0 a=ebeQFi2P/qHVC0Yw9JDJ4g==:117 a=PO7r1zJSAAAA:8 a=m8NKqXhiB1EA:10 a=kj9zAlcOel0A:10 a=JzwRw_2MAAAA:8 a=ElNGO9NIuDsA:10 a=QbVCLMTJzu0JbE78RuEA:9 a=CjuIK1q_8ugA:10 Cc: freebsd-numerics@FreeBSD.org, Filipe Maia , Steve Kargl 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: Fri, 06 Dec 2013 01:55:28 -0000 On Fri, 6 Dec 2013, Bruce Evans wrote: > On Thu, 5 Dec 2013, Steve Kargl wrote: > ... >> If we again look at the code from __ieee754_lgamma_r(), we see >> that sin_pi() is called if ix < 0x43300000, so by the time we >> arrive at the 'if(ix<0x43300000)' statement we already know that >> the condition is true. > > No, only for negative integers. hx<0 classifies negative values, and > then ix>=0x43300000 classifies numbers that are so large negative that > they must be integers, and the result is a sort of pole error. We > are just filtering out this case, perhaps as an optimization. Oops, sin_pi() is only called for negative integers, so your change seems to be correct. Just add a comment about the limited domain of sin_pi() (it already has one saying that "x is assumed negative". With the limited range, we can improve more things. floor() can be replaced by adding and subtracting 0x1p52 or 0x1.8p52 like we do for trig and exp functions, followed by an adjustment to get floor() if necessary. We can use irint(z) instead of bit magic or (int)z to extract the parity bit of z z == y case. Note that fdlibm can't use (int)z since this can overflow, and overflow is often not benign in practice. irint(z) can overflow too, but it is easy to specify irint() as having benign behaviour. Bruce