Date: Sun, 9 Jun 2013 22:58:34 -0700 From: David Schultz <das@FreeBSD.ORG> To: Bruce Evans <brde@optusnet.com.au> Cc: freebsd-numerics@FreeBSD.org, Steve Kargl <sgk@troutmask.apl.washington.edu> Subject: Re: Implementation for coshl. Message-ID: <20130610055834.GA68643@zim.MIT.EDU> In-Reply-To: <20130610110740.V24058@besplex.bde.org> References: <20130610003645.GA16444@troutmask.apl.washington.edu> <20130610110740.V24058@besplex.bde.org>
next in thread | previous in thread | raw e-mail | index | archive | help
Awesome, thanks for working on this. On Mon, Jun 10, 2013, Bruce Evans wrote: > On Sun, 9 Jun 2013, Steve Kargl wrote: > > Arch | Interval | #calls | Time (s) | Max ULP | Compiler | Value > > -----------+---------------------+--------+----------+---------+----------+------- > > i386 [1] | [ 0.00: 0.35] | 100M | 15.0198 | 0.58583 | gcc | 1 > > i386 [1] | [ 0.35: 24.00] | 100M | 15.1858 | 1.01504 | gcc | 2 > > i386 [1] | [ 24.00:11356.52] | 100M | 12.9591 | 0.51198 | gcc | 3 > > i386 [1] | [11356.52:11357.22] | 100M | 13.3328 | 1.90988 | gcc | 4 > > -----------+---------------------+--------+----------+---------+----------+------- > > Quite large errors, unfortunately much the same as in double precision. [...] > > Bruce and I exchanged emails a long time ago about possible ways > > to reduce the ulp in this range by either computer exp(x) with > > extra precision or using a table with cosh(x) = cosh(x_i) * cosh(d) > > + sinh(x_i) * sinh(d) with d = x - x_i. I tried the latter with > > disappointing results. This is fine for now. We might eventually want to improve the accuracy of the float, double, and long double implementations, but that ought to come later. It's much easier to start with float and double when making improvements in the algorithm. > > The former would require a refactoring of > > s_expl.c into a kernel __kernel_expl(x, hi, lo). I have no plans on > > pursuing this at the this time. > > But you need this soon for __ldexp_exp() and __ldexp_cexp(), which are > needed for hyperbolic functions (the large args already fixed for float > and double precision) and for cexp() and trig and hyperbolic complex > functions. It is much easier to implement these using a kernel. I do > this only for float precision. That's fine for now. We need a better k_exp, too. I believe Bruce was working on one. I have a few comments on the code... mostly a subset of Bruce's comments: - You can clear the sign bit on x after the Inf/NaN check. - I think the way you handle the tiny-x case is an improvement over the double version. However, float-to-int conversion is very slow on x86, so you might instead try something like RETURNI(1.0L + tiny), where tiny is volatile and e.g., 0x1p-1000. That will also return the correct result in all rounding modes. - Please check whether volatile hacks are needed to induce the compiler to do the right thing. I'm increasingly finding that these are needed with clang. - I agree with Bruce about the thresholds not needing to be exact. In fact, if you try to use exact thresholds, the threshold might round the wrong way and give a subtle off-by-one-ulp error. Often just looking at the exponent is sufficient.
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20130610055834.GA68643>