Mercurial > hg > octave-kai > gnulib-hg
annotate tests/test-fma2.h @ 17415:6550127da196
argp: typo fix
* lib/argp-help.c: Typo in comment.
| author | Alexandre Duret-Lutz <adl@lrde.epita.fr> |
|---|---|
| date | Fri, 17 May 2013 19:01:14 +0200 |
| parents | e542fd46ad6f |
| children |
| rev | line source |
|---|---|
| 16039 | 1 /* Test of fused multiply-add. |
|
17249
e542fd46ad6f
maint: update all copyright year number ranges
Eric Blake <eblake@redhat.com>
parents:
16201
diff
changeset
|
2 Copyright (C) 2011-2013 Free Software Foundation, Inc. |
| 16039 | 3 |
| 4 This program is free software: you can redistribute it and/or modify | |
| 5 it under the terms of the GNU General Public License as published by | |
| 6 the Free Software Foundation; either version 3 of the License, or | |
| 7 (at your option) any later version. | |
| 8 | |
| 9 This program is distributed in the hope that it will be useful, | |
| 10 but WITHOUT ANY WARRANTY; without even the implied warranty of | |
| 11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
| 12 GNU General Public License for more details. | |
| 13 | |
| 14 You should have received a copy of the GNU General Public License | |
| 15 along with this program. If not, see <http://www.gnu.org/licenses/>. */ | |
| 16 | |
| 17 /* Written by Bruno Haible <bruno@clisp.org>, 2011. */ | |
| 18 | |
| 19 /* Returns 2^e as a DOUBLE. */ | |
| 20 #define POW2(e) \ | |
| 21 LDEXP (L_(1.0), e) | |
| 22 | |
| 23 /* One could define XE_RANGE and YE_RANGE to 5 or 60, but this would slow down | |
| 24 the test without the expectation of catching more bugs. */ | |
| 25 #define XE_RANGE 0 | |
| 26 #define YE_RANGE 0 | |
| 27 | |
| 28 /* Define to 1 if you want to allow the behaviour of the 'double-double' | |
| 29 implementation of 'long double' (seen on IRIX 6.5 and Linux/PowerPC). | |
| 30 This floating-point type does not follow IEEE 754. */ | |
| 31 #if MANT_BIT == LDBL_MANT_BIT && LDBL_MANT_BIT == 2 * DBL_MANT_BIT | |
| 32 # define FORGIVE_DOUBLEDOUBLE_BUG 1 | |
| 33 #else | |
| 34 # define FORGIVE_DOUBLEDOUBLE_BUG 0 | |
| 35 #endif | |
| 36 | |
| 37 /* Subnormal numbers appear to not work as expected on IRIX 6.5. */ | |
| 38 #ifdef __sgi | |
| 39 # define MIN_SUBNORMAL_EXP (MIN_EXP - 1) | |
| 40 #else | |
| 41 # define MIN_SUBNORMAL_EXP (MIN_EXP - MANT_BIT) | |
| 42 #endif | |
| 43 | |
| 44 /* Check rounding behaviour. */ | |
| 45 | |
| 46 static void | |
| 47 test_function (DOUBLE (*my_fma) (DOUBLE, DOUBLE, DOUBLE)) | |
| 48 { | |
| 49 /* Array mapping n to (-1)^n. */ | |
| 50 static const DOUBLE pow_m1[] = | |
| 51 { | |
| 52 L_(1.0), - L_(1.0), L_(1.0), - L_(1.0), | |
| 53 L_(1.0), - L_(1.0), L_(1.0), - L_(1.0) | |
| 54 }; | |
| 55 | |
| 56 /* A product x * y that consists of two bits. */ | |
| 57 { | |
|
16110
146b5e66e142
fma tests: Avoid shadowing local variables.
Bruno Haible <bruno@clisp.org>
parents:
16044
diff
changeset
|
58 volatile DOUBLE x; |
|
146b5e66e142
fma tests: Avoid shadowing local variables.
Bruno Haible <bruno@clisp.org>
parents:
16044
diff
changeset
|
59 volatile DOUBLE y; |
|
146b5e66e142
fma tests: Avoid shadowing local variables.
Bruno Haible <bruno@clisp.org>
parents:
16044
diff
changeset
|
60 volatile DOUBLE z; |
|
146b5e66e142
fma tests: Avoid shadowing local variables.
Bruno Haible <bruno@clisp.org>
parents:
16044
diff
changeset
|
61 volatile DOUBLE result; |
|
146b5e66e142
fma tests: Avoid shadowing local variables.
Bruno Haible <bruno@clisp.org>
parents:
16044
diff
changeset
|
62 volatile DOUBLE expected; |
| 16039 | 63 int xs; |
| 64 int xe; | |
| 65 int ys; | |
| 66 int ye; | |
| 67 int ze; | |
| 68 DOUBLE sign; | |
| 69 | |
| 70 for (xs = 0; xs < 2; xs++) | |
| 71 for (xe = -XE_RANGE; xe <= XE_RANGE; xe++) | |
| 72 { | |
| 73 x = pow_m1[xs] * POW2 (xe); /* (-1)^xs * 2^xe */ | |
| 74 | |
| 75 for (ys = 0; ys < 2; ys++) | |
| 76 for (ye = -YE_RANGE; ye <= YE_RANGE; ye++) | |
| 77 { | |
| 78 y = pow_m1[ys] * POW2 (ye); /* (-1)^ys * 2^ye */ | |
| 79 | |
| 80 sign = pow_m1[xs + ys]; | |
| 81 | |
| 82 /* Test addition (same signs). */ | |
| 83 for (ze = MIN_EXP - MANT_BIT; ze <= MAX_EXP - 1;) | |
| 84 { | |
| 85 z = sign * POW2 (ze); /* (-1)^(xs+ys) * 2^ze */ | |
| 86 result = my_fma (x, y, z); | |
| 87 if (xe + ye >= ze + MANT_BIT) | |
| 88 expected = sign * POW2 (xe + ye); | |
| 89 else if (xe + ye > ze - MANT_BIT) | |
| 90 expected = sign * (POW2 (xe + ye) + POW2 (ze)); | |
| 91 else | |
| 92 expected = z; | |
| 93 ASSERT (result == expected); | |
| 94 | |
| 95 ze++; | |
| 96 /* Shortcut some values of ze, to speed up the test. */ | |
| 97 if (ze == MIN_EXP + MANT_BIT) | |
| 98 ze = - 2 * MANT_BIT - 1; | |
| 99 else if (ze == 2 * MANT_BIT) | |
| 100 ze = MAX_EXP - MANT_BIT - 1; | |
| 101 } | |
| 102 | |
| 103 /* Test subtraction (opposite signs). */ | |
| 104 for (ze = MIN_EXP - MANT_BIT; ze <= MAX_EXP - 1;) | |
| 105 { | |
| 106 z = - sign * POW2 (ze); /* (-1)^(xs+ys+1) * 2^ze */ | |
| 107 result = my_fma (x, y, z); | |
| 108 if (xe + ye > ze + MANT_BIT) | |
| 109 expected = sign * POW2 (xe + ye); | |
| 110 else if (xe + ye >= ze) | |
| 111 expected = sign * (POW2 (xe + ye) - POW2 (ze)); | |
| 112 else if (xe + ye > ze - 1 - MANT_BIT) | |
| 113 expected = - sign * (POW2 (ze) - POW2 (xe + ye)); | |
| 114 else | |
| 115 expected = z; | |
| 116 ASSERT (result == expected); | |
| 117 | |
| 118 ze++; | |
| 119 /* Shortcut some values of ze, to speed up the test. */ | |
| 120 if (ze == MIN_EXP + MANT_BIT) | |
| 121 ze = - 2 * MANT_BIT - 1; | |
| 122 else if (ze == 2 * MANT_BIT) | |
| 123 ze = MAX_EXP - MANT_BIT - 1; | |
| 124 } | |
| 125 } | |
| 126 } | |
| 127 } | |
| 128 /* A product x * y that consists of three bits. */ | |
| 129 { | |
|
16110
146b5e66e142
fma tests: Avoid shadowing local variables.
Bruno Haible <bruno@clisp.org>
parents:
16044
diff
changeset
|
130 volatile DOUBLE x; |
|
146b5e66e142
fma tests: Avoid shadowing local variables.
Bruno Haible <bruno@clisp.org>
parents:
16044
diff
changeset
|
131 volatile DOUBLE y; |
|
146b5e66e142
fma tests: Avoid shadowing local variables.
Bruno Haible <bruno@clisp.org>
parents:
16044
diff
changeset
|
132 volatile DOUBLE z; |
|
146b5e66e142
fma tests: Avoid shadowing local variables.
Bruno Haible <bruno@clisp.org>
parents:
16044
diff
changeset
|
133 volatile DOUBLE result; |
|
146b5e66e142
fma tests: Avoid shadowing local variables.
Bruno Haible <bruno@clisp.org>
parents:
16044
diff
changeset
|
134 volatile DOUBLE expected; |
| 16039 | 135 int i; |
| 136 int xs; | |
| 137 int xe; | |
| 138 int ys; | |
| 139 int ye; | |
| 140 int ze; | |
| 141 DOUBLE sign; | |
| 142 | |
| 143 for (i = 1; i <= MANT_BIT - 1; i++) | |
| 144 for (xs = 0; xs < 2; xs++) | |
| 145 for (xe = -XE_RANGE; xe <= XE_RANGE; xe++) | |
| 146 { | |
| 147 x = /* (-1)^xs * (2^xe + 2^(xe-i)) */ | |
| 148 pow_m1[xs] * (POW2 (xe) + POW2 (xe - i)); | |
| 149 | |
| 150 for (ys = 0; ys < 2; ys++) | |
| 151 for (ye = -YE_RANGE; ye <= YE_RANGE; ye++) | |
| 152 { | |
| 153 y = /* (-1)^ys * (2^ye + 2^(ye-i)) */ | |
| 154 pow_m1[ys] * (POW2 (ye) + POW2 (ye - i)); | |
| 155 | |
| 156 sign = pow_m1[xs + ys]; | |
| 157 | |
| 158 /* The exact value of x * y is | |
| 159 (-1)^(xs+ys) * (2^(xe+ye) + 2^(xe+ye-i+1) + 2^(xe+ye-2*i)) */ | |
| 160 | |
| 161 /* Test addition (same signs). */ | |
| 162 for (ze = MIN_SUBNORMAL_EXP; ze <= MAX_EXP - 1;) | |
| 163 { | |
| 164 z = sign * POW2 (ze); /* (-1)^(xs+ys) * 2^ze */ | |
| 165 result = my_fma (x, y, z); | |
| 166 if (FORGIVE_DOUBLEDOUBLE_BUG) | |
| 167 if ((xe + ye > ze | |
| 168 && xe + ye < ze + MANT_BIT | |
| 169 && i == DBL_MANT_BIT) | |
| 170 || (xe + ye == ze + DBL_MANT_BIT && i == DBL_MANT_BIT + 1) | |
| 171 || (xe + ye == ze + MANT_BIT - 1 && i == 1)) | |
| 172 goto skip1; | |
| 173 if (xe + ye > ze + MANT_BIT) | |
| 174 { | |
| 175 if (2 * i > MANT_BIT) | |
| 176 expected = | |
| 177 sign * (POW2 (xe + ye) | |
| 178 + POW2 (xe + ye - i + 1)); | |
| 179 else if (2 * i == MANT_BIT) | |
| 16044 | 180 expected = |
| 181 sign * (POW2 (xe + ye) | |
| 182 + POW2 (xe + ye - i + 1) | |
| 183 + POW2 (xe + ye - MANT_BIT + 1)); | |
| 16039 | 184 else |
| 185 expected = | |
| 186 sign * (POW2 (xe + ye) | |
| 187 + POW2 (xe + ye - i + 1) | |
| 188 + POW2 (xe + ye - 2 * i)); | |
| 189 } | |
| 190 else if (xe + ye == ze + MANT_BIT) | |
| 191 { | |
| 192 if (2 * i >= MANT_BIT) | |
| 16044 | 193 expected = |
| 194 sign * (POW2 (xe + ye) | |
| 195 + POW2 (xe + ye - i + 1) | |
| 196 + POW2 (xe + ye - MANT_BIT + 1)); | |
| 16039 | 197 else if (2 * i == MANT_BIT - 1) |
| 198 /* round-to-even rounds up */ | |
| 199 expected = | |
| 200 sign * (POW2 (xe + ye) | |
| 201 + POW2 (xe + ye - i + 1) | |
| 202 + POW2 (xe + ye - 2 * i + 1)); | |
| 203 else | |
| 204 expected = | |
| 205 sign * (POW2 (xe + ye) | |
| 206 + POW2 (xe + ye - i + 1) | |
| 207 + POW2 (xe + ye - 2 * i)); | |
| 208 } | |
| 209 else if (xe + ye > ze - MANT_BIT + 2 * i) | |
| 16044 | 210 expected = |
| 211 sign * (POW2 (ze) | |
| 212 + POW2 (xe + ye) | |
| 213 + POW2 (xe + ye - i + 1) | |
| 214 + POW2 (xe + ye - 2 * i)); | |
| 16039 | 215 else if (xe + ye >= ze - MANT_BIT + i) |
| 216 expected = | |
| 217 sign * (POW2 (ze) | |
| 218 + POW2 (xe + ye) | |
| 219 + POW2 (xe + ye - i + 1)); | |
| 220 else if (xe + ye == ze - MANT_BIT + i - 1) | |
| 221 { | |
| 222 if (i == 1) | |
| 223 expected = | |
| 224 sign * (POW2 (ze) + POW2 (ze - MANT_BIT + 1)); | |
| 225 else | |
| 226 expected = | |
| 227 sign * (POW2 (ze) | |
| 228 + POW2 (xe + ye) | |
| 229 + POW2 (ze - MANT_BIT + 1)); | |
| 230 } | |
| 231 else if (xe + ye >= ze - MANT_BIT + 1) | |
| 232 expected = sign * (POW2 (ze) + POW2 (xe + ye)); | |
| 233 else if (xe + ye == ze - MANT_BIT) | |
| 234 expected = | |
| 235 sign * (POW2 (ze) + POW2 (ze - MANT_BIT + 1)); | |
| 236 else if (xe + ye == ze - MANT_BIT - 1) | |
| 237 { | |
| 238 if (i == 1) | |
| 16044 | 239 expected = |
| 240 sign * (POW2 (ze) + POW2 (ze - MANT_BIT + 1)); | |
| 16039 | 241 else |
| 242 expected = z; | |
| 243 } | |
| 244 else | |
| 245 expected = z; | |
| 246 ASSERT (result == expected); | |
| 247 | |
| 248 skip1: | |
| 249 ze++; | |
| 250 /* Shortcut some values of ze, to speed up the test. */ | |
| 251 if (ze == MIN_EXP + MANT_BIT) | |
| 252 ze = - 2 * MANT_BIT - 1; | |
| 253 else if (ze == 2 * MANT_BIT) | |
| 254 ze = MAX_EXP - MANT_BIT - 1; | |
| 255 } | |
| 256 | |
| 257 /* Test subtraction (opposite signs). */ | |
| 258 if (i > 1) | |
| 259 for (ze = MIN_SUBNORMAL_EXP; ze <= MAX_EXP - 1;) | |
| 260 { | |
| 261 z = - sign * POW2 (ze); /* (-1)^(xs+ys+1) * 2^ze */ | |
| 262 result = my_fma (x, y, z); | |
| 263 if (FORGIVE_DOUBLEDOUBLE_BUG) | |
| 264 if ((xe + ye == ze && i == MANT_BIT - 1) | |
| 265 || (xe + ye > ze | |
| 266 && xe + ye <= ze + DBL_MANT_BIT - 1 | |
| 267 && i == DBL_MANT_BIT + 1) | |
| 268 || (xe + ye >= ze + DBL_MANT_BIT - 1 | |
| 269 && xe + ye < ze + MANT_BIT | |
| 270 && xe + ye == ze + i - 1) | |
| 271 || (xe + ye > ze + DBL_MANT_BIT | |
| 272 && xe + ye < ze + MANT_BIT | |
| 273 && i == DBL_MANT_BIT)) | |
| 274 goto skip2; | |
| 275 if (xe + ye == ze) | |
| 276 { | |
| 277 /* maximal extinction */ | |
| 278 expected = | |
| 279 sign * (POW2 (xe + ye - i + 1) | |
| 280 + POW2 (xe + ye - 2 * i)); | |
| 281 } | |
| 282 else if (xe + ye == ze - 1) | |
| 283 { | |
| 284 /* significant extinction */ | |
| 285 if (2 * i > MANT_BIT) | |
| 286 expected = | |
| 287 sign * (- POW2 (xe + ye) | |
| 288 + POW2 (xe + ye - i + 1)); | |
| 289 else | |
| 16044 | 290 expected = |
| 291 sign * (- POW2 (xe + ye) | |
| 292 + POW2 (xe + ye - i + 1) | |
| 293 + POW2 (xe + ye - 2 * i)); | |
| 16039 | 294 } |
| 295 else if (xe + ye > ze + MANT_BIT) | |
| 296 { | |
| 297 if (2 * i >= MANT_BIT) | |
| 298 expected = | |
| 299 sign * (POW2 (xe + ye) | |
| 300 + POW2 (xe + ye - i + 1)); | |
| 301 else | |
| 302 expected = | |
| 303 sign * (POW2 (xe + ye) | |
| 304 + POW2 (xe + ye - i + 1) | |
| 305 + POW2 (xe + ye - 2 * i)); | |
| 306 } | |
| 307 else if (xe + ye == ze + MANT_BIT) | |
| 308 { | |
| 309 if (2 * i >= MANT_BIT) | |
| 310 expected = | |
| 311 sign * (POW2 (xe + ye) | |
| 312 + POW2 (xe + ye - i + 1)); | |
| 313 else if (2 * i == MANT_BIT - 1) | |
| 314 /* round-to-even rounds down */ | |
| 315 expected = | |
| 316 sign * (POW2 (xe + ye) | |
| 317 + POW2 (xe + ye - i + 1)); | |
| 318 else | |
| 319 /* round-to-even rounds up */ | |
| 320 expected = | |
| 321 sign * (POW2 (xe + ye) | |
| 322 + POW2 (xe + ye - i + 1) | |
| 323 + POW2 (xe + ye - 2 * i)); | |
| 324 } | |
| 325 else if (xe + ye >= ze - MANT_BIT + 2 * i) | |
| 16044 | 326 expected = |
| 327 sign * (- POW2 (ze) | |
| 328 + POW2 (xe + ye) | |
| 329 + POW2 (xe + ye - i + 1) | |
| 330 + POW2 (xe + ye - 2 * i)); | |
| 16039 | 331 else if (xe + ye >= ze - MANT_BIT + i - 1) |
| 332 expected = | |
| 333 sign * (- POW2 (ze) | |
| 334 + POW2 (xe + ye) | |
| 335 + POW2 (xe + ye - i + 1)); | |
| 336 else if (xe + ye == ze - MANT_BIT + i - 2) | |
| 16044 | 337 expected = |
| 338 sign * (- POW2 (ze) | |
| 339 + POW2 (xe + ye) | |
| 340 + POW2 (ze - MANT_BIT)); | |
| 16039 | 341 else if (xe + ye >= ze - MANT_BIT) |
| 342 expected = | |
| 343 sign * (- POW2 (ze) | |
| 344 + POW2 (xe + ye)); | |
| 345 else if (xe + ye == ze - MANT_BIT - 1) | |
| 16044 | 346 expected = |
| 347 sign * (- POW2 (ze) | |
| 348 + POW2 (ze - MANT_BIT)); | |
| 16039 | 349 else |
| 350 expected = z; | |
| 351 ASSERT (result == expected); | |
| 352 | |
| 353 skip2: | |
| 354 ze++; | |
| 355 /* Shortcut some values of ze, to speed up the test. */ | |
| 356 if (ze == MIN_EXP + MANT_BIT) | |
| 357 ze = - 2 * MANT_BIT - 1; | |
| 358 else if (ze == 2 * MANT_BIT) | |
| 359 ze = MAX_EXP - MANT_BIT - 1; | |
| 360 } | |
| 361 } | |
| 362 } | |
| 363 } | |
| 364 /* TODO: Similar tests with | |
| 365 x = (-1)^xs * (2^xe - 2^(xe-i)), y = (-1)^ys * (2^ye - 2^(ye-i)) */ | |
| 366 /* A product x * y that consists of one segment of bits (or, if you prefer, | |
| 367 two bits, one with positive weight and one with negative weight). */ | |
| 368 { | |
|
16110
146b5e66e142
fma tests: Avoid shadowing local variables.
Bruno Haible <bruno@clisp.org>
parents:
16044
diff
changeset
|
369 volatile DOUBLE x; |
|
146b5e66e142
fma tests: Avoid shadowing local variables.
Bruno Haible <bruno@clisp.org>
parents:
16044
diff
changeset
|
370 volatile DOUBLE y; |
|
146b5e66e142
fma tests: Avoid shadowing local variables.
Bruno Haible <bruno@clisp.org>
parents:
16044
diff
changeset
|
371 volatile DOUBLE z; |
|
146b5e66e142
fma tests: Avoid shadowing local variables.
Bruno Haible <bruno@clisp.org>
parents:
16044
diff
changeset
|
372 volatile DOUBLE result; |
|
146b5e66e142
fma tests: Avoid shadowing local variables.
Bruno Haible <bruno@clisp.org>
parents:
16044
diff
changeset
|
373 volatile DOUBLE expected; |
| 16039 | 374 int i; |
| 375 int xs; | |
| 376 int xe; | |
| 377 int ys; | |
| 378 int ye; | |
| 379 int ze; | |
| 380 DOUBLE sign; | |
| 381 | |
| 382 for (i = 1; i <= MANT_BIT - 1; i++) | |
| 383 for (xs = 0; xs < 2; xs++) | |
| 384 for (xe = -XE_RANGE; xe <= XE_RANGE; xe++) | |
| 385 { | |
| 386 x = /* (-1)^xs * (2^xe + 2^(xe-i)) */ | |
| 387 pow_m1[xs] * (POW2 (xe) + POW2 (xe - i)); | |
| 388 | |
| 389 for (ys = 0; ys < 2; ys++) | |
| 390 for (ye = -YE_RANGE; ye <= YE_RANGE; ye++) | |
| 391 { | |
| 392 y = /* (-1)^ys * (2^ye - 2^(ye-i)) */ | |
| 393 pow_m1[ys] * (POW2 (ye) - POW2 (ye - i)); | |
| 394 | |
| 395 sign = pow_m1[xs + ys]; | |
| 396 | |
| 397 /* The exact value of x * y is | |
| 398 (-1)^(xs+ys) * (2^(xe+ye) - 2^(xe+ye-2*i)) */ | |
| 399 | |
| 400 /* Test addition (same signs). */ | |
| 401 for (ze = MIN_EXP - MANT_BIT; ze <= MAX_EXP - 1;) | |
| 402 { | |
| 403 z = sign * POW2 (ze); /* (-1)^(xs+ys) * 2^ze */ | |
| 404 result = my_fma (x, y, z); | |
| 405 if (FORGIVE_DOUBLEDOUBLE_BUG) | |
| 406 if ((xe + ye == ze + MANT_BIT && i > DBL_MANT_BIT) | |
| 407 || (xe + ye < ze + MANT_BIT | |
| 408 && xe + ye >= ze | |
| 409 && i == DBL_MANT_BIT) | |
| 410 || (xe + ye < ze | |
| 411 && xe + ye == ze - MANT_BIT + 2 * i)) | |
| 412 goto skip3; | |
| 413 if (xe + ye > ze + MANT_BIT + 1) | |
| 414 { | |
| 415 if (2 * i > MANT_BIT) | |
| 416 expected = sign * POW2 (xe + ye); | |
| 417 else | |
| 418 expected = | |
| 419 sign * (POW2 (xe + ye) | |
| 420 - POW2 (xe + ye - 2 * i)); | |
| 421 } | |
| 422 else if (xe + ye == ze + MANT_BIT + 1) | |
| 423 { | |
| 424 if (2 * i >= MANT_BIT) | |
| 425 expected = sign * POW2 (xe + ye); | |
| 426 else | |
| 427 expected = | |
| 428 sign * (POW2 (xe + ye) | |
| 429 - POW2 (xe + ye - 2 * i)); | |
| 430 } | |
| 431 else if (xe + ye >= ze - MANT_BIT + 2 * i) | |
| 432 { | |
| 433 if (2 * i > MANT_BIT) | |
| 434 expected = | |
| 435 sign * (POW2 (xe + ye) + POW2 (ze)); | |
| 436 else | |
| 437 expected = | |
| 438 sign * (POW2 (xe + ye) | |
| 439 - POW2 (xe + ye - 2 * i) | |
| 440 + POW2 (ze)); | |
| 441 } | |
| 442 else if (xe + ye >= ze - MANT_BIT + 1) | |
| 443 expected = | |
| 444 sign * (POW2 (ze) + POW2 (xe + ye)); | |
| 445 else | |
| 446 expected = z; | |
| 447 ASSERT (result == expected); | |
| 448 | |
| 449 skip3: | |
| 450 ze++; | |
| 451 /* Shortcut some values of ze, to speed up the test. */ | |
| 452 if (ze == MIN_EXP + MANT_BIT) | |
| 453 ze = - 2 * MANT_BIT - 1; | |
| 454 else if (ze == 2 * MANT_BIT) | |
| 455 ze = MAX_EXP - MANT_BIT - 1; | |
| 456 } | |
| 457 | |
| 458 /* Test subtraction (opposite signs). */ | |
| 459 if (i > 1) | |
| 460 for (ze = MIN_SUBNORMAL_EXP; ze <= MAX_EXP - 1;) | |
| 461 { | |
| 462 z = - sign * POW2 (ze); /* (-1)^(xs+ys+1) * 2^ze */ | |
| 463 result = my_fma (x, y, z); | |
| 464 if (FORGIVE_DOUBLEDOUBLE_BUG) | |
| 465 if (xe + ye > ze | |
| 466 && xe + ye < ze + DBL_MANT_BIT | |
| 467 && xe + ye == ze + 2 * i - LDBL_MANT_BIT) | |
| 468 goto skip4; | |
| 469 if (xe + ye == ze) | |
| 470 { | |
| 471 /* maximal extinction */ | |
| 472 expected = sign * - POW2 (xe + ye - 2 * i); | |
| 473 } | |
| 474 else if (xe + ye > ze + MANT_BIT + 1) | |
| 475 { | |
| 476 if (2 * i > MANT_BIT + 1) | |
| 477 expected = sign * POW2 (xe + ye); | |
| 478 else if (2 * i == MANT_BIT + 1) | |
| 16044 | 479 expected = |
| 480 sign * (POW2 (xe + ye) | |
| 481 - POW2 (xe + ye - MANT_BIT)); | |
| 16039 | 482 else |
| 483 expected = | |
| 484 sign * (POW2 (xe + ye) | |
| 485 - POW2 (xe + ye - 2 * i)); | |
| 486 } | |
| 487 else if (xe + ye == ze + MANT_BIT + 1) | |
| 488 { | |
| 489 if (2 * i > MANT_BIT) | |
| 16044 | 490 expected = |
| 491 sign * (POW2 (xe + ye) | |
| 492 - POW2 (xe + ye - MANT_BIT)); | |
| 16039 | 493 else if (2 * i == MANT_BIT) |
| 494 expected = | |
| 495 sign * (POW2 (xe + ye) | |
| 496 - POW2 (xe + ye - MANT_BIT + 1)); | |
| 497 else | |
| 498 expected = | |
| 499 sign * (POW2 (xe + ye) | |
| 500 - POW2 (xe + ye - 2 * i)); | |
| 501 } | |
| 502 else if (xe + ye == ze + MANT_BIT) | |
| 503 { | |
| 504 if (2 * i > MANT_BIT + 1) | |
| 505 expected = | |
| 506 sign * (POW2 (xe + ye) | |
| 507 - POW2 (xe + ye - MANT_BIT)); | |
| 508 else if (2 * i == MANT_BIT + 1) | |
| 509 expected = | |
| 510 sign * (POW2 (xe + ye) | |
| 511 - POW2 (xe + ye - MANT_BIT + 1)); | |
| 512 else | |
| 513 expected = | |
| 514 sign * (POW2 (xe + ye) | |
| 515 - POW2 (ze) | |
| 516 - POW2 (xe + ye - 2 * i)); | |
| 517 } | |
| 518 else if (xe + ye > ze - MANT_BIT + 2 * i) | |
| 519 { | |
| 520 if (2 * i > MANT_BIT) | |
| 521 expected = | |
| 522 sign * (POW2 (xe + ye) - POW2 (ze)); | |
| 523 else | |
| 524 expected = | |
| 525 sign * (POW2 (xe + ye) | |
| 526 - POW2 (ze) | |
| 527 - POW2 (xe + ye - 2 * i)); | |
| 528 } | |
| 529 else if (xe + ye == ze - MANT_BIT + 2 * i) | |
| 530 expected = | |
| 531 sign * (- POW2 (ze) | |
| 532 + POW2 (xe + ye) | |
| 533 - POW2 (xe + ye - 2 * i)); | |
| 16044 | 534 else if (xe + ye >= ze - MANT_BIT) |
| 16039 | 535 expected = sign * (- POW2 (ze) + POW2 (xe + ye)); |
| 536 else | |
| 537 expected = z; | |
| 538 ASSERT (result == expected); | |
| 539 | |
| 540 skip4: | |
| 541 ze++; | |
| 542 /* Shortcut some values of ze, to speed up the test. */ | |
| 543 if (ze == MIN_EXP + MANT_BIT) | |
| 544 ze = - 2 * MANT_BIT - 1; | |
| 545 else if (ze == 2 * MANT_BIT) | |
| 546 ze = MAX_EXP - MANT_BIT - 1; | |
| 547 } | |
| 548 } | |
| 549 } | |
| 550 } | |
| 551 /* TODO: Tests with denormalized results. */ | |
| 552 /* Tests with temporary overflow. */ | |
| 553 { | |
| 554 volatile DOUBLE x = POW2 (MAX_EXP - 1); | |
| 555 volatile DOUBLE y = POW2 (MAX_EXP - 1); | |
| 556 volatile DOUBLE z = - INFINITY; | |
| 557 volatile DOUBLE result = my_fma (x, y, z); | |
| 558 ASSERT (result == - INFINITY); | |
| 559 } | |
| 560 { | |
| 561 volatile DOUBLE x = POW2 (MAX_EXP - 1); /* 2^(MAX_EXP-1) */ | |
| 562 volatile DOUBLE y = L_(2.0); /* 2^1 */ | |
| 563 volatile DOUBLE z = /* -(2^MAX_EXP - 2^(MAX_EXP-MANT_BIT)) */ | |
| 564 - LDEXP (POW2 (MAX_EXP - 1) - POW2 (MAX_EXP - MANT_BIT - 1), 1); | |
| 565 volatile DOUBLE result = my_fma (x, y, z); | |
| 566 if (!FORGIVE_DOUBLEDOUBLE_BUG) | |
| 567 ASSERT (result == POW2 (MAX_EXP - MANT_BIT)); | |
| 568 } | |
| 569 { | |
| 570 volatile DOUBLE x = POW2 (MAX_EXP - 1); /* 2^(MAX_EXP-1) */ | |
| 571 volatile DOUBLE y = L_(3.0); /* 3 */ | |
| 572 volatile DOUBLE z = - LDEXP (L_(5.0), MAX_EXP - 3); /* -5*2^(MAX_EXP-3) */ | |
| 573 volatile DOUBLE result = my_fma (x, y, z); | |
| 574 ASSERT (result == LDEXP (L_(7.0), MAX_EXP - 3)); | |
| 575 } | |
| 576 } |