1 /* generic inverse of a function. 2 3 Copyright 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4 Contributed by the AriC and Caramel projects, INRIA. 5 6 This file is part of the GNU MPFR Library. 7 8 The GNU MPFR Library is free software; you can redistribute it and/or modify 9 it under the terms of the GNU Lesser General Public License as published by 10 the Free Software Foundation; either version 3 of the License, or (at your 11 option) any later version. 12 13 The GNU MPFR Library is distributed in the hope that it will be useful, but 14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16 License for more details. 17 18 You should have received a copy of the GNU Lesser General Public License 19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23 #define MPFR_NEED_LONGLONG_H 24 #include "mpfr-impl.h" 25 26 #ifndef ACTION_SPECIAL 27 #define ACTION_SPECIAL 28 #endif 29 30 #ifndef ACTION_TINY 31 #define ACTION_TINY 32 #endif 33 34 /* example of use: 35 #define FUNCTION mpfr_sec 36 #define INVERSE mpfr_cos 37 #define ACTION_NAN(y) do { MPFR_SET_NAN(y); MPFR_RET_NAN; } while (1) 38 #define ACTION_INF(y) do { MPFR_SET_NAN(y); MPFR_RET_NAN; } while (1) 39 #define ACTION_ZERO(y) return mpfr_set_ui (y, 1, MPFR_RNDN) 40 #include "gen_inverse.h" 41 */ 42 43 int 44 FUNCTION (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) 45 { 46 mpfr_prec_t precy; /* target precision */ 47 mpfr_prec_t m; /* working precision */ 48 mpfr_t z; /* temporary variable to store INVERSE(x) */ 49 int inexact; /* inexact flag */ 50 MPFR_ZIV_DECL (loop); 51 MPFR_SAVE_EXPO_DECL (expo); 52 53 if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(x))) 54 { 55 if (MPFR_IS_NAN(x)) 56 ACTION_NAN(y); 57 else if (MPFR_IS_INF(x)) 58 ACTION_INF(y); 59 else /* x = 0 */ 60 ACTION_ZERO(y,x); 61 } 62 63 /* x is neither NaN, Inf nor zero */ 64 MPFR_SAVE_EXPO_MARK (expo); 65 ACTION_TINY (y, x, rnd_mode); /* special case for very small input x */ 66 precy = MPFR_PREC(y); 67 m = precy + MPFR_INT_CEIL_LOG2 (precy) + 3; 68 mpfr_init2 (z, m); 69 70 MPFR_ZIV_INIT (loop, m); 71 for(;;) 72 { 73 MPFR_BLOCK_DECL (flags); 74 75 MPFR_BLOCK (flags, INVERSE (z, x, MPFR_RNDZ)); /* error k_u < 1 ulp */ 76 /* FIXME: the following assumes that if an overflow happens with 77 MPFR_EMAX_MAX, then necessarily an underflow happens with 78 __gmpfr_emin */ 79 if (MPFR_OVERFLOW (flags)) 80 { 81 int s = MPFR_SIGN(z); 82 MPFR_ZIV_FREE (loop); 83 mpfr_clear (z); 84 MPFR_SAVE_EXPO_FREE (expo); 85 return mpfr_underflow (y, (rnd_mode == MPFR_RNDN) ? 86 MPFR_RNDZ : rnd_mode, s); 87 } 88 mpfr_ui_div (z, 1, z, MPFR_RNDN); 89 /* the error is less than c_w + 2*c_u*k_u (see algorithms.tex), 90 where c_w = 1/2, c_u = 1 since z was rounded toward zero, 91 thus 1/2 + 2 < 4 */ 92 if (MPFR_LIKELY (MPFR_CAN_ROUND (z, m - 2, precy, rnd_mode))) 93 break; 94 ACTION_SPECIAL; 95 MPFR_ZIV_NEXT (loop, m); 96 mpfr_set_prec (z, m); 97 } 98 MPFR_ZIV_FREE (loop); 99 100 inexact = mpfr_set (y, z, rnd_mode); 101 mpfr_clear (z); 102 103 end: 104 MPFR_SAVE_EXPO_FREE (expo); 105 return mpfr_check_range (y, inexact, rnd_mode); 106 } 107