1 /* mpfr_sqr -- Floating square
3 Copyright 2004, 2005, 2006, 2007, 2008, 2009, 2010 Free Software Foundation, Inc.
4 Contributed by the Arenaire and Cacao projects, INRIA.
6 This file is part of the GNU MPFR Library.
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.
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.
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. */
23 #include "mpfr-impl.h"
26 mpfr_sqr (mpfr_ptr a, mpfr_srcptr b, mpfr_rnd_t rnd_mode)
34 MPFR_TMP_DECL(marker);
36 MPFR_LOG_FUNC (("x[%#R]=%R rnd=%d", b, b, rnd_mode),
37 ("y[%#R]=%R inexact=%d", a, a, inexact));
39 /* deal with special cases */
40 if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(b)))
51 ( MPFR_ASSERTD(MPFR_IS_ZERO(b)), MPFR_SET_ZERO(a) );
54 ax = 2 * MPFR_GET_EXP (b);
57 MPFR_ASSERTD (2 * bq > bq); /* PREC_MAX is /2 so no integer overflow */
59 bn = MPFR_LIMB_SIZE(b); /* number of limbs of b */
60 tn = 1 + (2 * bq - 1) / GMP_NUMB_BITS; /* number of limbs of square,
63 MPFR_TMP_MARK(marker);
64 tmp = (mp_limb_t *) MPFR_TMP_ALLOC((size_t) 2 * bn * BYTES_PER_MP_LIMB);
66 /* Multiplies the mantissa in temporary allocated space */
67 mpn_sqr_n (tmp, MPFR_MANT(b), bn);
70 /* now tmp[0]..tmp[2*bn-1] contains the product of both mantissa,
71 with tmp[2*bn-1]>=2^(GMP_NUMB_BITS-2) */
72 b1 >>= GMP_NUMB_BITS - 1; /* msb from the product */
74 /* if the mantissas of b and c are uniformly distributed in ]1/2, 1],
75 then their product is in ]1/4, 1/2] with probability 2*ln(2)-1 ~ 0.386
76 and in [1/2, 1] with probability 2-2*ln(2) ~ 0.614 */
77 tmp += 2 * bn - tn; /* +0 or +1 */
78 if (MPFR_UNLIKELY(b1 == 0))
79 mpn_lshift (tmp, tmp, tn, 1); /* tn <= k, so no stack corruption */
81 cc = mpfr_round_raw (MPFR_MANT (a), tmp, 2 * bq, 0,
82 MPFR_PREC (a), rnd_mode, &inexact);
83 /* cc = 1 ==> result is a power of two */
84 if (MPFR_UNLIKELY(cc))
85 MPFR_MANT(a)[MPFR_LIMB_SIZE(a)-1] = MPFR_LIMB_HIGHBIT;
87 MPFR_TMP_FREE(marker);
89 mpfr_exp_t ax2 = ax + (mpfr_exp_t) (b1 - 1 + cc);
90 if (MPFR_UNLIKELY( ax2 > __gmpfr_emax))
91 return mpfr_overflow (a, rnd_mode, MPFR_SIGN_POS);
92 if (MPFR_UNLIKELY( ax2 < __gmpfr_emin))
94 /* In the rounding to the nearest mode, if the exponent of the exact
95 result (i.e. before rounding, i.e. without taking cc into account)
96 is < __gmpfr_emin - 1 or the exact result is a power of 2 (i.e. if
97 both arguments are powers of 2), then round to zero. */
98 if (rnd_mode == MPFR_RNDN &&
99 (ax + (mpfr_exp_t) b1 < __gmpfr_emin || mpfr_powerof2_raw (b)))
100 rnd_mode = MPFR_RNDZ;
101 return mpfr_underflow (a, rnd_mode, MPFR_SIGN_POS);
103 MPFR_SET_EXP (a, ax2);