summaryrefslogtreecommitdiff
path: root/mpfr/sqrt.c
diff options
context:
space:
mode:
Diffstat (limited to 'mpfr/sqrt.c')
-rw-r--r--mpfr/sqrt.c123
1 files changed, 68 insertions, 55 deletions
diff --git a/mpfr/sqrt.c b/mpfr/sqrt.c
index 8489581bd..839f769bd 100644
--- a/mpfr/sqrt.c
+++ b/mpfr/sqrt.c
@@ -1,20 +1,20 @@
/* mpfr_sqrt -- square root of a floating-point number
-Copyright (C) 1999, 2001 Free Software Foundation.
+Copyright (C) 1999, 2001 Free Software Foundation, Inc.
This file is part of the MPFR Library.
The MPFR Library is free software; you can redistribute it and/or modify
-it under the terms of the GNU Library General Public License as published by
-the Free Software Foundation; either version 2 of the License, or (at your
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 2.1 of the License, or (at your
option) any later version.
The MPFR Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
-or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
+or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
-You should have received a copy of the GNU Library General Public License
+You should have received a copy of the GNU Lesser General Public License
along with the MPFR Library; see the file COPYING.LIB. If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
@@ -29,14 +29,7 @@ MA 02111-1307, USA. */
/* #define DEBUG */
int
-#if __STDC__
mpfr_sqrt (mpfr_ptr r, mpfr_srcptr u, mp_rnd_t rnd_mode)
-#else
-mpfr_sqrt (r, u, rnd_mode)
- mpfr_ptr r;
- mpfr_srcptr u;
- mp_rnd_t rnd_mode;
-#endif
{
mp_ptr up, rp, tmp, remp;
mp_size_t usize, rrsize;
@@ -44,25 +37,25 @@ mpfr_sqrt (r, u, rnd_mode)
mp_size_t prec, err;
mp_limb_t q_limb;
long rw, nw, k;
- int exact = 0, t;
+ int inexact = 0, t;
unsigned long cc = 0;
char can_round = 0;
TMP_DECL(marker0);
if (MPFR_IS_NAN(u)) {
MPFR_SET_NAN(r);
- return 1;
+ return 1; /* NaN is always inexact */
}
if (MPFR_SIGN(u) < 0) {
if (MPFR_IS_INF(u) || MPFR_NOTZERO(u)) {
MPFR_SET_NAN(r);
- return 1;
+ return 1; /* NaN is always inexact */
}
else { /* sqrt(-0) = -0 */
MPFR_SET_ZERO(r);
if (MPFR_SIGN(r) > 0) MPFR_CHANGE_SIGN(r);
- return 0;
+ return 0; /* zero is exact */
}
}
@@ -72,7 +65,7 @@ mpfr_sqrt (r, u, rnd_mode)
if (MPFR_IS_INF(u))
{
MPFR_SET_INF(r);
- return 1;
+ return 0; /* infinity is exact */
}
MPFR_CLEAR_INF(r);
@@ -84,11 +77,10 @@ mpfr_sqrt (r, u, rnd_mode)
MPFR_EXP(r) = 0;
rsize = (prec-1)/BITS_PER_MP_LIMB + 1;
MPN_ZERO(MPFR_MANT(r), rsize);
- return 0;
+ return 0; /* zero is exact */
}
up = MPFR_MANT(u);
-
usize = (MPFR_PREC(u) - 1)/BITS_PER_MP_LIMB + 1;
#ifdef DEBUG
@@ -121,7 +113,7 @@ mpfr_sqrt (r, u, rnd_mode)
{
up = TMP_ALLOC((usize + 1)*BYTES_PER_MP_LIMB);
if (mpn_rshift(up + 1, MPFR_MANT(u), usize, 1))
- up [0] = ((mp_limb_t) 1) << (BITS_PER_MP_LIMB - 1);
+ up [0] = MP_LIMB_T_HIGHBIT;
else up[0] = 0;
usize++;
}
@@ -164,14 +156,15 @@ mpfr_sqrt (r, u, rnd_mode)
q_limb = mpn_sqrtrem (rp, remp, tmp, rsize);
#ifdef DEBUG
- printf("The result is : \n");
- printf("sqrt : ");
- for(k = rrsize - 1; k >= 0; k--) { printf("%lu ", rp[k]); }
- printf("(q_limb = %lu)\n", q_limb);
+ printf ("The result is : \n");
+ printf ("sqrt : ");
+ for (k = rrsize - 1; k >= 0; k--)
+ printf ("%lu ", rp[k]);
+ printf ("(inexact = %lu)\n", q_limb);
#endif
- can_round = (mpfr_can_round_raw(rp, rrsize, 1, err,
- GMP_RNDZ, rnd_mode, MPFR_PREC(r)));
+ can_round = mpfr_can_round_raw(rp, rrsize, 1, err,
+ GMP_RNDZ, rnd_mode, MPFR_PREC(r));
/* If we used all the limbs of both the dividend and the divisor,
then we have the correct RNDZ rounding */
@@ -186,67 +179,87 @@ mpfr_sqrt (r, u, rnd_mode)
}
while (!can_round && (rsize < 2*usize)
&& (rsize += 2) && (rrsize ++));
-
+#ifdef DEBUG
+ printf ("can_round = %d\n", can_round);
+#endif
/* This part may be deplaced upper to avoid a few mpfr_can_round_raw */
/* when the square root is exact. It is however very unprobable that */
/* it would improve the behaviour of the present code on average. */
- if (!q_limb) /* possibly exact */
+ if (!q_limb) /* the sqrtrem call was exact, possible exact square root */
{
/* if we have taken into account the whole of up */
for (k = usize - rsize - 1; k >= 0; k ++)
if (up[k]) break;
- if (k < 0) { exact = 1; goto fin; }
+ if (k < 0)
+ goto fin; /* exact square root ==> inexact = 0 */
}
if (can_round)
{
- cc = mpfr_round_raw(rp, rp, err, 0, MPFR_PREC(r), rnd_mode);
+ cc = mpfr_round_raw (rp, rp, err, 0, MPFR_PREC(r), rnd_mode, &inexact);
+ if (!inexact) /* exact high part: inexact flag depends from remainder */
+ inexact = -q_limb;
rrsize = (MPFR_PREC(r) - 1)/BITS_PER_MP_LIMB + 1;
}
else
/* Use the return value of sqrtrem to decide of the rounding */
/* Note that at this point the sqrt has been computed */
- /* EXACTLY. If rounding = GMP_RNDZ, do nothing [comes from */
- /* the fact that the exact square root can end with a bunch of ones, */
- /* and in that case we indeed cannot round if we do not know that */
- /* the computation was exact. */
+ /* EXACTLY. */
switch (rnd_mode)
{
case GMP_RNDZ :
- case GMP_RNDD : break;
+ case GMP_RNDD :
+ inexact = -1; /* result is truncated */
+ break;
case GMP_RNDN :
/* Not in the situation ...0 111111 */
rw = (MPFR_PREC(r) + 1) & (BITS_PER_MP_LIMB - 1);
- if (rw) { rw = BITS_PER_MP_LIMB - rw; nw = 0; } else nw = 1;
+ if (rw)
+ {
+ rw = BITS_PER_MP_LIMB - rw;
+ nw = 0;
+ }
+ else
+ nw = 1;
if ((rp[nw] >> rw) & 1 && /* Not 0111111111 */
(q_limb || /* Nonzero remainder */
(rw ? (rp[nw] >> (rw + 1)) & 1 :
(rp[nw] >> (BITS_PER_MP_LIMB - 1)) & 1))) /* or even rounding */
- cc = mpn_add_1(rp + nw, rp + nw, rrsize, ((mp_limb_t)1) << rw);
+ {
+ cc = mpn_add_1 (rp + nw, rp + nw, rrsize, MP_LIMB_T_ONE << rw);
+ inexact = 1;
+ }
+ else
+ inexact = -1;
break;
- case GMP_RNDU :
- if (q_limb)
- {
- t = MPFR_PREC(r) & (BITS_PER_MP_LIMB - 1);
- if (t)
- {
- cc = mpn_add_1(rp, rp, rrsize, 1 << (BITS_PER_MP_LIMB - t));
- }
- else
- cc = mpn_add_1 (rp, rp, rrsize, 1);
- }
+ case GMP_RNDU:
+ /* we should arrive here only when the result is inexact,
+ i.e. either q_limb > 0 (the remainder from mpn_sqrtrem is non-zero)
+ or up[0..usize-rsize-1] is non zero, thus we have to add one
+ ulp, and inexact = 1 */
+ inexact = 1;
+ t = MPFR_PREC(r) & (BITS_PER_MP_LIMB - 1);
+ rsize = (MPFR_PREC(r) - 1)/BITS_PER_MP_LIMB + 1;
+ if (t)
+ cc = mpn_add_1 (rp + rrsize - rsize, rp + rrsize - rsize, rsize,
+ MP_LIMB_T_ONE << (BITS_PER_MP_LIMB - t));
+ else
+ cc = mpn_add_1 (rp + rrsize - rsize, rp + rrsize - rsize, rsize,
+ MP_LIMB_T_ONE);
}
- if (cc) {
- mpn_rshift(rp, rp, rrsize, 1);
- rp[rrsize-1] |= (mp_limb_t) 1 << (BITS_PER_MP_LIMB-1);
- MPFR_EXP(r)++;
- }
+ if (cc)
+ {
+ /* Is a shift necessary here? Isn't the result 1.0000...? */
+ mpn_rshift (rp, rp, rrsize, 1);
+ rp[rrsize-1] |= MP_LIMB_T_HIGHBIT;
+ MPFR_EXP(r)++;
+ }
fin:
rsize = rrsize;
@@ -254,11 +267,11 @@ mpfr_sqrt (r, u, rnd_mode)
MPN_COPY(MPFR_MANT(r), rp + rsize - rrsize, rrsize);
if (MPFR_PREC(r) & (BITS_PER_MP_LIMB - 1))
- MPFR_MANT(r) [0] &= ~(((mp_limb_t)1 << (BITS_PER_MP_LIMB -
+ MPFR_MANT(r) [0] &= ~((MP_LIMB_T_ONE << (BITS_PER_MP_LIMB -
(MPFR_PREC(r) & (BITS_PER_MP_LIMB - 1)))) - 1) ;
TMP_FREE (marker);
}
TMP_FREE(marker0);
- return exact;
+ return inexact;
}
View Lorry Controller Status