← Back to branch summary

~ubuntu-branches/ubuntu/intrepid/ecl/intrepid

« back to all changes in this revision

Viewing changes to src/gmp/mpfr/div_ui.c

  • Committer: Bazaar Package Importer
  • Author(s): Peter Van Eynde
  • Date: 2006-05-17 02:46:26 UTC
  • Revision ID: james.westby@ubuntu.com-20060517024626-lljr08ftv9g9vefl
Tags: upstream-0.9h-20060510
Import upstream version 0.9h-20060510

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/gmp/mpfr/div_ui.c
 
1
/* mpfr_div_ui -- divide a floating-point number by a machine integer
 
2
 
 
3
Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
 
4
 
 
5
This file is part of the MPFR Library.
 
6
 
 
7
The MPFR Library is free software; you can redistribute it and/or modify
 
8
it under the terms of the GNU Lesser General Public License as published by
 
9
the Free Software Foundation; either version 2.1 of the License, or (at your
 
10
option) any later version.
 
11
 
 
12
The MPFR Library is distributed in the hope that it will be useful, but
 
13
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 
14
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 
15
License for more details.
 
16
 
 
17
You should have received a copy of the GNU Lesser General Public License
 
18
along with the MPFR Library; see the file COPYING.LIB.  If not, write to
 
19
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 
20
MA 02111-1307, USA. */
 
21
 
 
22
#include <stdio.h>
 
23
#include "gmp.h"
 
24
#include "gmp-impl.h"
 
25
#include "longlong.h"
 
26
#include "mpfr.h"
 
27
#include "mpfr-impl.h"
 
28
 
 
29
/* returns 0 if result exact, non-zero otherwise */
 
30
int
 
31
mpfr_div_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mp_rnd_t rnd_mode)
 
