From owner-freebsd-numerics@FreeBSD.ORG Sat Sep 13 19:51:54 2014 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 06DFB391 for ; Sat, 13 Sep 2014 19:51:54 +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)) (Client CN "troutmask.apl.washington.edu", Issuer "troutmask.apl.washington.edu" (not verified)) by mx1.freebsd.org (Postfix) with ESMTPS id D78F5184 for ; Sat, 13 Sep 2014 19:51:53 +0000 (UTC) Received: from troutmask.apl.washington.edu (localhost.apl.washington.edu [127.0.0.1]) by troutmask.apl.washington.edu (8.14.9/8.14.9) with ESMTP id s8DJprNR002429 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-GCM-SHA384 bits=256 verify=NO) for ; Sat, 13 Sep 2014 12:51:53 -0700 (PDT) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.14.9/8.14.9/Submit) id s8DJprNV002428 for freebsd-numerics@freebsd.org; Sat, 13 Sep 2014 12:51:53 -0700 (PDT) (envelope-from sgk) Date: Sat, 13 Sep 2014 12:51:53 -0700 From: Steve Kargl To: freebsd-numerics@freebsd.org Subject: Re: lgammal[_r] patch Message-ID: <20140913195153.GB2323@troutmask.apl.washington.edu> References: <20140913193634.GA2323@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20140913193634.GA2323@troutmask.apl.washington.edu> User-Agent: Mutt/1.5.23 (2014-03-12) X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.18-1 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: Sat, 13 Sep 2014 19:51:54 -0000 On Sat, Sep 13, 2014 at 12:36:34PM -0700, Steve Kargl wrote: > Following my .sig is a patch that contains the ld80 and ld128 > implementations of lgamma[_r](). It also contains an update > to lgammaf_r() and lgamma_r(). Specifically, it has > Yes, a follow-up to my own post. The polynomial approximation for the "t" coefficients has been changed from whatever occurs in lgamma_r() to one that actually works in lgamma[fl]_r(). To explain, I tried to match the approximation in lgamma_r() in the other functions in tc-0.23 <= x <= tc+0.27, but my attempts were fraught with instabilities. Specifically, I tried approximations of the form 1. f(x) = (lgamma(x) - lgamma(tc)) / (x - tc)**2 = t(x - tc) 2. f(x) = (lgamma(x) - lgamma(tc)) = (x - tc)**2 * t(x - tc) 3. f(x) = (lgamma(x) - lgamma(tc)) = t(x - tc) For 1., there were stability issue for x -> tc. For 2., the special Remes algorithm I developed would loose an extrema. I speculate, but haven't pursued, that multiple zeros occurring in the rhs at x = tc is the root of the problem. For 3., I had sensitive issues near the endpoints of the domain. So, I adjusted the domain to [tc-0.24, tc+0.28]. This yielded the coefficients found in the code. -- steve