Annotation of OpenXM_contrib/gmp/mpfr/tests/texp.c, Revision 1.1
1.1 ! ohara 1: /* Test file for mpfr_exp.
! 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: int check3 _PROTO((double, mp_rnd_t, double));
! 33: int check_large _PROTO((double, int, mp_rnd_t));
! 34: int check_worst_case _PROTO((double, double));
! 35: int check_worst_cases _PROTO((void));
! 36: void compare_exp2_exp3 _PROTO((int));
! 37:
! 38: int maxu=0;
! 39:
! 40: #define check(d, r) check3(d, r, 0.0)
! 41:
! 42: /* returns the number of ulp of error */
! 43: int
! 44: check3 (double d, mp_rnd_t rnd, double e)
! 45: {
! 46: mpfr_t x, y; double f; int u=0, ck=0;
! 47:
! 48: mpfr_init2(x, 53); mpfr_init2(y, 53);
! 49: #ifdef MPFR_HAVE_FESETROUND
! 50: mpfr_set_machine_rnd_mode(rnd);
! 51: #endif
! 52: if (e==0.0) e = exp(d); else ck=1; /* really check */
! 53: mpfr_set_d(x, d, rnd);
! 54: mpfr_exp(y, x, rnd);
! 55: f = mpfr_get_d1 (y);
! 56: if (f != e && (!isnan(f) || !isnan(e))) {
! 57: u = ulp(e, f);
! 58: if (u<0) {
! 59: if (u == (mp_limb_t) 1 << (mp_bits_per_limb-1)) u += 1;
! 60: u=-u;
! 61: }
! 62: if (u!=0) {
! 63: if (ck) {
! 64: printf("mpfr_exp failed for x=%1.20e, rnd=%s\n", d,
! 65: mpfr_print_rnd_mode(rnd));
! 66: printf("expected result is %1.20e, got %1.20e, dif=%d ulp\n",e,f,u);
! 67: exit(1);
! 68: }
! 69: else if (u>maxu) {
! 70: maxu=u;
! 71: printf("mpfr_exp differs from libm.a for x=%1.20e, rnd=%s\n",d,
! 72: mpfr_print_rnd_mode(rnd));
! 73: printf("libm.a gave %1.20e, mpfr_exp got %1.20e, dif=%d ulp\n",e,f,u);
! 74: }
! 75: }
! 76: }
! 77: mpfr_clear(x); mpfr_clear(y);
! 78: return u;
! 79: }
! 80:
! 81: /* computes n bits of exp(d) */
! 82: int
! 83: check_large (double d, int n, mp_rnd_t rnd)
! 84: {
! 85: mpfr_t x; mpfr_t y;
! 86:
! 87: mpfr_init2(x, n); mpfr_init2(y, n);
! 88: if (d==0.0) { /* try exp(Pi*sqrt(163)/3)-640320 */
! 89: mpfr_set_d(x, 163.0, rnd);
! 90: mpfr_sqrt(x, x, rnd);
! 91: mpfr_const_pi(y, rnd);
! 92: mpfr_mul(x, x, y, rnd);
! 93: mpfr_div_ui(x, x, 3, rnd);
! 94: }
! 95: else mpfr_set_d(x, d, rnd);
! 96: mpfr_exp (y, x, rnd);
! 97: if (d==0.0) {
! 98: mpfr_set_d(x, 640320.0, rnd);
! 99: mpfr_sub(y, y, x, rnd);
! 100: printf("exp(Pi*sqrt(163)/3)-640320=");
! 101: }
! 102: else printf("exp(%1.20e)=",d);
! 103: mpfr_out_str(stdout, 10, 0, y, rnd);
! 104: putchar('\n');
! 105: printf(" ="); mpfr_print_binary(y); putchar('\n');
! 106: if (n==53) printf(" =%1.20e\n", mpfr_get_d1 (y));
! 107:
! 108: mpfr_clear(x); mpfr_clear(y);
! 109: return 0;
! 110: }
! 111:
! 112: /* expx is the value of exp(X) rounded towards -infinity */
! 113: int
! 114: check_worst_case (double X, double expx)
! 115: {
! 116: mpfr_t x, y;
! 117:
! 118: mpfr_init2(x, 53); mpfr_init2(y, 53);
! 119: mpfr_set_d(x, X, GMP_RNDN);
! 120: mpfr_exp(y, x, GMP_RNDD);
! 121: if (mpfr_get_d1 (y) != expx) {
! 122: fprintf(stderr, "exp(x) rounded towards -infinity is wrong\n"); exit(1);
! 123: }
! 124: mpfr_exp(x, x, GMP_RNDN);
! 125: mpfr_set_d(x, X, GMP_RNDN);
! 126: mpfr_exp(x, x, GMP_RNDU);
! 127: mpfr_add_one_ulp(y, GMP_RNDN);
! 128: if (mpfr_cmp(x,y)) {
! 129: fprintf(stderr, "exp(x) rounded towards +infinity is wrong\n"); exit(1);
! 130: }
! 131: mpfr_clear(x); mpfr_clear(y);
! 132: return 0;
! 133: }
! 134:
! 135: /* worst cases communicated by Jean-Michel Muller and Vincent Lefevre */
! 136: int
! 137: check_worst_cases (void)
! 138: {
! 139: mpfr_t x; mpfr_t y;
! 140:
! 141: mpfr_init(x);
! 142: mpfr_set_prec (x, 53);
! 143:
! 144: check_worst_case(4.44089209850062517562e-16, 1.00000000000000022204);
! 145: check_worst_case(6.39488462184069720009e-14, 1.00000000000006372680);
! 146: check_worst_case(1.84741111297455401935e-12, 1.00000000000184718907);
! 147: check_worst_case(1.76177628026265550074e-10, 1.00000000017617751702);
! 148: check3(1.76177628026265550074e-10, GMP_RNDN, 1.00000000017617773906);
! 149: check_worst_case(7.54175277499595900852e-10, 1.00000000075417516676);
! 150: check3(7.54175277499595900852e-10, GMP_RNDN, 1.00000000075417538881);
! 151: /* bug found by Vincent Lefe`vre on December 8, 1999 */
! 152: check3(-5.42410311287441459172e+02, GMP_RNDN, 2.7176584868845723e-236);
! 153: /* further cases communicated by Vincent Lefe`vre on January 27, 2000 */
! 154: check3(-1.32920285897904911589e-10, GMP_RNDN, 0.999999999867079769622);
! 155: check3(-1.44037948245738330735e-10, GMP_RNDN, 0.9999999998559621072757);
! 156: check3(-1.66795910430705305937e-10, GMP_RNDZ, 0.9999999998332040895832);
! 157: check3(-1.64310953745426656203e-10, GMP_RNDN, 0.9999999998356891017792);
! 158: check3(-1.38323574826034659172e-10, GMP_RNDZ, 0.9999999998616764251835);
! 159: check3(-1.23621668465115401498e-10, GMP_RNDZ, 0.9999999998763783315425);
! 160:
! 161: mpfr_set_prec (x, 601);
! 162: mpfr_set_str_raw (x, "0.1000100010110110101110100101000100001110000100000100010100001110110111000010010110000111010010001011110010011101111111011101010001100110111100100001101101000111111011010010011001001100110111110010010010101010100011110110010010101111000111110011111110101101100111101100001000110000000111010100001111000000011101000011111101010011010010110101101010100010000000001001000111111111011011010011010100101101111101000101100011101111000110111010010100011001100000010001111011110110111101011011000100011000010100110101001101001111110110001111101000110010011101100100101000001010011011010010110100001101110100100E0");
! 163: mpfr_init2 (y, 601);
! 164: mpfr_exp_2 (y, x, GMP_RNDD);
! 165: mpfr_exp3 (x, x, GMP_RNDD);
! 166: if (mpfr_cmp (x, y)) {
! 167: fprintf (stderr, "mpfr_exp_2 and mpfr_exp3 for prec=601\n");
! 168: exit (1);
! 169: }
! 170:
! 171: mpfr_clear (x);
! 172: mpfr_clear (y);
! 173: return 0;
! 174: }
! 175:
! 176: void
! 177: compare_exp2_exp3 (int n)
! 178: {
! 179: mpfr_t x, y, z; int prec; mp_rnd_t rnd;
! 180:
! 181: mpfr_init (x);
! 182: mpfr_init (y);
! 183: mpfr_init (z);
! 184: for (prec=20; prec<=n; prec++)
! 185: {
! 186: mpfr_set_prec (x, prec);
! 187: mpfr_set_prec (y, prec);
! 188: mpfr_set_prec (z, prec);
! 189: mpfr_random (x);
! 190: rnd = LONG_RAND() % 4;
! 191: mpfr_exp_2 (y, x, rnd);
! 192: mpfr_exp3 (z, x, rnd);
! 193: if (mpfr_cmp (y,z))
! 194: {
! 195: printf ("mpfr_exp_2 and mpfr_exp3 disagree for rnd=%s and\nx=",
! 196: mpfr_print_rnd_mode (rnd));
! 197: mpfr_print_binary (x);
! 198: putchar ('\n');
! 199: printf ("mpfr_exp_2 gives ");
! 200: mpfr_print_binary (y);
! 201: putchar ('\n');
! 202: printf ("mpfr_exp3 gives ");
! 203: mpfr_print_binary (z);
! 204: putchar ('\n');
! 205: exit (1);
! 206: }
! 207: }
! 208: mpfr_clear (x);
! 209: mpfr_clear (y);
! 210: mpfr_clear (z);
! 211: }
! 212:
! 213: #define TEST_FUNCTION mpfr_exp
! 214: #include "tgeneric.c"
! 215:
! 216: int
! 217: main (int argc, char *argv[])
! 218: {
! 219: #ifdef MPFR_HAVE_FESETROUND
! 220: int i, N, s=0, e, maxe=0;
! 221: double lo, hi;
! 222: #endif
! 223: double d;
! 224:
! 225: test_generic (2, 100, 100);
! 226:
! 227: if (argc == 4)
! 228: {
! 229: check_large (atof(argv[1]), atoi(argv[2]), atoi(argv[3]));
! 230: exit(1);
! 231: }
! 232:
! 233: compare_exp2_exp3(500);
! 234: check_worst_cases();
! 235: check3(0.0, GMP_RNDU, 1.0);
! 236: check3(-8.88024741073346941839e-17, GMP_RNDU, 1.0);
! 237: check3(8.70772839244701057915e-01, GMP_RNDN, 2.38875626491680437269);
! 238: check3(1.0, GMP_RNDN, 2.71828182845904509080);
! 239: check3(-3.42135637628104173534e-07, GMP_RNDZ, 0.999999657864420798958);
! 240: /* worst case for argument reduction, very near from 5*log(2),
! 241: thanks to Jean-Michel Muller
! 242: */
! 243: check3(3.4657359027997265421, GMP_RNDN, 32.0);
! 244: check3(3.4657359027997265421, GMP_RNDU, 32.0);
! 245: check3(3.4657359027997265421, GMP_RNDD, 31.999999999999996447);
! 246: check3(2.26523754332090625496e+01, GMP_RNDD, 6.8833785261699581146e9);
! 247: check3(1.31478962104089092122e+01, GMP_RNDZ, 5.12930793917860137299e+05);
! 248: check3(4.25637507920002378103e-01, GMP_RNDU, 1.53056585656161181497e+00);
! 249: check3(6.26551618962329307459e-16, GMP_RNDU, 1.00000000000000066613e+00);
! 250: check3(-3.35589513871216568383e-03, GMP_RNDD, 9.96649729583626853291e-01);
! 251: check3(1.95151388850007272424e+01, GMP_RNDU, 2.98756340674767792225e+08);
! 252: check3(2.45045953503350730784e+01, GMP_RNDN, 4.38743344916128387451e+10);
! 253: check3(2.58165606081678085104e+01, GMP_RNDD, 1.62925781879432281494e+11);
! 254: check3(-2.36539020084338638128e+01, GMP_RNDZ, 5.33630792749924762447e-11);
! 255: check3(2.39211946135858077866e+01, GMP_RNDU, 2.44817704330214385986e+10);
! 256: check3(-2.78190533055889162029e+01, GMP_RNDZ, 8.2858803483596879512e-13);
! 257: check3(2.64028186174889789584e+01, GMP_RNDD, 2.9281844652878973388e11);
! 258: check3(2.92086338843268329413e+01, GMP_RNDZ, 4.8433797301907177734e12);
! 259: check3(-2.46355324071459982349e+01, GMP_RNDZ, 1.9995129297760994791e-11);
! 260: check3(-2.23509444608605427618e+01, GMP_RNDZ, 1.9638492867489702307e-10);
! 261: check3(-2.41175390197331687148e+01, GMP_RNDD, 3.3564940885530624592e-11);
! 262: check3(2.46363885231578088053e+01, GMP_RNDU, 5.0055014282693267822e10);
! 263: d = 7819821913254249.0 / 70368744177664.0;
! 264: check3(d, GMP_RNDN, 1.8262572323517295459e48);
! 265: check3(-3.56196340354684821250e+02, GMP_RNDN, 2.0225297096141478156e-155);
! 266: check3(6.59678273772710895173e+02, GMP_RNDU, 3.1234469273830195529e286);
! 267: check3(5.13772529701934331570e+02, GMP_RNDD, 1.3445427121297197752e223);
! 268: check3(3.57430211008718345056e+02, GMP_RNDD, 1.6981197246857298443e155);
! 269: check3(3.82001814471465536371e+02, GMP_RNDU, 7.9667300591087367805e165);
! 270: check3(5.92396038219384422518e+02, GMP_RNDD, 1.880747529554661989e257);
! 271: check3(-5.02678550462488090034e+02, GMP_RNDU, 4.8919201895446217839e-219);
! 272: check3(5.30015757134837031117e+02, GMP_RNDD, 1.5237672861171573939e230);
! 273: check3(5.16239362447650933063e+02, GMP_RNDZ, 1.5845518406744492105e224);
! 274: check3(6.00812634798592370977e-01, GMP_RNDN, 1.823600119339019443);
! 275: #ifdef MPFR_HAVE_FESETROUND
! 276: SEED_RAND (time(NULL));
! 277: N = (argc==1) ? 0 : atoi(argv[1]);
! 278: lo = (argc>=3) ? atof(argv[2]) : -7.083964185e2;
! 279: hi = (argc>=4) ? atof(argv[3]) : 7.097827129e2;
! 280: for (i=0;i<N;i++) {
! 281: /* select d such that exp(d) can be represented as a normalized
! 282: machine double-precision number,
! 283: i.e. 2^(-1022) <= exp(d) <= 2^(1023)*(2-2^(-52)) */
! 284: d = lo + (hi-lo)*DBL_RAND();
! 285: e = check(d, LONG_RAND() % 4);
! 286: s += e;
! 287: if (e>maxe) maxe=e;
! 288: }
! 289: if (N) printf("mean error=%1.2e max error=%d\n", (double)s/(double)N,maxe);
! 290: #endif
! 291:
! 292: return 0;
! 293: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to texp.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