/* Test file for mpfr_log. Copyright (C) 1999 PolKA project, Inria Lorraine and Loria This file is part of the MPFR Library. The MPFR Library is free software; you can redistribute it and/or modify it under the terms of the GNU Library General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. The MPFR Library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public License for more details. You should have received a copy of the GNU Library General Public License along with the MPFR Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ #include #include #include #include #include "gmp.h" #include "mpfr.h" #include "mpfr-impl.h" extern int isnan(); double drand_log() { double d; long int *i; i = (long int*) &d; do { i[0] = lrand48(); i[1] = lrand48(); } while ((d<1e-153)||(d>1e153)); /* to avoid underflow or overflow in double calculus in sqrt(u*v) */ return d; } #define check2(a,rnd,res) check1(a,rnd,res,1) #define check(a,r) check2(a,r,0.0) int check1(double a, unsigned char rnd_mode, double res1, int ck) { mpfr_t ta, tres; double res2; /* ck=1 iff res1 is certified correct */ #ifdef TEST mpfr_set_machine_rnd_mode(rnd_mode); #endif if (ck==0 && res1==0.0) res1=log(a); mpfr_init2(ta, 53); mpfr_init2(tres, 53); mpfr_set_d(ta, a, GMP_RNDN); mpfr_log(tres, ta, rnd_mode); res2=mpfr_get_d(tres); mpfr_clear(ta); mpfr_clear(tres); if (res1!=res2 && (!isnan(res1) || !isnan(res2))) { if (ck) { printf("mpfr_log failed for a=%1.20e, rnd_mode=%s\n", a, mpfr_print_rnd_mode(rnd_mode)); printf("correct result is %1.20e\n mpfr_log gives %1.20e (%d ulp)\n",res1,res2,ulp(res1,res2)); exit(1); } else { printf("mpfr_log differs from libm.a for a=%1.20e, rnd_mode=%s\n", a, mpfr_print_rnd_mode(rnd_mode)); printf(" double calculus gives %1.20e\n mpfr_log gives %1.20e (%d ulp)\n",res1,res2,ulp(res1,res2)); } } if (!isnan(res1) || !isnan(res2)) return ulp(res1,res2); else return 0; } void check3(double d, unsigned long prec, unsigned char rnd) { mpfr_t x, y; mpfr_init2(x, prec); mpfr_init2(y, prec); mpfr_set_d(x, d, rnd); mpfr_log(y, x, rnd); mpfr_out_str(stdout, 10, 0, y, rnd); putchar('\n'); mpfr_print_raw(y); putchar('\n'); mpfr_clear(x); mpfr_clear(y); } void check4(int N) { int i, max=-1, sum=0, cur; double d; for(i=0;i max) max=cur; sum+=cur; } d=(double)sum / (double)N; printf("max error : %i \t mean error : %f (in ulps)\n",max,d); } void slave(int N, int p) { int i; double d; mpfr_t ta, tres; mpfr_init2(ta, 53); mpfr_init2(tres, p); for(i=0;i