t-get_d.c revision 1.1.1.2
1/* Test mpq_get_d and mpq_set_d 2 3Copyright 1991, 1993, 1994, 1996, 2000, 2001, 2002, 2003, 2012, 2013 Free 4Software Foundation, Inc. 5 6This file is part of the GNU MP Library test suite. 7 8The GNU MP Library test suite is free software; you can redistribute it 9and/or modify it under the terms of the GNU General Public License as 10published by the Free Software Foundation; either version 3 of the License, 11or (at your option) any later version. 12 13The GNU MP Library test suite is distributed in the hope that it will be 14useful, but WITHOUT ANY WARRANTY; without even the implied warranty of 15MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 16Public License for more details. 17 18You should have received a copy of the GNU General Public License along with 19the GNU MP Library test suite. 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) || defined (__vax__) && SIZE > 4 35#undef SIZE 36#define SIZE 4 37#define EPSIZE 3 38#else 39#define EPSIZE SIZE 40#endif 41 42void dump (mpq_t); 43 44void 45check_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 119double 120my_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#define MAXEXP 500 157 158#if defined (__vax) || defined (__vax__) 159#undef MAXEXP 160#define MAXEXP 30 161#endif 162 163void 164check_random (int argc, char **argv) 165{ 166 gmp_randstate_ptr rands = RANDS; 167 168 double d; 169 mpq_t q; 170 mpz_t a, t; 171 int exp; 172 173 int test, reps = 100000; 174 175 if (argc == 2) 176 reps = 100 * atoi (argv[1]); 177 178 mpq_init (q); 179 mpz_init (a); 180 mpz_init (t); 181 182 for (test = 0; test < reps; test++) 183 { 184 mpz_rrandomb (a, rands, 53); 185 mpz_urandomb (t, rands, 32); 186 exp = mpz_get_ui (t) % (2*MAXEXP) - MAXEXP; 187 188 d = my_ldexp (mpz_get_d (a), exp); 189 mpq_set_d (q, d); 190 /* Check that n/d = a * 2^exp, or 191 d*a 2^{exp} = n */ 192 mpz_mul (t, a, mpq_denref (q)); 193 if (exp > 0) 194 mpz_mul_2exp (t, t, exp); 195 else 196 { 197 if (!mpz_divisible_2exp_p (t, -exp)) 198 goto fail; 199 mpz_div_2exp (t, t, -exp); 200 } 201 if (mpz_cmp (t, mpq_numref (q)) != 0) 202 { 203 fail: 204 printf ("ERROR (check_random test %d): bad mpq_set_d results\n", test); 205 printf ("%.16g\n", d); 206 gmp_printf ("%Qd\n", q); 207 abort (); 208 } 209 } 210 mpq_clear (q); 211 mpz_clear (t); 212 mpz_clear (a); 213} 214 215void 216dump (mpq_t x) 217{ 218 mpz_out_str (stdout, 10, mpq_numref (x)); 219 printf ("/"); 220 mpz_out_str (stdout, 10, mpq_denref (x)); 221 printf ("\n"); 222} 223 224/* Check various values 2^n and 1/2^n. */ 225void 226check_onebit (void) 227{ 228 static const long data[] = { 229 -3*GMP_NUMB_BITS-1, -3*GMP_NUMB_BITS, -3*GMP_NUMB_BITS+1, 230 -2*GMP_NUMB_BITS-1, -2*GMP_NUMB_BITS, -2*GMP_NUMB_BITS+1, 231 -GMP_NUMB_BITS-1, -GMP_NUMB_BITS, -GMP_NUMB_BITS+1, 232 -5, -2, -1, 0, 1, 2, 5, 233 GMP_NUMB_BITS-1, GMP_NUMB_BITS, GMP_NUMB_BITS+1, 234 2*GMP_NUMB_BITS-1, 2*GMP_NUMB_BITS, 2*GMP_NUMB_BITS+1, 235 3*GMP_NUMB_BITS-1, 3*GMP_NUMB_BITS, 3*GMP_NUMB_BITS+1, 236 }; 237 238 int i, neg; 239 long exp, l; 240 mpq_t q; 241 double got, want; 242 243 mpq_init (q); 244 245 for (i = 0; i < numberof (data); i++) 246 { 247 exp = data[i]; 248 249 mpq_set_ui (q, 1L, 1L); 250 if (exp >= 0) 251 mpq_mul_2exp (q, q, exp); 252 else 253 mpq_div_2exp (q, q, -exp); 254 255 want = 1.0; 256 for (l = 0; l < exp; l++) 257 want *= 2.0; 258 for (l = 0; l > exp; l--) 259 want /= 2.0; 260 261 for (neg = 0; neg <= 1; neg++) 262 { 263 if (neg) 264 { 265 mpq_neg (q, q); 266 want = -want; 267 } 268 269 got = mpq_get_d (q); 270 271 if (got != want) 272 { 273 printf ("mpq_get_d wrong on %s2**%ld\n", neg ? "-" : "", exp); 274 mpq_trace (" q ", q); 275 d_trace (" want ", want); 276 d_trace (" got ", got); 277 abort(); 278 } 279 } 280 } 281 mpq_clear (q); 282} 283 284int 285main (int argc, char **argv) 286{ 287 tests_start (); 288 289 check_onebit (); 290 check_monotonic (argc, argv); 291 check_random (argc, argv); 292 293 tests_end (); 294 exit (0); 295} 296