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