Annotation of OpenXM_contrib/gmp/mpfr/tests/tadd.c, Revision 1.1
1.1 ! ohara 1: /* Test file for mpfr_add and mpfr_sub.
! 2:
! 3: Copyright 1999, 2000, 2001, 2002 Free Software Foundation.
! 4:
! 5: This file is part of the MPFR Library.
! 6:
! 7: The MPFR Library is free software; you can redistribute it and/or modify
! 8: it under the terms of the GNU Lesser General Public License as published by
! 9: the Free Software Foundation; either version 2.1 of the License, or (at your
! 10: option) any later version.
! 11:
! 12: The MPFR Library is distributed in the hope that it will be useful, but
! 13: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
! 14: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
! 15: License for more details.
! 16:
! 17: You should have received a copy of the GNU Lesser General Public License
! 18: along with the MPFR Library; see the file COPYING.LIB. If not, write to
! 19: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
! 20: MA 02111-1307, USA. */
! 21:
! 22: #define N 100000
! 23:
! 24: #include <stdio.h>
! 25: #include <stdlib.h>
! 26: #include <time.h>
! 27: #include "gmp.h"
! 28: #include "mpfr.h"
! 29: #include "mpfr-impl.h"
! 30: #include "mpfr-test.h"
! 31:
! 32: void checknan _PROTO((double, double, mp_rnd_t, unsigned int, unsigned int, unsigned int));
! 33: void check3 _PROTO((double, double, mp_rnd_t));
! 34: void check4 _PROTO((double, double, mp_rnd_t));
! 35: void check5 _PROTO((double, mp_rnd_t));
! 36: void check2 _PROTO((double, int, double, int, int, int));
! 37: void check2a _PROTO((double, int, double, int, int, int, char *));
! 38: void check64 _PROTO((void));
! 39: void check_same _PROTO((void));
! 40: void check_case_1b _PROTO((void));
! 41: void check_case_2 _PROTO((void));
! 42: void check_inexact _PROTO((void));
! 43:
! 44:
! 45: /* Parameter "z1" of check() used to be last in the argument list, but that
! 46: tickled a bug in 32-bit sparc gcc 2.95.2. A "double" in that position is
! 47: passed on the stack at an address which is 4mod8, but the generated code
! 48: didn't take into account that alignment, resulting in bus errors. The
! 49: easiest workaround is to move it to the start of the arg list (where it's
! 50: passed in registers), this macro does that. FIXME: Change the actual
! 51: calls to check(), rather than using a macro. */
! 52:
! 53: #define check(x,y,rnd_mode,px,py,pz,z1) _check(x,y,z1,rnd_mode,px,py,pz)
! 54:
! 55: /* checks that x+y gives the same results in double
! 56: and with mpfr with 53 bits of precision */
! 57: void
! 58: _check (double x, double y, double z1, mp_rnd_t rnd_mode, unsigned int px,
! 59: unsigned int py, unsigned int pz)
! 60: {
! 61: double z2; mpfr_t xx,yy,zz; int cert=0;
! 62:
! 63: mpfr_init2(xx, px);
! 64: mpfr_init2(yy, py);
! 65: mpfr_init2(zz, pz);
! 66: mpfr_set_d(xx, x, rnd_mode);
! 67: mpfr_set_d(yy, y, rnd_mode);
! 68: mpfr_add(zz, xx, yy, rnd_mode);
! 69: #ifdef MPFR_HAVE_FESETROUND
! 70: mpfr_set_machine_rnd_mode(rnd_mode);
! 71: if (px==53 && py==53 && pz==53) cert=1;
! 72: #endif
! 73: if (z1==0.0) z1=x+y; else cert=1;
! 74: z2 = mpfr_get_d1 (zz);
! 75: mpfr_set_d (yy, z2, GMP_RNDN);
! 76: if (!mpfr_cmp (zz, yy) && cert && z1!=z2 && !(isnan(z1) && isnan(z2))) {
! 77: printf("expected sum is %1.20e, got %1.20e\n",z1,z2);
! 78: printf("mpfr_add failed for x=%1.20e y=%1.20e with rnd_mode=%s\n",
! 79: x, y, mpfr_print_rnd_mode(rnd_mode));
! 80: exit(1);
! 81: }
! 82: mpfr_clear(xx); mpfr_clear(yy); mpfr_clear(zz);
! 83: }
! 84:
! 85: void
! 86: checknan (double x, double y, mp_rnd_t rnd_mode, unsigned int px,
! 87: unsigned int py, unsigned int pz)
! 88: {
! 89: double z2; mpfr_t xx, yy, zz;
! 90:
! 91: mpfr_init2(xx, px);
! 92: mpfr_init2(yy, py);
! 93: mpfr_init2(zz, pz);
! 94: mpfr_set_d(xx, x, rnd_mode);
! 95: mpfr_set_d(yy, y, rnd_mode);
! 96: mpfr_add(zz, xx, yy, rnd_mode);
! 97: #ifdef MPFR_HAVE_FESETROUND
! 98: mpfr_set_machine_rnd_mode(rnd_mode);
! 99: #endif
! 100: if (MPFR_IS_NAN(zz) == 0) { printf("Error, not an MPFR_NAN for xx = %1.20e, y = %1.20e\n", x, y); exit(1); }
! 101: z2 = mpfr_get_d1 (zz);
! 102: if (!isnan(z2)) { printf("Error, not a NaN after conversion, xx = %1.20e yy = %1.20e, got %1.20e\n", x, y, z2); exit(1); }
! 103:
! 104: mpfr_clear(xx); mpfr_clear(yy); mpfr_clear(zz);
! 105: }
! 106:
! 107: #ifdef MPFR_HAVE_FESETROUND
! 108: /* idem than check for mpfr_add(x, x, y) */
! 109: void
! 110: check3 (double x, double y, mp_rnd_t rnd_mode)
! 111: {
! 112: double z1,z2; mpfr_t xx,yy; int neg;
! 113:
! 114: neg = LONG_RAND() % 2;
! 115: mpfr_init2(xx, 53);
! 116: mpfr_init2(yy, 53);
! 117: mpfr_set_d(xx, x, rnd_mode);
! 118: mpfr_set_d(yy, y, rnd_mode);
! 119: if (neg) mpfr_sub(xx, xx, yy, rnd_mode);
! 120: else mpfr_add(xx, xx, yy, rnd_mode);
! 121: mpfr_set_machine_rnd_mode(rnd_mode);
! 122: z1 = (neg) ? x-y : x+y;
! 123: z2 = mpfr_get_d1 (xx);
! 124: mpfr_set_d (yy, z2, GMP_RNDN);
! 125: if (!mpfr_cmp (xx, yy) && z1!=z2 && !(isnan(z1) && isnan(z2))) {
! 126: printf("expected result is %1.20e, got %1.20e\n",z1,z2);
! 127: printf("mpfr_%s(x,x,y) failed for x=%1.20e y=%1.20e with rnd_mode=%u\n",
! 128: (neg) ? "sub" : "add",x,y,rnd_mode);
! 129: exit(1);
! 130: }
! 131: mpfr_clear(xx); mpfr_clear(yy);
! 132: }
! 133:
! 134: /* idem than check for mpfr_add(x, y, x) */
! 135: void
! 136: check4 (double x, double y, mp_rnd_t rnd_mode)
! 137: {
! 138: double z1, z2;
! 139: mpfr_t xx, yy;
! 140: int neg;
! 141:
! 142: neg = LONG_RAND() % 2;
! 143: mpfr_init2(xx, 53);
! 144: mpfr_init2(yy, 53);
! 145: mpfr_set_d(xx, x, rnd_mode);
! 146: mpfr_set_d(yy, y, rnd_mode);
! 147: if (neg) mpfr_sub(xx, yy, xx, rnd_mode);
! 148: else mpfr_add(xx, yy, xx, rnd_mode);
! 149: mpfr_set_machine_rnd_mode(rnd_mode);
! 150: z1 = (neg) ? y-x : x+y;
! 151: z2 = mpfr_get_d1 (xx);
! 152: mpfr_set_d (yy, z2, GMP_RNDN);
! 153: /* check that xx is representable as a double and no overflow occurred */
! 154: if ((mpfr_cmp (xx, yy) == 0) && (z1 != z2)) {
! 155: printf("expected result is %1.20e, got %1.20e\n", z1, z2);
! 156: printf("mpfr_%s(x,y,x) failed for x=%1.20e y=%1.20e with rnd_mode=%s\n",
! 157: (neg) ? "sub" : "add", x, y, mpfr_print_rnd_mode(rnd_mode));
! 158: exit(1);
! 159: }
! 160: mpfr_clear(xx); mpfr_clear(yy);
! 161: }
! 162:
! 163: /* idem than check for mpfr_add(x, x, x) */
! 164: void
! 165: check5 (double x, mp_rnd_t rnd_mode)
! 166: {
! 167: double z1,z2; mpfr_t xx, yy; int neg;
! 168:
! 169: mpfr_init2(xx, 53);
! 170: mpfr_init2(yy, 53);
! 171: neg = LONG_RAND() % 2;
! 172: mpfr_set_d(xx, x, rnd_mode);
! 173: if (neg) mpfr_sub(xx, xx, xx, rnd_mode);
! 174: else mpfr_add(xx, xx, xx, rnd_mode);
! 175: mpfr_set_machine_rnd_mode(rnd_mode);
! 176: z1 = (neg) ? x-x : x+x;
! 177: z2 = mpfr_get_d1 (xx);
! 178: mpfr_set_d (yy, z2, GMP_RNDN);
! 179: /* check NaNs first since mpfr_cmp does not like them */
! 180: if (!(isnan(z1) && isnan(z2)) && !mpfr_cmp (xx, yy) && z1!=z2)
! 181: {
! 182: printf ("expected result is %1.20e, got %1.20e\n",z1,z2);
! 183: printf ("mpfr_%s(x,x,x) failed for x=%1.20e with rnd_mode=%s\n",
! 184: (neg) ? "sub" : "add", x, mpfr_print_rnd_mode (rnd_mode));
! 185: exit (1);
! 186: }
! 187: mpfr_clear(xx);
! 188: mpfr_clear(yy);
! 189: }
! 190:
! 191: void
! 192: check2 (double x, int px, double y, int py, int pz, mp_rnd_t rnd_mode)
! 193: {
! 194: mpfr_t xx, yy, zz; double z,z2; int u;
! 195:
! 196: mpfr_init2(xx,px); mpfr_init2(yy,py); mpfr_init2(zz,pz);
! 197: mpfr_set_d(xx, x, rnd_mode);
! 198: mpfr_set_d(yy, y, rnd_mode);
! 199: mpfr_add(zz, xx, yy, rnd_mode);
! 200: mpfr_set_machine_rnd_mode(rnd_mode);
! 201: z = x+y; z2=mpfr_get_d1 (zz); u=ulp(z,z2);
! 202: /* one ulp difference is possible due to composed rounding */
! 203: if (px>=53 && py>=53 && pz>=53 && ABS(u)>1) {
! 204: printf("x=%1.20e,%d y=%1.20e,%d pz=%d,rnd=%s\n",
! 205: x,px,y,py,pz,mpfr_print_rnd_mode(rnd_mode));
! 206: printf("got %1.20e\n",z2);
! 207: printf("result should be %1.20e (diff=%d ulp)\n",z,u);
! 208: mpfr_set_d(zz, z, rnd_mode);
! 209: printf("i.e."); mpfr_print_binary(zz); putchar('\n');
! 210: exit(1); }
! 211: mpfr_clear(xx); mpfr_clear(yy); mpfr_clear(zz);
! 212: }
! 213: #endif
! 214:
! 215: void
! 216: check2a (double x, int px, double y, int py, int pz, mp_rnd_t rnd_mode,
! 217: char *res)
! 218: {
! 219: mpfr_t xx, yy, zz;
! 220:
! 221: mpfr_init2(xx,px); mpfr_init2(yy,py); mpfr_init2(zz,pz);
! 222: mpfr_set_d(xx, x, rnd_mode);
! 223: mpfr_set_d(yy, y, rnd_mode);
! 224: mpfr_add(zz, xx, yy, rnd_mode);
! 225: mpfr_set_prec(xx, pz);
! 226: mpfr_set_str(xx, res, 16, GMP_RNDN);
! 227: if (mpfr_cmp(xx, zz)) {
! 228: printf("x=%1.20e,%d y=%1.20e,%d pz=%d,rnd=%s\n",
! 229: x,px,y,py,pz,mpfr_print_rnd_mode(rnd_mode));
! 230: printf("got "); mpfr_print_binary(zz); putchar('\n');
! 231: printf("instead of "); mpfr_print_binary(xx); putchar('\n');
! 232: exit(1);
! 233: }
! 234: mpfr_clear(xx); mpfr_clear(yy); mpfr_clear(zz);
! 235: }
! 236:
! 237: void
! 238: check64 (void)
! 239: {
! 240: mpfr_t x, t, u;
! 241:
! 242: mpfr_init (x);
! 243: mpfr_init (t);
! 244: mpfr_init (u);
! 245:
! 246: mpfr_set_prec (x, 29);
! 247: mpfr_set_str_raw (x, "1.1101001000101111011010010110e-3");
! 248: mpfr_set_prec (t, 58);
! 249: mpfr_set_str_raw (t, "0.11100010011111001001100110010111110110011000000100101E-1");
! 250: mpfr_set_prec (u, 29);
! 251: mpfr_add (u, x, t, GMP_RNDD);
! 252: mpfr_set_str_raw (t, "1.0101011100001000011100111110e-1");
! 253: if (mpfr_cmp (u, t))
! 254: {
! 255: fprintf (stderr, "mpfr_add(u, x, t) failed for prec(x)=29, prec(t)=58\n");
! 256: printf ("expected "); mpfr_out_str (stdout, 2, 29, t, GMP_RNDN);
! 257: putchar ('\n');
! 258: printf ("got "); mpfr_out_str (stdout, 2, 29, u, GMP_RNDN);
! 259: putchar ('\n');
! 260: exit(1);
! 261: }
! 262:
! 263: mpfr_set_prec (x, 4);
! 264: mpfr_set_str_raw (x, "-1.0E-2");
! 265: mpfr_set_prec (t, 2);
! 266: mpfr_set_str_raw (t, "-1.1e-2");
! 267: mpfr_set_prec (u, 2);
! 268: mpfr_add (u, x, t, GMP_RNDN);
! 269: if (MPFR_MANT(u)[0] << 2)
! 270: {
! 271: fprintf (stderr, "result not normalized for prec=2\n");
! 272: mpfr_print_binary (u); putchar ('\n');
! 273: exit (1);
! 274: }
! 275: mpfr_set_str_raw (t, "-1.0e-1");
! 276: if (mpfr_cmp (u, t))
! 277: {
! 278: fprintf (stderr, "mpfr_add(u, x, t) failed for prec(x)=4, prec(t)=2\n");
! 279: printf ("expected -1.0e-1\n");
! 280: printf ("got "); mpfr_out_str (stdout, 2, 4, u, GMP_RNDN);
! 281: putchar ('\n');
! 282: exit (1);
! 283: }
! 284:
! 285: mpfr_set_prec (x, 8);
! 286: mpfr_set_str_raw (x, "-0.10011010"); /* -77/128 */
! 287: mpfr_set_prec (t, 4);
! 288: mpfr_set_str_raw (t, "-1.110e-5"); /* -7/128 */
! 289: mpfr_set_prec (u, 4);
! 290: mpfr_add (u, x, t, GMP_RNDN); /* should give -5/8 */
! 291: mpfr_set_str_raw (t, "-1.010e-1");
! 292: if (mpfr_cmp (u, t)) {
! 293: fprintf (stderr, "mpfr_add(u, x, t) failed for prec(x)=8, prec(t)=4\n");
! 294: printf ("expected -1.010e-1\n");
! 295: printf ("got "); mpfr_out_str (stdout, 2, 4, u, GMP_RNDN);
! 296: putchar ('\n');
! 297: exit (1);
! 298: }
! 299:
! 300: mpfr_set_prec (x, 112); mpfr_set_prec (t, 98); mpfr_set_prec (u, 54);
! 301: mpfr_set_str_raw (x, "-0.11111100100000000011000011100000101101010001000111E-401");
! 302: mpfr_set_str_raw (t, "0.10110000100100000101101100011111111011101000111000101E-464");
! 303: mpfr_add (u, x, t, GMP_RNDN);
! 304: if (mpfr_cmp (u, x)) {
! 305: fprintf (stderr, "mpfr_add(u, x, t) failed for prec(x)=112, prec(t)=98\n");
! 306: exit (1);
! 307: }
! 308:
! 309: mpfr_set_prec (x, 92); mpfr_set_prec (t, 86); mpfr_set_prec (u, 53);
! 310: mpfr_set_d (x, -5.03525136761487735093e-74, GMP_RNDN);
! 311: mpfr_set_d (t, 8.51539046314262304109e-91, GMP_RNDN);
! 312: mpfr_add (u, x, t, GMP_RNDN);
! 313: if (mpfr_get_d1 (u) != -5.0352513676148773509283672e-74) {
! 314: fprintf (stderr, "mpfr_add(u, x, t) failed for prec(x)=92, prec(t)=86\n");
! 315: exit (1);
! 316: }
! 317:
! 318: mpfr_set_prec(x, 53); mpfr_set_prec(t, 76); mpfr_set_prec(u, 76);
! 319: mpfr_set_str_raw(x, "-0.10010010001001011011110000000000001010011011011110001E-32");
! 320: mpfr_set_str_raw(t, "-0.1011000101110010000101111111011111010001110011110111100110101011110010011111");
! 321: mpfr_sub(u, x, t, GMP_RNDU);
! 322: mpfr_set_str_raw(t, "0.1011000101110010000101111111011100111111101010011011110110101011101000000100");
! 323: if (mpfr_cmp(u,t)) {
! 324: printf("expect "); mpfr_print_binary(t); putchar('\n');
! 325: fprintf (stderr, "mpfr_add failed for precisions 53-76\n"); exit(1);
! 326: }
! 327: mpfr_set_prec(x, 53); mpfr_set_prec(t, 108); mpfr_set_prec(u, 108);
! 328: mpfr_set_str_raw(x, "-0.10010010001001011011110000000000001010011011011110001E-32");
! 329: mpfr_set_str_raw(t, "-0.101100010111001000010111111101111101000111001111011110011010101111001001111000111011001110011000000000111111");
! 330: mpfr_sub(u, x, t, GMP_RNDU);
! 331: mpfr_set_str_raw(t, "0.101100010111001000010111111101110011111110101001101111011010101110100000001011000010101110011000000000111111");
! 332: if (mpfr_cmp(u,t)) {
! 333: printf("expect "); mpfr_print_binary(t); putchar('\n');
! 334: fprintf(stderr, "mpfr_add failed for precisions 53-108\n"); exit(1);
! 335: }
! 336: mpfr_set_prec(x, 97); mpfr_set_prec(t, 97); mpfr_set_prec(u, 97);
! 337: mpfr_set_str_raw(x, "0.1111101100001000000001011000110111101000001011111000100001000101010100011111110010000000000000000E-39");
! 338: mpfr_set_ui(t, 1, GMP_RNDN);
! 339: mpfr_add(u, x, t, GMP_RNDN);
! 340: mpfr_set_str_raw(x, "0.1000000000000000000000000000000000000000111110110000100000000101100011011110100000101111100010001E1");
! 341: if (mpfr_cmp(u,x)) {
! 342: fprintf(stderr, "mpfr_add failed for precision 97\n"); exit(1);
! 343: }
! 344: mpfr_set_prec(x, 128); mpfr_set_prec(t, 128); mpfr_set_prec(u, 128);
! 345: mpfr_set_str_raw(x, "0.10101011111001001010111011001000101100111101000000111111111011010100001100011101010001010111111101111010100110111111100101100010E-4");
! 346: mpfr_set(t, x, GMP_RNDN);
! 347: mpfr_sub(u, x, t, GMP_RNDN);
! 348: mpfr_set_prec(x, 96); mpfr_set_prec(t, 96); mpfr_set_prec(u, 96);
! 349: mpfr_set_str_raw(x, "0.111000000001110100111100110101101001001010010011010011100111100011010100011001010011011011000010E-4");
! 350: mpfr_set(t, x, GMP_RNDN);
! 351: mpfr_sub(u, x, t, GMP_RNDN);
! 352: mpfr_set_prec(x, 85); mpfr_set_prec(t, 85); mpfr_set_prec(u, 85);
! 353: mpfr_set_str_raw(x, "0.1111101110100110110110100010101011101001100010100011110110110010010011101100101111100E-4");
! 354: mpfr_set_str_raw(t, "0.1111101110100110110110100010101001001000011000111000011101100101110100001110101010110E-4");
! 355: mpfr_sub(u, x, t, GMP_RNDU);
! 356: mpfr_sub(x, x, t, GMP_RNDU);
! 357: if (mpfr_cmp(x, u) != 0) {
! 358: printf("Error in mpfr_sub: u=x-t and x=x-t give different results\n");
! 359: exit(1);
! 360: }
! 361: if ((MPFR_MANT(u)[(MPFR_PREC(u)-1)/mp_bits_per_limb] &
! 362: ((mp_limb_t)1<<(mp_bits_per_limb-1)))==0) {
! 363: printf("Error in mpfr_sub: result is not msb-normalized (1)\n"); exit(1);
! 364: }
! 365: mpfr_set_prec(x, 65); mpfr_set_prec(t, 65); mpfr_set_prec(u, 65);
! 366: mpfr_set_str_raw(x, "0.10011010101000110101010000000011001001001110001011101011111011101E623");
! 367: mpfr_set_str_raw(t, "0.10011010101000110101010000000011001001001110001011101011111011100E623");
! 368: mpfr_sub(u, x, t, GMP_RNDU);
! 369: if (mpfr_get_d1 (u) != 9.4349060620538533806e167) { /* 2^558 */
! 370: printf("Error (1) in mpfr_sub\n"); exit(1);
! 371: }
! 372:
! 373: mpfr_set_prec(x, 64); mpfr_set_prec(t, 64); mpfr_set_prec(u, 64);
! 374: mpfr_set_str_raw(x, "0.1000011110101111011110111111000011101011101111101101101100000100E-220");
! 375: mpfr_set_str_raw(t, "0.1000011110101111011110111111000011101011101111101101010011111101E-220");
! 376: mpfr_add(u, x, t, GMP_RNDU);
! 377: if ((MPFR_MANT(u)[0] & 1) != 1) {
! 378: printf("error in mpfr_add with rnd_mode=GMP_RNDU\n");
! 379: printf("b= "); mpfr_print_binary(x); putchar('\n');
! 380: printf("c= "); mpfr_print_binary(t); putchar('\n');
! 381: printf("b+c="); mpfr_print_binary(u); putchar('\n');
! 382: exit(1);
! 383: }
! 384:
! 385: /* bug found by Norbert Mueller, 14 Sep 2000 */
! 386: mpfr_set_prec(x, 56); mpfr_set_prec(t, 83); mpfr_set_prec(u, 10);
! 387: mpfr_set_str_raw(x, "0.10001001011011001111101100110100000101111010010111010111E-7");
! 388: mpfr_set_str_raw(t, "0.10001001011011001111101100110100000101111010010111010111000000000111110110110000100E-7");
! 389: mpfr_sub(u, x, t, GMP_RNDU);
! 390:
! 391: /* array bound write found by Norbert Mueller, 26 Sep 2000 */
! 392: mpfr_set_prec(x, 109); mpfr_set_prec(t, 153); mpfr_set_prec(u, 95);
! 393: mpfr_set_str_raw(x,"0.1001010000101011101100111000110001111111111111111111111111111111111111111111111111111111111111100000000000000E33");
! 394: mpfr_set_str_raw(t,"-0.100101000010101110110011100011000000000000000000000000000000000000000000000000000000000000000000000000000000000000000011100101101000000100100001100110111E33");
! 395: mpfr_add(u, x, t, GMP_RNDN);
! 396:
! 397: /* array bound writes found by Norbert Mueller, 27 Sep 2000 */
! 398: mpfr_set_prec(x, 106); mpfr_set_prec(t, 53); mpfr_set_prec(u, 23);
! 399: mpfr_set_str_raw(x, "-0.1000011110101111111001010001000100001011000000000000000000000000000000000000000000000000000000000000000000E-59");
! 400: mpfr_set_str_raw(t, "-0.10000111101011111110010100010001101100011100110100000E-59");
! 401: mpfr_sub(u, x, t, GMP_RNDN);
! 402: mpfr_set_prec(x, 177); mpfr_set_prec(t, 217); mpfr_set_prec(u, 160);
! 403: mpfr_set_str_raw(x, "-0.111010001011010000111001001010010000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E35");
! 404: mpfr_set_str_raw(t, "0.1110100010110100001110010010100100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000111011010011100001111001E35");
! 405: mpfr_add(u, x, t, GMP_RNDN);
! 406: mpfr_set_prec(x, 214); mpfr_set_prec(t, 278); mpfr_set_prec(u, 207);
! 407: mpfr_set_str_raw(x, "0.1000100110100110101101101101000000010000100111000001001110001000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E66");
! 408: mpfr_set_str_raw(t, "-0.10001001101001101011011011010000000100001001110000010011100010000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001111011111001001100011E66");
! 409: mpfr_add(u, x, t, GMP_RNDN);
! 410: mpfr_set_prec(x, 32); mpfr_set_prec(t, 247); mpfr_set_prec(u, 223);
! 411: mpfr_set_str_raw(x, "0.10000000000000000000000000000000E1");
! 412: mpfr_set_str_raw(t, "0.1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111100000110001110100000100011110000101110110011101110100110110111111011010111100100000000000000000000000000E0");
! 413: mpfr_sub(u, x, t, GMP_RNDN);
! 414: if ((MPFR_MANT(u)[(MPFR_PREC(u)-1)/mp_bits_per_limb] &
! 415: ((mp_limb_t)1<<(mp_bits_per_limb-1)))==0) {
! 416: printf("Error in mpfr_sub: result is not msb-normalized (2)\n"); exit(1);
! 417: }
! 418:
! 419: /* bug found by Nathalie Revol, 21 March 2001 */
! 420: mpfr_set_prec (x, 65);
! 421: mpfr_set_prec (t, 65);
! 422: mpfr_set_prec (u, 65);
! 423: mpfr_set_str_raw (x, "0.11100100101101001100111011111111110001101001000011101001001010010E-35");
! 424: mpfr_set_str_raw (t, "0.10000000000000000000000000000000000001110010010110100110011110000E1");
! 425: mpfr_sub (u, t, x, GMP_RNDU);
! 426: if ((MPFR_MANT(u)[(MPFR_PREC(u)-1)/mp_bits_per_limb] &
! 427: ((mp_limb_t)1<<(mp_bits_per_limb-1)))==0) {
! 428: fprintf(stderr, "Error in mpfr_sub: result is not msb-normalized (3)\n");
! 429: exit (1);
! 430: }
! 431:
! 432: /* bug found by Fabrice Rouillier, 27 Mar 2001 */
! 433: mpfr_set_prec (x, 107);
! 434: mpfr_set_prec (t, 107);
! 435: mpfr_set_prec (u, 107);
! 436: mpfr_set_str_raw (x, "0.10111001001111010010001000000010111111011011011101000001001000101000000000000000000000000000000000000000000E315");
! 437: mpfr_set_str_raw (t, "0.10000000000000000000000000000000000101110100100101110110000001100101011111001000011101111100100100111011000E350");
! 438: mpfr_sub (u, x, t, GMP_RNDU);
! 439: if ((MPFR_MANT(u)[(MPFR_PREC(u)-1)/mp_bits_per_limb] &
! 440: ((mp_limb_t)1<<(mp_bits_per_limb-1)))==0) {
! 441: fprintf(stderr, "Error in mpfr_sub: result is not msb-normalized (4)\n");
! 442: exit (1);
! 443: }
! 444:
! 445: /* checks that NaN flag is correctly reset */
! 446: mpfr_set_d (t, 1.0, GMP_RNDN);
! 447: mpfr_set_d (u, 1.0, GMP_RNDN);
! 448: mpfr_set_nan (x);
! 449: mpfr_add (x, t, u, GMP_RNDN);
! 450: if (mpfr_cmp_ui (x, 2)) {
! 451: fprintf (stderr, "Error in mpfr_add: 1+1 gives %e\n", mpfr_get_d1 (x));
! 452: exit (1);
! 453: }
! 454:
! 455: mpfr_clear(x); mpfr_clear(t); mpfr_clear(u);
! 456: }
! 457:
! 458: /* check case when c does not overlap with a, but both b and c count
! 459: for rounding */
! 460: void
! 461: check_case_1b (void)
! 462: {
! 463: mpfr_t a, b, c;
! 464: unsigned int prec_a, prec_b, prec_c, dif;
! 465:
! 466: mpfr_init (a);
! 467: mpfr_init (b);
! 468: mpfr_init (c);
! 469:
! 470: for (prec_a = 2; prec_a <= 64; prec_a++)
! 471: {
! 472: mpfr_set_prec (a, prec_a);
! 473: for (prec_b = prec_a + 2; prec_b <= 64; prec_b++)
! 474: {
! 475: dif = prec_b - prec_a;
! 476: mpfr_set_prec (b, prec_b);
! 477: /* b = 1 - 2^(-prec_a) + 2^(-prec_b) */
! 478: mpfr_set_ui (b, 1, GMP_RNDN);
! 479: mpfr_div_2exp (b, b, dif, GMP_RNDN);
! 480: mpfr_sub_ui (b, b, 1, GMP_RNDN);
! 481: mpfr_div_2exp (b, b, prec_a, GMP_RNDN);
! 482: mpfr_add_ui (b, b, 1, GMP_RNDN);
! 483: for (prec_c = dif; prec_c <= 64; prec_c++)
! 484: {
! 485: /* c = 2^(-prec_a) - 2^(-prec_b) */
! 486: mpfr_set_prec (c, prec_c);
! 487: mpfr_set_si (c, -1, GMP_RNDN);
! 488: mpfr_div_2exp (c, c, dif, GMP_RNDN);
! 489: mpfr_add_ui (c, c, 1, GMP_RNDN);
! 490: mpfr_div_2exp (c, c, prec_a, GMP_RNDN);
! 491: mpfr_add (a, b, c, GMP_RNDN);
! 492: if (mpfr_cmp_ui (a, 1) != 0)
! 493: {
! 494: fprintf (stderr, "case (1b) failed for prec_a=%u, prec_b=%u, prec_c=%u\n", prec_a, prec_b, prec_c);
! 495: printf("b="); mpfr_print_binary(b); putchar('\n');
! 496: printf("c="); mpfr_print_binary(c); putchar('\n');
! 497: printf("a="); mpfr_print_binary(a); putchar('\n');
! 498: exit (1);
! 499: }
! 500: }
! 501: }
! 502: }
! 503:
! 504: mpfr_clear (a);
! 505: mpfr_clear (b);
! 506: mpfr_clear (c);
! 507: }
! 508:
! 509: /* check case when c overlaps with a */
! 510: void
! 511: check_case_2 (void)
! 512: {
! 513: mpfr_t a, b, c, d;
! 514:
! 515: mpfr_init2 (a, 300);
! 516: mpfr_init2 (b, 800);
! 517: mpfr_init2 (c, 500);
! 518: mpfr_init2 (d, 800);
! 519:
! 520: mpfr_set_str_raw(a, "1E110"); /* a = 2^110 */
! 521: mpfr_set_str_raw(b, "1E900"); /* b = 2^900 */
! 522: mpfr_set_str_raw(c, "1E500"); /* c = 2^500 */
! 523: mpfr_add(c, c, a, GMP_RNDZ); /* c = 2^500 + 2^110 */
! 524: mpfr_sub(d, b, c, GMP_RNDZ); /* d = 2^900 - 2^500 - 2^110 */
! 525: mpfr_add(b, b, c, GMP_RNDZ); /* b = 2^900 + 2^500 + 2^110 */
! 526: mpfr_add(a, b, d, GMP_RNDZ); /* a = 2^901 */
! 527: if (mpfr_cmp_ui_2exp (a, 1, 901))
! 528: {
! 529: fprintf (stderr, "b + d fails for b=2^900+2^500+2^110, d=2^900-2^500-2^110\n");
! 530: fprintf (stderr, "expected 1.0e901, got ");
! 531: mpfr_out_str (stderr, 2, 0, a, GMP_RNDN);
! 532: fprintf (stderr, "\n");
! 533: exit (1);
! 534: }
! 535:
! 536: mpfr_clear (a);
! 537: mpfr_clear (b);
! 538: mpfr_clear (c);
! 539: mpfr_clear (d);
! 540: }
! 541:
! 542: /* checks when source and destination are equal */
! 543: void
! 544: check_same (void)
! 545: {
! 546: mpfr_t x;
! 547:
! 548: mpfr_init(x); mpfr_set_d(x, 1.0, GMP_RNDZ);
! 549: mpfr_add(x, x, x, GMP_RNDZ);
! 550: if (mpfr_get_d1 (x) != 2.0) {
! 551: printf("Error when all 3 operands are equal\n"); exit(1);
! 552: }
! 553: mpfr_clear(x);
! 554: }
! 555:
! 556: #define check53(x, y, r, z) check(x, y, r, 53, 53, 53, z)
! 557: #define check53nan(x, y, r) checknan(x, y, r, 53, 53, 53);
! 558:
! 559: #define MAX_PREC 100
! 560:
! 561: void
! 562: check_inexact (void)
! 563: {
! 564: mpfr_t x, y, z, u;
! 565: mp_prec_t px, py, pu, pz;
! 566: int inexact, cmp;
! 567: mp_rnd_t rnd;
! 568:
! 569: mpfr_init (x);
! 570: mpfr_init (y);
! 571: mpfr_init (z);
! 572: mpfr_init (u);
! 573:
! 574: mpfr_set_prec (x, 2);
! 575: mpfr_set_str_raw (x, "0.1E-4");
! 576: mpfr_set_prec (u, 33);
! 577: mpfr_set_str_raw (u, "0.101110100101101100000000111100000E-1");
! 578: mpfr_set_prec (y, 31);
! 579: if ((inexact = mpfr_add (y, x, u, GMP_RNDN)))
! 580: {
! 581: fprintf (stderr, "Wrong inexact flag (2): expected 0, got %d\n", inexact);
! 582: exit (1);
! 583: }
! 584:
! 585: mpfr_set_prec (x, 2);
! 586: mpfr_set_str_raw (x, "0.1E-4");
! 587: mpfr_set_prec (u, 33);
! 588: mpfr_set_str_raw (u, "0.101110100101101100000000111100000E-1");
! 589: mpfr_set_prec (y, 28);
! 590: if ((inexact = mpfr_add (y, x, u, GMP_RNDN)))
! 591: {
! 592: fprintf (stderr, "Wrong inexact flag (1): expected 0, got %d\n", inexact);
! 593: exit (1);
! 594: }
! 595:
! 596: for (px=2; px<MAX_PREC; px++)
! 597: {
! 598: mpfr_set_prec (x, px);
! 599: mpfr_random (x);
! 600: for (pu=2; pu<MAX_PREC; pu++)
! 601: {
! 602: mpfr_set_prec (u, pu);
! 603: mpfr_random (u);
! 604: for (py=2; py<MAX_PREC; py++)
! 605: {
! 606: mpfr_set_prec (y, py);
! 607: pz = (mpfr_cmp_abs (x, u) >= 0) ? MPFR_EXP(x)-MPFR_EXP(u)
! 608: : MPFR_EXP(u)-MPFR_EXP(x);
! 609: /* x + u is exactly representable with precision
! 610: abs(EXP(x)-EXP(u)) + max(prec(x), prec(u)) + 1 */
! 611: pz = pz + MAX(MPFR_PREC(x), MPFR_PREC(u)) + 1;
! 612: mpfr_set_prec (z, pz);
! 613: rnd = LONG_RAND () % 4;
! 614: if (mpfr_add (z, x, u, rnd))
! 615: {
! 616: fprintf (stderr, "z <- x + u should be exact\n");
! 617: printf ("x="); mpfr_print_binary (x); putchar ('\n');
! 618: printf ("u="); mpfr_print_binary (u); putchar ('\n');
! 619: printf ("z="); mpfr_print_binary (z); putchar ('\n');
! 620: exit (1);
! 621: }
! 622: for (rnd=0; rnd<4; rnd++)
! 623: {
! 624: inexact = mpfr_add (y, x, u, rnd);
! 625: cmp = mpfr_cmp (y, z);
! 626: if (((inexact == 0) && (cmp != 0)) ||
! 627: ((inexact > 0) && (cmp <= 0)) ||
! 628: ((inexact < 0) && (cmp >= 0)))
! 629: {
! 630: fprintf (stderr, "Wrong inexact flag for rnd=%s\n",
! 631: mpfr_print_rnd_mode(rnd));
! 632: printf ("expected %d, got %d\n", cmp, inexact);
! 633: printf ("x="); mpfr_print_binary (x); putchar ('\n');
! 634: printf ("u="); mpfr_print_binary (u); putchar ('\n');
! 635: printf ("y= "); mpfr_print_binary (y); putchar ('\n');
! 636: printf ("x+u="); mpfr_print_binary (z); putchar ('\n');
! 637: exit (1);
! 638: }
! 639: }
! 640: }
! 641: }
! 642: }
! 643:
! 644: mpfr_clear (x);
! 645: mpfr_clear (y);
! 646: mpfr_clear (z);
! 647: mpfr_clear (u);
! 648: }
! 649:
! 650: int
! 651: main (int argc, char *argv[])
! 652: {
! 653: #ifdef MPFR_HAVE_FESETROUND
! 654: int prec, rnd_mode;
! 655: int rnd;
! 656: double y;
! 657: #endif
! 658: double x;
! 659: int i;
! 660:
! 661: mpfr_test_init ();
! 662: check_inexact ();
! 663: check_case_1b ();
! 664: check_case_2 ();
! 665: check64();
! 666: check(293607738.0, 1.9967571564050541e-5, GMP_RNDU, 64, 53, 53,
! 667: 2.9360773800002003e8);
! 668: check(880524.0, -2.0769715792901673e-5, GMP_RNDN, 64, 53, 53,
! 669: 8.8052399997923023e5);
! 670: check(1196426492.0, -1.4218093058435347e-3, GMP_RNDN, 64, 53, 53,
! 671: 1.1964264919985781e9);
! 672: check(982013018.0, -8.941829477291838e-7, GMP_RNDN, 64, 53, 53,
! 673: 9.8201301799999905e8);
! 674: check(1092583421.0, 1.0880649218158844e9, GMP_RNDN, 64, 53, 53,
! 675: 2.1806483428158846e9);
! 676: check(1.8476886419022969e-6, 961494401.0, GMP_RNDN, 53, 64, 53,
! 677: 9.6149440100000179e8);
! 678: check(-2.3222118418069868e5, 1229318102.0, GMP_RNDN, 53, 64, 53,
! 679: 1.2290858808158193e9);
! 680: check(-3.0399171300395734e-6, 874924868.0, GMP_RNDN, 53, 64, 53,
! 681: 8.749248679999969e8);
! 682: check(9.064246624706179e1, 663787413.0, GMP_RNDN, 53, 64, 53,
! 683: 6.6378750364246619e8);
! 684: check(-1.0954322421551264e2, 281806592.0, GMP_RNDD, 53, 64, 53,
! 685: 2.8180648245677572e8);
! 686: check(5.9836930386056659e-8, 1016217213.0, GMP_RNDN, 53, 64, 53,
! 687: 1.0162172130000001e9);
! 688: check(-1.2772161928500301e-7, 1237734238.0, GMP_RNDN, 53, 64, 53,
! 689: 1.2377342379999998e9);
! 690: check(-4.567291988483277e8, 1262857194.0, GMP_RNDN, 53, 64, 53,
! 691: 8.0612799515167236e8);
! 692: check(4.7719471752925262e7, 196089880.0, GMP_RNDN, 53, 53, 53,
! 693: 2.4380935175292528e8);
! 694: check(4.7719471752925262e7, 196089880.0, GMP_RNDN, 53, 64, 53,
! 695: 2.4380935175292528e8);
! 696: check(-1.716113812768534e-140, 1271212614.0, GMP_RNDZ, 53, 64, 53,
! 697: 1.2712126139999998e9);
! 698: check(-1.2927455200185474e-50, 1675676122.0, GMP_RNDD, 53, 64, 53,
! 699: 1.6756761219999998e9);
! 700: check53(1.22191250737771397120e+20, 948002822.0, GMP_RNDN,
! 701: 122191250738719408128.0);
! 702: check53(9966027674114492.0, 1780341389094537.0, GMP_RNDN,
! 703: 11746369063209028.0);
! 704: check53(2.99280481918991653800e+272, 5.34637717585790933424e+271, GMP_RNDN,
! 705: 3.5274425367757071711e272);
! 706: check_same();
! 707: check53(6.14384195492641560499e-02, -6.14384195401037683237e-02, GMP_RNDU,
! 708: 9.1603877261370314499e-12);
! 709: check53(1.16809465359248765399e+196, 7.92883212101990665259e+196, GMP_RNDU,
! 710: 9.0969267746123943065e196);
! 711: check53(3.14553393112021279444e-67, 3.14553401015952024126e-67, GMP_RNDU,
! 712: 6.2910679412797336946e-67);
! 713:
! 714: SEED_RAND (time(NULL));
! 715: check53(5.43885304644369509058e+185,-1.87427265794105342763e-57,GMP_RNDN,
! 716: 5.4388530464436950905e185);
! 717: check53(5.43885304644369509058e+185,-1.87427265794105342763e-57, GMP_RNDZ,
! 718: 5.4388530464436944867e185);
! 719: check53(5.43885304644369509058e+185,-1.87427265794105342763e-57, GMP_RNDU,
! 720: 5.4388530464436950905e185);
! 721: check53(5.43885304644369509058e+185,-1.87427265794105342763e-57, GMP_RNDD,
! 722: 5.4388530464436944867e185);
! 723: check2a(6.85523243386777784171e+107,187,-2.78148588123699111146e+48,87,178,
! 724: GMP_RNDD, "4.ab980a5cb9407ffffffffffffffffffffffffffffffe@89");
! 725: check2a(-1.21510626304662318398e+145,70,1.21367733647758957118e+145,65,61,
! 726: GMP_RNDD, "-1.2bfad031d94@118");
! 727: check2a(2.73028857032080744543e+155,83,-1.16446121423113355603e+163,59,125,
! 728: GMP_RNDZ, "-3.3c42dee09703d0639a6@135");
! 729: check2a(-4.38589520019641698848e+78,155,-1.09923643769309483415e+72,15,159,
! 730: GMP_RNDD, "-2.5e09955c663d@65");
! 731: check2a(-1.49963910666191123860e+265,76,-2.30915090591874527520e-191,8,75,
! 732: GMP_RNDZ, "-1.dc3ec027da54e@220");
! 733: check2a(3.25471707846623300604e-160,81,-7.93846654265839958715e-274,58,54,
! 734: GMP_RNDN, "4.936a52bc17254@-133");
! 735: check2a(5.17945380930936917508e+112,119,1.11369077158813567738e+108,15,150,
! 736: GMP_RNDZ, "5.62661692498ec@93");
! 737: check2a(-2.66910493504493276454e-52,117,1.61188644159592323415e-52,61,68,
! 738: GMP_RNDZ, "-a.204acdd25d788@-44");
! 739: check2a(-1.87427265794105342764e-57,175,1.76570844587489516446e+190,2,115,
! 740: GMP_RNDZ, "b.fffffffffffffffffffffffffffe@157");
! 741: check2a(-1.15706375390780417299e-135,94,-1.07455137477117851576e-129,66,111,
! 742: GMP_RNDU, "-b.eae2643497ff6286b@-108");
! 743: check2a(-1.15706375390780417299e-135,94,-1.07455137477117851576e-129,66,111,
! 744: GMP_RNDD, "-b.eae2643497ff6286b@-108");
! 745: check2a(-3.31624349995221499866e-22,107,-8.20150212714204839621e+156,79,99,
! 746: GMP_RNDD, "-2.63b22b55697e8000000000008@130");
! 747: x = -5943982715394951.0; for (i=0; i<446; i++) x *= 2.0;
! 748: check2a(x, 63, 1.77607317509426332389e+73, 64, 64, GMP_RNDN,
! 749: "-5.4781549356e1c@124");
! 750: check2a(4.49465557237618783128e+53,108,-2.45103927353799477871e+48,60,105,
! 751: GMP_RNDN, "4.b14f230f909dc803e@44");
! 752: check2a(2.26531902208967707071e+168,99,-2.67795218510613988524e+168,67,94,
! 753: GMP_RNDU, "-1.bfd7ff2647098@139");
! 754: check2a(-8.88471912490080158206e+253,79,-7.84488427404526918825e+124,95,53,
! 755: GMP_RNDD, "-c.1e533b8d835@210");
! 756: check2a(-2.18548638152863831959e-125,61,-1.22788940592412363564e+226,71,54,
! 757: GMP_RNDN, "-8.4b0f99ffa3b58@187");
! 758: check2a(-7.94156823309993162569e+77,74,-5.26820160805275124882e+80,99,101,
! 759: GMP_RNDD, "-1.1cc90f11d6af26f4@67");
! 760: check2a(-3.85170653452493859064e+189,62,2.18827389706660314090e+158,94,106,
! 761: GMP_RNDD, "-3.753ac0935b701ffffffffffffd@157");
! 762: check2a(1.07966151149311101950e+46,88,1.13198076934147457400e+46,67,53,
! 763: GMP_RNDN, "3.dfbc152dd4368@38");
! 764: check2a(3.36768223223409657622e+209,55,-9.61624007357265441884e+219,113,53,
! 765: GMP_RNDN, "-6.cf7217a451388@182");
! 766: check2a(-6.47376909368979326475e+159,111,5.11127211034490340501e+159,99,62,
! 767: GMP_RNDD, "-1.8cf3aadf537c@132");
! 768: check2a(-4.95229483271607845549e+220,110,-6.06992115033276044773e+213,109,55,
! 769: GMP_RNDN, "-2.3129f1f63b31b@183");
! 770: check2a(-6.47376909368979326475e+159,74,5.11127211034490340501e+159,111,75,
! 771: GMP_RNDU, "-1.8cf3aadf537c@132");
! 772: check2a(2.26531902208967707070e+168,99,-2.67795218510613988525e+168,67,94,
! 773: GMP_RNDU, "-1.bfd7ff2647098@139");
! 774: check2a(-2.28886326552077479586e-188,67,3.41419438647157839320e-177,60,110,
! 775: GMP_RNDU, "3.75740b4fe8f17f90258907@-147");
! 776: check2a(-2.66910493504493276454e-52,117,1.61188644159592323415e-52,61,68,
! 777: GMP_RNDZ, "-a.204acdd25d788@-44");
! 778: check2a(2.90983392714730768886e+50,101,2.31299792168440591870e+50,74,105,
! 779: GMP_RNDZ, "1.655c53ff5719c8@42");
! 780: check2a(2.72046257722708717791e+243,97,-1.62158447436486437113e+243,83,96,
! 781: GMP_RNDN, "a.4cc63e002d2e8@201");
! 782: /* Checking double precision (53 bits) */
! 783: check53(-8.22183238641455905806e-19, 7.42227178769761587878e-19, GMP_RNDD,
! 784: -7.9956059871694317927e-20);
! 785: check53(5.82106394662028628236e+234, -5.21514064202368477230e+89, GMP_RNDD,
! 786: 5.8210639466202855763e234);
! 787: check53(5.72931679569871602371e+122, -5.72886070363264321230e+122, GMP_RNDN,
! 788: 4.5609206607281141508e118);
! 789: check53(-5.09937369394650450820e+238, 2.70203299854862982387e+250, GMP_RNDD,
! 790: 2.7020329985435301323e250);
! 791: check53(-2.96695924472363684394e+27, 1.22842938251111500000e+16, GMP_RNDD,
! 792: -2.96695924471135255027e27);
! 793: check53(1.74693641655743793422e-227, -7.71776956366861843469e-229, GMP_RNDN,
! 794: 1.669758720920751867e-227);
! 795: x = -7883040437021647.0; for (i=0; i<468; i++) x = x / 2.0;
! 796: check53(-1.03432206392780011159e-125, 1.30127034799251347548e-133, GMP_RNDN,
! 797: x);
! 798: check53(1.05824655795525779205e+71, -1.06022698059744327881e+71, GMP_RNDZ,
! 799: -1.9804226421854867632e68);
! 800: check53(-5.84204911040921732219e+240, 7.26658169050749590763e+240, GMP_RNDD,
! 801: 1.4245325800982785854e240);
! 802: check53(1.00944884131046636376e+221, 2.33809162651471520268e+215, GMP_RNDN,
! 803: 1.0094511794020929787e221);
! 804: x = 7045852550057985.0; for (i=0; i<986; i++) x = x / 2.0;
! 805: check53(4.29232078932667367325e-278, x, GMP_RNDU,
! 806: 4.2933981418314132787e-278);
! 807: check53(5.27584773801377058681e-80, 8.91207657803547196421e-91, GMP_RNDN,
! 808: 5.2758477381028917269e-80);
! 809: check53(2.99280481918991653800e+272, 5.34637717585790933424e+271, GMP_RNDN,
! 810: 3.5274425367757071711e272);
! 811: check53(4.67302514390488041733e-184, 2.18321376145645689945e-190, GMP_RNDN,
! 812: 4.6730273271186420541e-184);
! 813: check53(5.57294120336300389254e+71, 2.60596167942024924040e+65, GMP_RNDZ,
! 814: 5.5729438093246831053e71);
! 815: check53(6.6052588496951015469e24, 4938448004894539.0, GMP_RNDU,
! 816: 6.6052588546335505068e24);
! 817: check53(1.23056185051606761523e-190, 1.64589756643433857138e-181, GMP_RNDU,
! 818: 1.6458975676649006598e-181);
! 819: check53(2.93231171510175981584e-280, 3.26266919161341483877e-273, GMP_RNDU,
! 820: 3.2626694848445867288e-273);
! 821: check53(5.76707395945001907217e-58, 4.74752971449827687074e-51, GMP_RNDD,
! 822: 4.747530291205672325e-51);
! 823: check53(277363943109.0, 11.0, GMP_RNDN, 277363943120.0);
! 824: #if 0 /* disabled since it seems silly to use denorms *
! 825: /* test denormalized numbers too */
! 826: check53(8.06294740693074521573e-310, 6.95250701071929654575e-310, GMP_RNDU,
! 827: 1.5015454417650041761e-309);
! 828: #endif
! 829: #ifdef HAVE_INFS
! 830: /* the following check double overflow */
! 831: check53(6.27557402141211962228e+307, 1.32141396570101687757e+308,
! 832: GMP_RNDZ, DBL_POS_INF);
! 833: check53(DBL_POS_INF, 6.95250701071929654575e-310, GMP_RNDU, DBL_POS_INF);
! 834: check53(DBL_NEG_INF, 6.95250701071929654575e-310, GMP_RNDU, DBL_NEG_INF);
! 835: check53(6.95250701071929654575e-310, DBL_POS_INF, GMP_RNDU, DBL_POS_INF);
! 836: check53(6.95250701071929654575e-310, DBL_NEG_INF, GMP_RNDU, DBL_NEG_INF);
! 837: check53nan (DBL_POS_INF, DBL_NEG_INF, GMP_RNDN);
! 838: #endif
! 839: check53(1.44791789689198883921e-140, -1.90982880222349071284e-121,
! 840: GMP_RNDN, -1.90982880222349071e-121);
! 841:
! 842:
! 843: /* tests for particular cases (Vincent Lefevre, 22 Aug 2001) */
! 844: check53(9007199254740992.0, 1.0, GMP_RNDN, 9007199254740992.0);
! 845: check53(9007199254740994.0, 1.0, GMP_RNDN, 9007199254740996.0);
! 846: check53(9007199254740992.0, -1.0, GMP_RNDN, 9007199254740991.0);
! 847: check53(9007199254740994.0, -1.0, GMP_RNDN, 9007199254740992.0);
! 848: check53(9007199254740996.0, -1.0, GMP_RNDN, 9007199254740996.0);
! 849:
! 850: #ifdef MPFR_HAVE_FESETROUND
! 851: prec = (argc<2) ? 53 : atoi(argv[1]);
! 852: rnd_mode = (argc<3) ? -1 : atoi(argv[2]);
! 853: /* Comparing to double precision using machine arithmetic */
! 854: for (i=0;i<N;i++) {
! 855: x = drand();
! 856: y = drand();
! 857: if (ABS(x)>2.2e-307 && ABS(y)>2.2e-307 && x+y<1.7e+308 && x+y>-1.7e308) {
! 858: /* avoid denormalized numbers and overflows */
! 859: rnd = (rnd_mode==-1) ? LONG_RAND()%4 : rnd_mode;
! 860: check(x, y, rnd, prec, prec, prec, 0.0);
! 861: }
! 862: }
! 863: /* tests with random precisions */
! 864: for (i=0;i<N;i++) {
! 865: int px, py, pz;
! 866: px = 53 + (LONG_RAND() % 64);
! 867: py = 53 + (LONG_RAND() % 64);
! 868: pz = 53 + (LONG_RAND() % 64);
! 869: rnd_mode = LONG_RAND() % 4;
! 870: do { x = drand(); } while (isnan(x));
! 871: do { y = drand(); } while (isnan(y));
! 872: check2 (x, px, y, py, pz, rnd_mode);
! 873: }
! 874: /* Checking mpfr_add(x, x, y) with prec=53 */
! 875: for (i=0;i<N;i++) {
! 876: x = drand();
! 877: y = drand();
! 878: if (ABS(x)>2.2e-307 && ABS(y)>2.2e-307 && x+y<1.7e+308 && x+y>-1.7e308) {
! 879: /* avoid denormalized numbers and overflows */
! 880: rnd = (rnd_mode==-1) ? LONG_RAND()%4 : rnd_mode;
! 881: check3(x, y, rnd);
! 882: }
! 883: }
! 884: /* Checking mpfr_add(x, y, x) with prec=53 */
! 885: for (i=0;i<N;i++) {
! 886: x = drand();
! 887: y = drand();
! 888: if (ABS(x)>2.2e-307 && ABS(y)>2.2e-307 && x+y<1.7e+308 && x+y>-1.7e308) {
! 889: /* avoid denormalized numbers and overflows */
! 890: rnd = (rnd_mode==-1) ? LONG_RAND()%4 : rnd_mode;
! 891: check4(x, y, rnd);
! 892: }
! 893: }
! 894: /* Checking mpfr_add(x, x, x) with prec=53 */
! 895: for (i=0;i<N;i++) {
! 896: do { x = drand(); } while ((ABS(x)<2.2e-307) || (ABS(x)>0.8e308));
! 897: /* avoid denormalized numbers and overflows */
! 898: rnd = (rnd_mode==-1) ? LONG_RAND()%4 : rnd_mode;
! 899: check5(x, rnd);
! 900: }
! 901: #endif
! 902:
! 903: return 0;
! 904: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to tadd.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr / tests