This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Re: correct rounding or not?
- From: Wilco Dijkstra <Wilco dot Dijkstra at arm dot com>
- To: paul zimmermann <Paul dot Zimmermann at inria dot fr>
- Cc: "libc-alpha at sourceware dot org" <libc-alpha at sourceware dot org>
- Date: Fri, 31 Jan 2020 14:02:48 +0000
- Subject: Re: correct rounding or not?
- Arc-authentication-results: i=1; mx.microsoft.com 1; spf=pass smtp.mailfrom=arm.com; dmarc=pass action=none header.from=arm.com; dkim=pass header.d=arm.com; arc=none
- Arc-message-signature: i=1; a=rsa-sha256; c=relaxed/relaxed; d=microsoft.com; s=arcselector9901; h=From:Date:Subject:Message-ID:Content-Type:MIME-Version:X-MS-Exchange-SenderADCheck; bh=epnZxxPi5cCMYjs8tACSdnGTSi6yZucehl4PUo+cpyM=; b=XcNJGc/EQaD9qRCcguB3/Urk2sUrQf0X/Lm2Dh3ZtVT7YPpcQ2Pc3SQGhFcWu7+GiV9+Mmkf96dKlSFVyCg2TcPBL5DbyHuuDbexDEJYzrMoy9/vjISXHMFDiFnLLYKt5EFEA7P2RlatMMtInaP3wiZ8zOmXluFCWEl2Ruc04wGITrC59WzkKa0lqLnZehad5ziY7V0QSKzSmJ8YHp+cebFpuOHEtUpni0/Cjz6Omj37TgvLMC0jIfFxN1zZhh/5Va1s4djqkY3ayy+4mX1vakCU9cC5xIr58cUsRY6Fv8XeMrLULyM8t3JIYaOh0Yuttj/8PH1GooWQ/4/I8xOT2g==
- Arc-seal: i=1; a=rsa-sha256; s=arcselector9901; d=microsoft.com; cv=none; b=I9dz/beNa9wBy2ohvBN3014W/puyaVGrnulnSCI8x5KWpaFcZ0qhtd3Y6YyrgewVKmVnJHHEIS4Yb0Rku1ry7Rl3mRcgaVm+EXmcqIwRdAoQPsGBFYv65EeSuex7RHCetkfSlmM/tz7jEw4hE7amM8OCKoldL59eAZYwwLXRc0zw7AtRgsg2oUfrKIqF8wnfDDJW1hkfSG1aTmaDaOjpS9gKwBRLXze5gfxyi5GXNBCH14fH5tWzboDX1fgeVSQ5jY9dKgzRGdwNirnMKbWU5xFr+N/H4KR6W4ORvTnG3V2IgjWML8KlLKP+ML8YxNTH1kSajg2yPKaLA7pHE2htrg==
- Original-authentication-results: spf=none (sender IP is ) smtp.mailfrom=Wilco dot Dijkstra at arm dot com;
- References: <AM5PR0801MB2035BC46C0EE0E135316E21183040@AM5PR0801MB2035.eurprd08.prod.outlook.com>,<mw7e17j4p5.fsf@tomate.loria.fr>
Hi Paul,
> for [1] I don't know what "it computes *the* rounded value of sin(x)" means.
> The header of the file says that the error is ~0.55ulp, that would be clearer.
> Maybe the original libultim source contained an error bound for the fast path,
> but unfortunately I do not have it any more.
Well we don't have good error bounds for sin/cos, the ULP estimates are based
on the original code with the slow paths.
> However, in revision e4d8276, I see the largest 'corr' factor is 1.07,
> which corresponds to a relative error of 0.07/2^55, thus to a maximum
> error of 0.5175 ulps. But maybe the code was changed since then.
__sin used to have this:
else if (k < 0x3feb6000)
{
res = do_sin (x, 0, &cor);
retval = (res == res + 1.096 * cor) ? res : slow1 (x);
>From this I estimated the worst-case ULP as < 0.55. It could well be a bit
better in reality.
> Note also that we can still find "correctly rounded" in e_atan2.c, s_atan.c,
> and s_tan.c.
Yes, those still have slow paths which give correctly rounded results.
> /tmp/mpcheck-double --num=14400000 --seed=1
That looks like a useful tool - we did something similar for random testing
using long double math to test double functions.
> Also, "make bench" on my machine gives the following results:
function & fast path & slow path \\ \hline
acos & 168 & 160894 \\ % 235 after patch [2]
asin & 160 & 163083 \\ % 216 after patch [2]
atan & 127 & 29389 \\ % 144 bits
sin & 92 & 198 \\ % 768 bits
cos & 115 & 211 \\ % 768 bits
tan & 189 & 172531 \\ % 768 bits
exp & 25 & 26 \\ % 144 bits
log & 20 & \\
pow & 35 & 67 \\ % 240 bits
> Note that atan and tan still have slow "slow paths".
Yes the 1000x slowdown to try to round correctly is completely ridiculous!
It's also interesting to see that the new exp/log worst case is 5x faster than
sin/cos fast case...
> Also I wonder, now that the "slow path" for asin takes only 216 cycles,
> whereas the fast path takes 160 cycles, what would you obtain if you use
> the slow path only?
The fast path is likely already more than accurate enough. It's unlikely to be
worth it though - rewriting these functions from scratch is the way forward.
Cheers,
Wilco