
Date: Tue, 1 May 2012 10:44:38 0400 From: Rich Felker <dalias@...ifal.cx> To: musl@...ts.openwall.com Subject: Re: math todo On Tue, May 01, 2012 at 11:14:58AM +0200, Szabolcs Nagy wrote: > yes, actually > fesetround(FE_UPWARD); > exp(inf) = 0x1p1074; // 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, 0x1p1074 is the correctly rounded answer in roundsupward 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 binarysearch 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<y) rep_x++; modulo a little handling for sign and ±0? > 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: 0x1p1074 arg1: 0x0p+0 want: 0x0p+0 got: 0x1p1074 > 99 ulp want: 0 got: 1 round: m atan arg0: 0x1p1074 arg1: 0x0p+0 want: 0x0p+0 got: 0x1p1074 > 99 ulp want: 0 got: 1 round: z atan arg0: 0x1p1074 arg1: 0x0p+0 want: 0x0p+0 got: 0x1p1074 > 99 ulp want: 8000000000000000 got: 8000000000000001 round: z atan arg0: 0x1p1074 arg1: 0x0p+0 want: 0x0p+0 got: 0x1p1074 IMO this is 1ulp not >99ulp. For denormals ulp is not x * 0x1p52 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: 0x1p1074 > 99 ulp want: 8000000000000000 got: 8000000000000001 round: p expm1 arg0: 0x1p1074 arg1: 0x0p+0 want: 0x0p+0 got: 0x1p1074 > 99 ulp want: 1 got: 0 round: m expm1 arg0: 0x1p1074 arg1: 0x0p+0 want: 0x1p1074 got: 0x0p+0 > 99 ulp want: 1 got: 0 round: z expm1 arg0: 0x1p1074 arg1: 0x0p+0 want: 0x1p1074 got: 0x0p+0 > 99 ulp want: 8000000000000001 got: 8000000000000000 round: z log1p arg0: 0x1p1074 arg1: 0x0p+0 want: 0x1p1074 got: 0x0p+0 > 99 ulp want: 8000000000000000 got: 8000000000000001 round: p sin arg0: 0x1p1074 arg1: 0x0p+0 want: 0x0p+0 got: 0x1p1074 > 99 ulp want: 0 got: 1 round: m sin arg0: 0x1p1074 arg1: 0x0p+0 want: 0x0p+0 got: 0x1p1074 > 99 ulp want: 0 got: 1 round: z sin arg0: 0x1p1074 arg1: 0x0p+0 want: 0x0p+0 got: 0x1p1074 > 99 ulp want: 8000000000000000 got: 8000000000000001 round: z sin arg0: 0x1p1074 arg1: 0x0p+0 want: 0x0p+0 got: 0x1p1074 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: 0x1p1074 arg1: 0x0p+0 want: 0x1p1074 got: 0x0p+0 > 99 ulp want: 8000000000000000 got: 8000000000000001 round: m atan2 arg0: 0x1p1022 arg1: 0x1.fffffffffffffp+1023 want: 0x0p+0 got: 0x1p1074 > 99 ulp want: 8000000000000000 got: 8000000000000001 round: m atan2 arg0: 0x1p1022 arg1: 0x1.fffffffffffffp+1023 want: 0x0p+0 got: 0x1p1074 > 99 ulp want: 0 got: 1 round: p exp arg0: inf arg1: 0x0p+0 want: 0x0p+0 got: 0x1p1074 > 99 ulp want: 0 got: 1 round: p exp arg0: 0x1.75p+9 arg1: 0x0p+0 want: 0x0p+0 got: 0x1p1074 > 99 ulp want: 0 got: 1 round: p exp arg0: 0x1.c9c8p+13 arg1: 0x0p+0 want: 0x0p+0 got: 0x1p1074 And these. Of course in some cases the result is still nonconformant 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
Powered by blists  more mailing lists
Confused about mailing lists and their use? Read about mailing lists on Wikipedia and check out these guidelines on proper formatting of your messages.