1/* Test mpq_get_d and mpq_set_d 2 3Copyright 1991, 1993, 1994, 1996, 2000, 2001, 2002, 2003 Free Software 4Foundation, Inc. 5 6This file is part of the GNU MP Library. 7 8The GNU MP Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MP Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along 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 42void dump __GMP_PROTO ((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 156void 157check_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 192void 193dump (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. */ 202void 203check_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 261int 262main (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} 273