[BACK]Return to tlog.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / gmp / mpfr / tests

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>