Annotation of OpenXM_contrib/gmp/mpfr/tests/tmul.c, Revision 1.1
1.1 ! ohara 1: /* Test file for mpfr_mul.
! 2:
! 3: Copyright 1999, 2001, 2002 Free Software Foundation, Inc.
! 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 <stdio.h>
! 23: #include <stdlib.h>
! 24: #include <time.h>
! 25: #include "gmp.h"
! 26: #include "gmp-impl.h"
! 27: #include "mpfr.h"
! 28: #include "mpfr-impl.h"
! 29: #include "mpfr-test.h"
! 30:
! 31: void check53 _PROTO((double, double, mp_rnd_t, double));
! 32: void check24 _PROTO((float, float, mp_rnd_t, float));
! 33: void check_float _PROTO((void));
! 34: void check_sign _PROTO((void));
! 35: void check_exact _PROTO((void));
! 36: void check_max _PROTO((void));
! 37: void check_min _PROTO((void));
! 38:
! 39:
! 40: /* Workaround for sparc gcc 2.95.x bug, see notes in tadd.c. */
! 41: #define check(x,y,rnd_mode,px,py,pz,res) _check(x,y,res,rnd_mode,px,py,pz)
! 42:
! 43: /* checks that x*y gives the same results in double
! 44: and with mpfr with 53 bits of precision */
! 45: void
! 46: _check (double x, double y, double res, mp_rnd_t rnd_mode, unsigned int px,
! 47: unsigned int py, unsigned int pz)
! 48: {
! 49: double z1, z2; mpfr_t xx, yy, zz;
! 50:
! 51: mpfr_init2 (xx, px);
! 52: mpfr_init2 (yy, py);
! 53: mpfr_init2 (zz, pz);
! 54: mpfr_set_d(xx, x, rnd_mode);
! 55: mpfr_set_d(yy, y, rnd_mode);
! 56: mpfr_mul(zz, xx, yy, rnd_mode);
! 57: #ifdef MPFR_HAVE_FESETROUND
! 58: mpfr_set_machine_rnd_mode(rnd_mode);
! 59: #endif
! 60: z1 = (res==0.0) ? x*y : res;
! 61: z2 = mpfr_get_d1 (zz);
! 62: if (z1!=z2 && (z1>=MINNORM || z1<=-MINNORM)) {
! 63: printf("mpfr_mul ");
! 64: if (res==0.0) printf("differs from libm.a"); else printf("failed");
! 65: printf(" for x=%1.20e y=%1.20e with rnd_mode=%s\n", x, y,
! 66: mpfr_print_rnd_mode(rnd_mode));
! 67: printf("libm.a gives %1.20e, mpfr_mul gives %1.20e\n", z1, z2);
! 68: if (res!=0.0) exit(1);
! 69: }
! 70: mpfr_clear(xx); mpfr_clear(yy); mpfr_clear(zz);
! 71: }
! 72:
! 73: void
! 74: check53 (double x, double y, mp_rnd_t rnd_mode, double z1)
! 75: {
! 76: double z2; mpfr_t xx, yy, zz;
! 77:
! 78: mpfr_init2 (xx, 53);
! 79: mpfr_init2 (yy, 53);
! 80: mpfr_init2 (zz, 53);
! 81: mpfr_set_d (xx, x, rnd_mode);
! 82: mpfr_set_d (yy, y, rnd_mode);
! 83: mpfr_mul (zz, xx, yy, rnd_mode);
! 84: z2 = mpfr_get_d1 (zz);
! 85: if (z1!=z2 && (!isnan(z1) || !isnan(z2))) {
! 86: printf("mpfr_mul failed for x=%1.20e y=%1.20e with rnd_mode=%s\n",
! 87: x, y, mpfr_print_rnd_mode(rnd_mode));
! 88: printf("libm.a gives %1.20e, mpfr_mul gives %1.20e\n", z1, z2);
! 89: exit(1);
! 90: }
! 91: mpfr_clear(xx); mpfr_clear(yy); mpfr_clear(zz);
! 92: }
! 93:
! 94: /* checks that x*y gives the same results in double
! 95: and with mpfr with 24 bits of precision */
! 96: void
! 97: check24 (float x, float y, mp_rnd_t rnd_mode, float z1)
! 98: {
! 99: float z2;
! 100: mpfr_t xx, yy, zz;
! 101:
! 102: mpfr_init2 (xx, 24);
! 103: mpfr_init2 (yy, 24);
! 104: mpfr_init2 (zz, 24);
! 105: mpfr_set_d (xx, x, rnd_mode);
! 106: mpfr_set_d (yy, y, rnd_mode);
! 107: mpfr_mul (zz, xx, yy, rnd_mode);
! 108: z2 = (float) mpfr_get_d1 (zz);
! 109: if (z1 != z2)
! 110: {
! 111: fprintf (stderr, "mpfr_mul failed for x=%1.0f y=%1.0f with prec=24 and"
! 112: "rnd=%s\n", x, y, mpfr_print_rnd_mode(rnd_mode));
! 113: fprintf (stderr, "libm.a gives %.10e, mpfr_mul gives %.10e\n", z1, z2);
! 114: exit (1);
! 115: }
! 116: mpfr_clear(xx);
! 117: mpfr_clear(yy);
! 118: mpfr_clear(zz);
! 119: }
! 120:
! 121: /* the following examples come from the paper "Number-theoretic Test
! 122: Generation for Directed Rounding" from Michael Parks, Table 1 */
! 123: void
! 124: check_float (void)
! 125: {
! 126: float z;
! 127:
! 128: check24(8388609.0, 8388609.0, GMP_RNDN, 70368760954880.0);
! 129: check24(16777213.0, 8388609.0, GMP_RNDN, 140737479966720.0);
! 130: check24(8388611.0, 8388609.0, GMP_RNDN, 70368777732096.0);
! 131: check24(12582911.0, 8388610.0, GMP_RNDN, 105553133043712.0);
! 132: check24(12582914.0, 8388610.0, GMP_RNDN, 105553158209536.0);
! 133: check24(13981013.0, 8388611.0, GMP_RNDN, 117281279442944.0);
! 134: check24(11184811.0, 8388611.0, GMP_RNDN, 93825028587520.0);
! 135: check24(11184810.0, 8388611.0, GMP_RNDN, 93825020198912.0);
! 136: check24(13981014.0, 8388611.0, GMP_RNDN, 117281287831552.0);
! 137:
! 138: check24(8388609.0, 8388609.0, GMP_RNDZ, 70368760954880.0);
! 139: check24(16777213.0, 8388609.0, GMP_RNDZ, 140737471578112.0);
! 140: z = 70368777732096.0;
! 141: check24(8388611.0, 8388609.0, GMP_RNDZ, z);
! 142: check24(12582911.0, 8388610.0, GMP_RNDZ, 105553124655104.0);
! 143: check24(12582914.0, 8388610.0, GMP_RNDZ, 105553158209536.0);
! 144: check24(13981013.0, 8388611.0, GMP_RNDZ, 117281271054336.0);
! 145: check24(11184811.0, 8388611.0, GMP_RNDZ, 93825028587520.0);
! 146: check24(11184810.0, 8388611.0, GMP_RNDZ, 93825011810304.0);
! 147: check24(13981014.0, 8388611.0, GMP_RNDZ, 117281287831552.0);
! 148:
! 149: check24(8388609.0, 8388609.0, GMP_RNDU, 70368769343488.0);
! 150: check24(16777213.0, 8388609.0, GMP_RNDU, 140737479966720.0);
! 151: check24(8388611.0, 8388609.0, GMP_RNDU, 70368786120704.0);
! 152: check24(12582911.0, 8388610.0, GMP_RNDU, 105553133043712.0);
! 153: check24(12582914.0, 8388610.0, GMP_RNDU, 105553166598144.0);
! 154: check24(13981013.0, 8388611.0, GMP_RNDU, 117281279442944.0);
! 155: check24(11184811.0, 8388611.0, GMP_RNDU, 93825036976128.0);
! 156: check24(11184810.0, 8388611.0, GMP_RNDU, 93825020198912.0);
! 157: check24(13981014.0, 8388611.0, GMP_RNDU, 117281296220160.0);
! 158:
! 159: check24(8388609.0, 8388609.0, GMP_RNDD, 70368760954880.0);
! 160: check24(16777213.0, 8388609.0, GMP_RNDD, 140737471578112.0);
! 161: check24(8388611.0, 8388609.0, GMP_RNDD, 70368777732096.0);
! 162: check24(12582911.0, 8388610.0, GMP_RNDD, 105553124655104.0);
! 163: check24(12582914.0, 8388610.0, GMP_RNDD, 105553158209536.0);
! 164: check24(13981013.0, 8388611.0, GMP_RNDD, 117281271054336.0);
! 165: check24(11184811.0, 8388611.0, GMP_RNDD, 93825028587520.0);
! 166: check24(11184810.0, 8388611.0, GMP_RNDD, 93825011810304.0);
! 167: check24(13981014.0, 8388611.0, GMP_RNDD, 117281287831552.0);
! 168: }
! 169:
! 170: /* check sign of result */
! 171: void
! 172: check_sign (void)
! 173: {
! 174: mpfr_t a, b;
! 175:
! 176: mpfr_init2 (a, 53);
! 177: mpfr_init2 (b, 53);
! 178: mpfr_set_d(a, -1.0, GMP_RNDN);
! 179: mpfr_set_d(b, 2.0, GMP_RNDN);
! 180: mpfr_mul(a, b, b, GMP_RNDN);
! 181: if (mpfr_get_d1 (a) != 4.0) {
! 182: fprintf(stderr,"2.0*2.0 gives %1.20e\n", mpfr_get_d1 (a)); exit(1);
! 183: }
! 184: mpfr_clear(a); mpfr_clear(b);
! 185: }
! 186:
! 187: /* checks that the inexact return value is correct */
! 188: void
! 189: check_exact (void)
! 190: {
! 191: mpfr_t a, b, c, d;
! 192: mp_prec_t prec;
! 193: int i, inexact;
! 194: mp_rnd_t rnd;
! 195:
! 196: mpfr_init (a);
! 197: mpfr_init (b);
! 198: mpfr_init (c);
! 199: mpfr_init (d);
! 200:
! 201: mpfr_set_prec (a, 17);
! 202: mpfr_set_prec (b, 17);
! 203: mpfr_set_prec (c, 32);
! 204: mpfr_set_str_raw (a, "1.1000111011000100e-1");
! 205: mpfr_set_str_raw (b, "1.0010001111100111e-1");
! 206: if (mpfr_mul (c, a, b, GMP_RNDZ))
! 207: {
! 208: fprintf (stderr, "wrong return value (1)\n");
! 209: exit (1);
! 210: }
! 211:
! 212: for (prec = 2; prec < 100; prec++)
! 213: {
! 214: mpfr_set_prec (a, prec);
! 215: mpfr_set_prec (b, prec);
! 216: mpfr_set_prec (c, 2 * prec - 2);
! 217: mpfr_set_prec (d, 2 * prec);
! 218: for (i = 0; i < 1000; i++)
! 219: {
! 220: mpfr_random (a);
! 221: mpfr_random (b);
! 222: rnd = LONG_RAND() % 4;
! 223: inexact = mpfr_mul (c, a, b, rnd);
! 224: if (mpfr_mul (d, a, b, rnd)) /* should be always exact */
! 225: {
! 226: fprintf (stderr, "unexpected inexact return value\n");
! 227: exit (1);
! 228: }
! 229: if ((inexact == 0) && mpfr_cmp (c, d))
! 230: {
! 231: fprintf (stderr, "inexact=0 but results differ\n");
! 232: exit (1);
! 233: }
! 234: else if (inexact && (mpfr_cmp (c, d) == 0))
! 235: {
! 236: fprintf (stderr, "inexact!=0 but results agree\n");
! 237: fprintf (stderr, "prec=%u rnd=%s a=", (unsigned int) prec,
! 238: mpfr_print_rnd_mode (rnd));
! 239: mpfr_out_str (stderr, 2, 0, a, rnd);
! 240: fprintf (stderr, "\nb=");
! 241: mpfr_out_str (stderr, 2, 0, b, rnd);
! 242: fprintf (stderr, "\nc=");
! 243: mpfr_out_str (stderr, 2, 0, c, rnd);
! 244: fprintf (stderr, "\nd=");
! 245: mpfr_out_str (stderr, 2, 0, d, rnd);
! 246: fprintf (stderr, "\n");
! 247: exit (1);
! 248: }
! 249: }
! 250: }
! 251:
! 252: mpfr_clear (a);
! 253: mpfr_clear (b);
! 254: mpfr_clear (c);
! 255: mpfr_clear (d);
! 256: }
! 257:
! 258: void
! 259: check_max(void)
! 260: {
! 261: mpfr_t xx, yy, zz;
! 262:
! 263: mpfr_init2(xx, 4);
! 264: mpfr_init2(yy, 4);
! 265: mpfr_init2(zz, 4);
! 266: mpfr_set_d(xx, 11.0/16, GMP_RNDN);
! 267: mpfr_mul_2si(xx, xx, MPFR_EMAX_DEFAULT/2, GMP_RNDN);
! 268: mpfr_set_d(yy, 11.0/16, GMP_RNDN);
! 269: mpfr_mul_2si(yy, yy, MPFR_EMAX_DEFAULT - MPFR_EMAX_DEFAULT/2 + 1, GMP_RNDN);
! 270: mpfr_clear_flags();
! 271: mpfr_mul(zz, xx, yy, GMP_RNDU);
! 272: if (!(mpfr_overflow_p() && MPFR_IS_INF(zz)))
! 273: {
! 274: printf("check_max failed (should be an overflow)\n");
! 275: exit(1);
! 276: }
! 277:
! 278: mpfr_clear_flags();
! 279: mpfr_mul(zz, xx, yy, GMP_RNDD);
! 280: if (mpfr_overflow_p() || MPFR_IS_INF(zz))
! 281: {
! 282: printf("check_max failed (should NOT be an overflow)\n");
! 283: exit(1);
! 284: }
! 285: mpfr_set_d(xx, 15.0/16, GMP_RNDN);
! 286: mpfr_mul_2si(xx, xx, MPFR_EMAX_DEFAULT, GMP_RNDN);
! 287: if (!(MPFR_IS_FP(xx) && MPFR_IS_FP(zz)))
! 288: {
! 289: printf("check_max failed (internal error)\n");
! 290: exit(1);
! 291: }
! 292: if (mpfr_cmp(xx, zz) != 0)
! 293: {
! 294: printf("check_max failed: got ");
! 295: mpfr_out_str(stdout, 2, 0, zz, GMP_RNDZ);
! 296: printf(" instead of ");
! 297: mpfr_out_str(stdout, 2, 0, xx, GMP_RNDZ);
! 298: printf("\n");
! 299: exit(1);
! 300: }
! 301:
! 302: mpfr_clear(xx);
! 303: mpfr_clear(yy);
! 304: mpfr_clear(zz);
! 305: }
! 306:
! 307: void
! 308: check_min(void)
! 309: {
! 310: mpfr_t xx, yy, zz;
! 311:
! 312: mpfr_init2(xx, 4);
! 313: mpfr_init2(yy, 4);
! 314: mpfr_init2(zz, 3);
! 315: mpfr_set_d(xx, 15.0/16, GMP_RNDN);
! 316: mpfr_mul_2si(xx, xx, MPFR_EMIN_DEFAULT/2, GMP_RNDN);
! 317: mpfr_set_d(yy, 15.0/16, GMP_RNDN);
! 318: mpfr_mul_2si(yy, yy, MPFR_EMIN_DEFAULT - MPFR_EMIN_DEFAULT/2 - 1, GMP_RNDN);
! 319: mpfr_mul(zz, xx, yy, GMP_RNDD);
! 320: if (mpfr_sgn(zz) != 0)
! 321: {
! 322: printf("check_min failed: got ");
! 323: mpfr_out_str(stdout, 2, 0, zz, GMP_RNDZ);
! 324: printf(" instead of 0\n");
! 325: exit(1);
! 326: }
! 327:
! 328: mpfr_mul(zz, xx, yy, GMP_RNDU);
! 329: mpfr_set_d(xx, 0.5, GMP_RNDN);
! 330: mpfr_mul_2si(xx, xx, MPFR_EMIN_DEFAULT, GMP_RNDN);
! 331: if (mpfr_sgn(xx) <= 0)
! 332: {
! 333: printf("check_min failed (internal error)\n");
! 334: exit(1);
! 335: }
! 336: if (mpfr_cmp(xx, zz) != 0)
! 337: {
! 338: printf("check_min failed: got ");
! 339: mpfr_out_str(stdout, 2, 0, zz, GMP_RNDZ);
! 340: printf(" instead of ");
! 341: mpfr_out_str(stdout, 2, 0, xx, GMP_RNDZ);
! 342: printf("\n");
! 343: exit(1);
! 344: }
! 345:
! 346: mpfr_clear(xx);
! 347: mpfr_clear(yy);
! 348: mpfr_clear(zz);
! 349: }
! 350:
! 351: int
! 352: main (int argc, char *argv[])
! 353: {
! 354: #ifdef MPFR_HAVE_FESETROUND
! 355: double x, y, z;
! 356: int i, prec, rnd_mode;
! 357:
! 358: mpfr_test_init ();
! 359: #endif
! 360:
! 361: check_exact ();
! 362: check_float ();
! 363: #ifdef HAVE_INFS
! 364: check53 (0.0, DBL_POS_INF, GMP_RNDN, DBL_NAN);
! 365: check53(1.0, DBL_POS_INF, GMP_RNDN, DBL_POS_INF);
! 366: check53(-1.0, DBL_POS_INF, GMP_RNDN, DBL_NEG_INF);
! 367: check53(DBL_NAN, 0.0, GMP_RNDN, DBL_NAN);
! 368: check53(1.0, DBL_NAN, GMP_RNDN, DBL_NAN);
! 369: #endif
! 370: check53(6.9314718055994530941514e-1, 0.0, GMP_RNDZ, 0.0);
! 371: check53(0.0, 6.9314718055994530941514e-1, GMP_RNDZ, 0.0);
! 372: check_sign();
! 373: check53(-4.165000000e4, -0.00004801920768307322868063274915, GMP_RNDN, 2.0);
! 374: check53(2.71331408349172961467e-08, -6.72658901114033715233e-165, GMP_RNDZ,
! 375: -1.8251348697787782844e-172);
! 376: check53(0.31869277231188065, 0.88642843322303122, GMP_RNDZ,
! 377: 2.8249833483992453642e-1);
! 378: check(8.47622108205396074254e-01, 3.24039313247872939883e-01, GMP_RNDU,
! 379: 28, 45, 2, 0.375);
! 380: check(2.63978122803639081440e-01, 6.8378615379333496093e-1, GMP_RNDN,
! 381: 34, 23, 31, 0.180504585267044603);
! 382: check(1.0, 0.11835170935876249132, GMP_RNDU, 6, 41, 36, 0.1183517093595583);
! 383: check53(67108865.0, 134217729.0, GMP_RNDN, 9.007199456067584e15);
! 384: check(1.37399642157394197284e-01, 2.28877275604219221350e-01, GMP_RNDN,
! 385: 49, 15, 32, 0.0314472340833162888);
! 386: check(4.03160720978664954828e-01, 5.85483042917246621073e-01, GMP_RNDZ,
! 387: 51, 22, 32, 0.2360436821472831);
! 388: check(3.90798504668055102229e-14, 9.85394674650308388664e-04, GMP_RNDN,
! 389: 46, 22, 12, 0.385027296503914762e-16);
! 390: check(4.58687081072827851358e-01, 2.20543551472118792844e-01, GMP_RNDN,
! 391: 49, 3, 2, 0.09375);
! 392: check_max();
! 393: check_min();
! 394: #ifdef MPFR_HAVE_FESETROUND
! 395: SEED_RAND (time(NULL));
! 396: prec = (argc<2) ? 53 : atoi(argv[1]);
! 397: rnd_mode = (argc<3) ? -1 : atoi(argv[2]);
! 398: for (i=0;i<1000000;) {
! 399: x = drand();
! 400: y = drand();
! 401: z = x*y; if (z<0) z=-z;
! 402: if (z<1e+308 && z>1e-308) /* don't test overflow/underflow for now */
! 403: { i++;
! 404: check(x, y, (rnd_mode==-1) ? LONG_RAND()%4 : rnd_mode,
! 405: prec, prec, prec, 0.0);
! 406: }
! 407: }
! 408: #endif
! 409:
! 410: return 0;
! 411: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to tmul.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