![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to rint.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr|
Return to rint.c CVS log | Up to [local] / OpenXM_contrib / gmp / mpfr |
1.1 ohara 1: /* mpfr_rint -- Round to an 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 "gmp.h"
23: #include "gmp-impl.h"
24: #include "mpfr.h"
25: #include "mpfr-impl.h"
26:
27: int
28: mpfr_rint (mpfr_ptr r, mpfr_srcptr u, mp_rnd_t rnd_mode)
29: {
30: int sign;
31: mp_exp_t exp;
32:
33: if (MPFR_IS_NAN(u))
34: {
35: MPFR_SET_NAN(r);
36: MPFR_RET_NAN;
37: }
38:
39: MPFR_CLEAR_NAN(r);
40: MPFR_SET_SAME_SIGN(r, u);
41:
42: if (MPFR_IS_INF(u))
43: {
44: MPFR_SET_INF(r);
45: MPFR_RET(0); /* infinity is exact */
46: }
47:
48: MPFR_CLEAR_INF(r);
49:
50: if (MPFR_IS_ZERO(u))
51: {
52: MPFR_SET_ZERO(r);
53: MPFR_RET(0); /* zero is exact */
54: }
55:
56: sign = MPFR_SIGN(u);
57: exp = MPFR_EXP(u);
58:
59: /* Single out the case where |u| < 1. */
60: if (exp <= 0) /* 0 < |u| < 1 */
61: {
62: if ((rnd_mode == GMP_RNDD && sign < 0) ||
63: (rnd_mode == GMP_RNDU && sign > 0) ||
64: (rnd_mode == GMP_RNDN && exp == 0))
65: {
66: mp_limb_t *rp;
67: mp_size_t rm;
68:
69: rp = MPFR_MANT(r);
70: rm = (MPFR_PREC(r) - 1) / BITS_PER_MP_LIMB;
71: rp[rm] = GMP_LIMB_HIGHBIT;
72: MPN_ZERO(rp, rm);
73: MPFR_EXP(r) = 1; /* |r| = 1 */
74: MPFR_RET(sign > 0 ? 2 : -2);
75: }
76: else
77: {
78: MPFR_SET_ZERO(r); /* r = 0 */
79: MPFR_RET(sign > 0 ? -2 : 2);
80: }
81: }
82: else /* exp > 0, |u| >= 1 */
83: {
84: mp_limb_t *up, *rp;
85: mp_size_t un, rn, ui;
86: int sh, idiff;
87: int uflags;
88:
89: /*
90: * uflags will contain:
91: * _ 0 if u is an integer representable in r,
92: * _ 1 if u is an integer not representable in r,
93: * _ 2 if u is not an integer.
94: */
95:
96: up = MPFR_MANT(u);
97: rp = MPFR_MANT(r);
98:
99: un = MPFR_ESIZE(u);
100: rn = MPFR_ESIZE(r);
101: sh = rn * BITS_PER_MP_LIMB - MPFR_PREC(r);
102:
103: MPFR_EXP(r) = exp;
104:
105: if ((exp - 1) / BITS_PER_MP_LIMB >= un)
106: {
107: ui = un;
108: idiff = 0;
109: uflags = 0; /* u is an integer */
110: }
111: else
112: {
113: mp_size_t uj;
114:
115: ui = (exp - 1) / BITS_PER_MP_LIMB + 1; /* #limbs of the int part */
116: uj = un - ui; /* lowest limb of the integer part */
117: idiff = exp % BITS_PER_MP_LIMB; /* #int-part bits in up[uj] or 0 */
118:
119: uflags = idiff == 0 || (up[uj] << idiff) == 0 ? 0 : 2;
120: if (uflags == 0)
121: while (uj > 0)
122: if (up[--uj] != 0)
123: {
124: uflags = 2;
125: break;
126: }
127: }
128:
129: if (ui > rn)
130: {
131: MPFR_ASSERTN(rp != up && un > rn);
132: MPN_COPY(rp, up + (un - rn), rn);
133: /* In the rounding to the nearest mode, if the rounding bit
134: is 0, change the rounding mode to GMP_RNDZ. */
135: if (rnd_mode == GMP_RNDN &&
136: ((sh != 0 && (rp[0] & (MP_LIMB_T_ONE << (sh - 1))) == 0) ||
137: (sh == 0 && (up[un - rn - 1] & GMP_LIMB_HIGHBIT) == 0)))
138: rnd_mode = GMP_RNDZ;
139: if (uflags == 0)
140: { /* u is an integer; determine if it is representable */
141: if (sh != 0 && rp[0] << (BITS_PER_MP_LIMB - sh) != 0)
142: uflags = 1; /* u is not representable in r */
143: else
144: {
145: mp_size_t i;
146: for (i = un - rn - 1; i >= 0; i--)
147: if (up[i] != 0)
148: {
149: uflags = 1; /* u is not representable in r */
150: break;
151: }
152: }
153: }
154: if (sh != 0)
155: rp[0] &= MP_LIMB_T_MAX << sh;
156: }
157: else
158: {
159: mp_size_t uj, rj;
160: int ush;
161:
162: uj = un - ui; /* lowest limb of the integer part in u */
163: rj = rn - ui; /* lowest limb of the integer part in r */
164:
165: MPN_ZERO(rp, rj);
166:
167: if (rp != up)
168: MPN_COPY(rp + rj, up + uj, ui);
169:
170: rp += rj;
171: rn = ui;
172:
173: ush = idiff == 0 ? 0 : BITS_PER_MP_LIMB - idiff;
174: if (rj == 0 && ush < sh)
175: {
176: /* In the rounding to the nearest mode, if the rounding bit
177: is 0, change the rounding mode to GMP_RNDZ. */
178: if (rnd_mode == GMP_RNDN &&
179: (rp[0] & (MP_LIMB_T_ONE << (sh - 1))) == 0)
180: rnd_mode = GMP_RNDZ; /* rounding bit is 0 */
181: if (uflags == 0)
182: { /* u is an integer; determine if it is representable */
183: mp_limb_t mask;
184: mask = (MP_LIMB_T_ONE << sh) - (MP_LIMB_T_ONE << ush);
185: if ((rp[0] & mask) != 0)
186: uflags = 1; /* u is not representable in r */
187: }
188: }
189: else
190: {
191: sh = ush;
192: if (rnd_mode == GMP_RNDN &&
193: ((ush != 0 &&
194: (up[uj] & (MP_LIMB_T_ONE << (ush - 1))) == 0) ||
195: (ush == 0 &&
196: (uj == 0 || (up[uj - 1] & GMP_LIMB_HIGHBIT) == 0))))
197: rnd_mode = GMP_RNDZ; /* rounding bit is 0 */
198: }
199: if (sh != 0)
200: rp[0] &= MP_LIMB_T_MAX << sh;
201: }
202:
203: if (uflags == 0)
204: MPFR_RET(0);
205:
206: /* Note: if rnd_mode == GMP_RNDN, then round away from 0 (if
207: the rounding bit was 0 and rnd_mode == GMP_RNDN, rnd_mode
208: has been changed to GMP_RNDZ). */
209:
210: if ((rnd_mode == GMP_RNDN ||
211: (rnd_mode == GMP_RNDD && sign < 0) ||
212: (rnd_mode == GMP_RNDU && sign > 0))
213: && mpn_add_1(rp, rp, rn, MP_LIMB_T_ONE << sh))
214: {
215: if (exp == __mpfr_emax)
216: return mpfr_set_overflow(r, rnd_mode, MPFR_SIGN(r)) >= 0 ?
217: uflags : -uflags;
218: else
219: {
220: MPFR_EXP(r)++;
221: rp[rn-1] = GMP_LIMB_HIGHBIT;
222: }
223: }
224:
225: MPFR_RET(rnd_mode == GMP_RNDU ||
226: (rnd_mode == GMP_RNDZ && sign < 0) ||
227: (rnd_mode == GMP_RNDN && sign > 0) ? uflags : -uflags);
228: } /* exp > 0, |u| >= 1 */
229: }
230:
231: #undef mpfr_round
232:
233: int
234: mpfr_round (mpfr_ptr r, mpfr_srcptr u)
235: {
236: return mpfr_rint(r, u, GMP_RNDN);
237: }
238:
239: #undef mpfr_trunc
240:
241: int
242: mpfr_trunc (mpfr_ptr r, mpfr_srcptr u)
243: {
244: return mpfr_rint(r, u, GMP_RNDZ);
245: }
246:
247: #undef mpfr_ceil
248:
249: int
250: mpfr_ceil (mpfr_ptr r, mpfr_srcptr u)
251: {
252: return mpfr_rint(r, u, GMP_RNDU);
253: }
254:
255: #undef mpfr_floor
256:
257: int
258: mpfr_floor (mpfr_ptr r, mpfr_srcptr u)
259: {
260: return mpfr_rint(r, u, GMP_RNDD);
261: }