From mboxrd@z Thu Jan 1 00:00:00 1970 X-Msuck: nntp://news.gmane.org/gmane.linux.lib.musl.general/782 Path: news.gmane.org!not-for-mail From: Rich Felker Newsgroups: gmane.linux.lib.musl.general Subject: Re: math todo Date: Tue, 1 May 2012 10:44:38 -0400 Message-ID: <20120501144438.GL14673@brightrain.aerifal.cx> References: <20120501000503.GI16237@port70.net> <20120501011921.GK14673@brightrain.aerifal.cx> <20120501091458.GJ16237@port70.net> Reply-To: musl@lists.openwall.com NNTP-Posting-Host: plane.gmane.org Mime-Version: 1.0 Content-Type: text/plain; charset=utf-8 Content-Transfer-Encoding: 8bit X-Trace: dough.gmane.org 1335883281 20313 80.91.229.3 (1 May 2012 14:41:21 GMT) X-Complaints-To: usenet@dough.gmane.org NNTP-Posting-Date: Tue, 1 May 2012 14:41:21 +0000 (UTC) To: musl@lists.openwall.com Original-X-From: musl-return-783-gllmg-musl=m.gmane.org@lists.openwall.com Tue May 01 16:41:21 2012 Return-path: Envelope-to: gllmg-musl@plane.gmane.org Original-Received: from mother.openwall.net ([195.42.179.200]) by plane.gmane.org with smtp (Exim 4.69) (envelope-from ) id 1SPEGW-0001LU-Vj for gllmg-musl@plane.gmane.org; Tue, 01 May 2012 16:41:21 +0200 Original-Received: (qmail 27707 invoked by uid 550); 1 May 2012 14:41:20 -0000 Mailing-List: contact musl-help@lists.openwall.com; run by ezmlm Precedence: bulk List-Post: List-Help: List-Unsubscribe: List-Subscribe: Original-Received: (qmail 27697 invoked from network); 1 May 2012 14:41:19 -0000 Content-Disposition: inline In-Reply-To: <20120501091458.GJ16237@port70.net> User-Agent: Mutt/1.5.21 (2010-09-15) Xref: news.gmane.org gmane.linux.lib.musl.general:782 Archived-At: On Tue, May 01, 2012 at 11:14:58AM +0200, Szabolcs Nagy wrote: > yes, actually > fesetround(FE_UPWARD); > exp(-inf) = 0x1p-1074; // smallest subnormal > instead of 0, which is 1ulp depending on your definition of ulp :) This is incorrect (because the result is exact, the error is ==1ulp, and IEEE requires <1ulp). But for the large finite negative arguments, 0x1p-1074 is the correctly rounded answer in rounds-upward mode. > > BTW under asm, we may also want to switch to the faster acos > > implementation. > ok i'll add it > (for double it is known to work, but i havent tested it for > the long double case) Even if we can't switch ld we could still switch it for the float/double versions and keep our original asm for ld80. > > Until we have a better one, couldn't we just use sqrt() as a first > > approximation then use Newton's method or a binary-search for the > > correct long double result? Or just return sqrt() since the only arch > > that doesn't have sqrtl asm yet is the one where ld==double. > yes that's a good idea Already committed the latter idea. > > > tgamma, tgammaf > > > > Well we have these but they're inaccurate for some inputs. > ah yes i added dummy versions, but they are not usable for serious work Good to know. > > > (long double bessel) > > > nextafterf on ld64 > > > > Eh? > actually it's nexttowardf(float x, long double y) > it's not implemented when long double == double > it cannot be just a simple wrapper around nextafter or nextafterf All of the nextafter stuff looks overly complicated. Shouldn't these functions all just be (essentially): if (x>y) rep_x--; else if (x 99 ulp want: 0 got: 8000000000000000 round: m acos arg0: 0x1p+0 arg1: 0x0p+0 want: 0x0p+0 got: -0x0p+0 > 99 ulp want: 8000000000000000 got: 8000000000000001 round: p atan arg0: -0x1p-1074 arg1: 0x0p+0 want: -0x0p+0 got: -0x1p-1074 > 99 ulp want: 0 got: 1 round: m atan arg0: 0x1p-1074 arg1: 0x0p+0 want: 0x0p+0 got: 0x1p-1074 > 99 ulp want: 0 got: 1 round: z atan arg0: 0x1p-1074 arg1: 0x0p+0 want: 0x0p+0 got: 0x1p-1074 > 99 ulp want: 8000000000000000 got: 8000000000000001 round: z atan arg0: -0x1p-1074 arg1: 0x0p+0 want: -0x0p+0 got: -0x1p-1074 IMO this is 1ulp not >99ulp. For denormals ulp is not x * 0x1p-52 but the actual last place of significance. > 99 ulp want: 0 got: 1 round: p exp arg0: -inf arg1: 0x0p+0 want: 0x0p+0 got: 0x1p-1074 > 99 ulp want: 8000000000000000 got: 8000000000000001 round: p expm1 arg0: -0x1p-1074 arg1: 0x0p+0 want: -0x0p+0 got: -0x1p-1074 > 99 ulp want: 1 got: 0 round: m expm1 arg0: 0x1p-1074 arg1: 0x0p+0 want: 0x1p-1074 got: 0x0p+0 > 99 ulp want: 1 got: 0 round: z expm1 arg0: 0x1p-1074 arg1: 0x0p+0 want: 0x1p-1074 got: 0x0p+0 > 99 ulp want: 8000000000000001 got: 8000000000000000 round: z log1p arg0: -0x1p-1074 arg1: 0x0p+0 want: -0x1p-1074 got: -0x0p+0 > 99 ulp want: 8000000000000000 got: 8000000000000001 round: p sin arg0: -0x1p-1074 arg1: 0x0p+0 want: -0x0p+0 got: -0x1p-1074 > 99 ulp want: 0 got: 1 round: m sin arg0: 0x1p-1074 arg1: 0x0p+0 want: 0x0p+0 got: 0x1p-1074 > 99 ulp want: 0 got: 1 round: z sin arg0: 0x1p-1074 arg1: 0x0p+0 want: 0x0p+0 got: 0x1p-1074 > 99 ulp want: 8000000000000000 got: 8000000000000001 round: z sin arg0: -0x1p-1074 arg1: 0x0p+0 want: -0x0p+0 got: -0x1p-1074 Same for these. > 99 ulp want: 0 got: 8000000000000000 round: m acos arg0: 0x1p+0 arg1: 0x0p+0 want: 0x0p+0 got: -0x0p+0 > 99 ulp want: 8000000000000001 got: 8000000000000000 round: p asin arg0: -0x1p-1074 arg1: 0x0p+0 want: -0x1p-1074 got: -0x0p+0 > 99 ulp want: 8000000000000000 got: 8000000000000001 round: m atan2 arg0: -0x1p-1022 arg1: 0x1.fffffffffffffp+1023 want: -0x0p+0 got: -0x1p-1074 > 99 ulp want: 8000000000000000 got: 8000000000000001 round: m atan2 arg0: -0x1p-1022 arg1: 0x1.fffffffffffffp+1023 want: -0x0p+0 got: -0x1p-1074 > 99 ulp want: 0 got: 1 round: p exp arg0: -inf arg1: 0x0p+0 want: 0x0p+0 got: 0x1p-1074 > 99 ulp want: 0 got: 1 round: p exp arg0: -0x1.75p+9 arg1: 0x0p+0 want: 0x0p+0 got: 0x1p-1074 > 99 ulp want: 0 got: 1 round: p exp arg0: -0x1.c9c8p+13 arg1: 0x0p+0 want: 0x0p+0 got: 0x1p-1074 And these. Of course in some cases the result is still non-conformant for other reasons like being outside the range of the function, or ==1ulp error (when correct answer is exact) rather than <1ulp. > all: 69734 fail: 12787 failbad: 107 failepic: 72 Like your choice of naming. :-) Rich