1257483Sian/* mpfr_sqr -- Floating square 2257483Sian 3257483SianCopyright 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4257483SianContributed by the AriC and Caramel projects, INRIA. 5257483Sian 6257483SianThis file is part of the GNU MPFR Library. 7257483Sian 8257483SianThe GNU MPFR Library is free software; you can redistribute it and/or modify 9257483Sianit under the terms of the GNU Lesser General Public License as published by 10257483Sianthe Free Software Foundation; either version 3 of the License, or (at your 11261946Sianoption) any later version. 12257483Sian 13257483SianThe GNU MPFR Library is distributed in the hope that it will be useful, but 14257483SianWITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15262427Sianor FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16257483SianLicense for more details. 17290807Sgonzo 18323418SianYou should have received a copy of the GNU Lesser General Public License 19271550Sianalong with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20257483Sianhttp://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 21323418Sian51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22268835Sbr 23268977Sbr#include "mpfr-impl.h" 24277644Sbr 25277644Sbrint 26277644Sbrmpfr_sqr (mpfr_ptr a, mpfr_srcptr b, mpfr_rnd_t rnd_mode) 27257483Sian{ 28292571Sgonzo int cc, inexact; 29292571Sgonzo mpfr_exp_t ax; 30292571Sgonzo mp_limb_t *tmp; 31292574Sgonzo mp_limb_t b1; 32292574Sgonzo mpfr_prec_t bq; 33257483Sian mp_size_t bn, tn; 34257483Sian MPFR_TMP_DECL(marker); 35257483Sian 36314510Sian MPFR_LOG_FUNC 37257483Sian (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (b), mpfr_log_prec, b, rnd_mode), 38257483Sian ("y[%Pu]=%.*Rg inexact=%d", 39257483Sian mpfr_get_prec (a), mpfr_log_prec, a, inexact)); 40257483Sian 41257483Sian /* deal with special cases */ 42257483Sian if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(b))) 43257483Sian { 44257483Sian if (MPFR_IS_NAN(b)) 45257483Sian { 46257483Sian MPFR_SET_NAN(a); 47257483Sian MPFR_RET_NAN; 48257483Sian } 49257483Sian MPFR_SET_POS (a); 50257483Sian if (MPFR_IS_INF(b)) 51257483Sian MPFR_SET_INF(a); 52257483Sian else 53257483Sian ( MPFR_ASSERTD(MPFR_IS_ZERO(b)), MPFR_SET_ZERO(a) ); 54257483Sian MPFR_RET(0); 55257483Sian } 56277644Sbr ax = 2 * MPFR_GET_EXP (b); 57277644Sbr bq = MPFR_PREC(b); 58277644Sbr 59277644Sbr MPFR_ASSERTN (2 * (mpfr_uprec_t) bq <= MPFR_PREC_MAX); 60277644Sbr 61277644Sbr bn = MPFR_LIMB_SIZE (b); /* number of limbs of b */ 62277644Sbr tn = MPFR_PREC2LIMBS (2 * bq); /* number of limbs of square, 63277644Sbr 2*bn or 2*bn-1 */ 64277644Sbr 65277644Sbr if (MPFR_UNLIKELY(bn > MPFR_SQR_THRESHOLD)) 66277644Sbr return mpfr_mul (a, b, b, rnd_mode); 67277644Sbr 68277644Sbr MPFR_TMP_MARK(marker); 69277644Sbr tmp = MPFR_TMP_LIMBS_ALLOC (2 * bn); 70277644Sbr 71277644Sbr /* Multiplies the mantissa in temporary allocated space */ 72 mpn_sqr_n (tmp, MPFR_MANT(b), bn); 73 b1 = tmp[2 * bn - 1]; 74 75 /* now tmp[0]..tmp[2*bn-1] contains the product of both mantissa, 76 with tmp[2*bn-1]>=2^(GMP_NUMB_BITS-2) */ 77 b1 >>= GMP_NUMB_BITS - 1; /* msb from the product */ 78 79 /* if the mantissas of b and c are uniformly distributed in ]1/2, 1], 80 then their product is in ]1/4, 1/2] with probability 2*ln(2)-1 ~ 0.386 81 and in [1/2, 1] with probability 2-2*ln(2) ~ 0.614 */ 82 tmp += 2 * bn - tn; /* +0 or +1 */ 83 if (MPFR_UNLIKELY(b1 == 0)) 84 mpn_lshift (tmp, tmp, tn, 1); /* tn <= k, so no stack corruption */ 85 86 cc = mpfr_round_raw (MPFR_MANT (a), tmp, 2 * bq, 0, 87 MPFR_PREC (a), rnd_mode, &inexact); 88 /* cc = 1 ==> result is a power of two */ 89 if (MPFR_UNLIKELY(cc)) 90 MPFR_MANT(a)[MPFR_LIMB_SIZE(a)-1] = MPFR_LIMB_HIGHBIT; 91 92 MPFR_TMP_FREE(marker); 93 { 94 mpfr_exp_t ax2 = ax + (mpfr_exp_t) (b1 - 1 + cc); 95 if (MPFR_UNLIKELY( ax2 > __gmpfr_emax)) 96 return mpfr_overflow (a, rnd_mode, MPFR_SIGN_POS); 97 if (MPFR_UNLIKELY( ax2 < __gmpfr_emin)) 98 { 99 /* In the rounding to the nearest mode, if the exponent of the exact 100 result (i.e. before rounding, i.e. without taking cc into account) 101 is < __gmpfr_emin - 1 or the exact result is a power of 2 (i.e. if 102 both arguments are powers of 2), then round to zero. */ 103 if (rnd_mode == MPFR_RNDN && 104 (ax + (mpfr_exp_t) b1 < __gmpfr_emin || mpfr_powerof2_raw (b))) 105 rnd_mode = MPFR_RNDZ; 106 return mpfr_underflow (a, rnd_mode, MPFR_SIGN_POS); 107 } 108 MPFR_SET_EXP (a, ax2); 109 MPFR_SET_POS (a); 110 } 111 MPFR_RET (inexact); 112} 113