Sourceware Bugzilla – Attachment 6014 Details for
Bug 13304
fma, fmaf, fmal produce wrong results
Home
|
New
|
Browse
|
Search
|
[?]
|
Reports
|
Requests
|
Help
|
New Account
|
Log In
[x]
|
Forgot Password
Login:
[x]
test suite for fma(), version 2
fma-testsuite.c (text/x-csrc), 22.12 KB, created by
Bruno Haible
on 2011-10-17 22:39:55 UTC
(
hide
)
Description:
test suite for fma(), version 2
Filename:
MIME Type:
Creator:
Bruno Haible
Created:
2011-10-17 22:39:55 UTC
Size:
22.12 KB
patch
obsolete
>/* Test of fma(). > Copyright (C) 2011 Free Software Foundation, Inc. > > This program is free software: you can redistribute it and/or modify it > under the terms of the GNU Lesser General Public License as published > by the Free Software Foundation; either version 3 of the License, or > (at your option) any later version. > > This program is distributed in the hope that it will be useful, > but WITHOUT ANY WARRANTY; without even the implied warranty of > MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU > Lesser General Public License for more details. > > You should have received a copy of the GNU Lesser General Public License > along with this program. If not, see <http://www.gnu.org/licenses/>. */ > >/* Written by Bruno Haible <bruno@clisp.org>, 2011. */ > >#include <math.h> >#include <float.h> >#include <stdio.h> >#include <stdlib.h> > >#define DOUBLE double >#define ISNAN isnand >#define LDEXP ldexp >#define MIN_EXP DBL_MIN_EXP >#define MAX_EXP DBL_MAX_EXP >#define MANT_BIT DBL_MANT_DIG >#define L_(literal) literal > >#define ASSERT(expr) \ > do \ > { \ > if (!(expr)) \ > { \ > fprintf (stderr, "%s:%d: assertion failed\n", \ > __FILE__, __LINE__); \ > fflush (stderr); \ > abort (); \ > } \ > } \ > while (0) > >/* Returns 2^e as a DOUBLE. */ >#define POW2(e) \ > LDEXP (L_(1.0), e) > >#define XE_RANGE 2 >#define YE_RANGE 2 > >static void >test_function (DOUBLE (*my_fma) (DOUBLE, DOUBLE, DOUBLE)) >{ > /* Array mapping n to (-1)^n. */ > static const DOUBLE pow_m1[] = > { > L_(1.0), - L_(1.0), L_(1.0), - L_(1.0), > L_(1.0), - L_(1.0), L_(1.0), - L_(1.0) > }; > volatile DOUBLE x; > volatile DOUBLE y; > volatile DOUBLE z; > volatile DOUBLE result; > volatile DOUBLE expected; > > /* A product x * y that consists of two bits. */ > { > int xs; > int xe; > int ys; > int ye; > int ze; > DOUBLE sign; > > for (xs = 0; xs < 2; xs++) > for (xe = -XE_RANGE; xe <= XE_RANGE; xe++) > { > x = pow_m1[xs] * POW2 (xe); /* (-1)^xs * 2^xe */ > > for (ys = 0; ys < 2; ys++) > for (ye = -YE_RANGE; ye <= YE_RANGE; ye++) > { > y = pow_m1[ys] * POW2 (ye); /* (-1)^ys * 2^ye */ > > sign = pow_m1[xs + ys]; > > /* Test addition (same signs). */ > for (ze = MIN_EXP - MANT_BIT; ze <= MAX_EXP - 1;) > { > z = sign * POW2 (ze); /* (-1)^(xs+ys) * 2^ze */ > result = my_fma (x, y, z); > if (xe + ye >= ze + MANT_BIT) > expected = sign * POW2 (xe + ye); > else if (xe + ye > ze - MANT_BIT) > expected = sign * (POW2 (xe + ye) + POW2 (ze)); > else > expected = z; > ASSERT (result == expected); > > ze++; > /* Shortcut some values of ze, to speed up the test. */ > if (ze == MIN_EXP + MANT_BIT) > ze = - 2 * MANT_BIT - 1; > else if (ze == 2 * MANT_BIT) > ze = MAX_EXP - MANT_BIT - 1; > } > > /* Test subtraction (opposite signs). */ > for (ze = MIN_EXP - MANT_BIT; ze <= MAX_EXP - 1;) > { > z = - sign * POW2 (ze); /* (-1)^(xs+ys+1) * 2^ze */ > result = my_fma (x, y, z); > if (xe + ye > ze + MANT_BIT) > expected = sign * POW2 (xe + ye); > else if (xe + ye >= ze) > expected = sign * (POW2 (xe + ye) - POW2 (ze)); > else if (xe + ye > ze - 1 - MANT_BIT) > expected = - sign * (POW2 (ze) - POW2 (xe + ye)); > else > expected = z; > ASSERT (result == expected); > > ze++; > /* Shortcut some values of ze, to speed up the test. */ > if (ze == MIN_EXP + MANT_BIT) > ze = - 2 * MANT_BIT - 1; > else if (ze == 2 * MANT_BIT) > ze = MAX_EXP - MANT_BIT - 1; > } > } > } > } > /* A product x * y that consists of three bits. */ > { > int i; > int xs; > int xe; > int ys; > int ye; > int ze; > DOUBLE sign; > > for (i = 1; i <= MANT_BIT - 1; i++) > for (xs = 0; xs < 2; xs++) > for (xe = -XE_RANGE; xe <= XE_RANGE; xe++) > { > x = /* (-1)^xs * (2^xe + 2^(xe-i)) */ > pow_m1[xs] * (POW2 (xe) + POW2 (xe - i)); > > for (ys = 0; ys < 2; ys++) > for (ye = -YE_RANGE; ye <= YE_RANGE; ye++) > { > y = /* (-1)^ys * (2^ye + 2^(ye-i)) */ > pow_m1[ys] * (POW2 (ye) + POW2 (ye - i)); > > sign = pow_m1[xs + ys]; > > /* The exact value of x * y is > (-1)^(xs+ys) * (2^(xe+ye) + 2^(xe+ye-i+1) + 2^(xe+ye-2*i)) */ > > /* Test addition (same signs). */ > for (ze = MIN_EXP - MANT_BIT; ze <= MAX_EXP - 1;) > { > z = sign * POW2 (ze); /* (-1)^(xs+ys) * 2^ze */ > result = my_fma (x, y, z); > if (xe + ye > ze + MANT_BIT) > { > if (2 * i > MANT_BIT) > expected = > sign * (POW2 (xe + ye) > + POW2 (xe + ye - i + 1)); > else if (2 * i == MANT_BIT) > expected = > sign * (POW2 (xe + ye) > + POW2 (xe + ye - i + 1) > + POW2 (xe + ye - MANT_BIT + 1)); > else > expected = > sign * (POW2 (xe + ye) > + POW2 (xe + ye - i + 1) > + POW2 (xe + ye - 2 * i)); > } > else if (xe + ye == ze + MANT_BIT) > { > if (2 * i >= MANT_BIT) > expected = > sign * (POW2 (xe + ye) > + POW2 (xe + ye - i + 1) > + POW2 (xe + ye - MANT_BIT + 1)); > else if (2 * i == MANT_BIT - 1) > /* round-to-even rounds up */ > expected = > sign * (POW2 (xe + ye) > + POW2 (xe + ye - i + 1) > + POW2 (xe + ye - 2 * i + 1)); > else > expected = > sign * (POW2 (xe + ye) > + POW2 (xe + ye - i + 1) > + POW2 (xe + ye - 2 * i)); > } > else if (xe + ye > ze - MANT_BIT + 2 * i) > expected = > sign * (POW2 (ze) > + POW2 (xe + ye) > + POW2 (xe + ye - i + 1) > + POW2 (xe + ye - 2 * i)); > else if (xe + ye >= ze - MANT_BIT + i) > expected = > sign * (POW2 (ze) > + POW2 (xe + ye) > + POW2 (xe + ye - i + 1)); > else if (xe + ye == ze - MANT_BIT + i - 1) > expected = > sign * (POW2 (ze) > + POW2 (xe + ye) > + POW2 (ze - MANT_BIT + 1)); > else if (xe + ye >= ze - MANT_BIT + 1) > expected = sign * (POW2 (ze) + POW2 (xe + ye)); > else if (xe + ye == ze - MANT_BIT) > expected = > sign * (POW2 (ze) + POW2 (ze - MANT_BIT + 1)); > else if (xe + ye == ze - MANT_BIT - 1) > { > if (i == 1) > expected = > sign * (POW2 (ze) + POW2 (ze - MANT_BIT + 1)); > else > expected = z; > } > else > expected = z; > ASSERT (result == expected); > > ze++; > /* Shortcut some values of ze, to speed up the test. */ > if (ze == MIN_EXP + MANT_BIT) > ze = - 2 * MANT_BIT - 1; > else if (ze == 2 * MANT_BIT) > ze = MAX_EXP - MANT_BIT - 1; > } > > /* Test subtraction (opposite signs). */ > if (i > 1) > for (ze = MIN_EXP - MANT_BIT; ze <= MAX_EXP - 1;) > { > z = - sign * POW2 (ze); /* (-1)^(xs+ys+1) * 2^ze */ > result = my_fma (x, y, z); > if (xe + ye == ze) > /* maximal extinction */ > expected = > sign * (POW2 (xe + ye - i + 1) > + POW2 (xe + ye - 2 * i)); > else if (xe + ye == ze - 1) > { > /* significant extinction */ > if (2 * i > MANT_BIT) > expected = > sign * (- POW2 (xe + ye) > + POW2 (xe + ye - i + 1)); > else > expected = > sign * (- POW2 (xe + ye) > + POW2 (xe + ye - i + 1) > + POW2 (xe + ye - 2 * i)); > } > else if (xe + ye > ze + MANT_BIT) > { > if (2 * i >= MANT_BIT) > expected = > sign * (POW2 (xe + ye) > + POW2 (xe + ye - i + 1)); > else > expected = > sign * (POW2 (xe + ye) > + POW2 (xe + ye - i + 1) > + POW2 (xe + ye - 2 * i)); > } > else if (xe + ye == ze + MANT_BIT) > { > if (2 * i >= MANT_BIT) > expected = > sign * (POW2 (xe + ye) > + POW2 (xe + ye - i + 1)); > else if (2 * i == MANT_BIT - 1) > /* round-to-even rounds down */ > expected = > sign * (POW2 (xe + ye) > + POW2 (xe + ye - i + 1)); > else > /* round-to-even rounds up */ > expected = > sign * (POW2 (xe + ye) > + POW2 (xe + ye - i + 1) > + POW2 (xe + ye - 2 * i)); > } > else if (xe + ye >= ze - MANT_BIT + 2 * i) > expected = > sign * (- POW2 (ze) > + POW2 (xe + ye) > + POW2 (xe + ye - i + 1) > + POW2 (xe + ye - 2 * i)); > else if (xe + ye >= ze - MANT_BIT + i - 1) > expected = > sign * (- POW2 (ze) > + POW2 (xe + ye) > + POW2 (xe + ye - i + 1)); > else if (xe + ye == ze - MANT_BIT + i - 2) > expected = > sign * (- POW2 (ze) > + POW2 (xe + ye) > + POW2 (ze - MANT_BIT)); > else if (xe + ye >= ze - MANT_BIT) > expected = > sign * (- POW2 (ze) > + POW2 (xe + ye)); > else if (xe + ye == ze - MANT_BIT - 1) > expected = > sign * (- POW2 (ze) > + POW2 (ze - MANT_BIT)); > else > expected = z; > ASSERT (result == expected); > > ze++; > /* Shortcut some values of ze, to speed up the test. */ > if (ze == MIN_EXP + MANT_BIT) > ze = - 2 * MANT_BIT - 1; > else if (ze == 2 * MANT_BIT) > ze = MAX_EXP - MANT_BIT - 1; > } > } > } > } > /* TODO: Similar tests with > x = (-1)^xs * (2^xe - 2^(xe-i)), y = (-1)^ys * (2^ye - 2^(ye-i)) */ > /* A product x * y that consists of one segment of bits (or, if you prefer, > two bits, one with positive weight and one with negative weight). */ > { > int i; > int xs; > int xe; > int ys; > int ye; > int ze; > DOUBLE sign; > > for (i = 1; i <= MANT_BIT - 1; i++) > for (xs = 0; xs < 2; xs++) > for (xe = -XE_RANGE; xe <= XE_RANGE; xe++) > { > x = /* (-1)^xs * (2^xe + 2^(xe-i)) */ > pow_m1[xs] * (POW2 (xe) + POW2 (xe - i)); > > for (ys = 0; ys < 2; ys++) > for (ye = -YE_RANGE; ye <= YE_RANGE; ye++) > { > y = /* (-1)^ys * (2^ye - 2^(ye-i)) */ > pow_m1[ys] * (POW2 (ye) - POW2 (ye - i)); > > sign = pow_m1[xs + ys]; > > /* The exact value of x * y is > (-1)^(xs+ys) * (2^(xe+ye) - 2^(xe+ye-2*i)) */ > > /* Test addition (same signs). */ > for (ze = MIN_EXP - MANT_BIT; ze <= MAX_EXP - 1;) > { > z = sign * POW2 (ze); /* (-1)^(xs+ys) * 2^ze */ > result = my_fma (x, y, z); > if (xe + ye > ze + MANT_BIT + 1) > { > if (2 * i > MANT_BIT) > expected = sign * POW2 (xe + ye); > else > expected = > sign * (POW2 (xe + ye) > - POW2 (xe + ye - 2 * i)); > } > else if (xe + ye == ze + MANT_BIT + 1) > { > if (2 * i >= MANT_BIT) > expected = sign * POW2 (xe + ye); > else > expected = > sign * (POW2 (xe + ye) > - POW2 (xe + ye - 2 * i)); > } > else if (xe + ye >= ze - MANT_BIT + 2 * i) > { > if (2 * i > MANT_BIT) > expected = > sign * (POW2 (xe + ye) + POW2 (ze)); > else > expected = > sign * (POW2 (xe + ye) > - POW2 (xe + ye - 2 * i) > + POW2 (ze)); > } > else if (xe + ye >= ze - MANT_BIT + 1) > expected = > sign * (POW2 (ze) + POW2 (xe + ye)); > else > expected = z; > ASSERT (result == expected); > > ze++; > /* Shortcut some values of ze, to speed up the test. */ > if (ze == MIN_EXP + MANT_BIT) > ze = - 2 * MANT_BIT - 1; > else if (ze == 2 * MANT_BIT) > ze = MAX_EXP - MANT_BIT - 1; > } > > /* Test subtraction (opposite signs). */ > if (i > 1) > for (ze = MIN_EXP - MANT_BIT; ze <= MAX_EXP - 1;) > { > z = - sign * POW2 (ze); /* (-1)^(xs+ys+1) * 2^ze */ > result = my_fma (x, y, z); > if (xe + ye == ze) > /* maximal extinction */ > expected = sign * - POW2 (xe + ye - 2 * i); > else if (xe + ye > ze + MANT_BIT + 1) > { > if (2 * i > MANT_BIT + 1) > expected = sign * POW2 (xe + ye); > else if (2 * i == MANT_BIT + 1) > expected = > sign * (POW2 (xe + ye) > - POW2 (xe + ye - MANT_BIT)); > else > expected = > sign * (POW2 (xe + ye) > - POW2 (xe + ye - 2 * i)); > } > else if (xe + ye == ze + MANT_BIT + 1) > { > if (2 * i > MANT_BIT) > expected = > sign * (POW2 (xe + ye) > - POW2 (xe + ye - MANT_BIT)); > else if (2 * i == MANT_BIT) > expected = > sign * (POW2 (xe + ye) > - POW2 (xe + ye - MANT_BIT + 1)); > else > expected = > sign * (POW2 (xe + ye) > - POW2 (xe + ye - 2 * i)); > } > else if (xe + ye == ze + MANT_BIT) > { > if (2 * i > MANT_BIT + 1) > expected = > sign * (POW2 (xe + ye) > - POW2 (xe + ye - MANT_BIT)); > else if (2 * i == MANT_BIT + 1) > expected = > sign * (POW2 (xe + ye) > - POW2 (xe + ye - MANT_BIT + 1)); > else > expected = > sign * (POW2 (xe + ye) > - POW2 (ze) > - POW2 (xe + ye - 2 * i)); > } > else if (xe + ye > ze - MANT_BIT + 2 * i) > { > if (2 * i > MANT_BIT) > expected = > sign * (POW2 (xe + ye) - POW2 (ze)); > else > expected = > sign * (POW2 (xe + ye) > - POW2 (ze) > - POW2 (xe + ye - 2 * i)); > } > else if (xe + ye == ze - MANT_BIT + 2 * i) > expected = > sign * (- POW2 (ze) > + POW2 (xe + ye) > - POW2 (xe + ye - 2 * i)); > else if (xe + ye > ze - MANT_BIT) > expected = sign * (- POW2 (ze) + POW2 (xe + ye)); > else if (xe + ye == ze - MANT_BIT) > expected = sign * (- POW2 (ze) + POW2 (xe + ye)); > else > expected = z; > ASSERT (result == expected); > > ze++; > /* Shortcut some values of ze, to speed up the test. */ > if (ze == MIN_EXP + MANT_BIT) > ze = - 2 * MANT_BIT - 1; > else if (ze == 2 * MANT_BIT) > ze = MAX_EXP - MANT_BIT - 1; > } > } > } > } > /* TODO: Tests with denormalized results. */ > /* Tests with temporary overflow. */ > { > volatile DOUBLE x = POW2 (MAX_EXP - 1); > volatile DOUBLE y = POW2 (MAX_EXP - 1); > volatile DOUBLE z = - INFINITY; > volatile DOUBLE result = my_fma (x, y, z); > ASSERT (result == - INFINITY); > } > { > volatile DOUBLE x = POW2 (MAX_EXP - 1); /* 2^(MAX_EXP-1) */ > volatile DOUBLE y = L_(2.0); /* 2^1 */ > volatile DOUBLE z = /* -(2^MAX_EXP - 2^(MAX_EXP-MANT_BIT)) */ > - LDEXP (POW2 (MAX_EXP - 1) - POW2 (MAX_EXP - MANT_BIT - 1), 1); > volatile DOUBLE result = my_fma (x, y, z); > ASSERT (result == POW2 (MAX_EXP - MANT_BIT)); > } > { > volatile DOUBLE x = POW2 (MAX_EXP - 1); /* 2^(MAX_EXP-1) */ > volatile DOUBLE y = L_(3.0); /* 3 */ > volatile DOUBLE z = - LDEXP (L_(5.0), MAX_EXP - 3); /* -5*2^(MAX_EXP-3) */ > volatile DOUBLE result = my_fma (x, y, z); > ASSERT (result == LDEXP (L_(7.0), MAX_EXP - 3)); > } >} > >int >main () >{ > test_function (fma); > > return 0; >}
<b>You cannot view the attachment while viewing its details because your browser does not support IFRAMEs. <a href="attachment.cgi?id=6014">View the attachment on a separate page</a>.</b>
View Attachment As Raw
Actions:
View
Attachments on
bug 13304
:
5990
|
5991
|
5992
|
5993
|
5994
|
5995
|
5996
|
5997
|
5998
|
5999
|
6000
|
6001
|
6002
|
6003
|
6004
|
6005
|
6006
|
6008
| 6014 |
6015
|
6017
|
6657