Annotation of OpenXM_contrib/gmp/mpfr/tests/tlog.c, Revision 1.1
1.1 ! ohara 1: /* Test file for mpfr_log.
! 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 <math.h>
! 25: #include <time.h>
! 26: #include "gmp.h"
! 27: #include "mpfr.h"
! 28: #include "mpfr-test.h"
! 29:
! 30: double drand_log _PROTO((void));
! 31: int check1 _PROTO((double, mp_rnd_t, double, int, int));
! 32: void check3 _PROTO((double, unsigned long, mp_rnd_t));
! 33: void check4 _PROTO((int));
! 34: void slave _PROTO((int, int));
! 35: void check_worst_cases _PROTO((void));
! 36: void special _PROTO((void));
! 37:
! 38: #if (BITS_PER_LONGINT == 32)
! 39: #define INT32 long int
! 40: #else
! 41: #define INT32 int
! 42: #endif
! 43:
! 44:
! 45: double
! 46: drand_log (void)
! 47: {
! 48: double d; INT32 *i;
! 49:
! 50: i = (INT32*) &d;
! 51: do {
! 52: i[0] = LONG_RAND();
! 53: i[1] = LONG_RAND();
! 54: } while ((d<1e-153) || (d>1e153)); /* to avoid underflow or overflow
! 55: in double calculus in sqrt(u*v) */
! 56: return d;
! 57: }
! 58:
! 59: #define check2(a,rnd,res) check1(a,rnd,res,1,0)
! 60: #define check(a,r) check2(a,r,0.0)
! 61:
! 62: int
! 63: check1 (double a, mp_rnd_t rnd_mode, double res1, int ck, int max_ulp)
! 64: {
! 65: mpfr_t ta, tres;
! 66: double res2;
! 67: int diff=0;
! 68: /* ck=1 iff res1 is certified correct */
! 69:
! 70: #ifdef MPFR_HAVE_FESETROUND
! 71: mpfr_set_machine_rnd_mode(rnd_mode);
! 72: #endif
! 73: if (ck==0 && res1==0.0) res1=log(a);
! 74: mpfr_init2(ta, 53);
! 75: mpfr_init2(tres, 53);
! 76: mpfr_set_d(ta, a, GMP_RNDN);
! 77: mpfr_log(tres, ta, rnd_mode);
! 78: res2=mpfr_get_d1 (tres);
! 79: mpfr_clear(ta); mpfr_clear(tres);
! 80:
! 81: if (res1!=res2 && (!isnan(res1) || !isnan(res2))) {
! 82: diff = ulp(res1,res2);
! 83: if (ck) {
! 84: printf("mpfr_log failed for a=%1.20e, rnd_mode=%s\n", a,
! 85: mpfr_print_rnd_mode(rnd_mode));
! 86: printf("correct result is %1.20e\n mpfr_log gives %1.20e (%d ulp)\n",res1,res2,ulp(res1,res2));
! 87: exit(1);
! 88: }
! 89: else if (diff>max_ulp) {
! 90: printf("mpfr_log differs from libm.a for a=%1.20e, rnd_mode=%s\n", a,
! 91: mpfr_print_rnd_mode(rnd_mode));
! 92: printf(" double calculus gives %1.20e\n mpfr_log gives %1.20e (%d ulp)\n",res1,res2,ulp(res1,res2));
! 93: }
! 94: }
! 95: if (!isnan(res1) || !isnan(res2))
! 96: return diff;
! 97: else
! 98: return 0;
! 99: }
! 100:
! 101: void
! 102: check3 (double d, unsigned long prec, mp_rnd_t rnd)
! 103: {
! 104: mpfr_t x, y;
! 105:
! 106: mpfr_init2(x, prec); mpfr_init2(y, prec);
! 107: mpfr_set_d(x, d, rnd);
! 108: mpfr_log(y, x, rnd);
! 109: mpfr_out_str(stdout, 10, 0, y, rnd); putchar('\n');
! 110: mpfr_print_binary(y); putchar('\n');
! 111: mpfr_clear(x); mpfr_clear(y);
! 112: }
! 113:
! 114: void
! 115: check4 (int N)
! 116: {
! 117: int i, max=0, sum=0, cur;
! 118: double d;
! 119: mp_rnd_t rnd;
! 120:
! 121: for(i=0;i<N;i++) {
! 122: d = drand_log ();
! 123: rnd = LONG_RAND() % 4;
! 124: cur=check1 (d, rnd, 0.0, 0, max);
! 125: if (cur<0)
! 126: cur = -cur;
! 127: if (cur > max)
! 128: max=cur;
! 129: sum+=cur;
! 130: }
! 131: d=(double)sum / (double)N;
! 132: fprintf(stderr, "max error : %i \t mean error : %f (in ulps)\n",max,d);
! 133: }
! 134:
! 135: void
! 136: slave (int N, int p)
! 137: {
! 138: int i;
! 139: double d;
! 140: mpfr_t ta, tres;
! 141:
! 142: mpfr_init2(ta, 53);
! 143: mpfr_init2(tres, p);
! 144: for(i=0;i<N;i++) {
! 145: d = drand_log();
! 146: mpfr_set_d (ta, d, GMP_RNDN);
! 147: mpfr_log (tres, ta, LONG_RAND() % 4 );
! 148: }
! 149: mpfr_clear(ta); mpfr_clear(tres);
! 150: printf("fin\n");
! 151: }
! 152:
! 153: /* examples from Jean-Michel Muller and Vincent Lefevre
! 154: Cf http://www.ens-lyon.fr/~jmmuller/Intro-to-TMD.htm
! 155: */
! 156:
! 157: void
! 158: check_worst_cases (void)
! 159: {
! 160: check2(1.00089971802309629645, GMP_RNDD, 8.99313519443722736088e-04);
! 161: check2(1.00089971802309629645, GMP_RNDN, 8.99313519443722844508e-04);
! 162: check2(1.00089971802309629645, GMP_RNDU, 8.99313519443722844508e-04);
! 163:
! 164: check2(1.01979300812244555452, GMP_RNDD, 1.95996734891603630047e-02);
! 165: check2(1.01979300812244555452, GMP_RNDN, 1.95996734891603664741e-02);
! 166: check2(1.01979300812244555452, GMP_RNDU, 1.95996734891603664741e-02);
! 167:
! 168: check2(1.02900871924604464525, GMP_RNDD, 2.85959303301472726744e-02);
! 169: check2(1.02900871924604464525, GMP_RNDN, 2.85959303301472761438e-02);
! 170: check2(1.02900871924604464525, GMP_RNDU, 2.85959303301472761438e-02);
! 171:
! 172: check2(1.27832870030418943585, GMP_RNDD, 2.45553521871417795852e-01);
! 173: check2(1.27832870030418943585, GMP_RNDN, 2.45553521871417823608e-01);
! 174: check2(1.27832870030418943585, GMP_RNDU, 2.45553521871417823608e-01);
! 175:
! 176: check2(1.31706530746788241792, GMP_RNDD, 2.75406009586277422674e-01);
! 177: check2(1.31706530746788241792, GMP_RNDN, 2.75406009586277478185e-01);
! 178: check2(1.31706530746788241792, GMP_RNDU, 2.75406009586277478185e-01);
! 179:
! 180: check2(1.47116981099449883885, GMP_RNDD, 3.86057874110010412760e-01);
! 181: check2(1.47116981099449883885, GMP_RNDN, 3.86057874110010412760e-01);
! 182: check2(1.47116981099449883885, GMP_RNDU, 3.86057874110010468272e-01);
! 183:
! 184: check2(1.58405446812987782401, GMP_RNDD, 4.59987679246663727639e-01);
! 185: check2(1.58405446812987782401, GMP_RNDN, 4.59987679246663783150e-01);
! 186: check2(1.58405446812987782401, GMP_RNDU, 4.59987679246663783150e-01);
! 187:
! 188: check2(1.67192331263391547047, GMP_RNDD, 5.13974647961076613889e-01);
! 189: check2(1.67192331263391547047, GMP_RNDN, 5.13974647961076724911e-01);
! 190: check2(1.67192331263391547047, GMP_RNDU, 5.13974647961076724911e-01);
! 191:
! 192: check2(1.71101198068990645318, GMP_RNDD, 5.37084997042120315669e-01);
! 193: check2(1.71101198068990645318, GMP_RNDN, 5.37084997042120315669e-01);
! 194: check2(1.71101198068990645318, GMP_RNDU, 5.37084997042120426691e-01);
! 195:
! 196: check2(1.72634853551388700588, GMP_RNDD, 5.46008504786553605648e-01);
! 197: check2(1.72634853551388700588, GMP_RNDN, 5.46008504786553716670e-01);
! 198: check2(1.72634853551388700588, GMP_RNDU, 5.46008504786553716670e-01);
! 199:
! 200: check2(2.00028876593004323325, GMP_RNDD, 6.93291553102749702475e-01);
! 201: check2(2.00028876593004323325, GMP_RNDN, 6.93291553102749813497e-01);
! 202: check2(2.00028876593004323325, GMP_RNDU, 6.93291553102749813497e-01);
! 203:
! 204: check2(6.27593230200363105808, GMP_RNDD, 1.83672204800630312072);
! 205: check2(6.27593230200363105808, GMP_RNDN, 1.83672204800630334276);
! 206: check2(6.27593230200363105808, GMP_RNDU, 1.83672204800630334276);
! 207:
! 208: check2(7.47216682321367997588, GMP_RNDD, 2.01118502712453661729);
! 209: check2(7.47216682321367997588, GMP_RNDN, 2.01118502712453706138);
! 210: check2(7.47216682321367997588, GMP_RNDU, 2.01118502712453706138);
! 211:
! 212: check2(9.34589857718275318632, GMP_RNDD, 2.23493759221664944903);
! 213: check2(9.34589857718275318632, GMP_RNDN, 2.23493759221664989312);
! 214: check2(9.34589857718275318632, GMP_RNDU, 2.23493759221664989312);
! 215:
! 216: check2(10.6856587560831854944, GMP_RNDD, 2.36890253928838445674);
! 217: check2(10.6856587560831854944, GMP_RNDN, 2.36890253928838445674);
! 218: check2(10.6856587560831854944, GMP_RNDU, 2.36890253928838490083);
! 219:
! 220: check2(12.4646345033981766903, GMP_RNDD, 2.52289539471636015122);
! 221: check2(12.4646345033981766903, GMP_RNDN, 2.52289539471636015122);
! 222: check2(12.4646345033981766903, GMP_RNDU, 2.52289539471636059531);
! 223:
! 224: check2(17.0953275851761752335, GMP_RNDD, 2.83880518553861849185);
! 225: check2(17.0953275851761752335, GMP_RNDN, 2.83880518553861893594);
! 226: check2(17.0953275851761752335, GMP_RNDU, 2.83880518553861893594);
! 227:
! 228: check2(19.8509496207496916043, GMP_RNDD, 2.98825184582516722998);
! 229: check2(19.8509496207496916043, GMP_RNDN, 2.98825184582516722998);
! 230: check2(19.8509496207496916043, GMP_RNDU, 2.98825184582516767406);
! 231:
! 232: check2(23.9512076062771335216, GMP_RNDD, 3.17601874455977206679);
! 233: check2(23.9512076062771335216, GMP_RNDN, 3.17601874455977206679);
! 234: check2(23.9512076062771335216, GMP_RNDU, 3.17601874455977251088);
! 235:
! 236: check2(428.315247165198229595, GMP_RNDD, 6.05985948325268264369);
! 237: check2(428.315247165198229595, GMP_RNDN, 6.05985948325268353187);
! 238: check2(428.315247165198229595, GMP_RNDU, 6.05985948325268353187);
! 239: }
! 240:
! 241: void
! 242: special (void)
! 243: {
! 244: mpfr_t x, y;
! 245:
! 246: mpfr_init2 (x, 53);
! 247: mpfr_init2 (y, 53);
! 248: mpfr_set_ui (x, 3, GMP_RNDD);
! 249: mpfr_log (y, x, GMP_RNDD);
! 250: if (mpfr_get_d1 (y) != 1.09861228866810956) {
! 251: fprintf (stderr, "Error in mpfr_log(3) for GMP_RNDD\n");
! 252: exit (1);
! 253: }
! 254:
! 255: /* check large precision */
! 256: mpfr_set_prec (x, 3322);
! 257: mpfr_set_prec (y, 3322);
! 258: mpfr_set_d (x, 3.0, GMP_RNDN);
! 259: mpfr_sqrt (x, x, GMP_RNDN);
! 260: mpfr_log (y, x, GMP_RNDN);
! 261:
! 262: mpfr_clear (x);
! 263: mpfr_clear (y);
! 264: }
! 265:
! 266: #define TEST_FUNCTION mpfr_log
! 267: #include "tgeneric.c"
! 268:
! 269: int
! 270: main (int argc, char *argv[])
! 271: {
! 272: int N = 0;
! 273: double d;
! 274:
! 275: SEED_RAND (time(NULL));
! 276: if (argc==4) { /* tlog x prec rnd */
! 277: check3(atof(argv[1]), atoi(argv[2]), atoi(argv[3]));
! 278: return 0;
! 279: }
! 280:
! 281: if (argc==3) { /* tlog N p : N calculus with precision p*/
! 282: printf("Doing %d random tests in %d precision\n",atoi(argv[1]),atoi(argv[2]));
! 283: slave(atoi(argv[1]),atoi(argv[2]));
! 284: return 0;
! 285: }
! 286:
! 287: if (argc==2) { /* tlog N: N tests with random double's */
! 288: N=atoi(argv[1]);
! 289: printf("Doing %d random tests in double precision\n", N);
! 290: check4(N);
! 291: }
! 292: else {
! 293: special ();
! 294: check_worst_cases();
! 295:
! 296: check2(1.01979300812244555452, GMP_RNDN, 1.95996734891603664741e-02);
! 297: check2(10.0,GMP_RNDU,2.30258509299404590110e+00);
! 298: check2(6.0,GMP_RNDU,1.79175946922805517936);
! 299: check2(1.0,GMP_RNDZ,0.0);
! 300: check2(62.0,GMP_RNDU,4.12713438504509166905);
! 301: check2(0.5,GMP_RNDZ,-6.93147180559945286226e-01);
! 302: check2(3.0,GMP_RNDZ,1.09861228866810956006e+00);
! 303: check2(234375765.0,GMP_RNDU,1.92724362186836231104e+01);
! 304: check2(8.0,GMP_RNDZ,2.07944154167983574765e+00);
! 305: check2(44.0,GMP_RNDU,3.78418963391826146392e+00);
! 306: check2(1.01979300812244555452, GMP_RNDN, 1.95996734891603664741e-02);
! 307: /* bugs found by Vincent Lefe`vre */
! 308: d = -4723773766428415.0 / 1180591620717411303424.0;
! 309: check2(0.99999599881598921769, GMP_RNDN, d);
! 310: check2(9.99995576063808955247e-01, GMP_RNDZ, -4.42394597667932383816e-06);
! 311: check2(9.99993687357856209097e-01, GMP_RNDN, -6.31266206860017342601e-06);
! 312: check2(9.99995223520736886691e-01, GMP_RNDN, -4.77649067052670982220e-06);
! 313: check2(9.99993025794720935551e-01, GMP_RNDN, -6.97422959894716163837e-06);
! 314: check2(9.99987549017837484833e-01, GMP_RNDN, -1.24510596766369924330e-05);
! 315: check2(9.99985901426543311032e-01, GMP_RNDN, -1.40986728425098585229e-05);
! 316: d = -8232353813100321.0 / 590295810358705651712.0;
! 317: check2(9.99986053947420794330e-01, GMP_RNDN, d);
! 318: check2(9.99971938247442126979e-01, GMP_RNDN, -2.80621462962173414790e-05);
! 319: /* other bugs found by Vincent Lefe`vre */
! 320: check2(1.18615436389927785905e+77, GMP_RNDN, 1.77469768607706015473e+02);
! 321: check2(9.48868723578399476187e+77, GMP_RNDZ, 1.79549152432275803903e+02);
! 322: check2(2.31822210096938820854e+89, GMP_RNDN, 2.05770873832573869322e+02);
! 323: /* further bugs found by Vincent Lefe`vre */
! 324: check2(9.99999989485669482647e-01, GMP_RNDZ, -1.05143305726283042331e-08);
! 325: check2(9.99999989237970177136e-01, GMP_RNDZ, -1.07620298807745377934e-08);
! 326: check2(9.99999989239339082125e-01, GMP_RNDN, -1.07606609757704445430e-08);
! 327:
! 328: check2(7.3890560989306504,GMP_RNDU,2.0000000000000004); /* exp(2.0) */
! 329: check2(7.3890560989306495,GMP_RNDU,2.0); /* exp(2.0) */
! 330: check2(7.53428236571286402512e+34,GMP_RNDZ,8.03073567492226345621e+01);
! 331: check2(6.18784121531737948160e+19,GMP_RNDZ,4.55717030391710693493e+01);
! 332: check2(1.02560267603047283735e+00,GMP_RNDD,2.52804164149448735987e-02);
! 333: check2(7.53428236571286402512e+34,GMP_RNDZ,8.03073567492226345621e+01);
! 334: d = 497773706319601.0 / 4398046511104.0;
! 335: check2(1.42470900831881198052e+49, GMP_RNDZ, d);
! 336:
! 337: check2(1.08013816255293777466e+11,GMP_RNDN,2.54055249841782604392e+01);
! 338: check2(6.72783635300509015581e-37,GMP_RNDU,-8.32893948416799503320e+01);
! 339: check2(2.25904918906057891180e-52,GMP_RNDU,-1.18919480823735682406e+02);
! 340: check2(1.48901209246462951085e+00,GMP_RNDD,3.98112874867437460668e-01);
! 341: check2(1.70322470467612341327e-01,GMP_RNDN,-1.77006175364294615626);
! 342: check2(1.94572026316065240791e+01,GMP_RNDD,2.96821731676437838842);
! 343: check2(4.01419512207026418764e+04,GMP_RNDD,1.06001772315501128218e+01);
! 344: check2(9.47077365236487591672e-04,GMP_RNDZ,-6.96212977303956748187e+00);
! 345: check2(3.95906157687589643802e-109,GMP_RNDD,-2.49605768114704119399e+02);
! 346: check2(2.73874914516503004113e-02,GMP_RNDD,-3.59766888618655977794e+00);
! 347: check2(9.18989072589566467669e-17,GMP_RNDZ,-3.69258425351464083519e+01);
! 348: check2(dbl(2830750724514701.0,131),GMP_RNDZ,dbl(1111664301085491.0,-43));
! 349: check2(1.74827399630587801934e-23,GMP_RNDZ,-5.24008281254547156891e+01);
! 350: check2(4.35302958401482307665e+22,GMP_RNDD,5.21277441046519527390e+01);
! 351: check2(9.70791868689332915209e+00,GMP_RNDD,2.27294191194272210410e+00);
! 352: check2(2.22183639799464011100e-01,GMP_RNDN,-1.50425103275253957413e+00);
! 353: check2(2.27313466156682375540e+00,GMP_RNDD,8.21159787095675608448e-01);
! 354: check2(6.58057413965851156767e-01,GMP_RNDZ,-4.18463096196088235600e-01);
! 355: d = 7107588635148285.0 / 70368744177664.0;
! 356: check2 (7.34302197248998461006e+43, GMP_RNDZ, d);
! 357: check2(6.09969788341579732815e+00,GMP_RNDD,1.80823924264386204363e+00);
! 358: }
! 359:
! 360: test_generic (2, 100, 40);
! 361:
! 362: return 0;
! 363: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>