1/* Test file for mpfr_agm. 2 3Copyright 1999, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc. 4Contributed by the Arenaire and Cacao 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 20http://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 <stdio.h> 24#include <stdlib.h> 25 26#include "mpfr-test.h" 27 28#define check(a,b,r) check4(a,b,r,0.0) 29 30static void 31check4 (const char *as, const char *bs, mpfr_rnd_t rnd_mode, 32 const char *res, int inex) 33{ 34 mpfr_t ta, tb, tc, tres; 35 mpfr_exp_t emin, emax; 36 int i; 37 38 emin = mpfr_get_emin (); 39 emax = mpfr_get_emax (); 40 41 mpfr_inits2 (53, ta, tb, tc, tres, (mpfr_ptr) 0); 42 43 for (i = 0; i <= 2; i++) 44 { 45 unsigned int expflags, newflags; 46 int inex2; 47 48 mpfr_set_str1 (ta, as); 49 mpfr_set_str1 (tb, bs); 50 mpfr_set_str1 (tc, res); 51 52 if (i > 0) 53 { 54 mpfr_exp_t ea, eb, ec, e0; 55 56 set_emin (MPFR_EMIN_MIN); 57 set_emax (MPFR_EMAX_MAX); 58 59 ea = mpfr_get_exp (ta); 60 eb = mpfr_get_exp (tb); 61 ec = mpfr_get_exp (tc); 62 63 e0 = i == 1 ? __gmpfr_emin : __gmpfr_emax; 64 if ((i == 1 && ea < eb) || (i == 2 && ea > eb)) 65 { 66 mpfr_set_exp (ta, e0); 67 mpfr_set_exp (tb, e0 + (eb - ea)); 68 mpfr_set_exp (tc, e0 + (ec - ea)); 69 } 70 else 71 { 72 mpfr_set_exp (ta, e0 + (ea - eb)); 73 mpfr_set_exp (tb, e0); 74 mpfr_set_exp (tc, e0 + (ec - eb)); 75 } 76 } 77 78 __gmpfr_flags = expflags = 79 (randlimb () & 1) ? MPFR_FLAGS_ALL ^ MPFR_FLAGS_ERANGE : 0; 80 inex2 = mpfr_agm (tres, ta, tb, rnd_mode); 81 newflags = __gmpfr_flags; 82 expflags |= MPFR_FLAGS_INEXACT; 83 84 if (SIGN (inex2) != inex || newflags != expflags || 85 ! mpfr_equal_p (tres, tc)) 86 { 87 printf ("mpfr_agm failed in rnd_mode=%s for\n", 88 mpfr_print_rnd_mode (rnd_mode)); 89 printf (" a = "); 90 mpfr_out_str (stdout, 10, 0, ta, MPFR_RNDN); 91 printf ("\n"); 92 printf (" b = "); 93 mpfr_out_str (stdout, 10, 0, tb, MPFR_RNDN); 94 printf ("\n"); 95 printf ("expected inex = %d, flags = %u,\n" 96 " ", inex, expflags); 97 mpfr_dump (tc); 98 printf ("got inex = %d, flags = %u,\n" 99 " ", inex2, newflags); 100 mpfr_dump (tres); 101 exit (1); 102 } 103 104 set_emin (emin); 105 set_emax (emax); 106 } 107 108 mpfr_clears (ta, tb, tc, tres, (mpfr_ptr) 0); 109} 110 111static void 112check_large (void) 113{ 114 mpfr_t a, b, agm; 115 int inex; 116 117 mpfr_init2 (a, 82); 118 mpfr_init2 (b, 82); 119 mpfr_init2 (agm, 82); 120 121 mpfr_set_ui (a, 1, MPFR_RNDN); 122 mpfr_set_str_binary (b, "0.1111101100001000000001011000110111101000001011111000100001000101010100011111110010E-39"); 123 mpfr_agm (agm, a, b, MPFR_RNDN); 124 mpfr_set_str_binary (a, "0.1110001000111101101010101010101101001010001001001011100101111011110101111001111100E-4"); 125 if (mpfr_cmp (agm, a)) 126 { 127 printf ("mpfr_agm failed for precision 82\n"); 128 exit (1); 129 } 130 131 /* problem found by Damien Fischer <damien@maths.usyd.edu.au> 4 Aug 2003: 132 produced a division by zero exception */ 133 mpfr_set_prec (a, 268); 134 mpfr_set_prec (b, 268); 135 mpfr_set_prec (agm, 268); 136 mpfr_set_str (a, "703.93543315330225238487276503953366664991725089988315253092140138947103394917006", 10, MPFR_RNDN); 137 mpfr_set_str (b, "703.93543315330225238487279020523738740563816490895994499256063816906728642622316", 10, MPFR_RNDN); 138 mpfr_agm (agm, a, b, MPFR_RNDN); 139 140 mpfr_set_prec (a, 18); 141 mpfr_set_prec (b, 70); 142 mpfr_set_prec (agm, 67); 143 mpfr_set_str_binary (a, "0.111001111100101000e8"); 144 mpfr_set_str_binary (b, "0.1101110111100100010100110000010111011011011100110100111001010100100001e10"); 145 inex = mpfr_agm (agm, a, b, MPFR_RNDN); 146 mpfr_set_str_binary (b, "0.1111110010011101101100010101011011010010010000001010100011000110011e9"); 147 if (mpfr_cmp (agm, b)) 148 { 149 printf ("Error in mpfr_agm (1)\n"); 150 exit (1); 151 } 152 if (inex >= 0) 153 { 154 printf ("Wrong flag for mpfr_agm (1)\n"); 155 exit (1); 156 } 157 158 /* test worst case: 9 consecutive ones after the last bit */ 159 mpfr_set_prec (a, 2); 160 mpfr_set_prec (b, 2); 161 mpfr_set_ui (a, 1, MPFR_RNDN); 162 mpfr_set_ui (b, 2, MPFR_RNDN); 163 mpfr_set_prec (agm, 904); 164 mpfr_agm (agm, a, b, MPFR_RNDZ); 165 166 mpfr_clear (a); 167 mpfr_clear (b); 168 mpfr_clear (agm); 169} 170 171static void 172check_nans (void) 173{ 174 mpfr_t x, y, m; 175 176 mpfr_init2 (x, 123L); 177 mpfr_init2 (y, 123L); 178 mpfr_init2 (m, 123L); 179 180 /* agm(1,nan) == nan */ 181 mpfr_set_ui (x, 1L, MPFR_RNDN); 182 mpfr_set_nan (y); 183 mpfr_agm (m, x, y, MPFR_RNDN); 184 MPFR_ASSERTN (mpfr_nan_p (m)); 185 186 /* agm(1,+inf) == +inf */ 187 mpfr_set_ui (x, 1L, MPFR_RNDN); 188 mpfr_set_inf (y, 1); 189 mpfr_agm (m, x, y, MPFR_RNDN); 190 MPFR_ASSERTN (mpfr_inf_p (m)); 191 MPFR_ASSERTN (mpfr_sgn (m) > 0); 192 193 /* agm(+inf,+inf) == +inf */ 194 mpfr_set_inf (x, 1); 195 mpfr_set_inf (y, 1); 196 mpfr_agm (m, x, y, MPFR_RNDN); 197 MPFR_ASSERTN (mpfr_inf_p (m)); 198 MPFR_ASSERTN (mpfr_sgn (m) > 0); 199 200 /* agm(-inf,+inf) == nan */ 201 mpfr_set_inf (x, -1); 202 mpfr_set_inf (y, 1); 203 mpfr_agm (m, x, y, MPFR_RNDN); 204 MPFR_ASSERTN (mpfr_nan_p (m)); 205 206 /* agm(+0,+inf) == nan */ 207 mpfr_set_ui (x, 0, MPFR_RNDN); 208 mpfr_set_inf (y, 1); 209 mpfr_agm (m, x, y, MPFR_RNDN); 210 MPFR_ASSERTN (mpfr_nan_p (m)); 211 212 /* agm(+0,1) == +0 */ 213 mpfr_set_ui (x, 0, MPFR_RNDN); 214 mpfr_set_ui (y, 1, MPFR_RNDN); 215 mpfr_agm (m, x, y, MPFR_RNDN); 216 MPFR_ASSERTN (MPFR_IS_ZERO (m) && MPFR_IS_POS(m)); 217 218 /* agm(-0,1) == +0 */ 219 mpfr_set_ui (x, 0, MPFR_RNDN); 220 mpfr_neg (x, x, MPFR_RNDN); 221 mpfr_set_ui (y, 1, MPFR_RNDN); 222 mpfr_agm (m, x, y, MPFR_RNDN); 223 MPFR_ASSERTN (MPFR_IS_ZERO (m) && MPFR_IS_POS(m)); 224 225 /* agm(-0,+0) == +0 */ 226 mpfr_set_ui (x, 0, MPFR_RNDN); 227 mpfr_neg (x, x, MPFR_RNDN); 228 mpfr_set_ui (y, 0, MPFR_RNDN); 229 mpfr_agm (m, x, y, MPFR_RNDN); 230 MPFR_ASSERTN (MPFR_IS_ZERO (m) && MPFR_IS_POS(m)); 231 232 /* agm(1,1) == 1 */ 233 mpfr_set_ui (x, 1, MPFR_RNDN); 234 mpfr_set_ui (y, 1, MPFR_RNDN); 235 mpfr_agm (m, x, y, MPFR_RNDN); 236 MPFR_ASSERTN (mpfr_cmp_ui (m ,1) == 0); 237 238 /* agm(-1,-2) == NaN */ 239 mpfr_set_si (x, -1, MPFR_RNDN); 240 mpfr_set_si (y, -2, MPFR_RNDN); 241 mpfr_agm (m, x, y, MPFR_RNDN); 242 MPFR_ASSERTN (mpfr_nan_p (m)); 243 244 mpfr_clear (x); 245 mpfr_clear (y); 246 mpfr_clear (m); 247} 248 249#define TEST_FUNCTION mpfr_agm 250#define TWO_ARGS 251#define TEST_RANDOM_POS 4 252#define TEST_RANDOM_POS2 4 253#include "tgeneric.c" 254 255int 256main (int argc, char* argv[]) 257{ 258 tests_start_mpfr (); 259 260 check_nans (); 261 262 check_large (); 263 check4 ("2.0", "1.0", MPFR_RNDN, "1.456791031046906869", -1); 264 check4 ("6.0", "4.0", MPFR_RNDN, "4.949360872472608925", 1); 265 check4 ("62.0", "61.0", MPFR_RNDN, "61.498983718845075902", -1); 266 check4 ("0.5", "1.0", MPFR_RNDN, "0.72839551552345343459", -1); 267 check4 ("1.0", "2.0", MPFR_RNDN, "1.456791031046906869", -1); 268 check4 ("234375765.0", "234375000.0", MPFR_RNDN, "234375382.49984394025", 1); 269 check4 ("8.0", "1.0", MPFR_RNDU, "3.61575617759736274873", 1); 270 check4 ("1.0", "44.0", MPFR_RNDU, "13.3658354512981243907", 1); 271 check4 ("1.0", "3.7252902984619140625e-9", MPFR_RNDU, 272 "0.07553933569711989657765", 1); 273 test_generic (2, 300, 17); 274 275 tests_end_mpfr (); 276 return 0; 277} 278