tests/mpq/t-get_d.c

   1 /* Test mpq_get_d and mpq_set_d
   2
   3 Copyright 1991, 1993, 1994, 1996, 2000, 2001, 2002, 2003 Free Software
   4 Foundation, Inc.
   5
   6 This file is part of the GNU MP Library.
   7
   8 The GNU MP Library is free software; you can redistribute it and/or modify
   9 it under the terms of the GNU Lesser General Public License as published by
  10 the Free Software Foundation; either version 3 of the License, or (at your
  11 option) any later version.
  12
  13 The GNU MP Library is distributed in the hope that it will be useful, but
  14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  15 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  16 License for more details.
  17
  18 You should have received a copy of the GNU Lesser General Public License
  19 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  20
  21 #include <stdio.h>
  22 #include <stdlib.h>
  23
  24 #include "gmp.h"
  25 #include "gmp-impl.h"
  26 #include "tests.h"
  27
  28 #ifndef SIZE
  29 #define SIZE 8
  30 #endif
  31
  32 /* VAX D floats only have an 8 bit signed exponent, so anything 2^128 or
  33    bigger will overflow, that being 4 limbs. */
  34 #if defined (__vax__) && SIZE > 4
  35 #undef SIZE
  36 #define SIZE 4
  37 #define EPSIZE 3
  38 #else
  39 #define EPSIZE SIZE
  40 #endif
  41
  42 void dump __GMP_PROTO ((mpq_t));
  43
  44 void
  45 check_monotonic (int argc, char **argv)
  46 {
  47   mpq_t a;
  48   mp_size_t size;
  49   int reps = 100;
  50   int i, j;
  51   double last_d, new_d;
  52   mpq_t qlast_d, qnew_d;
  53   mpq_t eps;
  54
  55   if (argc == 2)
  56      reps = atoi (argv[1]);
  57
  58   /* The idea here is to test the monotonousness of mpq_get_d by adding
  59      numbers to the numerator and denominator.  */
  60
  61   mpq_init (a);
  62   mpq_init (eps);
  63   mpq_init (qlast_d);
  64   mpq_init (qnew_d);
  65
  66   for (i = 0; i < reps; i++)
  67     {
  68       size = urandom () % SIZE - SIZE/2;
  69       mpz_random2 (mpq_numref (a), size);
  70       do
  71         {
  72           size = urandom () % SIZE - SIZE/2;
  73           mpz_random2 (mpq_denref (a), size);
  74         }
  75       while (mpz_cmp_ui (mpq_denref (a), 0) == 0);
  76
  77       mpq_canonicalize (a);
  78
  79       last_d = mpq_get_d (a);
  80       mpq_set_d (qlast_d, last_d);
  81       for (j = 0; j < 10; j++)
  82         {
  83           size = urandom () % EPSIZE + 1;
  84           mpz_random2 (mpq_numref (eps), size);
  85           size = urandom () % EPSIZE + 1;
  86           mpz_random2 (mpq_denref (eps), size);
  87           mpq_canonicalize (eps);
  88
  89           mpq_add (a, a, eps);
  90           mpq_canonicalize (a);
  91           new_d = mpq_get_d (a);
  92           if (last_d > new_d)
  93             {
  94               printf ("\nERROR (test %d/%d): bad mpq_get_d results\n", i, j);
  95               printf ("last: %.16g\n", last_d);
  96               printf (" new: %.16g\n", new_d); dump (a);
  97               abort ();
  98             }
  99           mpq_set_d (qnew_d, new_d);
 100           MPQ_CHECK_FORMAT (qnew_d);
 101           if (mpq_cmp (qlast_d, qnew_d) > 0)
 102             {
 103               printf ("ERROR (test %d/%d): bad mpq_set_d results\n", i, j);
 104               printf ("last: %.16g\n", last_d); dump (qlast_d);
 105               printf (" new: %.16g\n", new_d); dump (qnew_d);
 106               abort ();
 107             }
 108           last_d = new_d;
 109           mpq_set (qlast_d, qnew_d);
 110         }
 111     }
 112
 113   mpq_clear (a);
 114   mpq_clear (eps);
 115   mpq_clear (qlast_d);
 116   mpq_clear (qnew_d);
 117 }
 118
 119 double
 120 my_ldexp (double d, int e)
 121 {
 122   for (;;)
 123     {
 124       if (e > 0)
 125         {
 126           if (e >= 16)
 127             {
 128               d *= 65536.0;
 129               e -= 16;
 130             }
 131           else
 132             {
 133               d *= 2.0;
 134               e -= 1;
 135             }
 136         }
 137       else if (e < 0)
 138         {
 139
 140           if (e <= -16)
 141             {
 142               d /= 65536.0;
 143               e += 16;
 144             }
 145           else
 146             {
 147               d /= 2.0;
 148               e += 1;
 149             }
 150         }
 151       else
 152         return d;
 153     }
 154 }
 155
 156 void
 157 check_random (int argc, char **argv)
 158 {
 159   double d, d2, nd, dd;
 160   mpq_t q;
 161   mp_limb_t rp[LIMBS_PER_DOUBLE + 1];
 162   int test, reps = 100000;
 163   int i;
 164
 165   if (argc == 2)
 166      reps = 100 * atoi (argv[1]);
 167
 168   mpq_init (q);
 169
 170   for (test = 0; test < reps; test++)
 171     {
 172       mpn_random2 (rp, LIMBS_PER_DOUBLE + 1);
 173       d = 0.0;
 174       for (i = LIMBS_PER_DOUBLE - 1; i >= 0; i--)
 175         d = d * MP_BASE_AS_DOUBLE + rp[i];
 176       d = my_ldexp (d, (int) (rp[LIMBS_PER_DOUBLE] % 1000) - 500);
 177       mpq_set_d (q, d);
 178       nd = mpz_get_d (mpq_numref (q));
 179       dd = mpz_get_d (mpq_denref (q));
 180       d2 = nd / dd;
 181       if (d != d2)
 182         {
 183           printf ("ERROR (check_random test %d): bad mpq_set_d results\n", test);
 184           printf ("%.16g\n", d);
 185           printf ("%.16g\n", d2);
 186           abort ();
 187         }
 188     }
 189   mpq_clear (q);
 190 }
 191
 192 void
 193 dump (mpq_t x)
 194 {
 195   mpz_out_str (stdout, 10, mpq_numref (x));
 196   printf ("/");
 197   mpz_out_str (stdout, 10, mpq_denref (x));
 198   printf ("\n");
 199 }
 200
 201 /* Check various values 2^n and 1/2^n. */
 202 void
 203 check_onebit (void)
 204 {
 205   static const long data[] = {
 206     -3*GMP_NUMB_BITS-1, -3*GMP_NUMB_BITS, -3*GMP_NUMB_BITS+1,
 207     -2*GMP_NUMB_BITS-1, -2*GMP_NUMB_BITS, -2*GMP_NUMB_BITS+1,
 208     -GMP_NUMB_BITS-1, -GMP_NUMB_BITS, -GMP_NUMB_BITS+1,
 209     -5, -2, -1, 0, 1, 2, 5,
 210     GMP_NUMB_BITS-1, GMP_NUMB_BITS, GMP_NUMB_BITS+1,
 211     2*GMP_NUMB_BITS-1, 2*GMP_NUMB_BITS, 2*GMP_NUMB_BITS+1,
 212     3*GMP_NUMB_BITS-1, 3*GMP_NUMB_BITS, 3*GMP_NUMB_BITS+1,
 213   };
 214
 215   int     i, neg;
 216   long    exp, l;
 217   mpq_t   q;
 218   double  got, want;
 219
 220   mpq_init (q);
 221
 222   for (i = 0; i < numberof (data); i++)
 223     {
 224       exp = data[i];
 225
 226       mpq_set_ui (q, 1L, 1L);
 227       if (exp >= 0)
 228         mpq_mul_2exp (q, q, exp);
 229       else
 230         mpq_div_2exp (q, q, -exp);
 231
 232       want = 1.0;
 233       for (l = 0; l < exp; l++)
 234         want *= 2.0;
 235       for (l = 0; l > exp; l--)
 236         want /= 2.0;
 237
 238       for (neg = 0; neg <= 1; neg++)
 239         {
 240           if (neg)
 241             {
 242               mpq_neg (q, q);
 243               want = -want;
 244             }
 245
 246           got = mpq_get_d (q);
 247
 248           if (got != want)
 249             {
 250               printf    ("mpq_get_d wrong on %s2**%ld\n", neg ? "-" : "", exp);
 251               mpq_trace ("   q    ", q);
 252               d_trace   ("   want ", want);
 253               d_trace   ("   got  ", got);
 254               abort();
 255             }
 256         }
 257     }
 258   mpq_clear (q);
 259 }
 260
 261 int
 262 main (int argc, char **argv)
 263 {
 264   tests_start ();
 265
 266   check_onebit ();
 267   check_monotonic (argc, argv);
 268   check_random (argc, argv);
 269
 270   tests_end ();
 271   exit (0);
 272 }