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