1/*
2
3Copyright 2012, 2013, 2018, 2020 Free Software Foundation, Inc.
4
5This file is part of the GNU MP Library test suite.
6
7The GNU MP Library test suite is free software; you can redistribute it
8and/or modify it under the terms of the GNU General Public License as
9published by the Free Software Foundation; either version 3 of the License,
10or (at your option) any later version.
11
12The GNU MP Library test suite is distributed in the hope that it will be
13useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
14MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
15Public License for more details.
16
17You should have received a copy of the GNU General Public License along with
18the GNU MP Library test suite.  If not, see https://www.gnu.org/licenses/.  */
19
20#include <limits.h>
21#include <math.h>
22#include <stdlib.h>
23#include <stdio.h>
24#include <string.h>
25#include <float.h>
26
27#include "testutils.h"
28
29#define GMP_LIMB_BITS (sizeof(mp_limb_t) * CHAR_BIT)
30
31mp_bitcnt_t
32mpz_mantissasizeinbits (const mpz_t z)
33{
34  return ! mpz_cmp_ui (z, 0) ? 0 :
35    mpz_sizeinbase (z, 2) - mpz_scan1 (z, 0);
36}
37
38#if defined(DBL_MANT_DIG) && FLT_RADIX == 2
39int
40mpz_get_d_exact_p (const mpz_t z)
41{
42  return mpz_mantissasizeinbits (z) <= DBL_MANT_DIG;
43}
44#define HAVE_EXACT_P 1
45#endif
46
47#define COUNT 10000
48
49void
50test_matissa (void)
51{
52  mpz_t x, y;
53  int i, c;
54
55  mpz_init (x);
56  mpz_init (y);
57
58  mini_urandomb (y, 4);
59  c = i = mpz_get_ui (y);
60
61  do {
62    double d;
63    int cmp;
64
65    mpz_setbit (x, c);
66    d = mpz_get_d (x);
67    mpz_set_d (y, d);
68    if (mpz_cmp_d (y, d) != 0)
69      {
70	fprintf (stderr, "mpz_cmp_d (y, d) failed:\n"
71		 "d = %.20g\n"
72		 "i = %i\n"
73		 "c = %i\n",
74		 d, i, c);
75	abort ();
76      }
77
78    cmp = mpz_cmp (x, y);
79
80#if defined(HAVE_EXACT_P)
81    if ((mpz_get_d_exact_p (x) != 0) != (cmp == 0))
82      {
83	fprintf (stderr, "Not all bits converted:\n"
84		 "d = %.20g\n"
85		 "i = %i\n"
86		 "c = %i\n",
87		 d, i, c);
88	abort ();
89      }
90#endif
91
92    if (cmp < 0)
93      {
94	fprintf (stderr, "mpz_get_d failed:\n"
95		 "d = %.20g\n"
96		 "i = %i\n"
97		 "c = %i\n",
98		 d, i, c);
99	abort ();
100      }
101    else if (cmp > 0)
102      {
103	if (mpz_cmp_d (x, d) <= 0)
104	  {
105	    fprintf (stderr, "mpz_cmp_d (x, d) failed:\n"
106		     "d = %.20g\n"
107		     "i = %i\n"
108		     "c = %i\n",
109		     d, i, c);
110	    abort ();
111	  }
112	break;
113      }
114    ++c;
115  } while (1);
116
117  mpz_clear (x);
118  mpz_clear (y);
119}
120
121#ifndef M_PI
122#define M_PI 3.141592653589793238462643383279502884
123#endif
124
125static const struct
126{
127  double d;
128  const char *s;
129} values[] = {
130  { 0.0, "0" },
131  { 0.3, "0" },
132  { -0.3, "0" },
133  { M_PI, "3" },
134  { M_PI*1e15, "b29430a256d21" },
135  { -M_PI*1e15, "-b29430a256d21" },
136  /* 17 * 2^{200} =
137     27317946752402834684213355569799764242877450894307478200123392 */
138  {0.2731794675240283468421335556979976424288e62,
139    "1100000000000000000000000000000000000000000000000000" },
140  { 0.0, NULL }
141};
142
143void
144testmain (int argc, char **argv)
145{
146  unsigned i;
147  mpz_t x;
148
149  for (i = 0; values[i].s; i++)
150    {
151      char *s;
152      mpz_init_set_d (x, values[i].d);
153      s = mpz_get_str (NULL, 16, x);
154      if (strcmp (s, values[i].s) != 0)
155	{
156	  fprintf (stderr, "mpz_set_d failed:\n"
157		   "d = %.20g\n"
158		   "s = %s\n"
159		   "r = %s\n",
160		   values[i].d, s, values[i].s);
161	  abort ();
162	}
163      testfree (s);
164      mpz_clear (x);
165    }
166
167  mpz_init (x);
168
169  for (i = 0; i < COUNT; i++)
170    {
171      /* Use volatile, to avoid extended precision in floating point
172	 registers, e.g., on m68k and 80387. */
173      volatile double d, f;
174      unsigned long m;
175      int e;
176
177      mini_rrandomb (x, GMP_LIMB_BITS);
178      m = mpz_get_ui (x);
179      mini_urandomb (x, 8);
180      e = mpz_get_ui (x) - 100;
181
182      d = ldexp ((double) m, e);
183      mpz_set_d (x, d);
184      f = mpz_get_d (x);
185      if (f != floor (d))
186	{
187	  fprintf (stderr, "mpz_set_d/mpz_get_d failed:\n");
188	  goto dumperror;
189	}
190      if ((f == d) ? (mpz_cmp_d (x, d) != 0) : (mpz_cmp_d (x, d) >= 0))
191	{
192	  fprintf (stderr, "mpz_cmp_d (x, d) failed:\n");
193	  goto dumperror;
194	}
195      f = d + 1.0;
196      if (f > d && ! (mpz_cmp_d (x, f) < 0))
197	{
198	  fprintf (stderr, "mpz_cmp_d (x, f) failed:\n");
199	  goto dumperror;
200	}
201
202      d = - d;
203
204      mpz_set_d (x, d);
205      f = mpz_get_d (x);
206      if (f != ceil (d))
207	{
208	  fprintf (stderr, "mpz_set_d/mpz_get_d failed:\n");
209	dumperror:
210	  dump ("x", x);
211	  fprintf (stderr, "m = %lx, e = %i\n", m, e);
212	  fprintf (stderr, "d = %.15g\n", d);
213	  fprintf (stderr, "f = %.15g\n", f);
214	  fprintf (stderr, "f - d = %.5g\n", f - d);
215	  abort ();
216	}
217      if ((f == d) ? (mpz_cmp_d (x, d) != 0) : (mpz_cmp_d (x, d) <= 0))
218	{
219	  fprintf (stderr, "mpz_cmp_d (x, d) failed:\n");
220	  goto dumperror;
221	}
222      f = d - 1.0;
223      if (f < d && ! (mpz_cmp_d (x, f) > 0))
224	{
225	  fprintf (stderr, "mpz_cmp_d (x, f) failed:\n");
226	  goto dumperror;
227	}
228    }
229
230  mpz_clear (x);
231  test_matissa();
232}
233