Date: Tue, 29 May 2018 10:32:24 -0700 From: Steve Kargl <sgk@troutmask.apl.washington.edu> To: Adhemerval Zanella <adhemerval.zanella@linaro.org> Cc: Konstantin Belousov <kib@freebsd.org>, freebsd-hackers@freebsd.org, emaste@freebsd.org Subject: Re: Code with apache-2 on /usr/src Message-ID: <20180529173224.GA96547@troutmask.apl.washington.edu> In-Reply-To: <05943b3c-e2c6-4c03-93d9-5c2553e5865a@linaro.org> References: <20180528190444.GE3789@kib.kiev.ua> <f9f10762-651d-d2f2-c46f-6960b9a69705@linaro.org> <20180528193506.GA76705@troutmask.apl.washington.edu> <1c09023e-9bf5-d23a-dedc-1c4f4706bbde@linaro.org> <20180528202117.GA77184@troutmask.apl.washington.edu> <72101038-9e89-3f23-ab67-1c97b2a89803@linaro.org> <20180528210907.GA77475@troutmask.apl.washington.edu> <b79b4bc0-c584-1888-3207-9a7b640989fc@linaro.org> <20180528221819.GA77894@troutmask.apl.washington.edu> <05943b3c-e2c6-4c03-93d9-5c2553e5865a@linaro.org>
next in thread | previous in thread | raw e-mail | index | archive | help
On Tue, May 29, 2018 at 09:37:07AM -0300, Adhemerval Zanella wrote: > > > On 28/05/2018 19:18, Steve Kargl wrote: > > On Mon, May 28, 2018 at 06:12:13PM -0300, Adhemerval Zanella wrote: > >> > >>>> And is having a different algorithm for single and double prevision > >>>> a blocker for a future patch proposal? > >>> > >>> No. Given the comment in sinf.c that max ULP is 0.56072, I do note that > >>> the current implementation of sinf in lib/msun is more accurate (for > >>> interesting values of x). I also looked at single/s_sincosf.c. It is > >>> rather dubious to have 80+ digit numerical constants for a float, which > >>> at most has 9 relevant digits. > >>> > >> > >> Also keep in mind my initial idea is to propose patches only to expf, powf, > >> logf, expf2, and log2f. > > > > OK, so I peeked at expf. Comment claims max ulp of 0.502. > > Exhaustive testing for normal numbers in relevent range for > > the current implementation of expf(x) shows > > > > Interval tested: [-18,88.72] > > ULP: 0.90951, x = -5.19804668e+00f, /* 0xc0a65666 */ > > flt = 5.52735012e-03f, /* 0x3bb51ec6 */ > > dbl = 5.5273505437686398e-03, /* 0x3f76a3d8, 0xdd1aae8e */ > > > > But, then one looks at implementation details. msun's current > > implementation is written in terms of single precision; while > > the routine you're suggesting is written in terms of double_t. > > So, achieving 0.502 ULP is due to having 53-bits in intermediate > > results. It appears that the algorithm of the suggested code > > cannot easily be generalized to double and long double without > > implementing a multiple-precision routines. > > This is indeed true for the default implementation, although the same repo > has alternative implementation that uses only float for expf, powf, and > logf. However, as far as I could evaluated, the optimized expf and powf > single version does not yield any gain over current FreeBSD version, only > for the logf I see some gains. > > Do you see any issue about current approach of using intermediary double_t > for internal calculations? > No. The kernels for sinf and cosf (ie., k_sinf.c and k_cosf.c) use double for its intermediate computations. But, the main code in s_sin[fl].c and s_cos[f].c have the same internal structure: 1) Split argument into integer parts 2) Filter special values (+-inf, NaN) 3) Split into intervals a) for small x no range reduction is needed. b) do range reduction into [0,pi/4] 4) In (3a) deal with subnormal numbers with care to avoid spurious underflow. 5) In (3b), use polynomial approximations. Because the internal structure is similar for all precision, it makes maintenance easier. For maintenance and the importance of having the same structure, see the history of s_erff.c: https://svnweb.freebsd.org/base/head/lib/msun/src/s_erff.c?view=log > > Note, years ago, I submitted implementations for expf, exp, > > ld80/expl, ld128/expl, logf, log, ld80/logl, and ld128/logl > > based on papers by PTP Tang [1,2]. My versions for single > > and double precision were not adopted even though these had > > better accuracy. Either Bruce Evans improved or with Bruce's > > help I improved the ld80 and ld128 routines, which were added > > to msun. I know Bruce fixed minor issues with the single > > and double precision routines, but he has not submitted patches. > > > > 1. PTP Tang, "Table-driven implementation of the exponential > > function in IEEE floating-point arithmetic," ACM Trans. Math. > > Soft., 15, 144-157 (1989). > > > > 2. PTP Tang, "Table-driven implementation of the logarithm > > function in IEEE floating-point arithmetic," ACM Trans. Math. > > Soft., 16, 378-400 (1990). > > Thanks for the links, do you recall why exactly your implementations were > not adopted? Do you think a similar proposal based on the arm repo would > be also rejected? Mostly due to issues on my part. Bruce was/is the only person interested in reviewing patches to libm. At the time I submitted that code, his comments and suggestions could be characterized as drinking from a fire hose. When I had a commit bit, I finally gave up on the pursuit of perfect code and simply committed s_expl.c. Later, David Das committed s_logl.c. -- Steve
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20180529173224.GA96547>