Annotation of OpenXM_contrib/gmp/mpz/fdiv_q_2exp.c, Revision 1.1.1.1
1.1 maekawa 1: /* mpz_fdiv_q_2exp -- Divide an integer by 2**CNT. Round the quotient
2: towards -infinity.
3:
4: Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc.
5:
6: This file is part of the GNU MP Library.
7:
8: The GNU MP Library is free software; you can redistribute it and/or modify
9: it under the terms of the GNU Library General Public License as published by
10: the Free Software Foundation; either version 2 of the License, or (at your
11: option) any later version.
12:
13: The GNU MP 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 Library General Public
16: License for more details.
17:
18: You should have received a copy of the GNU Library General Public License
19: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
20: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
21: MA 02111-1307, USA. */
22:
23: #include "gmp.h"
24: #include "gmp-impl.h"
25:
26: void
27: #if __STDC__
28: mpz_fdiv_q_2exp (mpz_ptr w, mpz_srcptr u, unsigned long int cnt)
29: #else
30: mpz_fdiv_q_2exp (w, u, cnt)
31: mpz_ptr w;
32: mpz_srcptr u;
33: unsigned long int cnt;
34: #endif
35: {
36: mp_size_t usize = u->_mp_size;
37: mp_size_t wsize;
38: mp_size_t abs_usize = ABS (usize);
39: mp_size_t limb_cnt;
40: mp_ptr wp;
41: mp_limb_t round = 0;
42:
43: limb_cnt = cnt / BITS_PER_MP_LIMB;
44: wsize = abs_usize - limb_cnt;
45: if (wsize <= 0)
46: {
47: wp = w->_mp_d;
48: wsize = 0;
49: /* Set ROUND since we know we skip some non-zero words in this case.
50: Well, if U is zero, we don't, but then this will be taken care of
51: below, since rounding only really takes place for negative U. */
52: round = 1;
53: wp[0] = 1;
54: w->_mp_size = -(usize < 0);
55: return;
56: }
57: else
58: {
59: mp_size_t i;
60: mp_ptr up;
61:
62: /* Make sure there is enough space. We make an extra limb
63: here to account for possible rounding at the end. */
64: if (w->_mp_alloc < wsize + 1)
65: _mpz_realloc (w, wsize + 1);
66:
67: wp = w->_mp_d;
68: up = u->_mp_d;
69:
70: /* Set ROUND if we are about skip some non-zero limbs. */
71: for (i = 0; i < limb_cnt && round == 0; i++)
72: round = up[i];
73:
74: cnt %= BITS_PER_MP_LIMB;
75: if (cnt != 0)
76: {
77: round |= mpn_rshift (wp, up + limb_cnt, wsize, cnt);
78: wsize -= wp[wsize - 1] == 0;
79: }
80: else
81: {
82: MPN_COPY_INCR (wp, up + limb_cnt, wsize);
83: }
84: }
85:
86: if (usize < 0 && round != 0)
87: {
88: mp_limb_t cy;
89: cy = mpn_add_1 (wp, wp, wsize, 1);
90: wp[wsize] = cy;
91: wsize += cy;
92: }
93: w->_mp_size = usize >= 0 ? wsize : -wsize;
94: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>