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 } 273