From mboxrd@z Thu Jan 1 00:00:00 1970 Received: from eggs.gnu.org ([2001:4830:134:3::10]:56801) by lists.gnu.org with esmtp (Exim 4.71) (envelope-from ) id 1VghNC-0006ys-BK for qemu-devel@nongnu.org; Wed, 13 Nov 2013 15:49:20 -0500 Received: from Debian-exim by eggs.gnu.org with spam-scanned (Exim 4.71) (envelope-from ) id 1VghN6-0000X0-GG for qemu-devel@nongnu.org; Wed, 13 Nov 2013 15:49:14 -0500 Message-ID: <5283E5BE.1080809@gmail.com> Date: Wed, 13 Nov 2013 14:49:02 -0600 From: Tom Musta MIME-Version: 1.0 References: <1383769916-5582-1-git-send-email-tommusta@gmail.com> <1383769916-5582-13-git-send-email-tommusta@gmail.com> <527C221A.5000609@twiddle.net> <527C2297.7040407@twiddle.net> <527C2C92.7060006@twiddle.net> In-Reply-To: <527C2C92.7060006@twiddle.net> Content-Type: text/plain; charset=ISO-8859-1 Content-Transfer-Encoding: 7bit Subject: Re: [Qemu-devel] [PATCH 12/14] VSX Stage 4: Add Scalar SP Fused Multiply-Adds List-Id: List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , To: Richard Henderson , qemu-devel@nongnu.org Cc: qemu-ppc@nongnu.org On 11/7/2013 6:13 PM, Richard Henderson wrote: > On 11/08/2013 09:30 AM, Richard Henderson wrote: >> On 11/08/2013 09:28 AM, Richard Henderson wrote: >>> On 11/07/2013 06:31 AM, Tom Musta wrote: >>>> } \ >>>> + \ >>>> + if (r2sp) { \ >>>> + float32 tmp32 = float64_to_float32(xt_out.fld[i], \ >>>> + &env->fp_status); \ >>>> + xt_out.fld[i] = float32_to_float64(tmp32, &env->fp_status); \ >>>> + } \ >>>> + \ >>> >>> You can't get correct results for a single-precision fma from a >>> double-precision fma and merely rounding the results. >>> >>> See e.g. glibc's sysdeps/ieee754/dbl-64/s_fmaf.c. >> >> Blah, nevermind. That would be using separate add+mul in double-precision, not >> using a double-precision fma primitive. > > Hmph. I was right the first time. See > >> http://www.exploringbinary.com/double-rounding-errors-in-floating-point-conversions/ > > for example inputs that suffer from double-rounding. > > What's needed in each of the examples are infinite precision values containing > 55 bits. This is easy to accomplish with fma. > > Two 23-bit inputs can create a product with 46 significant bits. One can > append 23 more significant bits by choosing an exponent for the addend that > does not overlap the product. Thus one can create (almost) every intermediate > result with up to 69 consecutive bits (the exception being products without > factors that can be represented in 23-bits). > > I'm too lazy to decompose the examples therein to actual single-precision > inputs, but I'm certain it can be done. > > Thus you *do* need the round-to-zero plus inexact solution that glibc uses. > (Or to perform the fma in 128-bits and round once, but I think that would be > way more intrusive wrt the rest of the code, and more expensive than necessary.) I have reviewed the code and the spec and I cannot see a flaw. The sequence is effectively this: - float64_muladd - performs proper FMA for 64 bit numbers) - float64_to_float32 - converts to single precision, including proper rounding - float32_to_float64 The implementation of float64_muladd would seem to provide enough mantissa bits for proper handling of the case you describe. The only rounding occurs in the second step. I have also done quite a bit of random and targeted random testing using Power hardware to produce expected results. The targeted random tests followed your suggestion above: generate AxB + C where abs(exp(A) - exp(B)) = 23 and abs(exp(A) - exp(C)) = 46. Several million test patterns have been generated and played back through QEMU without any miscompares in the numerical results.