xref: /netbsd-src/external/lgpl3/gmp/dist/mini-gmp/tests/t-mpq_double.c (revision 9fd8799cb5ceb66c69f2eb1a6d26a1d587ba1f1e)
1 /* Test mpq_set_d.
2 
3 Copyright 2001-2003, 2005, 2013, 2018 Free Software Foundation, Inc.
4 
5 This file is part of the GNU MP Library test suite.
6 
7 The GNU MP Library test suite is free software; you can redistribute it
8 and/or modify it under the terms of the GNU General Public License as
9 published by the Free Software Foundation; either version 3 of the License,
10 or (at your option) any later version.
11 
12 The GNU MP Library test suite is distributed in the hope that it will be
13 useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
15 Public License for more details.
16 
17 You should have received a copy of the GNU General Public License along with
18 the GNU MP Library test suite.  If not, see https://www.gnu.org/licenses/.  */
19 
20 #include <math.h>
21 #include <float.h>
22 #include <limits.h>
23 
24 #include "testutils.h"
25 #include "../mini-mpq.h"
26 
27 #define COUNT 2000
28 
29 mp_bitcnt_t
30 mpz_mantissasizeinbits (const mpz_t z)
31 {
32   return ! mpz_cmp_ui (z, 0) ? 0 :
33     mpz_sizeinbase (z, 2) - mpz_scan1 (z, 0);
34 }
35 
36 int
37 mpz_abspow2_p (const mpz_t z)
38 {
39   return mpz_mantissasizeinbits (z) == 1;
40 }
41 
42 mp_bitcnt_t
43 mpq_mantissasizeinbits (const mpq_t q)
44 {
45   if (! mpz_abspow2_p (mpq_denref (q)))
46     return ~ (mp_bitcnt_t) 0;
47 
48   return mpz_mantissasizeinbits (mpq_numref (q));
49 }
50 
51 #if defined(DBL_MANT_DIG) && FLT_RADIX == 2
52 int
53 mpz_get_d_exact_p (const mpz_t z)
54 {
55   return mpz_mantissasizeinbits (z) <= DBL_MANT_DIG;
56 }
57 
58 int
59 mpq_get_d_exact_p (const mpq_t q)
60 {
61   /* return mpq_mantissasizeinbits (q) <= DBL_MANT_DIG; */
62   return
63     (mpz_sizeinbase (mpq_denref (q), 2) -
64      mpz_scan1 (mpq_denref (q), 0) == 1) &&
65     (mpz_sizeinbase (mpq_numref (q), 2) -
66      mpz_scan1 (mpq_numref (q), 0) <= DBL_MANT_DIG);
67   /* mpz_sizeinbase (zero, 2) - mpz_scan1 (zero, 0) == 2 */
68 }
69 #define HAVE_EXACT_P 1
70 #endif
71 
72 void
73 check_random (void)
74 {
75   unsigned i;
76   mpz_t x;
77   mpq_t y, z;
78 
79   mpz_init (x);
80   mpq_init (y);
81   mpq_init (z);
82 
83   for (i = 0; i < COUNT; i++)
84     {
85       /* Use volatile, to avoid extended precision in floating point
86 	 registers, e.g., on m68k and 80387. */
87       volatile double d, f;
88       unsigned long m;
89       int e, c;
90 
91       mini_rrandomb (x, CHAR_BIT * sizeof (unsigned long));
92       m = mpz_get_ui (x);
93       mini_urandomb (x, 8);
94       e = mpz_get_ui (x) - 128;
95 
96       d = ldexp ((double) m, e);
97       mpq_set_d (y, d);
98       f = mpq_get_d (y);
99       if (f != d)
100 	{
101 	  fprintf (stderr, "mpq_set_d/mpq_get_d failed:\n");
102 	  goto dumperror;
103 	}
104 
105       d = - d;
106       mpq_neg (y, y);
107 
108       mpq_set_d (z, d);
109       f = mpq_get_d (z);
110       if (f != d || !mpq_equal (y, z))
111 	{
112 	  fprintf (stderr, "mpq_set_d/mpq_get_d failed:\n");
113 	dumperror:
114 	  dump ("ny", mpq_numref (y));
115 	  dump ("dy", mpq_denref (y));
116 	  fprintf (stderr, "m = %lx, e = %i\n", m, e);
117 	  fprintf (stderr, "d = %.35g\n", d);
118 	  fprintf (stderr, "f = %.35g\n", f);
119 	  fprintf (stderr, "f - d = %.35g\n", f - d);
120 	  abort ();
121 	}
122 
123       mini_rrandomb (x, CHAR_BIT * sizeof (unsigned long));
124       m = mpz_get_ui (x);
125       mini_urandomb (x, 8);
126       e = mpz_get_ui (x) - 128;
127 
128       d = ldexp ((double) m, e);
129       mpq_set_d (y, d);
130 
131       mpq_add (y, y, z);
132       mpq_set_d (z, mpq_get_d (y));
133       f = mpq_get_d (z);
134       c = mpq_cmp (y, z);
135 
136 #if defined(HAVE_EXACT_P)
137       if (mpq_get_d_exact_p (y) ? c != 0 : (f > 0 ? c <= 0 : c >= 0))
138 #else
139       if (f > 0 ? c < 0 : c > 0)
140 #endif
141 	{
142 	  fprintf (stderr, "mpq_get_d/mpq_set_d failed: %i %i\n", i, c);
143 	  goto dumperror;
144 	}
145     }
146 
147   mpz_clear (x);
148   mpq_clear (y);
149   mpq_clear (z);
150 }
151 
152 
153 void
154 check_data (void)
155 {
156   static const struct {
157     double        y;
158     long int      n;
159     unsigned long d;
160   } data[] = {
161     {  0.0,  0, 1 },
162     {  1.0,  1, 1 },
163     { -1.0, -1, 1 },
164     { -1.5, -3, 2 },
165     {-1.25, -5, 4 },
166     {0.125,  1, 8 },
167 
168     {24685,24685,1},
169     {-9876,-9876,1},
170     {463.5,  927,2},
171 
172     {1234.5/8192,  2469, 16384 },
173     {-543.0/1024,  -543,  1024 },
174     {9876.5/ 512, 19753,  1024 },
175     {9753.0/ 128,  9753,   128 },
176     {-789.0/  32,  -789,    32 },
177     {4.580078125,  2345,   512 },
178   };
179 
180   mpq_t    x, r;
181   unsigned i;
182   double d;
183 
184   mpq_init (x);
185   mpq_init (r);
186 
187   for (i = 0; i < numberof (data); i++)
188     {
189       mpq_set_d (x, data[i].y);
190       mpq_set_si (r, data[i].n, data[i].d);
191       mpq_canonicalize (r);
192       if (!mpq_equal (x, r))
193 	{
194 	  fprintf (stderr, "mpq_set_d failed: %li / %lu != %g\n", data[i].n, data[i].d, data[i].y);
195 	  abort ();
196 	}
197       d = mpq_get_d (r);
198       if (d != data[i].y)
199 	{
200 	  fprintf (stderr, "mpq_get_d failed: %li / %lu != %g\n", data[i].n, data[i].d, data[i].y);
201 	  abort ();
202 	}
203     }
204 
205   mpq_clear (x);
206   mpq_clear (r);
207 }
208 
209 void
210 testmain (int argc, char *argv[])
211 {
212   check_data ();
213   check_random ();
214 }
215