summaryrefslogtreecommitdiff
path: root/round_p.c
blob: 85e6e3005d13444d73aae993eaf3c0926773f1d6 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
/* mpfr_round_p -- check if an approximation is roundable.

Copyright 2005 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 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 Lesser General Public
License for more details.

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., 51 Franklin Place, Fifth Floor, Boston,
MA 02110-1301, USA. */

#include "mpfr-impl.h"

/* Check against mpfr_can_round ? */
#ifdef WANT_ASSERT
# if WANT_ASSERT >= 2
int mpfr_round_p_2 (mp_limb_t *, mp_size_t, mp_exp_t, mp_prec_t);
int
mpfr_round_p (mp_limb_t *bp, mp_size_t bn, mp_exp_t err0, mp_prec_t prec)
{
  int i1, i2;

  i1 = mpfr_round_p_2 (bp, bn, err0, prec);
  i2 = mpfr_can_round_raw (bp, bn, MPFR_SIGN_POS, err0,
                           GMP_RNDN, GMP_RNDZ, prec);
  if (i1 != i2)
    {
      fprintf (stderr, "mpfr_round_p(%d) != mpfr_can_round(%d)!\n"
              "bn = %ld, err0 = %ld, prec = %ld\nbp = ",
              i1, i2, bn, err0, prec);
      gmp_fprintf (stderr, "%NX\n", bp, bn);
      MPFR_ASSERTN (0);
    }
  return i1;
}
# define mpfr_round_p mpfr_round_p_2
# endif
#endif

/*
 * Assuming {bp, bn} is an approximation of a non-singular number
 * with error at most equal to 2^(EXP(b)-err) (`err' bits of b are known)
 * of direction unknown, check if we can round b to zero with precision prec.
 */
int
mpfr_round_p (mp_limb_t *bp, mp_size_t bn, mp_exp_t err0, mp_prec_t prec)
{
  mp_prec_t err;
  mp_size_t k, n;
  mp_limb_t tmp, mask;
  int s;

  err = (mp_prec_t) bn * BITS_PER_MP_LIMB;
  if (MPFR_UNLIKELY (err0 <= 0 || (mpfr_uexp_t) err0 <= prec || prec >= err))
    return 0;  /* can't round */
  err = MIN (err, (mpfr_uexp_t) err0);

  k = prec / BITS_PER_MP_LIMB;
  s = BITS_PER_MP_LIMB - prec%BITS_PER_MP_LIMB;
  n = err / BITS_PER_MP_LIMB - k;

  MPFR_ASSERTD (n >= 0);
  MPFR_ASSERTD (bn > k);

  /* Check first limb */
  bp += bn-1-k;
  tmp = *bp--;
  mask = s == BITS_PER_MP_LIMB ? MP_LIMB_T_MAX : MPFR_LIMB_MASK (s);
  tmp &= mask;

  if (MPFR_LIKELY (n == 0))
    {
      /* prec and error are in the same limb */
      s = BITS_PER_MP_LIMB - err % BITS_PER_MP_LIMB;
      MPFR_ASSERTD (s < BITS_PER_MP_LIMB);
      tmp  >>= s;
      mask >>= s;
      return tmp != 0 && tmp != mask;
    }
  else if (MPFR_UNLIKELY (tmp == 0))
    {
      /* Check if all (n-1) limbs are 0 */
      while (--n)
        if (*bp-- != 0)
          return 1;
      /* Check if final error limb is 0 */
      s = BITS_PER_MP_LIMB - err % BITS_PER_MP_LIMB;
      if (s == BITS_PER_MP_LIMB)
        return 0;
      tmp = *bp >> s;
      return tmp != 0;
    }
  else if (MPFR_UNLIKELY (tmp == mask))
    {
      /* Check if all (n-1) limbs are 11111111111111111 */
      while (--n)
        if (*bp-- != MP_LIMB_T_MAX)
          return 1;
      /* Check if final error limb is 0 */
      s = BITS_PER_MP_LIMB - err % BITS_PER_MP_LIMB;
      if (s == BITS_PER_MP_LIMB)
        return 0;
      tmp = *bp >> s;
      return tmp != (MP_LIMB_T_MAX >> s);
    }
  else
    {
      /* First limb is different from 000000 or 1111111 */
      return 1;
    }
}