fma_test.c revision 251024
1193326Sed/*- 2193326Sed * Copyright (c) 2008 David Schultz <das@FreeBSD.org> 3193326Sed * All rights reserved. 4193326Sed * 5193326Sed * Redistribution and use in source and binary forms, with or without 6193326Sed * modification, are permitted provided that the following conditions 7193326Sed * are met: 8193326Sed * 1. Redistributions of source code must retain the above copyright 9193326Sed * notice, this list of conditions and the following disclaimer. 10193326Sed * 2. Redistributions in binary form must reproduce the above copyright 11193326Sed * notice, this list of conditions and the following disclaimer in the 12193326Sed * documentation and/or other materials provided with the distribution. 13193326Sed * 14193326Sed * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND 15193326Sed * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 16193326Sed * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 17193326Sed * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE 18193326Sed * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 19193326Sed * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 20193326Sed * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 21249423Sdim * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 22249423Sdim * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 23203955Srdivacky * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 24226633Sdim * SUCH DAMAGE. 25221345Sdim */ 26234353Sdim 27193326Sed/* 28193326Sed * Tests for fma{,f,l}(). 29249423Sdim */ 30201361Srdivacky 31249423Sdim#include <sys/cdefs.h> 32249423Sdim__FBSDID("$FreeBSD: head/tools/regression/lib/msun/test-fma.c 251024 2013-05-27 08:50:10Z das $"); 33193326Sed 34193326Sed#include <assert.h> 35249423Sdim#include <fenv.h> 36249423Sdim#include <float.h> 37249423Sdim#include <math.h> 38249423Sdim#include <stdio.h> 39249423Sdim 40193326Sed#define ALL_STD_EXCEPT (FE_DIVBYZERO | FE_INEXACT | FE_INVALID | \ 41193326Sed FE_OVERFLOW | FE_UNDERFLOW) 42249423Sdim 43210299Sed#pragma STDC FENV_ACCESS ON 44198092Srdivacky 45207619Srdivacky/* 46249423Sdim * Test that a function returns the correct value and sets the 47198092Srdivacky * exception flags correctly. The exceptmask specifies which 48198092Srdivacky * exceptions we should check. We need to be lenient for several 49198092Srdivacky * reasons, but mainly because on some architectures it's impossible 50198092Srdivacky * to raise FE_OVERFLOW without raising FE_INEXACT. 51200583Srdivacky * 52198092Srdivacky * These are macros instead of functions so that assert provides more 53198092Srdivacky * meaningful error messages. 54198092Srdivacky */ 55198092Srdivacky#define test(func, x, y, z, result, exceptmask, excepts) do { \ 56200583Srdivacky assert(feclearexcept(FE_ALL_EXCEPT) == 0); \ 57198092Srdivacky assert(fpequal((func)((x), (y), (z)), (result))); \ 58198092Srdivacky assert(((func), fetestexcept(exceptmask) == (excepts))); \ 59198092Srdivacky} while (0) 60198092Srdivacky 61200583Srdivacky#define testall(x, y, z, result, exceptmask, excepts) do { \ 62198092Srdivacky test(fma, (x), (y), (z), (double)(result), (exceptmask), (excepts)); \ 63198398Srdivacky test(fmaf, (x), (y), (z), (float)(result), (exceptmask), (excepts)); \ 64198398Srdivacky test(fmal, (x), (y), (z), (result), (exceptmask), (excepts)); \ 65198398Srdivacky} while (0) 66198092Srdivacky 67208600Srdivacky/* Test in all rounding modes. */ 68198092Srdivacky#define testrnd(func, x, y, z, rn, ru, rd, rz, exceptmask, excepts) do { \ 69198092Srdivacky fesetround(FE_TONEAREST); \ 70193326Sed test((func), (x), (y), (z), (rn), (exceptmask), (excepts)); \ 71193326Sed fesetround(FE_UPWARD); \ 72234353Sdim test((func), (x), (y), (z), (ru), (exceptmask), (excepts)); \ 73195341Sed fesetround(FE_DOWNWARD); \ 74198092Srdivacky test((func), (x), (y), (z), (rd), (exceptmask), (excepts)); \ 75201361Srdivacky fesetround(FE_TOWARDZERO); \ 76201361Srdivacky test((func), (x), (y), (z), (rz), (exceptmask), (excepts)); \ 77201361Srdivacky} while (0) 78201361Srdivacky 79195341Sed/* 80193326Sed * This is needed because clang constant-folds fma in ways that are incorrect 81195341Sed * in rounding modes other than FE_TONEAREST. 82201361Srdivacky */ 83193326Sedvolatile double one = 1.0; 84195341Sed 85198092Srdivacky/* 86201361Srdivacky * Determine whether x and y are equal, with two special rules: 87201361Srdivacky * +0.0 != -0.0 88201361Srdivacky * NaN == NaN 89193326Sed */ 90193326Sedint 91203955Srdivackyfpequal(long double x, long double y) 92203955Srdivacky{ 93193326Sed 94193326Sed return ((x == y && !signbit(x) == !signbit(y)) 95193326Sed || (isnan(x) && isnan(y))); 96193326Sed} 97193326Sed 98193326Sedstatic void 99193326Sedtest_zeroes(void) 100193326Sed{ 101193326Sed const int rd = (fegetround() == FE_DOWNWARD); 102239462Sdim 103193326Sed testall(0.0, 0.0, 0.0, 0.0, ALL_STD_EXCEPT, 0); 104234353Sdim testall(1.0, 0.0, 0.0, 0.0, ALL_STD_EXCEPT, 0); 105193326Sed testall(0.0, 1.0, 0.0, 0.0, ALL_STD_EXCEPT, 0); 106193326Sed testall(0.0, 0.0, 1.0, 1.0, ALL_STD_EXCEPT, 0); 107193326Sed 108193326Sed testall(-0.0, 0.0, 0.0, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0); 109193326Sed testall(0.0, -0.0, 0.0, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0); 110234353Sdim testall(-0.0, -0.0, 0.0, 0.0, ALL_STD_EXCEPT, 0); 111234353Sdim testall(0.0, 0.0, -0.0, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0); 112249423Sdim testall(-0.0, -0.0, -0.0, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0); 113234353Sdim 114193326Sed testall(-0.0, 0.0, -0.0, -0.0, ALL_STD_EXCEPT, 0); 115193326Sed testall(0.0, -0.0, -0.0, -0.0, ALL_STD_EXCEPT, 0); 116193326Sed 117193326Sed testall(-one, one, one, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0); 118193326Sed testall(one, -one, one, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0); 119193326Sed testall(-one, -one, -one, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0); 120193326Sed 121193326Sed switch (fegetround()) { 122193326Sed case FE_TONEAREST: 123193326Sed case FE_TOWARDZERO: 124193326Sed test(fmaf, -FLT_MIN, FLT_MIN, 0.0, -0.0, 125198398Srdivacky ALL_STD_EXCEPT, FE_INEXACT | FE_UNDERFLOW); 126198398Srdivacky test(fma, -DBL_MIN, DBL_MIN, 0.0, -0.0, 127198398Srdivacky ALL_STD_EXCEPT, FE_INEXACT | FE_UNDERFLOW); 128226633Sdim test(fmal, -LDBL_MIN, LDBL_MIN, 0.0, -0.0, 129198398Srdivacky ALL_STD_EXCEPT, FE_INEXACT | FE_UNDERFLOW); 130198398Srdivacky } 131198398Srdivacky} 132198398Srdivacky 133198398Srdivackystatic void 134198398Srdivackytest_infinities(void) 135198398Srdivacky{ 136198398Srdivacky 137198398Srdivacky testall(INFINITY, 1.0, -1.0, INFINITY, ALL_STD_EXCEPT, 0); 138198398Srdivacky testall(-1.0, INFINITY, 0.0, -INFINITY, ALL_STD_EXCEPT, 0); 139198398Srdivacky testall(0.0, 0.0, INFINITY, INFINITY, ALL_STD_EXCEPT, 0); 140198398Srdivacky testall(1.0, 1.0, INFINITY, INFINITY, ALL_STD_EXCEPT, 0); 141198398Srdivacky testall(1.0, 1.0, -INFINITY, -INFINITY, ALL_STD_EXCEPT, 0); 142198398Srdivacky 143198398Srdivacky testall(INFINITY, -INFINITY, 1.0, -INFINITY, ALL_STD_EXCEPT, 0); 144198398Srdivacky testall(INFINITY, INFINITY, 1.0, INFINITY, ALL_STD_EXCEPT, 0); 145226633Sdim testall(-INFINITY, -INFINITY, INFINITY, INFINITY, ALL_STD_EXCEPT, 0); 146208600Srdivacky 147193326Sed testall(0.0, INFINITY, 1.0, NAN, ALL_STD_EXCEPT, FE_INVALID); 148193326Sed testall(INFINITY, 0.0, -0.0, NAN, ALL_STD_EXCEPT, FE_INVALID); 149193326Sed 150193326Sed /* The invalid exception is optional in this case. */ 151193326Sed testall(INFINITY, 0.0, NAN, NAN, ALL_STD_EXCEPT & ~FE_INVALID, 0); 152193326Sed 153193326Sed testall(INFINITY, INFINITY, -INFINITY, NAN, 154249423Sdim ALL_STD_EXCEPT, FE_INVALID); 155193326Sed testall(-INFINITY, INFINITY, INFINITY, NAN, 156193326Sed ALL_STD_EXCEPT, FE_INVALID); 157249423Sdim testall(INFINITY, -1.0, INFINITY, NAN, 158193326Sed ALL_STD_EXCEPT, FE_INVALID); 159193326Sed 160249423Sdim test(fmaf, FLT_MAX, FLT_MAX, -INFINITY, -INFINITY, ALL_STD_EXCEPT, 0); 161249423Sdim test(fma, DBL_MAX, DBL_MAX, -INFINITY, -INFINITY, ALL_STD_EXCEPT, 0); 162249423Sdim test(fmal, LDBL_MAX, LDBL_MAX, -INFINITY, -INFINITY, 163249423Sdim ALL_STD_EXCEPT, 0); 164193326Sed test(fmaf, FLT_MAX, -FLT_MAX, INFINITY, INFINITY, ALL_STD_EXCEPT, 0); 165198092Srdivacky test(fma, DBL_MAX, -DBL_MAX, INFINITY, INFINITY, ALL_STD_EXCEPT, 0); 166193326Sed test(fmal, LDBL_MAX, -LDBL_MAX, INFINITY, INFINITY, 167198092Srdivacky ALL_STD_EXCEPT, 0); 168249423Sdim} 169198092Srdivacky 170198092Srdivackystatic void 171198092Srdivackytest_nans(void) 172249423Sdim{ 173198092Srdivacky 174249423Sdim testall(NAN, 0.0, 0.0, NAN, ALL_STD_EXCEPT, 0); 175198092Srdivacky testall(1.0, NAN, 1.0, NAN, ALL_STD_EXCEPT, 0); 176249423Sdim testall(1.0, -1.0, NAN, NAN, ALL_STD_EXCEPT, 0); 177198092Srdivacky testall(0.0, 0.0, NAN, NAN, ALL_STD_EXCEPT, 0); 178193326Sed testall(NAN, NAN, NAN, NAN, ALL_STD_EXCEPT, 0); 179193326Sed 180193326Sed /* x*y should not raise an inexact/overflow/underflow if z is NaN. */ 181193326Sed testall(M_PI, M_PI, NAN, NAN, ALL_STD_EXCEPT, 0); 182193326Sed test(fmaf, FLT_MIN, FLT_MIN, NAN, NAN, ALL_STD_EXCEPT, 0); 183193326Sed test(fma, DBL_MIN, DBL_MIN, NAN, NAN, ALL_STD_EXCEPT, 0); 184193326Sed test(fmal, LDBL_MIN, LDBL_MIN, NAN, NAN, ALL_STD_EXCEPT, 0); 185193326Sed test(fmaf, FLT_MAX, FLT_MAX, NAN, NAN, ALL_STD_EXCEPT, 0); 186193326Sed test(fma, DBL_MAX, DBL_MAX, NAN, NAN, ALL_STD_EXCEPT, 0); 187193326Sed test(fmal, LDBL_MAX, LDBL_MAX, NAN, NAN, ALL_STD_EXCEPT, 0); 188193326Sed} 189193326Sed 190234353Sdim/* 191234353Sdim * Tests for cases where z is very small compared to x*y. 192226633Sdim */ 193234353Sdimstatic void 194234353Sdimtest_small_z(void) 195226633Sdim{ 196201361Srdivacky 197201361Srdivacky /* x*y positive, z positive */ 198201361Srdivacky if (fegetround() == FE_UPWARD) { 199201361Srdivacky test(fmaf, one, one, 0x1.0p-100, 1.0 + FLT_EPSILON, 200201361Srdivacky ALL_STD_EXCEPT, FE_INEXACT); 201201361Srdivacky test(fma, one, one, 0x1.0p-200, 1.0 + DBL_EPSILON, 202201361Srdivacky ALL_STD_EXCEPT, FE_INEXACT); 203201361Srdivacky test(fmal, one, one, 0x1.0p-200, 1.0 + LDBL_EPSILON, 204201361Srdivacky ALL_STD_EXCEPT, FE_INEXACT); 205201361Srdivacky } else { 206201361Srdivacky testall(0x1.0p100, one, 0x1.0p-100, 0x1.0p100, 207201361Srdivacky ALL_STD_EXCEPT, FE_INEXACT); 208201361Srdivacky } 209201361Srdivacky 210234353Sdim /* x*y negative, z negative */ 211234353Sdim if (fegetround() == FE_DOWNWARD) { 212207619Srdivacky test(fmaf, -one, one, -0x1.0p-100, -(1.0 + FLT_EPSILON), 213207619Srdivacky ALL_STD_EXCEPT, FE_INEXACT); 214249423Sdim test(fma, -one, one, -0x1.0p-200, -(1.0 + DBL_EPSILON), 215249423Sdim ALL_STD_EXCEPT, FE_INEXACT); 216218893Sdim test(fmal, -one, one, -0x1.0p-200, -(1.0 + LDBL_EPSILON), 217249423Sdim ALL_STD_EXCEPT, FE_INEXACT); 218249423Sdim } else { 219249423Sdim testall(0x1.0p100, -one, -0x1.0p-100, -0x1.0p100, 220249423Sdim ALL_STD_EXCEPT, FE_INEXACT); 221218893Sdim } 222218893Sdim 223234353Sdim /* x*y positive, z negative */ 224249423Sdim if (fegetround() == FE_DOWNWARD || fegetround() == FE_TOWARDZERO) { 225234353Sdim test(fmaf, one, one, -0x1.0p-100, 1.0 - FLT_EPSILON / 2, 226218893Sdim ALL_STD_EXCEPT, FE_INEXACT); 227218893Sdim test(fma, one, one, -0x1.0p-200, 1.0 - DBL_EPSILON / 2, 228218893Sdim ALL_STD_EXCEPT, FE_INEXACT); 229218893Sdim test(fmal, one, one, -0x1.0p-200, 1.0 - LDBL_EPSILON / 2, 230249423Sdim ALL_STD_EXCEPT, FE_INEXACT); 231249423Sdim } else { 232249423Sdim testall(0x1.0p100, one, -0x1.0p-100, 0x1.0p100, 233249423Sdim ALL_STD_EXCEPT, FE_INEXACT); 234249423Sdim } 235249423Sdim 236221345Sdim /* x*y negative, z positive */ 237221345Sdim if (fegetround() == FE_UPWARD || fegetround() == FE_TOWARDZERO) { 238249423Sdim test(fmaf, -one, one, 0x1.0p-100, -1.0 + FLT_EPSILON / 2, 239249423Sdim ALL_STD_EXCEPT, FE_INEXACT); 240221345Sdim test(fma, -one, one, 0x1.0p-200, -1.0 + DBL_EPSILON / 2, 241249423Sdim ALL_STD_EXCEPT, FE_INEXACT); 242249423Sdim test(fmal, -one, one, 0x1.0p-200, -1.0 + LDBL_EPSILON / 2, 243249423Sdim ALL_STD_EXCEPT, FE_INEXACT); 244218893Sdim } else { 245195099Sed testall(-0x1.0p100, one, 0x1.0p-100, -0x1.0p100, 246195099Sed ALL_STD_EXCEPT, FE_INEXACT); 247234353Sdim } 248234353Sdim} 249234353Sdim 250234353Sdim/* 251234353Sdim * Tests for cases where z is very large compared to x*y. 252234353Sdim */ 253234353Sdimstatic void 254234353Sdimtest_big_z(void) 255195099Sed{ 256195099Sed 257195099Sed /* z positive, x*y positive */ 258198092Srdivacky if (fegetround() == FE_UPWARD) { 259203955Srdivacky test(fmaf, 0x1.0p-50, 0x1.0p-50, 1.0, 1.0 + FLT_EPSILON, 260210299Sed ALL_STD_EXCEPT, FE_INEXACT); 261193326Sed test(fma, 0x1.0p-100, 0x1.0p-100, 1.0, 1.0 + DBL_EPSILON, 262193326Sed ALL_STD_EXCEPT, FE_INEXACT); 263226633Sdim test(fmal, 0x1.0p-100, 0x1.0p-100, 1.0, 1.0 + LDBL_EPSILON, 264226633Sdim ALL_STD_EXCEPT, FE_INEXACT); 265207619Srdivacky } else { 266207619Srdivacky testall(-0x1.0p-50, -0x1.0p-50, 0x1.0p100, 0x1.0p100, 267207619Srdivacky ALL_STD_EXCEPT, FE_INEXACT); 268218893Sdim } 269218893Sdim 270218893Sdim /* z negative, x*y negative */ 271218893Sdim if (fegetround() == FE_DOWNWARD) { 272218893Sdim test(fmaf, -0x1.0p-50, 0x1.0p-50, -1.0, -(1.0 + FLT_EPSILON), 273218893Sdim ALL_STD_EXCEPT, FE_INEXACT); 274234353Sdim test(fma, -0x1.0p-100, 0x1.0p-100, -1.0, -(1.0 + DBL_EPSILON), 275218893Sdim ALL_STD_EXCEPT, FE_INEXACT); 276221345Sdim test(fmal, -0x1.0p-100, 0x1.0p-100, -1.0, -(1.0 + LDBL_EPSILON), 277221345Sdim ALL_STD_EXCEPT, FE_INEXACT); 278221345Sdim } else { 279221345Sdim testall(0x1.0p-50, -0x1.0p-50, -0x1.0p100, -0x1.0p100, 280221345Sdim ALL_STD_EXCEPT, FE_INEXACT); 281221345Sdim } 282221345Sdim 283221345Sdim /* z negative, x*y positive */ 284221345Sdim if (fegetround() == FE_UPWARD || fegetround() == FE_TOWARDZERO) { 285218893Sdim test(fmaf, -0x1.0p-50, -0x1.0p-50, -1.0, 286218893Sdim -1.0 + FLT_EPSILON / 2, ALL_STD_EXCEPT, FE_INEXACT); 287221345Sdim test(fma, -0x1.0p-100, -0x1.0p-100, -1.0, 288221345Sdim -1.0 + DBL_EPSILON / 2, ALL_STD_EXCEPT, FE_INEXACT); 289221345Sdim test(fmal, -0x1.0p-100, -0x1.0p-100, -1.0, 290221345Sdim -1.0 + LDBL_EPSILON / 2, ALL_STD_EXCEPT, FE_INEXACT); 291234353Sdim } else { 292234353Sdim testall(0x1.0p-50, 0x1.0p-50, -0x1.0p100, -0x1.0p100, 293218893Sdim ALL_STD_EXCEPT, FE_INEXACT); 294218893Sdim } 295221345Sdim 296221345Sdim /* z positive, x*y negative */ 297221345Sdim if (fegetround() == FE_DOWNWARD || fegetround() == FE_TOWARDZERO) { 298221345Sdim test(fmaf, 0x1.0p-50, -0x1.0p-50, 1.0, 1.0 - FLT_EPSILON / 2, 299234353Sdim ALL_STD_EXCEPT, FE_INEXACT); 300221345Sdim test(fma, 0x1.0p-100, -0x1.0p-100, 1.0, 1.0 - DBL_EPSILON / 2, 301221345Sdim ALL_STD_EXCEPT, FE_INEXACT); 302221345Sdim test(fmal, 0x1.0p-100, -0x1.0p-100, 1.0, 1.0 - LDBL_EPSILON / 2, 303218893Sdim ALL_STD_EXCEPT, FE_INEXACT); 304218893Sdim } else { 305218893Sdim testall(-0x1.0p-50, 0x1.0p-50, 0x1.0p100, 0x1.0p100, 306218893Sdim ALL_STD_EXCEPT, FE_INEXACT); 307234353Sdim } 308193326Sed} 309234353Sdim 310234353Sdimstatic void 311234353Sdimtest_accuracy(void) 312234353Sdim{ 313212904Sdim 314221345Sdim /* ilogb(x*y) - ilogb(z) = 20 */ 315221345Sdim testrnd(fmaf, -0x1.c139d8p-51, -0x1.600e7ap32, 0x1.26558cp-38, 316221345Sdim 0x1.34e48ap-18, 0x1.34e48cp-18, 0x1.34e48ap-18, 0x1.34e48ap-18, 317221345Sdim ALL_STD_EXCEPT, FE_INEXACT); 318221345Sdim testrnd(fma, -0x1.c139d7b84f1a3p-51, -0x1.600e7a2a16484p32, 319198092Srdivacky 0x1.26558cac31580p-38, 0x1.34e48a78aae97p-18, 320234353Sdim 0x1.34e48a78aae97p-18, 0x1.34e48a78aae96p-18, 321234353Sdim 0x1.34e48a78aae96p-18, ALL_STD_EXCEPT, FE_INEXACT); 322234353Sdim#if LDBL_MANT_DIG == 113 323234353Sdim testrnd(fmal, -0x1.c139d7b84f1a3079263afcc5bae3p-51L, 324234353Sdim -0x1.600e7a2a164840edbe2e7d301a72p32L, 325198092Srdivacky 0x1.26558cac315807eb07e448042101p-38L, 326234353Sdim 0x1.34e48a78aae96c76ed36077dd387p-18L, 327234353Sdim 0x1.34e48a78aae96c76ed36077dd388p-18L, 328234353Sdim 0x1.34e48a78aae96c76ed36077dd387p-18L, 329234353Sdim 0x1.34e48a78aae96c76ed36077dd387p-18L, 330234353Sdim ALL_STD_EXCEPT, FE_INEXACT); 331234353Sdim#elif LDBL_MANT_DIG == 64 332234353Sdim testrnd(fmal, -0x1.c139d7b84f1a307ap-51L, -0x1.600e7a2a164840eep32L, 333234353Sdim 0x1.26558cac315807ecp-38L, 0x1.34e48a78aae96c78p-18L, 334234353Sdim 0x1.34e48a78aae96c78p-18L, 0x1.34e48a78aae96c76p-18L, 335234353Sdim 0x1.34e48a78aae96c76p-18L, ALL_STD_EXCEPT, FE_INEXACT); 336234353Sdim#elif LDBL_MANT_DIG == 53 337234353Sdim testrnd(fmal, -0x1.c139d7b84f1a3p-51L, -0x1.600e7a2a16484p32L, 338234353Sdim 0x1.26558cac31580p-38L, 0x1.34e48a78aae97p-18L, 339234353Sdim 0x1.34e48a78aae97p-18L, 0x1.34e48a78aae96p-18L, 340234353Sdim 0x1.34e48a78aae96p-18L, ALL_STD_EXCEPT, FE_INEXACT); 341193326Sed#endif 342193326Sed 343234353Sdim /* ilogb(x*y) - ilogb(z) = -40 */ 344234353Sdim testrnd(fmaf, 0x1.98210ap53, 0x1.9556acp-24, 0x1.d87da4p70, 345234353Sdim 0x1.d87da4p70, 0x1.d87da6p70, 0x1.d87da4p70, 0x1.d87da4p70, 346198092Srdivacky ALL_STD_EXCEPT, FE_INEXACT); 347234353Sdim testrnd(fma, 0x1.98210ac83fe2bp53, 0x1.9556ac1475f0fp-24, 348234353Sdim 0x1.d87da3aafc60ep70, 0x1.d87da3aafda40p70, 349234353Sdim 0x1.d87da3aafda40p70, 0x1.d87da3aafda3fp70, 350234353Sdim 0x1.d87da3aafda3fp70, ALL_STD_EXCEPT, FE_INEXACT); 351234353Sdim#if LDBL_MANT_DIG == 113 352234353Sdim testrnd(fmal, 0x1.98210ac83fe2a8f65b6278b74cebp53L, 353234353Sdim 0x1.9556ac1475f0f28968b61d0de65ap-24L, 354234353Sdim 0x1.d87da3aafc60d830aa4c6d73b749p70L, 355212904Sdim 0x1.d87da3aafda3f36a69eb86488224p70L, 356212904Sdim 0x1.d87da3aafda3f36a69eb86488225p70L, 357212904Sdim 0x1.d87da3aafda3f36a69eb86488224p70L, 358208600Srdivacky 0x1.d87da3aafda3f36a69eb86488224p70L, 359212904Sdim ALL_STD_EXCEPT, FE_INEXACT); 360212904Sdim#elif LDBL_MANT_DIG == 64 361212904Sdim testrnd(fmal, 0x1.98210ac83fe2a8f6p53L, 0x1.9556ac1475f0f28ap-24L, 362212904Sdim 0x1.d87da3aafc60d83p70L, 0x1.d87da3aafda3f36ap70L, 363212904Sdim 0x1.d87da3aafda3f36ap70L, 0x1.d87da3aafda3f368p70L, 364198092Srdivacky 0x1.d87da3aafda3f368p70L, ALL_STD_EXCEPT, FE_INEXACT); 365198092Srdivacky#elif LDBL_MANT_DIG == 53 366198092Srdivacky testrnd(fmal, 0x1.98210ac83fe2bp53L, 0x1.9556ac1475f0fp-24L, 367198092Srdivacky 0x1.d87da3aafc60ep70L, 0x1.d87da3aafda40p70L, 368212904Sdim 0x1.d87da3aafda40p70L, 0x1.d87da3aafda3fp70L, 369212904Sdim 0x1.d87da3aafda3fp70L, ALL_STD_EXCEPT, FE_INEXACT); 370234353Sdim#endif 371212904Sdim 372212904Sdim /* ilogb(x*y) - ilogb(z) = 0 */ 373212904Sdim testrnd(fmaf, 0x1.31ad02p+100, 0x1.2fbf7ap-42, -0x1.c3e106p+58, 374212904Sdim -0x1.64c27cp+56, -0x1.64c27ap+56, -0x1.64c27cp+56, 375234353Sdim -0x1.64c27ap+56, ALL_STD_EXCEPT, FE_INEXACT); 376212904Sdim testrnd(fma, 0x1.31ad012ede8aap+100, 0x1.2fbf79c839067p-42, 377212904Sdim -0x1.c3e106929056ep+58, -0x1.64c282b970a5fp+56, 378234353Sdim -0x1.64c282b970a5ep+56, -0x1.64c282b970a5fp+56, 379234353Sdim -0x1.64c282b970a5ep+56, ALL_STD_EXCEPT, FE_INEXACT); 380234353Sdim#if LDBL_MANT_DIG == 113 381234353Sdim testrnd(fmal, 0x1.31ad012ede8aa282fa1c19376d16p+100L, 382234353Sdim 0x1.2fbf79c839066f0f5c68f6d2e814p-42L, 383234353Sdim -0x1.c3e106929056ec19de72bfe64215p+58L, 384218893Sdim -0x1.64c282b970a612598fc025ca8cddp+56L, 385208600Srdivacky -0x1.64c282b970a612598fc025ca8cddp+56L, 386208600Srdivacky -0x1.64c282b970a612598fc025ca8cdep+56L, 387234353Sdim -0x1.64c282b970a612598fc025ca8cddp+56L, 388234353Sdim ALL_STD_EXCEPT, FE_INEXACT); 389234353Sdim#elif LDBL_MANT_DIG == 64 390206084Srdivacky testrnd(fmal, 0x1.31ad012ede8aa4eap+100L, 0x1.2fbf79c839066aeap-42L, 391234353Sdim -0x1.c3e106929056e61p+58L, -0x1.64c282b970a60298p+56L, 392193326Sed -0x1.64c282b970a60298p+56L, -0x1.64c282b970a6029ap+56L, 393198092Srdivacky -0x1.64c282b970a60298p+56L, ALL_STD_EXCEPT, FE_INEXACT); 394208600Srdivacky#elif LDBL_MANT_DIG == 53 395208600Srdivacky testrnd(fmal, 0x1.31ad012ede8aap+100L, 0x1.2fbf79c839067p-42L, 396208600Srdivacky -0x1.c3e106929056ep+58L, -0x1.64c282b970a5fp+56L, 397208600Srdivacky -0x1.64c282b970a5ep+56L, -0x1.64c282b970a5fp+56L, 398234353Sdim -0x1.64c282b970a5ep+56L, ALL_STD_EXCEPT, FE_INEXACT); 399208600Srdivacky#endif 400208600Srdivacky 401234353Sdim /* x*y (rounded) ~= -z */ 402234353Sdim /* XXX spurious inexact exceptions */ 403201361Srdivacky testrnd(fmaf, 0x1.bbffeep-30, -0x1.1d164cp-74, 0x1.ee7296p-104, 404234353Sdim -0x1.c46ea8p-128, -0x1.c46ea8p-128, -0x1.c46ea8p-128, 405201361Srdivacky -0x1.c46ea8p-128, ALL_STD_EXCEPT & ~FE_INEXACT, 0); 406201361Srdivacky testrnd(fma, 0x1.bbffeea6fc7d6p-30, 0x1.1d164c6cbf078p-74, 407201361Srdivacky -0x1.ee72993aff948p-104, -0x1.71f72ac7d9d8p-159, 408234353Sdim -0x1.71f72ac7d9d8p-159, -0x1.71f72ac7d9d8p-159, 409201361Srdivacky -0x1.71f72ac7d9d8p-159, ALL_STD_EXCEPT & ~FE_INEXACT, 0); 410201361Srdivacky#if LDBL_MANT_DIG == 113 411234353Sdim testrnd(fmal, 0x1.bbffeea6fc7d65927d147f437675p-30L, 412234353Sdim 0x1.1d164c6cbf078b7a22607d1cd6a2p-74L, 413234353Sdim -0x1.ee72993aff94973876031bec0944p-104L, 414206084Srdivacky 0x1.64e086175b3a2adc36e607058814p-217L, 415234353Sdim 0x1.64e086175b3a2adc36e607058814p-217L, 416234353Sdim 0x1.64e086175b3a2adc36e607058814p-217L, 417234353Sdim 0x1.64e086175b3a2adc36e607058814p-217L, 418234353Sdim ALL_STD_EXCEPT & ~FE_INEXACT, 0); 419234353Sdim#elif LDBL_MANT_DIG == 64 420221345Sdim testrnd(fmal, 0x1.bbffeea6fc7d6592p-30L, 0x1.1d164c6cbf078b7ap-74L, 421193326Sed -0x1.ee72993aff949736p-104L, 0x1.af190e7a1ee6ad94p-168L, 422193326Sed 0x1.af190e7a1ee6ad94p-168L, 0x1.af190e7a1ee6ad94p-168L, 423234353Sdim 0x1.af190e7a1ee6ad94p-168L, ALL_STD_EXCEPT & ~FE_INEXACT, 0); 424221345Sdim#elif LDBL_MANT_DIG == 53 425221345Sdim testrnd(fmal, 0x1.bbffeea6fc7d6p-30L, 0x1.1d164c6cbf078p-74L, 426221345Sdim -0x1.ee72993aff948p-104L, -0x1.71f72ac7d9d8p-159L, 427221345Sdim -0x1.71f72ac7d9d8p-159L, -0x1.71f72ac7d9d8p-159L, 428193326Sed -0x1.71f72ac7d9d8p-159L, ALL_STD_EXCEPT & ~FE_INEXACT, 0); 429203955Srdivacky#endif 430203955Srdivacky} 431193326Sed 432193326Sedstatic void 433193326Sedtest_double_rounding(void) 434193326Sed{ 435193326Sed 436193326Sed /* 437234353Sdim * a = 0x1.8000000000001p0 438212904Sdim * b = 0x1.8000000000001p0 439212904Sdim * c = -0x0.0000000000000000000000000080...1p+1 440193326Sed * a * b = 0x1.2000000000001800000000000080p+1 441193326Sed * 442198092Srdivacky * The correct behavior is to round DOWN to 0x1.2000000000001p+1 in 443193326Sed * round-to-nearest mode. An implementation that computes a*b+c in 444198092Srdivacky * double+double precision, however, will get 0x1.20000000000018p+1, 445193326Sed * and then round UP. 446234353Sdim */ 447193326Sed fesetround(FE_TONEAREST); 448193326Sed test(fma, 0x1.8000000000001p0, 0x1.8000000000001p0, 449193326Sed -0x1.0000000000001p-104, 0x1.2000000000001p+1, 450193326Sed ALL_STD_EXCEPT, FE_INEXACT); 451198092Srdivacky fesetround(FE_DOWNWARD); 452193326Sed test(fma, 0x1.8000000000001p0, 0x1.8000000000001p0, 453193326Sed -0x1.0000000000001p-104, 0x1.2000000000001p+1, 454193326Sed ALL_STD_EXCEPT, FE_INEXACT); 455193326Sed fesetround(FE_UPWARD); 456198092Srdivacky test(fma, 0x1.8000000000001p0, 0x1.8000000000001p0, 457234353Sdim -0x1.0000000000001p-104, 0x1.2000000000002p+1, 458234353Sdim ALL_STD_EXCEPT, FE_INEXACT); 459249423Sdim 460234353Sdim fesetround(FE_TONEAREST); 461193326Sed test(fmaf, 0x1.800002p+0, 0x1.800002p+0, -0x1.000002p-46, 0x1.200002p+1, 462203955Srdivacky ALL_STD_EXCEPT, FE_INEXACT); 463210299Sed fesetround(FE_DOWNWARD); 464193326Sed test(fmaf, 0x1.800002p+0, 0x1.800002p+0, -0x1.000002p-46, 0x1.200002p+1, 465193326Sed ALL_STD_EXCEPT, FE_INEXACT); 466210299Sed fesetround(FE_UPWARD); 467210299Sed test(fmaf, 0x1.800002p+0, 0x1.800002p+0, -0x1.000002p-46, 0x1.200004p+1, 468210299Sed ALL_STD_EXCEPT, FE_INEXACT); 469219077Sdim 470221345Sdim fesetround(FE_TONEAREST); 471221345Sdim#if LDBL_MANT_DIG == 64 472221345Sdim test(fmal, 0x1.4p+0L, 0x1.0000000000000004p+0L, 0x1p-128L, 473221345Sdim 0x1.4000000000000006p+0L, ALL_STD_EXCEPT, FE_INEXACT); 474221345Sdim#elif LDBL_MANT_DIG == 113 475210299Sed test(fmal, 0x1.8000000000000000000000000001p+0L, 476221345Sdim 0x1.8000000000000000000000000001p+0L, 477210299Sed -0x1.0000000000000000000000000001p-224L, 478221345Sdim 0x1.2000000000000000000000000001p+1L, ALL_STD_EXCEPT, FE_INEXACT); 479221345Sdim#endif 480221345Sdim 481221345Sdim} 482210299Sed 483210299Sedint 484210299Sedmain(int argc, char *argv[]) 485219077Sdim{ 486221345Sdim int rmodes[] = { FE_TONEAREST, FE_UPWARD, FE_DOWNWARD, FE_TOWARDZERO }; 487221345Sdim int i; 488210299Sed 489210299Sed printf("1..19\n"); 490210299Sed 491210299Sed for (i = 0; i < 4; i++) { 492234353Sdim fesetround(rmodes[i]); 493210299Sed test_zeroes(); 494210299Sed printf("ok %d - fma zeroes\n", i + 1); 495243830Sdim } 496243830Sdim 497210299Sed for (i = 0; i < 4; i++) { 498210299Sed fesetround(rmodes[i]); 499198092Srdivacky test_infinities(); 500200583Srdivacky printf("ok %d - fma infinities\n", i + 5); 501198092Srdivacky } 502205219Srdivacky 503205219Srdivacky fesetround(FE_TONEAREST); 504210299Sed test_nans(); 505205219Srdivacky printf("ok 9 - fma NaNs\n"); 506205219Srdivacky 507198092Srdivacky for (i = 0; i < 4; i++) { 508205219Srdivacky fesetround(rmodes[i]); 509205219Srdivacky test_small_z(); 510221345Sdim printf("ok %d - fma small z\n", i + 10); 511221345Sdim } 512221345Sdim 513221345Sdim for (i = 0; i < 4; i++) { 514205219Srdivacky fesetround(rmodes[i]); 515205219Srdivacky test_big_z(); 516205219Srdivacky printf("ok %d - fma big z\n", i + 14); 517205219Srdivacky } 518198092Srdivacky 519198092Srdivacky fesetround(FE_TONEAREST); 520221345Sdim test_accuracy(); 521221345Sdim printf("ok 18 - fma accuracy\n"); 522221345Sdim 523221345Sdim test_double_rounding(); 524198092Srdivacky printf("ok 19 - fma double rounding\n"); 525198092Srdivacky 526205219Srdivacky /* 527205219Srdivacky * TODO: 528210299Sed * - Tests for subnormals 529205219Srdivacky * - Cancellation tests (e.g., z = (double)x*y, but x*y is inexact) 530205219Srdivacky */ 531205219Srdivacky 532205219Srdivacky return (0); 533210299Sed} 534205219Srdivacky