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