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

Annotation of OpenXM_contrib/gmp/mpfr/tests/tget_d.c, Revision 1.1

1.1     ! ohara       1: /* Test file for mpfr_get_d
        !             2:
        !             3: Copyright 1999, 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 <stdio.h>
        !            23: #include <stdlib.h>
        !            24: #include <float.h>
        !            25: #include "gmp.h"
        !            26: #include "gmp-impl.h"
        !            27: #include "mpfr.h"
        !            28: #include "mpfr-impl.h"
        !            29: #include "mpfr-test.h"
        !            30:
        !            31: int check_denorms _PROTO ((void));
        !            32:
        !            33: int
        !            34: check_denorms ()
        !            35: {
        !            36:   mp_rnd_t rnd_mode;
        !            37:   mpfr_t x;
        !            38:   double d, d2, dd, f;
        !            39:   int fail = 0, k, n;
        !            40:
        !            41:   mpfr_init2 (x, BITS_PER_MP_LIMB);
        !            42:
        !            43:   for (rnd_mode = 0; rnd_mode <= 3; rnd_mode++)
        !            44:     {
        !            45: #ifdef MPFR_HAVE_FESETROUND
        !            46:       mpfr_set_machine_rnd_mode (rnd_mode);
        !            47: #else
        !            48:       if (rnd_mode)
        !            49:         break; /* second iteration */
        !            50:       rnd_mode = GMP_RNDN;
        !            51: #endif
        !            52:       for (k = -17; k <= 17; k += 2)
        !            53:         {
        !            54:           d = k * DBL_MIN; /* k * 2^(-1022) */
        !            55:           f = 1.0;
        !            56:           mpfr_set_si (x, k, GMP_RNDN);
        !            57:           mpfr_div_2exp (x, x, 1022, GMP_RNDN); /* k * 2^(-1022) */
        !            58:           for (n = 0; n <= 58; n++)
        !            59:             {
        !            60:               d2 = d * f;
        !            61:               dd = mpfr_get_d (x, rnd_mode);
        !            62:               if (d2 != dd) /* should be k * 2^(-1022-n) for n < 53 */
        !            63:                 {
        !            64:                   fprintf (stderr,
        !            65:                            "Wrong result for %d * 2^(%d), rnd_mode %d\n",
        !            66:                            k, -1022-n, rnd_mode);
        !            67:                   fprintf (stderr, "got %.20e instead of %.20e\n", dd, d2);
        !            68:                   fail = 1;
        !            69:                 }
        !            70:               f *= 0.5;
        !            71:               mpfr_div_2exp (x, x, 1, GMP_RNDN);
        !            72:             }
        !            73:         }
        !            74:     }
        !            75:
        !            76:   mpfr_clear (x);
        !            77:   return fail;
        !            78: }
        !            79:
        !            80: int
        !            81: main (void)
        !            82: {
        !            83:
        !            84: #ifdef MPFR_HAVE_FESETROUND
        !            85:
        !            86:   mpfr_t half, x, y;
        !            87:   mp_rnd_t rnd_mode;
        !            88:
        !            89: #endif
        !            90:
        !            91:   mpfr_test_init ();
        !            92:
        !            93: #ifdef MPFR_HAVE_FESETROUND
        !            94:
        !            95:   mpfr_init2(half, 2);
        !            96:   mpfr_set_ui(half, 1, GMP_RNDZ);
        !            97:   mpfr_div_2ui(half, half, 1, GMP_RNDZ); /* has exponent 0 */
        !            98:
        !            99:   mpfr_init2(x, 128);
        !           100:   mpfr_init2(y, 128);
        !           101:
        !           102:   for (rnd_mode = 0; rnd_mode <= 3; rnd_mode++)
        !           103:     {
        !           104:       int i, j, si, sj;
        !           105:       double di, dj;
        !           106:
        !           107:       mpfr_set_machine_rnd_mode (rnd_mode);
        !           108:       for (i = 1, di = 0.25; i < 127; i++, di *= 0.5)
        !           109:         for (si = 0; si <= 1; si++)
        !           110:           {
        !           111:             mpfr_div_2ui (x, half, i, GMP_RNDZ);
        !           112:             (si ? mpfr_sub : mpfr_add)(x, half, x, GMP_RNDZ);
        !           113:             /* x = 1/2 +/- 1/2^(1+i) */
        !           114:             for (j = i+1, dj = di * 0.5; j < 128 && j < i+53; j++, dj *= 0.5)
        !           115:               for (sj = 0; sj <= 1; sj++)
        !           116:                 {
        !           117:                   double c, d, dd;
        !           118:                   int exp;
        !           119:                   char *f;
        !           120:
        !           121:                   mpfr_div_2ui (y, half, j, GMP_RNDZ);
        !           122:                   (sj ? mpfr_sub : mpfr_add)(y, x, y, GMP_RNDZ);
        !           123:                   /* y = 1/2 +/- 1/2^(1+i) +/- 1/2^(1+j) */
        !           124:                   exp = (LONG_RAND() % 47) - 23;
        !           125:                   mpfr_mul_2si (y, y, exp, GMP_RNDZ);
        !           126:                   if (mpfr_inexflag_p())
        !           127:                     {
        !           128:                       fprintf(stderr, "Error in tget_d: inexact flag for "
        !           129:                               "(i,si,j,sj,rnd,exp) = (%d,%d,%d,%d,%d,%d)\n",
        !           130:                               i, si, j, sj, rnd_mode, exp);
        !           131:                       exit(1);
        !           132:                     }
        !           133:                   dd = si != sj ? di - dj : di + dj;
        !           134:                   d = si ? 0.5 - dd : 0.5 + dd;
        !           135:                   if ((LONG_RAND() / 1024) & 1)
        !           136:                     {
        !           137:                       c = mpfr_get_d (y, rnd_mode);
        !           138:                       f = "mpfr_get_d";
        !           139:                     }
        !           140:                   else
        !           141:                     {
        !           142:                       exp = (LONG_RAND() % 47) - 23;
        !           143:                       c = mpfr_get_d3 (y, exp, rnd_mode);
        !           144:                       f = "mpfr_get_d3";
        !           145:                       if (si) /* then real d < 0.5 */
        !           146:                         d *= sj && i == 1 ? 4 : 2; /* normalize real d */
        !           147:                     }
        !           148:                   if (exp > 0)
        !           149:                     d *= 1 << exp;
        !           150:                   if (exp < 0)
        !           151:                     d /= 1 << -exp;
        !           152:                   if (c != d)
        !           153:                     {
        !           154:                       fprintf (stderr, "Error in tget_d (%s) for "
        !           155:                                "(i,si,j,sj,rnd,exp) = (%d,%d,%d,%d,%d,%d)\n"
        !           156:                                "got %.25Le instead of %.25Le\n"
        !           157:                                "Difference: %.19e\n",
        !           158:                                f, i, si, j, sj, rnd_mode, exp,
        !           159:                                (long double) c, (long double) d, d - c);
        !           160:                       exit (1);
        !           161:                     }
        !           162:                 }
        !           163:           }
        !           164:     }
        !           165:
        !           166:   mpfr_clear(half);
        !           167:   mpfr_clear(x);
        !           168:   mpfr_clear(y);
        !           169:
        !           170: #endif
        !           171:
        !           172:   if (check_denorms ())
        !           173:     exit (1);
        !           174:
        !           175:   return 0;
        !           176: }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>