Annotation of OpenXM_contrib/gmp/mpfr/tests/tsqrt.c, Revision 1.1
1.1 ! ohara 1: /* Test file for mpfr_sqrt.
! 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 <math.h>
! 23: #include <stdio.h>
! 24: #include <stdlib.h>
! 25: #include <time.h>
! 26: #include "gmp.h"
! 27: #include "gmp-impl.h"
! 28: #include "mpfr.h"
! 29: #include "mpfr-impl.h"
! 30: #include "mpfr-test.h"
! 31:
! 32: #define check(a,r) check3(a,r,-1.0)
! 33:
! 34: int maxulp=0;
! 35:
! 36: void check3 _PROTO((double, mp_rnd_t, double));
! 37: void check4 _PROTO((double, mp_rnd_t, char *));
! 38: void check24 _PROTO((float, mp_rnd_t, float));
! 39: void check_float _PROTO((void));
! 40: void special _PROTO((void));
! 41: void check_inexact _PROTO((mp_prec_t));
! 42: void check_nan _PROTO((void));
! 43:
! 44: void
! 45: check3 (double a, mp_rnd_t rnd_mode, double Q)
! 46: {
! 47: mpfr_t q; double Q2; int ck,u;
! 48:
! 49: ck = (Q!=-1.0); /* if ck=1, then Q is certified correct */
! 50: mpfr_init2(q, 53);
! 51: mpfr_set_d(q, a, rnd_mode);
! 52: #ifdef MPFR_HAVE_FESETROUND
! 53: mpfr_set_machine_rnd_mode(rnd_mode);
! 54: #endif
! 55: mpfr_sqrt(q, q, rnd_mode);
! 56: if (ck==0) Q = sqrt(a);
! 57: else {
! 58: if (Q != sqrt(a) && (!isnan(Q) || !isnan(sqrt(a)))) {
! 59: fprintf(stderr, "you've found a bug in your machine's sqrt for x=%1.20e\n", a);
! 60: mpfr_clear(q);
! 61: exit(1);
! 62:
! 63: }
! 64: }
! 65: Q2 = mpfr_get_d1 (q);
! 66: if (Q!=Q2 && (!isnan(Q) || !isnan(Q2))) {
! 67: u = ulp(Q2,Q);
! 68: if (ck) {
! 69: printf("mpfr_sqrt failed for a=%1.20e, rnd_mode=%s\n",
! 70: a, mpfr_print_rnd_mode(rnd_mode));
! 71: printf("expected sqrt is %1.20e, got %1.20e (%d ulp)\n",Q,Q2,u);
! 72: mpfr_clear(q);
! 73: exit(1);
! 74: }
! 75: else if (u>maxulp || u<-maxulp) {
! 76: maxulp = (u>maxulp) ? u : -u;
! 77: printf("libm.a differs from mpfr_sqrt for a=%1.20e, rnd_mode=%s\n",
! 78: a, mpfr_print_rnd_mode(rnd_mode));
! 79: printf("libm.a gives %1.20e, mpfr_sqrt gives %1.20e (%d ulp)\n",Q,Q2,u);
! 80: }
! 81: }
! 82: mpfr_clear(q);
! 83: }
! 84:
! 85: void
! 86: check4 (double a, mp_rnd_t rnd_mode, char *Q)
! 87: {
! 88: mpfr_t q, res;
! 89:
! 90: mpfr_init2(q, 53); mpfr_init2(res, 53);
! 91: mpfr_set_d(q, a, rnd_mode);
! 92: mpfr_sqrt(q, q, rnd_mode);
! 93: mpfr_set_str(res, Q, 16, GMP_RNDN);
! 94: if (mpfr_cmp(q, res)) {
! 95: printf("mpfr_sqrt failed for a=%1.20e, rnd_mode=%s\n",
! 96: a, mpfr_print_rnd_mode(rnd_mode));
! 97: printf("expected "); mpfr_print_binary(res); putchar('\n');
! 98: printf("got "); mpfr_print_binary(q); putchar('\n');
! 99: mpfr_clear(q); mpfr_clear(res);
! 100: exit(1);
! 101: }
! 102: mpfr_clear(res);
! 103: mpfr_clear(q);
! 104: }
! 105:
! 106: void
! 107: check24 (float a, mp_rnd_t rnd_mode, float Q)
! 108: {
! 109: mpfr_t q; float Q2;
! 110:
! 111: mpfr_init2(q, 24);
! 112: mpfr_set_d(q, a, rnd_mode);
! 113: mpfr_sqrt(q, q, rnd_mode);
! 114: Q2 = mpfr_get_d1 (q);
! 115: if (Q!=Q2) {
! 116: printf("mpfr_sqrt failed for a=%1.10e, prec=24, rnd_mode=%s\n",
! 117: a, mpfr_print_rnd_mode(rnd_mode));
! 118: printf("expected sqrt is %1.10e, got %1.10e\n",Q,Q2);
! 119: exit(1);
! 120: }
! 121: mpfr_clear(q);
! 122: }
! 123:
! 124: /* the following examples come from the paper "Number-theoretic Test
! 125: Generation for Directed Rounding" from Michael Parks, Table 3 */
! 126: void
! 127: check_float (void)
! 128: {
! 129: float b = 8388608.0; /* 2^23 */
! 130:
! 131: check24(b*8388610.0, GMP_RNDN, 8.388609e6);
! 132: check24(b*2.0*16777214.0, GMP_RNDN, 1.6777215e7);
! 133: check24(b*8388612.0, GMP_RNDN, 8.388610e6);
! 134: check24(b*2.0*16777212.0, GMP_RNDN, 1.6777214e7);
! 135: check24(b*11946704.0, GMP_RNDN, 1.0010805e7);
! 136: check24(b*14321479.0, GMP_RNDN, 1.0960715e7);
! 137: check24(b*2.0*13689673.0, GMP_RNDN, 1.5155019e7);
! 138: check24(b*8388614.0, GMP_RNDN, 8.388611e6);
! 139: check24(b*2.0*16777210.0, GMP_RNDN, 1.6777213e7);
! 140: check24(b*10873622.0, GMP_RNDN, 9.550631e6);
! 141:
! 142: check24(b*8388610.0, GMP_RNDZ, 8.388608e6);
! 143: check24(b*2.0*16777214.0, GMP_RNDZ, 1.6777214e7);
! 144: check24(b*8388612.0, GMP_RNDZ, 8.388609e6);
! 145: check24(b*2.0*16777212.0, GMP_RNDZ, 1.6777213e7);
! 146: check24(b*11946704.0, GMP_RNDZ, 1.0010805e7);
! 147: check24(b*14321479.0, GMP_RNDZ, 1.0960715e7);
! 148: check24(b*2.0*13689673.0, GMP_RNDZ, 1.5155019e7);
! 149: check24(b*8388614.0, GMP_RNDZ, 8.38861e6);
! 150: check24(b*2.0*16777210.0, GMP_RNDZ, 1.6777212e7);
! 151: check24(b*10873622.0, GMP_RNDZ, 9.550631e6);
! 152:
! 153: check24(b*8388610.0, GMP_RNDU, 8.388609e6);
! 154: check24(b*2.0*16777214.0, GMP_RNDU, 1.6777215e7);
! 155: check24(b*8388612.0, GMP_RNDU, 8.388610e6);
! 156: check24(b*2.0*16777212.0, GMP_RNDU, 1.6777214e7);
! 157: check24(b*11946704.0, GMP_RNDU, 1.0010806e7);
! 158: check24(b*14321479.0, GMP_RNDU, 1.0960716e7);
! 159: check24(b*2.0*13689673.0, GMP_RNDU, 1.515502e7);
! 160: check24(b*8388614.0, GMP_RNDU, 8.388611e6);
! 161: check24(b*2.0*16777210.0, GMP_RNDU, 1.6777213e7);
! 162: check24(b*10873622.0, GMP_RNDU, 9.550632e6);
! 163:
! 164: check24(b*8388610.0, GMP_RNDD, 8.388608e6);
! 165: check24(b*2.0*16777214.0, GMP_RNDD, 1.6777214e7);
! 166: check24(b*8388612.0, GMP_RNDD, 8.388609e6);
! 167: check24(b*2.0*16777212.0, GMP_RNDD, 1.6777213e7);
! 168: check24(b*11946704.0, GMP_RNDD, 1.0010805e7);
! 169: check24(b*14321479.0, GMP_RNDD, 1.0960715e7);
! 170: check24(b*2.0*13689673.0, GMP_RNDD, 1.5155019e7);
! 171: check24(b*8388614.0, GMP_RNDD, 8.38861e6);
! 172: check24(b*2.0*16777210.0, GMP_RNDD, 1.6777212e7);
! 173: check24(b*10873622.0, GMP_RNDD, 9.550631e6);
! 174: }
! 175:
! 176: void
! 177: special (void)
! 178: {
! 179: mpfr_t x, z;
! 180: int inexact;
! 181: mp_prec_t p;
! 182:
! 183: mpfr_init (x);
! 184: mpfr_init (z);
! 185:
! 186: mpfr_set_prec (x, 27);
! 187: mpfr_set_str_raw (x, "0.110100111010101000010001011");
! 188: if ((inexact = mpfr_sqrt (x, x, GMP_RNDZ)) >= 0)
! 189: {
! 190: fprintf (stderr, "Wrong inexact flag: expected -1, got %d\n", inexact);
! 191: exit (1);
! 192: }
! 193:
! 194: mpfr_set_prec (x, 2);
! 195: for (p=2; p<1000; p++)
! 196: {
! 197: mpfr_set_prec (z, p);
! 198: mpfr_set_ui (z, 1, GMP_RNDN);
! 199: mpfr_add_one_ulp (z, GMP_RNDN);
! 200: mpfr_sqrt (x, z, GMP_RNDU);
! 201: if (mpfr_get_d1 (x) != 1.5)
! 202: {
! 203: fprintf (stderr, "Error: sqrt(1+ulp(1), up) should give 1.5 (prec=%u)\n", (unsigned) p);
! 204: printf ("got "); mpfr_print_binary (x); putchar ('\n');
! 205: exit (1);
! 206: }
! 207: }
! 208:
! 209: /* check inexact flag */
! 210: mpfr_set_prec (x, 5);
! 211: mpfr_set_str_raw (x, "1.1001E-2");
! 212: if ((inexact = mpfr_sqrt (x, x, GMP_RNDN)))
! 213: {
! 214: fprintf (stderr, "Wrong inexact flag: expected 0, got %d\n", inexact);
! 215: exit (1);
! 216: }
! 217:
! 218: mpfr_set_prec (x, 2);
! 219: mpfr_set_prec (z, 2);
! 220:
! 221: /* checks the sign is correctly set */
! 222: mpfr_set_d (x, 1.0, GMP_RNDN);
! 223: mpfr_set_d (z, -1.0, GMP_RNDN);
! 224: mpfr_sqrt (z, x, GMP_RNDN);
! 225: if (mpfr_cmp_ui (z, 0) < 0) {
! 226: fprintf (stderr, "Error: square root of %e gives %e\n",
! 227: mpfr_get_d1 (x), mpfr_get_d1 (z));
! 228: exit (1);
! 229: }
! 230:
! 231: mpfr_set_prec (x, 192);
! 232: mpfr_set_prec (z, 160);
! 233: mpfr_set_str_raw (z, "0.1011010100000100100100100110011001011100100100000011000111011001011101101101110000110100001000100001100001011000E1");
! 234: mpfr_set_prec (x, 160);
! 235: mpfr_sqrt(x, z, GMP_RNDN);
! 236: mpfr_sqrt(z, x, GMP_RNDN);
! 237:
! 238: mpfr_clear (x);
! 239: mpfr_clear (z);
! 240: }
! 241:
! 242: void
! 243: check_inexact (mp_prec_t p)
! 244: {
! 245: mpfr_t x, y, z;
! 246: mp_rnd_t rnd;
! 247: int inexact, sign;
! 248:
! 249: mpfr_init2 (x, p);
! 250: mpfr_init2 (y, p);
! 251: mpfr_init2 (z, 2*p);
! 252: mpfr_random (x);
! 253: rnd = LONG_RAND() % 4;
! 254: inexact = mpfr_sqrt (y, x, rnd);
! 255: if (mpfr_mul (z, y, y, rnd)) /* exact since prec(z) = 2*prec(y) */
! 256: {
! 257: fprintf (stderr, "Error: multiplication should be exact\n");
! 258: exit (1);
! 259: }
! 260: mpfr_sub (z, z, x, rnd); /* exact also */
! 261: sign = mpfr_cmp_ui (z, 0);
! 262: if (((inexact == 0) && (sign)) ||
! 263: ((inexact > 0) && (sign <= 0)) ||
! 264: ((inexact < 0) && (sign >= 0)))
! 265: {
! 266: fprintf (stderr, "Error: wrong inexact flag, expected %d, got %d\n",
! 267: sign, inexact);
! 268: printf ("x=");
! 269: mpfr_print_binary (x);
! 270: printf (" rnd=%s\n", mpfr_print_rnd_mode (rnd));
! 271: printf ("y="); mpfr_print_binary (y); putchar ('\n');
! 272: exit (1);
! 273: }
! 274: mpfr_clear (x);
! 275: mpfr_clear (y);
! 276: mpfr_clear (z);
! 277: }
! 278:
! 279: void
! 280: check_nan (void)
! 281: {
! 282: mpfr_t x, got;
! 283:
! 284: mpfr_init2 (x, 100L);
! 285: mpfr_init2 (got, 100L);
! 286:
! 287: /* sqrt(NaN) == NaN */
! 288: MPFR_CLEAR_FLAGS (x);
! 289: MPFR_SET_NAN (x);
! 290: ASSERT_ALWAYS (mpfr_sqrt (got, x, GMP_RNDZ) == 0); /* exact */
! 291: ASSERT_ALWAYS (mpfr_nan_p (got));
! 292:
! 293: /* sqrt(-1) == NaN */
! 294: mpfr_set_si (x, -1L, GMP_RNDZ);
! 295: ASSERT_ALWAYS (mpfr_sqrt (got, x, GMP_RNDZ) == 0); /* exact */
! 296: ASSERT_ALWAYS (mpfr_nan_p (got));
! 297:
! 298: /* sqrt(+inf) == +inf */
! 299: MPFR_CLEAR_FLAGS (x);
! 300: MPFR_SET_INF (x);
! 301: MPFR_SET_POS (x);
! 302: ASSERT_ALWAYS (mpfr_sqrt (got, x, GMP_RNDZ) == 0); /* exact */
! 303: ASSERT_ALWAYS (mpfr_inf_p (got));
! 304:
! 305: /* sqrt(-inf) == NaN */
! 306: MPFR_CLEAR_FLAGS (x);
! 307: MPFR_SET_INF (x);
! 308: MPFR_SET_NEG (x);
! 309: ASSERT_ALWAYS (mpfr_sqrt (got, x, GMP_RNDZ) == 0); /* exact */
! 310: ASSERT_ALWAYS (mpfr_nan_p (got));
! 311:
! 312: /* sqrt(-0) == 0 */
! 313: mpfr_set_si (x, 0L, GMP_RNDZ);
! 314: MPFR_SET_NEG (x);
! 315: ASSERT_ALWAYS (mpfr_sqrt (got, x, GMP_RNDZ) == 0); /* exact */
! 316: ASSERT_ALWAYS (mpfr_number_p (got));
! 317: ASSERT_ALWAYS (mpfr_cmp_ui (got, 0L) == 0);
! 318:
! 319: mpfr_clear (x);
! 320: mpfr_clear (got);
! 321: }
! 322:
! 323: double five = 5.0;
! 324:
! 325: int
! 326: main (void)
! 327: {
! 328: double a;
! 329: mp_prec_t p;
! 330: int k;
! 331: #ifdef MPFR_HAVE_FESETROUND
! 332: int i;
! 333:
! 334: mpfr_test_init ();
! 335:
! 336: /* On Debian potato glibc 2.1.3-18, sqrt() doesn't seem to respect
! 337: fesetround. */
! 338: {
! 339: double a, b;
! 340: mpfr_set_machine_rnd_mode (GMP_RNDU);
! 341: a = sqrt (five);
! 342: mpfr_set_machine_rnd_mode (GMP_RNDD);
! 343: b = sqrt (five);
! 344: if (a == b)
! 345: {
! 346: printf ("Tests suppressed, mpfr_set_machine_rnd_mode doesn't affect sqrt()\n");
! 347: goto nogood;
! 348: }
! 349: }
! 350:
! 351: SEED_RAND (time(NULL));
! 352: for (i=0;i<100000;i++)
! 353: {
! 354: a = drand();
! 355: if (a < 0.0) a = -a; /* ensures a is positive */
! 356: check (a, LONG_RAND() % 4);
! 357: }
! 358: nogood:
! 359: #endif
! 360:
! 361: check_nan ();
! 362:
! 363: for (p=2; p<200; p++)
! 364: for (k=0; k<200; k++)
! 365: check_inexact (p);
! 366: special ();
! 367: check_float();
! 368: #ifdef HAVE_INFS
! 369: check3 (DBL_NAN, GMP_RNDN, DBL_NAN);
! 370: check3 (-1.0, GMP_RNDN, DBL_NAN);
! 371: check3 (DBL_POS_INF, GMP_RNDN, DBL_POS_INF);
! 372: check3 (DBL_NEG_INF, GMP_RNDN, DBL_NAN);
! 373: #endif
! 374: check3(-0.0, GMP_RNDN, 0.0);
! 375: check4(6.37983013646045901440e+32, GMP_RNDN, "5.9bc5036d09e0c@13");
! 376: check4(1.0, GMP_RNDN, "1");
! 377: check4(1.0, GMP_RNDZ, "1");
! 378: check4(3.725290298461914062500000e-9, GMP_RNDN, "4@-4");
! 379: check4(3.725290298461914062500000e-9, GMP_RNDZ, "4@-4");
! 380: a=1190456976439861.0;
! 381: check4(a, GMP_RNDZ, "2.0e7957873529a@6");
! 382: check4(1024.0*a, GMP_RNDZ, "4.1cf2af0e6a534@7");
! 383: /* the following examples are bugs in Cygnus compiler/system, found by
! 384: Fabrice Rouillier while porting mpfr to Windows */
! 385: check4(9.89438396044940256501e-134, GMP_RNDU, "8.7af7bf0ebbee@-56");
! 386: check4(7.86528588050363751914e+31, GMP_RNDZ, "1.f81fc40f32062@13");
! 387: check4(0.99999999999999988897, GMP_RNDN, "f.ffffffffffff8@-1");
! 388: check4(1.00000000000000022204, GMP_RNDN, "1");
! 389: /* the following examples come from the paper "Number-theoretic Test
! 390: Generation for Directed Rounding" from Michael Parks, Table 4 */
! 391: a = 4503599627370496.0; /* 2^52 */
! 392:
! 393: check4(a*2.0*8732221479794286.0, GMP_RNDN, "1.f81fc40f32063@13");
! 394: check4(a*8550954388695124.0, GMP_RNDN, "1.60c012a92fc65@13");
! 395: check4(a*7842344481681754.0, GMP_RNDN, "1.51d17526c7161@13");
! 396: check4(a*5935035262218600.0, GMP_RNDN, "1.25e19302f7e51@13");
! 397: check4(a*5039650445085418.0, GMP_RNDN, "1.0ecea7dd2ec3d@13");
! 398: check4(a*5039721545366078.0, GMP_RNDN, "1.0ecf250e8e921@13");
! 399: check4(a*8005963117781324.0, GMP_RNDN, "1.5552f3eedcf33@13");
! 400: check4(a*6703494707970582.0, GMP_RNDN, "1.3853ee10c9c99@13");
! 401: check4(a*8010323124937260.0, GMP_RNDN, "1.556abe212b56f@13");
! 402: check4(a*2.0*8010776873384260.0, GMP_RNDN, "1.e2d9a51977e6e@13");
! 403:
! 404: check4(a*2.0*8732221479794286.0, GMP_RNDZ, "1.f81fc40f32062@13");
! 405: check4(a*8550954388695124.0, GMP_RNDZ, "1.60c012a92fc64@13");
! 406: check4(a*7842344481681754.0, GMP_RNDZ, "1.51d17526c716@13");
! 407: check4(a*5935035262218600.0, GMP_RNDZ, "1.25e19302f7e5@13");
! 408: check4(a*5039650445085418.0, GMP_RNDZ, "1.0ecea7dd2ec3c@13");
! 409: check4(a*5039721545366078.0, GMP_RNDZ, "1.0ecf250e8e92@13");
! 410: check4(a*8005963117781324.0, GMP_RNDZ, "1.5552f3eedcf32@13");
! 411: check4(a*6703494707970582.0, GMP_RNDZ, "1.3853ee10c9c98@13");
! 412: check4(a*8010323124937260.0, GMP_RNDZ, "1.556abe212b56e@13");
! 413: check4(a*2.0*8010776873384260.0, GMP_RNDZ, "1.e2d9a51977e6d@13");
! 414:
! 415: check4(a*2.0*8732221479794286.0, GMP_RNDU, "1.f81fc40f32063@13");
! 416: check4(a*8550954388695124.0, GMP_RNDU, "1.60c012a92fc65@13");
! 417: check4(a*7842344481681754.0, GMP_RNDU, "1.51d17526c7161@13");
! 418: check4(a*5935035262218600.0, GMP_RNDU, "1.25e19302f7e51@13");
! 419: check4(a*5039650445085418.0, GMP_RNDU, "1.0ecea7dd2ec3d@13");
! 420: check4(a*5039721545366078.0, GMP_RNDU, "1.0ecf250e8e921@13");
! 421: check4(a*8005963117781324.0, GMP_RNDU, "1.5552f3eedcf33@13");
! 422: check4(a*6703494707970582.0, GMP_RNDU, "1.3853ee10c9c99@13");
! 423: check4(a*8010323124937260.0, GMP_RNDU, "1.556abe212b56f@13");
! 424: check4(a*2.0*8010776873384260.0, GMP_RNDU, "1.e2d9a51977e6e@13");
! 425:
! 426: check4(a*2.0*8732221479794286.0, GMP_RNDD, "1.f81fc40f32062@13");
! 427: check4(a*8550954388695124.0, GMP_RNDD, "1.60c012a92fc64@13");
! 428: check4(a*7842344481681754.0, GMP_RNDD, "1.51d17526c716@13");
! 429: check4(a*5935035262218600.0, GMP_RNDD, "1.25e19302f7e5@13");
! 430: check4(a*5039650445085418.0, GMP_RNDD, "1.0ecea7dd2ec3c@13");
! 431: check4(a*5039721545366078.0, GMP_RNDD, "1.0ecf250e8e92@13");
! 432: check4(a*8005963117781324.0, GMP_RNDD, "1.5552f3eedcf32@13");
! 433: check4(a*6703494707970582.0, GMP_RNDD, "1.3853ee10c9c98@13");
! 434: check4(a*8010323124937260.0, GMP_RNDD, "1.556abe212b56e@13");
! 435: check4(a*2.0*8010776873384260.0, GMP_RNDD, "1.e2d9a51977e6d@13");
! 436:
! 437: return 0;
! 438: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to tsqrt.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