From owner-freebsd-numerics@FreeBSD.ORG Fri Dec 6 02:08:08 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 645498E5 for ; Fri, 6 Dec 2013 02:08:08 +0000 (UTC) Received: from mail107.syd.optusnet.com.au (mail107.syd.optusnet.com.au [211.29.132.53]) by mx1.freebsd.org (Postfix) with ESMTP id 186DE1F01 for ; Fri, 6 Dec 2013 02:08:07 +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 mail107.syd.optusnet.com.au (Postfix) with ESMTPS id F02A0D4053E; Fri, 6 Dec 2013 13:07:58 +1100 (EST) Date: Fri, 6 Dec 2013 13:07:57 +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: <20131206114426.I1329@besplex.bde.org> Message-ID: <20131206130353.U1549@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> <20131206114426.I1329@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=bpB1Wiqi 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=NULW6ebW8qo2wbMJ_24A: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 02:08:08 -0000 On Fri, 6 Dec 2013, Bruce Evans wrote: > 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. Oops^2. s/easy/hard/. When irint(z) overflows on i386, it gives an overflow exception and a result of INT_MIN instead of benign truncation. To avoid these problems, it would have to extract the low bits like lgamma() does now, but that would slow down the usual case of a severely limited range. Bruce