From nobody Sat Dec 04 19:40:56 2021 X-Original-To: freebsd-current@mlmmj.nyi.freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2610:1c1:1:606c::19:1]) by mlmmj.nyi.freebsd.org (Postfix) with ESMTP id 2895A8D4384; Sat, 4 Dec 2021 19:41:11 +0000 (UTC) (envelope-from hps@selasky.org) Received: from mail.turbocat.net (turbocat.net [IPv6:2a01:4f8:c17:6c4b::2]) (using TLSv1.3 with cipher TLS_AES_256_GCM_SHA384 (256/256 bits) key-exchange X25519 server-signature RSA-PSS (4096 bits) server-digest SHA256) (Client did not present a certificate) by mx1.freebsd.org (Postfix) with ESMTPS id 4J60Rq09wtz3Gtl; Sat, 4 Dec 2021 19:41:10 +0000 (UTC) (envelope-from hps@selasky.org) Received: from [10.36.2.165] (unknown [178.17.145.105]) (using TLSv1.3 with cipher TLS_AES_128_GCM_SHA256 (128/128 bits) key-exchange X25519 server-signature RSA-PSS (2048 bits) server-digest SHA256) (No client certificate requested) by mail.turbocat.net (Postfix) with ESMTPSA id A502D2601B4; Sat, 4 Dec 2021 20:41:02 +0100 (CET) Message-ID: <46162af1-b63f-be10-ebdf-cd328dcfb6e2@selasky.org> Date: Sat, 4 Dec 2021 20:40:56 +0100 List-Id: Discussions about the use of FreeBSD-current List-Archive: https://lists.freebsd.org/archives/freebsd-current List-Help: List-Post: List-Subscribe: List-Unsubscribe: Sender: owner-freebsd-current@freebsd.org MIME-Version: 1.0 User-Agent: Mozilla/5.0 (X11; FreeBSD amd64; rv:91.0) Gecko/20100101 Thunderbird/91.3.0 Subject: Re: What to do about tgammal? Content-Language: en-US To: Steve Kargl , freebsd-hackers@freebsd.org, freebsd-current@freebsd.org References: <20211204185352.GA20452@troutmask.apl.washington.edu> From: Hans Petter Selasky In-Reply-To: <20211204185352.GA20452@troutmask.apl.washington.edu> Content-Type: text/plain; charset=UTF-8; format=flowed Content-Transfer-Encoding: 7bit X-Rspamd-Queue-Id: 4J60Rq09wtz3Gtl X-Spamd-Bar: ---- Authentication-Results: mx1.freebsd.org; none X-Spamd-Result: default: False [-4.00 / 15.00]; REPLY(-4.00)[] X-Spam: Yes X-ThisMailContainsUnwantedMimeParts: N On 12/4/21 19:53, Steve Kargl wrote: > What to do about tgammal? > > A long time ago (2013-09-06), theraven@ committed a kludge that mapped > several missing long double math functions to double math functions > (e.g., tanhl(x) was mapped to tanh(x)). Over the next few years, I > (along with bde and das reviews) provided Intel 80-bit (ld80) and IEEE > 128-bit (ld128) implementations for some of these functions; namely, > coshl(x), sinhl(x), tanhl(x), erfl(x), erfcl(x), and lgamma(x). The > last remaining function is tgammal(x). If one links a program that uses > tgammal(x) with libm, one sees > > /usr/local/bin/ld: fcn_list.o: in function `build_fcn_list': > fcn_list.c:(.text+0x7c4): warning: tgammal has lower than advertised > precision > > The warning is actually misleading. Not only does tgammal(x) have a > *MUCH* lower precision, it has a reduced domain. That is, tgammal(x) > produces +inf for x > 172 whereas tgammal(x) should produce a finite > result for values of x up to 1755 (or so). On amd64-*-freebsd, > testing 1000000 in the below intervals demonstrates pathetic accuracy. > > Current implmentation via imprecise.c > > Interval | Max ULP > -------------------+------------ > [6,171] | 1340542.2 > [1.0662,6] | 14293.3 > [1.01e-17,1.0661] | 3116.1 > [-1.9999,-1.0001] | 15330369.3 > -------------------+------------ > > Well, I finally have gotten around to removing theraven@'s last kludge > for FreeBSD on systems that support ld80. This is done with a straight > forward modification of the msun/bsdsrc code. The limitation on > domain is removed and the accuracy substantially improved. > > Interval | Max ULP > -------------------+---------- > [6,1755] | 8.457 > [1.0662,6] | 11.710 > [1.01e-17,1.0661] | 11.689 > [-1.9999,-1.0001] | 11.871 > -------------------+---------- > > My modifications leverage the fact that tgamma(x) (ie., double function) > uses extend arithmetic to do the computations (approximately 85 bits of > precision). To get the Max ULP below 1 (the desired upper limit), a few > minimax polynomials need to be determined and the mystery around a few > magic numbers need to be unraveled. > > Extending what I have done to an ld128 implementation requires much > more effort than I have time and energy to pursue. Someone with > interest in floating point math on ld128 system can provide an > implementation. > > So, is anyone interested in seeing a massive patch? > Hi, Do you need a implementation of tgamma() which is 100% correct, or a so-called speed-hack version of tgamma() which is almost correct? I've looked a bit into libm in FreeBSD and I see some functions are implemented so that they execute quickly, instead of producing exact results. Is this true? --HPS