Annotation of OpenXM_contrib/gmp/mpfr/tests/tsub.c, Revision 1.1
1.1 ! ohara 1: /* Test file for mpfr_sub.
! 2:
! 3: Copyright 2001 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: #include <math.h>
! 23: #include <stdio.h>
! 24: #include <stdlib.h>
! 25: #include "gmp.h"
! 26: #include "mpfr.h"
! 27: #include "mpfr-impl.h"
! 28: #include "mpfr-test.h"
! 29:
! 30: void check_diverse _PROTO((void));
! 31: void bug_ddefour _PROTO((void));
! 32: void check_two_sum _PROTO((mp_prec_t));
! 33: void check_inexact _PROTO((void));
! 34:
! 35: void
! 36: check_diverse (void)
! 37: {
! 38: mpfr_t x, y, z;
! 39: double res, got;
! 40: int inexact;
! 41:
! 42: mpfr_init (x);
! 43: mpfr_init (y);
! 44: mpfr_init (z);
! 45:
! 46: mpfr_set_prec (x, 32);
! 47: mpfr_set_prec (y, 63);
! 48: mpfr_set_prec (z, 63);
! 49: mpfr_set_str_raw (x, "0.101101111011011100100100100111E31");
! 50: mpfr_set_str_raw (y, "0.111110010010100100110101101010001001100101110001000101110111111E-1");
! 51: mpfr_sub (z, x, y, GMP_RNDN);
! 52: mpfr_set_str_raw (y, "0.1011011110110111001001001001101100000110110101101100101001011E31");
! 53: if (mpfr_cmp (z, y))
! 54: {
! 55: fprintf (stderr, "Error in mpfr_sub (5)\n");
! 56: printf ("expected "); mpfr_print_binary (y); putchar ('\n');
! 57: printf ("got "); mpfr_print_binary (z); putchar ('\n');
! 58: exit (1);
! 59: }
! 60:
! 61: mpfr_set_prec (y, 63);
! 62: mpfr_set_prec (z, 63);
! 63: mpfr_set_str_raw (y, "0.1011011110110111001001001001101100000110110101101100101001011E31");
! 64: mpfr_sub_ui (z, y, 1541116494, GMP_RNDN);
! 65: mpfr_set_str_raw (y, "-0.11111001001010010011010110101E-1");
! 66: if (mpfr_cmp (z, y))
! 67: {
! 68: fprintf (stderr, "Error in mpfr_sub (7)\n");
! 69: printf ("expected "); mpfr_print_binary (y); putchar ('\n');
! 70: printf ("got "); mpfr_print_binary (z); putchar ('\n');
! 71: exit (1);
! 72: }
! 73:
! 74: mpfr_set_prec (y, 63);
! 75: mpfr_set_prec (z, 63);
! 76: mpfr_set_str_raw (y, "0.1011011110110111001001001001101100000110110101101100101001011E31");
! 77: mpfr_sub_ui (z, y, 1541116494, GMP_RNDN);
! 78: mpfr_set_str_raw (y, "-0.11111001001010010011010110101E-1");
! 79: if (mpfr_cmp (z, y))
! 80: {
! 81: fprintf (stderr, "Error in mpfr_sub (6)\n");
! 82: printf ("expected "); mpfr_print_binary (y); putchar ('\n');
! 83: printf ("got "); mpfr_print_binary (z); putchar ('\n');
! 84: exit (1);
! 85: }
! 86:
! 87: mpfr_set_prec (x, 32);
! 88: mpfr_set_prec (y, 32);
! 89: mpfr_set_str_raw (x, "0.10110111101001110100100101111000E0");
! 90: mpfr_set_str_raw (y, "0.10001100100101000100110111000100E0");
! 91: if ((inexact = mpfr_sub (x, x, y, GMP_RNDN)))
! 92: {
! 93: fprintf (stderr, "Wrong inexact flag (2): got %d instead of 0\n",
! 94: inexact);
! 95: exit (1);
! 96: }
! 97:
! 98: mpfr_set_prec (x, 32);
! 99: mpfr_set_prec (y, 32);
! 100: mpfr_set_str_raw (x, "0.11111000110111011000100111011010E0");
! 101: mpfr_set_str_raw (y, "0.10011111101111000100001000000000E-3");
! 102: if ((inexact = mpfr_sub (x, x, y, GMP_RNDN)))
! 103: {
! 104: fprintf (stderr, "Wrong inexact flag (1): got %d instead of 0\n",
! 105: inexact);
! 106: exit (1);
! 107: }
! 108:
! 109: mpfr_set_prec (x, 33);
! 110: mpfr_set_prec (y, 33);
! 111: mpfr_set_ui (x, 1, GMP_RNDN);
! 112: mpfr_set_str_raw (y, "-0.1E-32");
! 113: mpfr_add (x, x, y, GMP_RNDN);
! 114: mpfr_set_str_raw (y, "0.111111111111111111111111111111111E0");
! 115: if (mpfr_cmp (x, y))
! 116: {
! 117: fprintf (stderr, "Error in mpfr_sub (1 - 1E-33) with prec=33\n");
! 118: printf ("Expected "); mpfr_print_binary (y); putchar ('\n');
! 119: printf ("got "); mpfr_print_binary (x); putchar ('\n');
! 120: exit (1);
! 121: }
! 122:
! 123: mpfr_set_prec (x, 32);
! 124: mpfr_set_prec (y, 33);
! 125: mpfr_set_ui (x, 1, GMP_RNDN);
! 126: mpfr_set_str_raw (y, "-0.1E-32");
! 127: mpfr_add (x, x, y, GMP_RNDN);
! 128: if (mpfr_cmp_ui (x, 1))
! 129: {
! 130: fprintf (stderr, "Error in mpfr_sub (1 - 1E-33) with prec=32\n");
! 131: printf ("Expected 1.0, got "); mpfr_print_binary (x); putchar ('\n');
! 132: exit (1);
! 133: }
! 134:
! 135: mpfr_set_prec (x, 65);
! 136: mpfr_set_prec (y, 65);
! 137: mpfr_set_prec (z, 64);
! 138: mpfr_set_str_raw (x, "1.1110111011110001110111011111111111101000011001011100101100101101");
! 139: mpfr_set_str_raw (y, "0.1110111011110001110111011111111111101000011001011100101100101100");
! 140: mpfr_sub (z, x, y, GMP_RNDZ);
! 141: if (mpfr_cmp_ui (z, 1))
! 142: {
! 143: fprintf (stderr, "Error in mpfr_sub (1)\n");
! 144: exit (1);
! 145: }
! 146: mpfr_sub (z, x, y, GMP_RNDU);
! 147: mpfr_set_str_raw (x, "1.000000000000000000000000000000000000000000000000000000000000001");
! 148: if (mpfr_cmp (z, x))
! 149: {
! 150: fprintf (stderr, "Error in mpfr_sub (2)\n");
! 151: exit (1);
! 152: }
! 153: mpfr_set_str_raw (x, "1.1110111011110001110111011111111111101000011001011100101100101101");
! 154: mpfr_sub (z, x, y, GMP_RNDN);
! 155: if (mpfr_cmp_ui (z, 1))
! 156: {
! 157: fprintf (stderr, "Error in mpfr_sub (3)\n");
! 158: exit (1);
! 159: }
! 160: mpfr_set_prec (x, 66);
! 161: mpfr_set_str_raw (x, "1.11101110111100011101110111111111111010000110010111001011001010111");
! 162: mpfr_sub (z, x, y, GMP_RNDN);
! 163: if (mpfr_cmp_ui (z, 1))
! 164: {
! 165: fprintf (stderr, "Error in mpfr_sub (4)\n");
! 166: exit (1);
! 167: }
! 168:
! 169: /* check in-place operations */
! 170: mpfr_set_d (x, 1.0, GMP_RNDN);
! 171: mpfr_sub (x, x, x, GMP_RNDN);
! 172: if (mpfr_get_d1 (x) != 0.0)
! 173: {
! 174: fprintf (stderr, "Error for mpfr_add (x, x, x, GMP_RNDN) with x=1.0\n");
! 175: exit (1);
! 176: }
! 177:
! 178: mpfr_set_prec (x, 53);
! 179: mpfr_set_prec (y, 53);
! 180: mpfr_set_prec (z, 53);
! 181: mpfr_set_d (x, 1.229318102e+09, GMP_RNDN);
! 182: mpfr_set_d (y, 2.32221184180698677665e+05, GMP_RNDN);
! 183: mpfr_sub (z, x, y, GMP_RNDN);
! 184: res = 1229085880.815819263458251953125;
! 185: got = mpfr_get_d1 (z);
! 186: if (got != res)
! 187: {
! 188: fprintf (stderr, "Error in mpfr_sub (1.22e9 - 2.32e5)\n");
! 189: printf ("expected %1.20e, got %1.20e\n", res, got);
! 190: exit (1);
! 191: }
! 192:
! 193: mpfr_set_prec (x, 112);
! 194: mpfr_set_prec (y, 98);
! 195: mpfr_set_prec (z, 54);
! 196: mpfr_set_str_raw (x, "0.11111100100000000011000011100000101101010001000111E-401");
! 197: mpfr_set_str_raw (y, "0.10110000100100000101101100011111111011101000111000101E-464");
! 198: mpfr_sub (z, x, y, GMP_RNDN);
! 199: if (mpfr_cmp (z, x)) {
! 200: fprintf (stderr, "mpfr_sub(z, x, y) failed for prec(x)=112, prec(y)=98\n");
! 201: printf ("expected "); mpfr_print_binary (x); putchar('\n');
! 202: printf ("got "); mpfr_print_binary (z); putchar('\n');
! 203: exit (1);
! 204: }
! 205:
! 206: mpfr_set_prec (x, 33);
! 207: mpfr_set_ui (x, 1, GMP_RNDN);
! 208: mpfr_div_2exp (x, x, 32, GMP_RNDN);
! 209: mpfr_sub_ui (x, x, 1, GMP_RNDN);
! 210:
! 211: mpfr_set_prec (x, 5);
! 212: mpfr_set_prec (y, 5);
! 213: mpfr_set_str_raw (x, "1e-12");
! 214: mpfr_set_ui (y, 1, GMP_RNDN);
! 215: mpfr_sub (x, y, x, GMP_RNDD);
! 216: mpfr_set_str_raw (y, "0.11111");
! 217: if (mpfr_cmp (x, y))
! 218: {
! 219: fprintf (stderr, "Error in mpfr_sub (x, y, x, GMP_RNDD) for x=2^(-12), y=1\n");
! 220: exit (1);
! 221: }
! 222:
! 223: mpfr_set_prec (x, 24);
! 224: mpfr_set_prec (y, 24);
! 225: mpfr_set_str_raw (x, "-0.100010000000000000000000E19");
! 226: mpfr_set_str_raw (y, "0.100000000000000000000100E15");
! 227: mpfr_add (x, x, y, GMP_RNDD);
! 228: mpfr_set_str_raw (y, "-0.1E19");
! 229: if (mpfr_cmp (x, y))
! 230: {
! 231: fprintf (stderr, "Error in mpfr_add (2)\n");
! 232: exit (1);
! 233: }
! 234:
! 235: mpfr_set_prec (x, 2);
! 236: mpfr_set_prec (y, 10);
! 237: mpfr_set_prec (z, 10);
! 238: mpfr_set_ui (y, 0, GMP_RNDN);
! 239: mpfr_set_str_raw (z, "0.100000000000000000000100E15");
! 240: if (mpfr_sub (x, y, z, GMP_RNDN) <= 0)
! 241: {
! 242: fprintf (stderr, "Wrong inexact flag in x=mpfr_sub(0,z) for prec(z)>prec(x)\n");
! 243: exit (1);
! 244: }
! 245: if (mpfr_sub (x, z, y, GMP_RNDN) >= 0)
! 246: {
! 247: fprintf (stderr, "Wrong inexact flag in x=mpfr_sub(z,0) for prec(z)>prec(x)\n");
! 248: exit (1);
! 249: }
! 250:
! 251: mpfr_clear (x);
! 252: mpfr_clear (y);
! 253: mpfr_clear (z);
! 254: }
! 255:
! 256: void
! 257: bug_ddefour(void)
! 258: {
! 259: mpfr_t ex, ex1, ex2, ex3, tot, tot1;
! 260:
! 261: mpfr_init2(ex, 53);
! 262: mpfr_init2(ex1, 53);
! 263: mpfr_init2(ex2, 53);
! 264: mpfr_init2(ex3, 53);
! 265: mpfr_init2(tot, 150);
! 266: mpfr_init2(tot1, 150);
! 267:
! 268: mpfr_set_ui( ex, 1, GMP_RNDN);
! 269: mpfr_mul_2exp( ex, ex, 906, GMP_RNDN);
! 270: mpfr_log( tot, ex, GMP_RNDN);
! 271: mpfr_set( ex1, tot, GMP_RNDN); /* ex1 = high(tot) */
! 272: mpfr_sub( ex2, tot, ex1, GMP_RNDN); /* ex2 = high(tot - ex1) */
! 273: mpfr_sub( tot1, tot, ex1, GMP_RNDN); /* tot1 = tot - ex1 */
! 274: mpfr_set( ex3, tot1, GMP_RNDN); /* ex3 = high(tot - ex1) */
! 275:
! 276: if (mpfr_cmp(ex2, ex3))
! 277: {
! 278: fprintf (stderr, "Error in ddefour test.\n");
! 279: printf ("ex2="); mpfr_print_binary (ex2); putchar ('\n');
! 280: printf ("ex3="); mpfr_print_binary (ex3); putchar ('\n');
! 281: exit (1);
! 282: }
! 283:
! 284: mpfr_clear (ex);
! 285: mpfr_clear (ex1);
! 286: mpfr_clear (ex2);
! 287: mpfr_clear (ex3);
! 288: mpfr_clear (tot);
! 289: mpfr_clear (tot1);
! 290: }
! 291:
! 292: /* if u = o(x-y), v = o(u-x), w = o(v+y), then x-y = u-w */
! 293: void
! 294: check_two_sum (mp_prec_t p)
! 295: {
! 296: mpfr_t x, y, u, v, w;
! 297: mp_rnd_t rnd;
! 298: int inexact;
! 299:
! 300: mpfr_init2 (x, p);
! 301: mpfr_init2 (y, p);
! 302: mpfr_init2 (u, p);
! 303: mpfr_init2 (v, p);
! 304: mpfr_init2 (w, p);
! 305: mpfr_random (x);
! 306: mpfr_random (y);
! 307: if (mpfr_cmp_abs (x, y) < 0)
! 308: mpfr_swap (x, y);
! 309: rnd = LONG_RAND() % 4;
! 310: rnd = GMP_RNDN;
! 311: inexact = mpfr_sub (u, x, y, GMP_RNDN);
! 312: mpfr_sub (v, u, x, GMP_RNDN);
! 313: mpfr_add (w, v, y, GMP_RNDN);
! 314: /* as u = (x-y) - w, we should have inexact and w of opposite signs */
! 315: if (((inexact == 0) && mpfr_cmp_ui (w, 0)) ||
! 316: ((inexact > 0) && (mpfr_cmp_ui (w, 0) <= 0)) ||
! 317: ((inexact < 0) && (mpfr_cmp_ui (w, 0) >= 0)))
! 318: {
! 319: fprintf (stderr, "Wrong inexact flag for prec=%u, rnd=%s\n", (unsigned)p,
! 320: mpfr_print_rnd_mode (rnd));
! 321: printf ("x="); mpfr_print_binary(x); putchar('\n');
! 322: printf ("y="); mpfr_print_binary(y); putchar('\n');
! 323: printf ("u="); mpfr_print_binary(u); putchar('\n');
! 324: printf ("v="); mpfr_print_binary(v); putchar('\n');
! 325: printf ("w="); mpfr_print_binary(w); putchar('\n');
! 326: printf ("inexact = %d\n", inexact);
! 327: exit (1);
! 328: }
! 329: mpfr_clear (x);
! 330: mpfr_clear (y);
! 331: mpfr_clear (u);
! 332: mpfr_clear (v);
! 333: mpfr_clear (w);
! 334: }
! 335:
! 336: #define MAX_PREC 100
! 337:
! 338: void
! 339: check_inexact (void)
! 340: {
! 341: mpfr_t x, y, z, u;
! 342: mp_prec_t px, py, pu, pz;
! 343: int inexact, cmp;
! 344: mp_rnd_t rnd;
! 345:
! 346: mpfr_init (x);
! 347: mpfr_init (y);
! 348: mpfr_init (z);
! 349: mpfr_init (u);
! 350:
! 351: for (px=2; px<MAX_PREC; px++)
! 352: {
! 353: mpfr_set_prec (x, px);
! 354: mpfr_random (x);
! 355: for (pu=2; pu<MAX_PREC; pu++)
! 356: {
! 357: mpfr_set_prec (u, pu);
! 358: mpfr_random (u);
! 359: for (py=2; py<MAX_PREC; py++)
! 360: {
! 361: mpfr_set_prec (y, py);
! 362: pz = (mpfr_cmp_abs (x, u) >= 0) ? MPFR_EXP(x)-MPFR_EXP(u)
! 363: : MPFR_EXP(u)-MPFR_EXP(x);
! 364: pz = pz + MAX(MPFR_PREC(x), MPFR_PREC(u));
! 365: mpfr_set_prec (z, pz);
! 366: rnd = LONG_RAND () % 4;
! 367: if (mpfr_sub (z, x, u, rnd))
! 368: {
! 369: fprintf (stderr, "z <- x - u should be exact\n");
! 370: exit (1);
! 371: }
! 372: for (rnd=0; rnd<4; rnd++)
! 373: {
! 374: inexact = mpfr_sub (y, x, u, rnd);
! 375: cmp = mpfr_cmp (y, z);
! 376: if (((inexact == 0) && (cmp != 0)) ||
! 377: ((inexact > 0) && (cmp <= 0)) ||
! 378: ((inexact < 0) && (cmp >= 0)))
! 379: {
! 380: fprintf (stderr, "Wrong inexact flag for rnd=%s\n",
! 381: mpfr_print_rnd_mode(rnd));
! 382: printf ("expected %d, got %d\n", cmp, inexact);
! 383: printf ("x="); mpfr_print_binary (x); putchar ('\n');
! 384: printf ("u="); mpfr_print_binary (u); putchar ('\n');
! 385: printf ("y= "); mpfr_print_binary (y); putchar ('\n');
! 386: printf ("x-u="); mpfr_print_binary (z); putchar ('\n');
! 387: exit (1);
! 388: }
! 389: }
! 390: }
! 391: }
! 392: }
! 393:
! 394: mpfr_clear (x);
! 395: mpfr_clear (y);
! 396: mpfr_clear (z);
! 397: mpfr_clear (u);
! 398: }
! 399:
! 400: int
! 401: main(void)
! 402: {
! 403: mp_prec_t p;
! 404: unsigned i;
! 405:
! 406: check_diverse ();
! 407: check_inexact ();
! 408: bug_ddefour ();
! 409:
! 410: for (p=2; p<200; p++)
! 411: for (i=0; i<200; i++)
! 412: check_two_sum (p);
! 413:
! 414: return 0;
! 415: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to tsub.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