1/* generic inverse of a function. 2 3Copyright 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4Contributed by the AriC and Caramel 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#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 43int 44FUNCTION (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