Annotation of OpenXM_contrib/gmp/mpfr/tests/tagm.c, Revision 1.1
1.1 ! ohara 1: /* Test file for mpfr_agm.
! 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 "mpfr.h"
! 28: #include "mpfr-test.h"
! 29:
! 30: double drand_agm _PROTO((void));
! 31: double dagm _PROTO((double, double));
! 32: void check4 _PROTO((double, double, mp_rnd_t, double));
! 33: void check_large _PROTO((void));
! 34: void slave _PROTO((int, int));
! 35:
! 36: double
! 37: drand_agm (void)
! 38: {
! 39: double d; long int *i;
! 40:
! 41: i = (long int*) &d;
! 42: do {
! 43: i[0] = LONG_RAND();
! 44: i[1] = LONG_RAND();
! 45: } while ((d<1e-153)||(d>1e153)); /* to avoid underflow or overflow
! 46: in double calculus in sqrt(u*v) */
! 47:
! 48: return d;
! 49: }
! 50:
! 51: double
! 52: dagm (double a, double b)
! 53: {
! 54: double u, v, tmpu, tmpv;
! 55:
! 56: if ((isnan(a))||(isnan(b)))
! 57: return a+b;
! 58:
! 59: tmpv=MAX(a,b);
! 60: tmpu=MIN(a,b);
! 61:
! 62: do
! 63: {
! 64: u=tmpu;
! 65: v=tmpv;
! 66: tmpu=sqrt(u*v);
! 67: tmpv=(u+v)/2.0;
! 68: }
! 69: while (!(((tmpu==u)&&(tmpv==v))||(ulp(u,v)==0)));
! 70:
! 71: /* printf("difference : %i ulp\n",ulp(u,v)); */
! 72: return u;
! 73: }
! 74:
! 75: #define check(a,b,r) check4(a,b,r,0.0)
! 76:
! 77: void
! 78: check4 (double a, double b, mp_rnd_t rnd_mode, double res1)
! 79: {
! 80: mpfr_t ta, tb, tres;
! 81: double res2;
! 82: int ck=0;
! 83:
! 84: mpfr_init2(ta, 53);
! 85: mpfr_init2(tb, 53);
! 86: mpfr_init2(tres, 53);
! 87:
! 88: mpfr_set_d(ta, a, rnd_mode);
! 89: mpfr_set_d(tb, b, rnd_mode);
! 90:
! 91: mpfr_agm(tres, ta, tb, rnd_mode);
! 92: #ifdef MPFR_HAVE_FESETROUND
! 93: mpfr_set_machine_rnd_mode(rnd_mode);
! 94: #endif
! 95:
! 96: if (res1==0.0) res1=dagm(a,b); else ck=1;
! 97: if (ck==0) printf("%1.20e\n", res1);
! 98: res2 = mpfr_get_d1 (tres);
! 99:
! 100: if (ck && res1!=res2 && (!isnan(res1) || !isnan(res2))) {
! 101: printf("mpfr_agm failed for a=%1.20e, b=%1.20e, rnd_mode=%d\n",a,b,rnd_mode);
! 102: printf("expected result is %1.20e, got %1.20e (%d ulp)\n",res1,res2,
! 103: ulp(res2,res1));
! 104: /*exit(1);*/
! 105: }
! 106: mpfr_clear(ta); mpfr_clear(tb); mpfr_clear(tres);
! 107: }
! 108:
! 109: void
! 110: check_large (void)
! 111: {
! 112: mpfr_t a, b, agm;
! 113:
! 114: mpfr_init2(a, 82); mpfr_init2(b, 82); mpfr_init2(agm, 82);
! 115: mpfr_set_ui(a, 1, GMP_RNDN);
! 116: mpfr_set_str_raw(b, "0.1111101100001000000001011000110111101000001011111000100001000101010100011111110010E-39");
! 117: mpfr_agm(agm, a, b, GMP_RNDN);
! 118: mpfr_set_str_raw(a, "0.1110001000111101101010101010101101001010001001001011100101111011110101111001111100E-4");
! 119: if (mpfr_cmp(agm, a)) {
! 120: fprintf(stderr, "mpfr_agm failed for precision 82\n"); exit(1);
! 121: }
! 122: mpfr_clear(a); mpfr_clear(b); mpfr_clear(agm);
! 123: }
! 124:
! 125: void
! 126: slave (int N, int p)
! 127: {
! 128: int i;
! 129: double a,b;
! 130: mpfr_t ta, tb, tres;
! 131:
! 132: mpfr_init2(ta, 53);
! 133: mpfr_init2(tb, 53);
! 134: mpfr_init2(tres, p);
! 135: for(i=0;i<N;i++) {
! 136: a = drand_agm();
! 137: b = drand_agm();
! 138: mpfr_set_d(ta, a, GMP_RNDN);
! 139: mpfr_set_d(tb, b, GMP_RNDN);
! 140: mpfr_agm(tres, ta, tb, LONG_RAND() % 4 );
! 141: }
! 142: mpfr_clear(ta); mpfr_clear(ta); mpfr_clear(tres);
! 143: printf("fin\n");
! 144: }
! 145:
! 146:
! 147: int
! 148: main (int argc, char* argv[])
! 149: {
! 150: int N;
! 151:
! 152: SEED_RAND (time(NULL));
! 153:
! 154: if (argc==3) { /* tagm N p : N calculus with precision p*/
! 155: printf("Doing %d random tests in %d precision\n",atoi(argv[1]),atoi(argv[2]));
! 156: slave(atoi(argv[1]),atoi(argv[2]));
! 157: return 0;
! 158: }
! 159:
! 160: if (argc==2) { /* tagm N: N tests with random double's */
! 161: int i;
! 162: double a,b;
! 163:
! 164: N = atoi(argv[1]);
! 165: for (i=0;i<N;i++) {
! 166: a = drand();
! 167: b = drand();
! 168: check(a, b, LONG_RAND() % 4);
! 169: }
! 170: return 0;
! 171: }
! 172: else {
! 173: check_large();
! 174: check4(2.0, 1.0, GMP_RNDN, 1.45679103104690677029);
! 175: check4(6.0, 4.0, GMP_RNDN, 4.94936087247260925182);
! 176: check4(62.0, 61.0, GMP_RNDN, 6.14989837188450749750e+01);
! 177: check4(0.5, 1.0, GMP_RNDN, 7.28395515523453385143e-01);
! 178: check4(1.0, 2.0, GMP_RNDN, 1.45679103104690677029);
! 179: check4(234375765.0, 234375000.0, GMP_RNDN, 2.3437538249984395504e8);
! 180: check4(8.0, 1.0, GMP_RNDU, 3.615756177597362786);
! 181: check4(1.0, 44.0, GMP_RNDU, 1.33658354512981247808e1);
! 182: check4(1.0, 3.7252902984619140625e-9, GMP_RNDU, 7.55393356971199025907e-02);
! 183: }
! 184:
! 185: /* TODO : tests des infinis dans tagm.c */
! 186:
! 187: return 0;
! 188: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to tagm.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