From: Rich Felker <dalias@aerifal.cx>
To: musl@lists.openwall.com
Subject: Re: math todo
Date: Tue, 1 May 2012 10:44:38 -0400 [thread overview]
Message-ID: <20120501144438.GL14673@brightrain.aerifal.cx> (raw)
In-Reply-To: <20120501091458.GJ16237@port70.net>
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<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: -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
next prev parent reply other threads:[~2012-05-01 14:44 UTC|newest]
Thread overview: 5+ messages / expand[flat|nested] mbox.gz Atom feed top
2012-05-01 0:05 Szabolcs Nagy
2012-05-01 1:19 ` Rich Felker
2012-05-01 9:14 ` Szabolcs Nagy
2012-05-01 14:44 ` Rich Felker [this message]
2012-05-01 18:33 ` Szabolcs Nagy
Reply instructions:
You may reply publicly to this message via plain-text email
using any one of the following methods:
* Save the following mbox file, import it into your mail client,
and reply-to-all from there: mbox
Avoid top-posting and favor interleaved quoting:
https://en.wikipedia.org/wiki/Posting_style#Interleaved_style
* Reply using the --to, --cc, and --in-reply-to
switches of git-send-email(1):
git send-email \
--in-reply-to=20120501144438.GL14673@brightrain.aerifal.cx \
--to=dalias@aerifal.cx \
--cc=musl@lists.openwall.com \
/path/to/YOUR_REPLY
https://kernel.org/pub/software/scm/git/docs/git-send-email.html
* If your mail client supports setting the In-Reply-To header
via mailto: links, try the mailto: link
Be sure your reply has a Subject: header at the top and a blank line
before the message body.
Code repositories for project(s) associated with this public inbox
https://git.vuxu.org/mirror/musl/
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for NNTP newsgroup(s).