32
{
 
33
  long int xn, yn, dif, sh, i;
 
34
  mp_limb_t *xp, *yp, *tmp, c, d;
 
35
  int inexact, middle = 1;
 
36
  TMP_DECL(marker);
 
37
 
 
38
  if (MPFR_IS_NAN(x))
 
39
    {
 
40
      MPFR_SET_NAN(y);
 
41
      MPFR_RET_NAN;
 
42
    }
 
43
 
 
44
  MPFR_CLEAR_NAN(y); /* clear NaN flag */
 
45
 
 
46
  if (MPFR_IS_INF(x)) 
 
47
    { 
 
48
      MPFR_SET_INF(y);
 
49
      MPFR_SET_SAME_SIGN(y, x);
 
50
      MPFR_RET(0);
 
51
    }
 
52
 
 
53
  if (u == 0)
 
54
    {
 
55
      if (MPFR_IS_ZERO(x)) /* 0/0 is NaN */
 
56
        {
 
57
          MPFR_SET_NAN(y);
 
58
          MPFR_RET_NAN;
 
59
        }
 
60
      else /* x/0 is Inf */
 
61
        {
 
62
          MPFR_SET_INF(y);
 
63
          MPFR_SET_SAME_SIGN(y, x);
 
64
          MPFR_RET(0);
 
65
        }
 
66
    }
 
67
 
 
68
  MPFR_CLEAR_INF(y);
 
69
  MPFR_SET_SAME_SIGN(y, x);
 
70
 
 
71
  if (MPFR_IS_ZERO(x))
 
72
    {
 
73
      MPFR_SET_ZERO(y);
 
74
      MPFR_RET(0);
 
75
    }
 
76
 
 
77
  TMP_MARK(marker);
 
78
  xn = (MPFR_PREC(x) - 1)/BITS_PER_MP_LIMB + 1;
 
79
  yn = (MPFR_PREC(y) - 1)/BITS_PER_MP_LIMB + 1;
 
80
 
 
81
  xp = MPFR_MANT(x);
 
82
  yp = MPFR_MANT(y);
 
83
  MPFR_EXP(y) = MPFR_EXP(x);
 
84
 
 
85
  dif = yn + 1 - xn;
 
86
 
 
87
  /* we need to store yn+1 = xn + dif limbs of the quotient */
 
88
  /* don't use tmp=yp since the mpn_lshift call below requires yp >= tmp+1 */
 
89
  tmp = TMP_ALLOC((yn + 1) * BYTES_PER_MP_LIMB);
 
90
 
 
91
  c = (mp_limb_t) u;
 
92
  if (dif >= 0)
 
93
    c = mpn_divrem_1 (tmp, dif, xp, xn, c); /* used all the dividend */
 
94
  else /* dif < 0 i.e. xn > yn, don't use the (-dif) low limbs from x */
 
95
    c = mpn_divrem_1 (tmp, 0, xp - dif, yn + 1, c);
 
96
 
 
97
  inexact = (c != 0);
 
98
  if (rnd_mode == GMP_RNDN)
 
99
    {
 
100
      if (2 * c < (mp_limb_t) u)
 
101
        middle = -1;
 
102
      else if (2 * c > (mp_limb_t) u)
 
103
        middle = 1;
 
104
      else
 
105
        middle = 0; /* exactly in the middle */
 
106
    }
 
107
  for (i=0; ((inexact == 0) || (middle == 0)) && (i < -dif); i++)
 
108
    if (xp[i])
 
109
      inexact = middle = 1; /* larger than middle */
 
110
 
 
111
  if (tmp[yn] == 0) /* high limb is zero */
 
112
    {
 
113
      tmp--;
 
114
      sh = 0;
 
115
      MPFR_EXP(y) -= BITS_PER_MP_LIMB;
 
116
    }
 
117
 
 
118
  /* now we have yn limbs starting from tmp[1], with tmp[yn]<>0 */
 
119
 
 
120
  /* shift left to normalize */
 
121
  count_leading_zeros(sh, tmp[yn]);
 
122
  if (sh)
 
123
    {
 
124
      mpn_lshift (yp, tmp + 1, yn, sh);
 
125
      yp[0] += tmp[0] >> (BITS_PER_MP_LIMB - sh);
 
126
      middle = middle || ((tmp[0] << sh) != 0);
 
127
      inexact = inexact || ((tmp[0] << sh) != 0);
 
128
      MPFR_EXP(y) -= sh; 
 
129
    }
 
130
  else
 
131
    MPN_COPY(yp, tmp + 1, yn);
 
132
 
 
133
  sh = yn * BITS_PER_MP_LIMB - MPFR_PREC(y);
 
134
  /* it remains sh bits in less significant limb of y */
 
135
 
 
136
  d = *yp & ((MP_LIMB_T_ONE << sh) - MP_LIMB_T_ONE);
 
137
  *yp ^= d; /* set to zero lowest sh bits */
 
138
 
 
139
  TMP_FREE(marker);
 
140
  if ((d == 0) && (inexact == 0))
 
141
    return 0; /* result is exact */
 
142
 
 
143
  switch (rnd_mode)
 
144
    {
 
145
    case GMP_RNDZ:
 
146
      MPFR_RET(-MPFR_SIGN(x)); /* result is inexact */
 
147
 
 
148
    case GMP_RNDU:
 
149
      if (MPFR_SIGN(y) > 0)
 
150
        mpfr_add_one_ulp (y, rnd_mode);
 
151
      MPFR_RET(1); /* result is inexact */
 
152
 
 
153
    case GMP_RNDD:
 
154
      if (MPFR_SIGN(y) < 0)
 
155
        mpfr_add_one_ulp (y, rnd_mode);
 
156
      MPFR_RET(-1); /* result is inexact */
 
157
 
 
158
    case GMP_RNDN:
 
159
      if (sh && d < (MP_LIMB_T_ONE << (sh - 1)))
 
160
        MPFR_RET(-MPFR_SIGN(x));
 
161
      else if (sh && d > (MP_LIMB_T_ONE << (sh - 1)))
 
162
        {
 
163
          mpfr_add_one_ulp (y, rnd_mode);
 
164
          MPFR_RET(MPFR_SIGN(x));
 
165
        }
 
166
    else /* sh = 0 or d = 1 << (sh-1) */
 
167
      {
 
168
        /* we are in the middle if:
 
169
           (a) sh > 0 and inexact == 0
 
170
           (b) sh=0 and middle=1
 
171
         */
 
172
        if ((sh && inexact) || (!sh && (middle > 0)) || (*yp & (MP_LIMB_T_ONE << sh)))
 
173
          {
 
174
            mpfr_add_one_ulp (y, rnd_mode);
 
175
            MPFR_RET(MPFR_SIGN(x));
 
176
          }
 
177
            else
 
178
              MPFR_RET(-MPFR_SIGN(x));
 
179
      }
 
180
    }
 
181
  MPFR_RET(inexact); /* should never go here */
 
182
}
 
183
 
 
184
 
 
185