Return to tsub_ui.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpfr / tests

Annotation of OpenXM_contrib/gmp/mpfr/tests/tsub_ui.c, Revision 1.1.1.1

1.1       ohara       1: /* Test file for mpfr_sub_ui
                      2:
                      3: Copyright 2000, 2001, 2002 Free Software Foundation.
                      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-impl.h"
                     29: #include "mpfr-test.h"
                     30:
                     31: void check_two_sum _PROTO ((mp_prec_t));
                     32: void check3 _PROTO ((double, unsigned long, mp_rnd_t, double));
                     33:
                     34: #define check(x,y,r) check3(x,y,r,0.0)
                     35:
                     36: /* checks that x+y gives the same results in double
                     37:    and with mpfr with 53 bits of precision */
                     38: void
                     39: check3 (double x, unsigned long y, mp_rnd_t rnd_mode, double z1)
                     40: {
                     41:   double z2;
                     42:   mpfr_t xx,zz;
                     43:
                     44:   mpfr_init(xx);
                     45:   mpfr_init(zz);
                     46:   mpfr_set_prec(xx, 53);
                     47:   mpfr_set_prec(zz, 53);
                     48:   mpfr_set_d(xx, x, rnd_mode);
                     49:   mpfr_sub_ui(zz, xx, y, rnd_mode);
                     50: #ifdef MPFR_HAVE_FESETROUND
                     51:   mpfr_set_machine_rnd_mode(rnd_mode);
                     52: #endif
                     53:   if (z1==0.0) z1 = x-y;
                     54:   z2 = mpfr_get_d1 (zz);
                     55:   if (z1!=z2 && !(isnan(z1) && isnan(z2))) {
                     56:     printf("expected sum is %1.20e, got %1.20e\n",z1,z2);
                     57:     printf("mpfr_sub_ui failed for x=%1.20e y=%lu with rnd_mode=%s\n",
                     58:           x, y, mpfr_print_rnd_mode(rnd_mode));
                     59:     exit(1);
                     60:   }
                     61:   mpfr_clear(xx);
                     62:   mpfr_clear(zz);
                     63: }
                     64:
                     65: /* FastTwoSum: if EXP(x) >= EXP(y), u = o(x+y), v = o(u-x), w = o(y-v),
                     66:                then x + y = u + w
                     67: thus if u = o(y-x), v = o(u+x), w = o(v-y), then y-x = u-w */
                     68: void
                     69: check_two_sum (mp_prec_t p)
                     70: {
                     71:   unsigned int x;
                     72:   mpfr_t y, u, v, w;
                     73:   mp_rnd_t rnd;
                     74:   int inexact;
                     75:
                     76:   mpfr_init2 (y, p);
                     77:   mpfr_init2 (u, p);
                     78:   mpfr_init2 (v, p);
                     79:   mpfr_init2 (w, p);
                     80:   do
                     81:     {
                     82:       x = LONG_RAND ();
                     83:     }
                     84:   while (x < 1);
                     85:   mpfr_random (y);
                     86:   rnd = LONG_RAND() % 4;
                     87:   rnd = GMP_RNDN;
                     88:   inexact = mpfr_sub_ui (u, y, x, GMP_RNDN);
                     89:   mpfr_add_ui (v, u, x, GMP_RNDN);
                     90:   mpfr_sub (w, v, y, GMP_RNDN);
                     91:   /* as u - (y-x) = w, we should have inexact and w of same sign */
                     92:   if (((inexact == 0) && mpfr_cmp_ui (w, 0)) ||
                     93:       ((inexact > 0) && (mpfr_cmp_ui (w, 0) <= 0)) ||
                     94:       ((inexact < 0) && (mpfr_cmp_ui (w, 0) >= 0)))
                     95:     {
                     96:       fprintf (stderr, "Wrong inexact flag for prec=%u, rnd=%s\n", (unsigned)p,
                     97:               mpfr_print_rnd_mode (rnd));
                     98:       printf ("x=%u\n", x);
                     99:       printf ("y="); mpfr_print_binary(y); putchar('\n');
                    100:       printf ("u="); mpfr_print_binary(u); putchar('\n');
                    101:       printf ("v="); mpfr_print_binary(v); putchar('\n');
                    102:       printf ("w="); mpfr_print_binary(w); putchar('\n');
                    103:       printf ("inexact = %d\n", inexact);
                    104:       exit (1);
                    105:     }
                    106:   mpfr_clear (y);
                    107:   mpfr_clear (u);
                    108:   mpfr_clear (v);
                    109:   mpfr_clear (w);
                    110: }
                    111:
                    112: int
                    113: main (int argc, char *argv[])
                    114: {
                    115:   mp_prec_t p;
                    116:   int k;
                    117: #ifdef MPFR_HAVE_FESETROUND
                    118:   double x; unsigned long y, N; int i,rnd_mode,rnd;
                    119:
                    120:   mpfr_test_init ();
                    121:
                    122:   SEED_RAND (time(NULL));
                    123:   N = (argc<2) ? 1000000 : atoi(argv[1]);
                    124:   rnd_mode = (argc<3) ? -1 : atoi(argv[2]);
                    125:   for (i=0;i<1000000;i++) {
                    126:     x = drand();
                    127:     y = LONG_RAND();
                    128:     if (ABS(x)>2.2e-307 && x+y<1.7e+308 && x+y>-1.7e308) {
                    129:       /* avoid denormalized numbers and overflows */
                    130:       rnd = (rnd_mode==-1) ? LONG_RAND()%4 : rnd_mode;
                    131:       check(x, y, rnd);
                    132:     }
                    133:   }
                    134: #endif
                    135:
                    136:   for (p=2; p<200; p++)
                    137:     for (k=0; k<200; k++)
                    138:       check_two_sum (p);
                    139:
                    140:   check3 (0.9999999999, 1, GMP_RNDN, -56295.0 / 562949953421312.0);
                    141: #ifdef HAVE_INFS
                    142:   check3 (DBL_NAN, 1, GMP_RNDN, DBL_NAN);
                    143:   check3 (DBL_POS_INF, 1, GMP_RNDN, DBL_POS_INF);
                    144:   check3 (DBL_NEG_INF, 1, GMP_RNDN, DBL_NEG_INF);
                    145: #endif
                    146:
                    147:   return 0;
                    148: }
                    149:
                    150: