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