1/* Test file for mpfr_ai.
2
3Copyright 2010-2023 Free Software Foundation, Inc.
4Contributed by the AriC and Caramba projects, INRIA.
5
6This file is part of the GNU MPFR Library.
7
8The GNU MPFR 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 MPFR 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 MPFR Library; see the file COPYING.LESSER.  If not, see
20https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22
23#include "mpfr-test.h"
24
25#define TEST_FUNCTION mpfr_ai
26#define TEST_RANDOM_EMIN -5
27#define TEST_RANDOM_EMAX 5
28#define REDUCE_EMAX 7 /* this is to avoid that test_generic() calls mpfr_ai
29                         with too large inputs. FIXME: remove this once
30                         mpfr_ai can handle large inputs */
31#include "tgeneric.c"
32
33static void
34check_large (void)
35{
36  mpfr_t x, y, z;
37
38  mpfr_init2 (x, 38);
39  mpfr_init2 (y, 110);
40  mpfr_init2 (z, 110);
41  mpfr_set_str_binary (x, "-1E8");
42  mpfr_ai (y, x, MPFR_RNDN);
43  mpfr_set_str_binary (z, "-10001110100001011111110001100011101100011100010000110100100101011111011100000101110101010010000000101110011111E-112");
44  if (mpfr_equal_p (y, z) == 0)
45    {
46      printf ("Error in mpfr_ai for x=-2^8\n");
47      exit (1);
48    }
49#if 0 /* disabled since mpfr_ai does not currently handle large arguments */
50  mpfr_set_str_binary (x, "-1E26");
51  mpfr_ai (y, x, MPFR_RNDN);
52  mpfr_set_str_binary (z, "-110001111100000011001010010101001101001011001011101011001010100100001110001101101101000010000011001000001011E-118");
53  if (mpfr_equal_p (y, z) == 0)
54    {
55      printf ("Error in mpfr_ai for x=-2^26\n");
56      exit (1);
57    }
58  mpfr_set_str_binary (x, "-0.11111111111111111111111111111111111111E1073741823");
59  mpfr_ai (y, x, MPFR_RNDN);
60  /* FIXME: compute the correctly rounded value we should get for Ai(x),
61     and check we get this value */
62#endif
63  mpfr_clear (x);
64  mpfr_clear (y);
65  mpfr_clear (z);
66}
67
68static void
69check_zero (void)
70{
71  mpfr_t x, y, r;
72
73  mpfr_init2 (x, 53);
74  mpfr_init2 (y, 53);
75  mpfr_init2 (r, 53);
76
77  mpfr_set_str_binary (r, "10110101110001100011110010110001001110001010110111E-51");
78
79  mpfr_set_ui (x, 0, MPFR_RNDN);
80  mpfr_ai (y, x, MPFR_RNDN);
81  if (mpfr_equal_p (y, r) == 0)
82    {
83      printf ("Error in mpfr_ai for x=0\n");
84      printf ("Expected "); mpfr_dump (r);
85      printf ("Got      "); mpfr_dump (y);
86      exit (1);
87    }
88  mpfr_clear (x);
89  mpfr_clear (y);
90  mpfr_clear (r);
91}
92
93static void
94bug20180107 (void)
95{
96  mpfr_t x, y, z;
97  mpfr_exp_t emin;
98  int inex;
99  mpfr_flags_t flags;
100
101  mpfr_init2 (x, 152);
102  mpfr_init2 (y, 11);
103  mpfr_init2 (z, 11);
104  mpfr_set_str_binary (x, "0.11010101100111000111001001010110101001100001011110101111000010100111011101011110000100111011101100100100001010000110100011001000111010010001110000011100E5");
105  emin = mpfr_get_emin ();
106  set_emin (-134);
107  mpfr_clear_flags ();
108  inex = mpfr_ai (y, x, MPFR_RNDA);
109  flags = __gmpfr_flags;
110  /* result should be 0.10011100000E-135 with unlimited exponent range,
111     and thus should be rounded to 0.1E-134 */
112  mpfr_set_str_binary (z, "0.1E-134");
113  MPFR_ASSERTN (mpfr_equal_p (y, z));
114  MPFR_ASSERTN (inex > 0);
115  MPFR_ASSERTN (flags == (MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT));
116
117  mpfr_set_prec (x, 2);
118  mpfr_set_str_binary (x, "0.11E7");
119  mpfr_set_prec (y, 2);
120  mpfr_clear_flags ();
121  inex = mpfr_ai (y, x, MPFR_RNDA);
122  flags = __gmpfr_flags;
123  /* result should be 1.0E-908 with unlimited exponent range,
124     and thus should be rounded to 0.1E-134 */
125  mpfr_set_str_binary (z, "0.1E-134");
126  MPFR_ASSERTN (mpfr_equal_p (y, z));
127  MPFR_ASSERTN (inex > 0);
128  MPFR_ASSERTN (flags == (MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT));
129
130  set_emin (emin);
131  mpfr_clear (x);
132  mpfr_clear (y);
133  mpfr_clear (z);
134}
135
136/* exercise mpfr_ai near m*2^e, for precision p */
137static void
138test_near_m2e (long m, mpfr_exp_t e, mpfr_prec_t pmax)
139{
140  mpfr_t x, xx, y, yy;
141  mpfr_prec_t p;
142  int inex;
143
144  mpfr_clear_flags ();
145
146  /* first determine the smallest precision for which m*2^e is exact */
147  for (p = MPFR_PREC_MIN; p <= pmax; p++)
148    {
149      mpfr_init2 (x, p);
150      inex = mpfr_set_si_2exp (x, m, e, MPFR_RNDN);
151      mpfr_clear (x);
152      if (inex == 0)
153        break;
154    }
155
156  mpfr_init2 (x, p);
157  inex = mpfr_set_si_2exp (x, m, e, MPFR_RNDN);
158  MPFR_ASSERTN(inex == 0);
159
160  for (; p <= pmax; p++)
161    {
162      mpfr_init2 (y, p);
163      mpfr_init2 (xx, p);
164      mpfr_init2 (yy, p);
165      mpfr_prec_round (x, p, MPFR_RNDN);
166      mpfr_ai (y, x, MPFR_RNDN);
167      while (1)
168        {
169          mpfr_set (xx, x, MPFR_RNDN);
170          mpfr_nextbelow (xx);
171          mpfr_ai (yy, xx, MPFR_RNDN);
172          if (mpfr_cmpabs (yy, y) >= 0)
173            break;
174          else
175            {
176              mpfr_set (x, xx, MPFR_RNDN);
177              mpfr_set (y, yy, MPFR_RNDN);
178            }
179        }
180      while (1)
181        {
182          mpfr_set (xx, x, MPFR_RNDN);
183          mpfr_nextabove (xx);
184          mpfr_ai (yy, xx, MPFR_RNDN);
185          if (mpfr_cmpabs (yy, y) >= 0)
186            break;
187          else
188            {
189              mpfr_set (x, xx, MPFR_RNDN);
190              mpfr_set (y, yy, MPFR_RNDN);
191            }
192        }
193      mpfr_clear (y);
194      mpfr_clear (xx);
195      mpfr_clear (yy);
196    }
197
198  mpfr_clear (x);
199
200  /* Since some tests don't really check that the result is not NaN... */
201  MPFR_ASSERTN (! mpfr_nanflag_p ());
202}
203
204/* example provided by Sylvain Chevillard, which exercises the case
205   wprec < err + 1, and thus correct_bits = 0, in src/ai.c */
206static void
207coverage (void)
208{
209  mpfr_t x, y;
210  int inex;
211
212  mpfr_init2 (x, 800);
213  mpfr_init2 (y, 20);
214  mpfr_set_str (x, "-2.3381074104597670384891972524467354406385401456723878524838544372136680027002836477821640417313293202847600938532659527752254668583598667448688987168197275409731526749911127480659996456283534915503672", 10, MPFR_RNDN);
215  inex = mpfr_ai (y, x, MPFR_RNDN);
216  MPFR_ASSERTN(inex < 0);
217  MPFR_ASSERTN(mpfr_cmp_ui_2exp (y, 593131, -682) == 0);
218
219  mpfr_clear (x);
220  mpfr_clear (y);
221}
222
223int
224main (int argc, char *argv[])
225{
226  tests_start_mpfr ();
227
228  coverage ();
229  test_near_m2e (-5, -1, 100); /* exercise near -2.5 */
230  test_near_m2e (-4, 0, 100); /* exercise near -4 */
231  test_near_m2e (-11, -1, 100); /* exercise near -5.5 */
232  test_near_m2e (-27, -2, 100); /* exercise near -6.75 */
233  test_near_m2e (-31, -2, 100); /* exercise near -7.75 */
234  test_near_m2e (-15, -1, 100); /* exercise near -7.5 */
235  bug20180107 ();
236  check_large ();
237  check_zero ();
238
239  test_generic (MPFR_PREC_MIN, 100, 5);
240
241  tests_end_mpfr ();
242  return 0;
243}
244