1 /* mpfr_sqrt -- square root of a floating-point number
3 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 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 2.1 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.LIB. If not, write to
20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
21 MA 02110-1301, USA. */
23 #include "mpfr-impl.h"
26 mpfr_sqrt (mpfr_ptr r
, mpfr_srcptr u
, mp_rnd_t rnd_mode
)
28 mp_size_t rsize
; /* number of limbs of r */
30 mp_size_t usize
; /* number of limbs of u */
31 mp_size_t tsize
; /* number of limbs of the sqrtrem remainder */
38 mp_limb_t sticky0
; /* truncated part of input */
39 mp_limb_t sticky1
; /* truncated part of rp[0] */
42 int sh
; /* number of extra bits in rp[0] */
43 int inexact
; /* return ternary flag */
45 MPFR_TMP_DECL(marker
);
47 MPFR_LOG_FUNC (("x[%#R]=%R rnd=%d", u
, u
, rnd_mode
),
48 ("y[%#R]=%R inexact=%d", r
, r
, inexact
));
50 if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(u
)))
57 else if (MPFR_IS_ZERO(u
))
60 MPFR_SET_SAME_SIGN(r
, u
);
62 MPFR_RET(0); /* zero is exact */
66 MPFR_ASSERTD(MPFR_IS_INF(u
));
67 /* sqrt(-Inf) = NAN */
78 if (MPFR_UNLIKELY(MPFR_IS_NEG(u
)))
86 rsize
= MPFR_LIMB_SIZE(r
); /* number of limbs of r */
87 rrsize
= rsize
+ rsize
;
88 usize
= MPFR_LIMB_SIZE(u
); /* number of limbs of u */
91 sticky0
= MPFR_LIMB_ZERO
; /* truncated part of input */
92 sticky1
= MPFR_LIMB_ZERO
; /* truncated part of rp[0] */
93 odd_exp
= (unsigned int) MPFR_GET_EXP (u
) & 1;
94 inexact
= -1; /* return ternary flag */
96 MPFR_TMP_MARK (marker
);
97 sp
= (mp_limb_t
*) MPFR_TMP_ALLOC (rrsize
* sizeof (mp_limb_t
));
99 /* copy the most significant limbs of u to {sp, rrsize} */
100 if (MPFR_LIKELY(usize
<= rrsize
)) /* in case r and u have the same precision,
101 we have indeed rrsize = 2 * usize */
109 sp
[k
- 1] = mpn_rshift (sp
+ k
, up
, usize
, 1);
111 sticky0
= mpn_rshift (sp
, up
, usize
, 1);
114 MPN_COPY (sp
+ rrsize
- usize
, up
, usize
);
116 else /* usize > rrsize: truncate the input */
120 sticky0
= mpn_rshift (sp
, up
+ k
, rrsize
, 1);
122 MPN_COPY (sp
, up
+ k
, rrsize
);
124 while (sticky0
== MPFR_LIMB_ZERO
&& l
!= 0)
128 /* sticky0 is non-zero iff the truncated part of the input is non-zero */
130 tsize
= mpn_sqrtrem (rp
, tp
= sp
, sp
, rrsize
);
134 while (sticky
== MPFR_LIMB_ZERO
&& l
!= 0)
137 /* truncated low bits of rp[0] */
138 MPFR_UNSIGNED_MINUS_MODULO(sh
,MPFR_PREC(r
));
139 sticky1
= rp
[0] & MPFR_LIMB_MASK(sh
);
142 sticky
= sticky
|| sticky1
;
144 expr
= (MPFR_GET_EXP(u
) + odd_exp
) / 2; /* exact */
146 if (rnd_mode
== GMP_RNDZ
|| rnd_mode
== GMP_RNDD
|| sticky
== MPFR_LIMB_ZERO
)
148 inexact
= (sticky
== MPFR_LIMB_ZERO
) ? 0 : -1;
151 else if (rnd_mode
== GMP_RNDN
)
153 /* if sh>0, the round bit is bit (sh-1) of sticky1
154 and the sticky bit is formed by the low sh-1 bits from
155 sticky1, together with {tp, tsize} and sticky0. */
158 if (sticky1
& (MPFR_LIMB_ONE
<< (sh
- 1)))
159 { /* round bit is set */
160 if (sticky1
== (MPFR_LIMB_ONE
<< (sh
- 1)) && tsize
== 0
166 else /* round bit is zero */
167 goto truncate
; /* with the default inexact=-1 */
171 /* if sh=0, we have to compare {tp, tsize} with {rp, rsize}:
172 if {tp, tsize} < {rp, rsize}: truncate
173 if {tp, tsize} > {rp, rsize}: round up
174 if {tp, tsize} = {rp, rsize}: compare the truncated part of the
178 if = 1/4: even rounding rule
179 Set inexact = -1 if truncate
180 inexact = 1 if add one ulp
181 inexact = 0 if even rounding rule
185 else if (tsize
> rsize
) /* FIXME: may happen? */
187 else /* tsize = rsize */
191 cmp
= mpn_cmp (tp
, rp
, rsize
);
194 else if (cmp
< 0 || sticky0
== MPFR_LIMB_ZERO
)
196 /* now tricky case {tp, tsize} = {rp, rsize} */
197 /* in case usize <= rrsize, the only case where sticky0 <> 0
198 is when the exponent of u is odd and usize = rrsize (k=0),
199 but in that case the truncated part is exactly 1/2, thus
201 If the exponent of u is odd, and up[k] is odd, the truncated
202 part is >= 1/2, so we round up too. */
203 else if (usize
<= rrsize
|| (odd_exp
&& (up
[k
] & MPFR_LIMB_ONE
)))
207 /* now usize > rrsize:
208 (a) if the exponent of u is even, the 1/4 bit is the
209 2nd most significant bit of up[k-1];
210 (b) if the exponent of u is odd, the 1/4 bit is the
211 1st most significant bit of up[k-1]; */
212 sticky1
= MPFR_LIMB_ONE
<< (BITS_PER_MP_LIMB
- 2 + odd_exp
);
213 if (up
[k
- 1] < sticky1
)
215 else if (up
[k
- 1] > sticky1
)
219 /* up[k - 1] == sticky1: consider low k-1 limbs */
220 while (--k
> 0 && up
[k
- 1] == MPFR_LIMB_ZERO
)
224 } /* end of case {tp, tsize} = {rp, rsize} */
225 } /* end of case tsize = rsize */
228 else if (inexact
== 1)
230 /* else go through even_rule */
233 else /* rnd_mode=GMP_RDNU, necessarily sticky <> 0, thus add 1 ulp */
236 even_rule
: /* has to set inexact */
237 inexact
= (rp
[0] & (MPFR_LIMB_ONE
<< sh
)) ? 1 : -1;
240 /* else go through add_one_ulp */
243 inexact
= 1; /* always here */
244 if (mpn_add_1 (rp
, rp
, rsize
, MPFR_LIMB_ONE
<< sh
))
247 rp
[rsize
- 1] = MPFR_LIMB_HIGHBIT
;
250 truncate
: /* inexact = 0 or -1 */
252 MPFR_ASSERTN (expr
>= MPFR_EMIN_MIN
&& expr
<= MPFR_EMAX_MAX
);
255 MPFR_TMP_FREE(marker
);
256 return mpfr_check_range (r
, inexact
, rnd_mode
);