Annotation of OpenXM_contrib/gmp/mpz/cfdiv_q_2exp.c, Revision 1.1.1.1
1.1 ohara 1: /* mpz_cdiv_q_2exp, mpz_fdiv_q_2exp -- quotient from mpz divided by 2^n.
2:
3: Copyright 1991, 1993, 1994, 1996, 1998, 1999, 2001, 2002 Free Software
4: 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 Lesser General Public License as published by
10: the Free Software Foundation; either version 2.1 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 Lesser General Public
16: License for more details.
17:
18: You should have received a copy of the GNU Lesser 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:
27: /* dir==1 for ceil, dir==-1 for floor */
28:
29: static void __gmpz_cfdiv_q_2exp _PROTO ((REGPARM_3_1 (mpz_ptr w, mpz_srcptr u, unsigned long cnt, int dir))) REGPARM_ATTR (1);
30: #define cfdiv_q_2exp(w,u,cnt,dir) __gmpz_cfdiv_q_2exp (REGPARM_3_1 (w,u,cnt,dir))
31:
32: static void
33: cfdiv_q_2exp (mpz_ptr w, mpz_srcptr u, unsigned long cnt, int dir)
34: {
35: mp_size_t wsize, usize, abs_usize, limb_cnt, i;
36: mp_srcptr up;
37: mp_ptr wp;
38: mp_limb_t round, rmask;
39:
40: usize = SIZ (u);
41: abs_usize = ABS (usize);
42: limb_cnt = cnt / GMP_NUMB_BITS;
43: wsize = abs_usize - limb_cnt;
44: if (wsize <= 0)
45: {
46: /* u < 2**cnt, so result 1, 0 or -1 according to rounding */
47: PTR(w)[0] = 1;
48: SIZ(w) = (usize == 0 || (usize ^ dir) < 0 ? 0 : dir);
49: return;
50: }
51:
52: /* +1 limb to allow for mpn_add_1 below */
53: MPZ_REALLOC (w, wsize+1);
54:
55: /* Check for rounding if direction matches u sign.
56: Set round if we're skipping non-zero limbs. */
57: up = PTR(u);
58: round = 0;
59: rmask = ((usize ^ dir) >= 0 ? MP_LIMB_T_MAX : 0);
60: if (rmask != 0)
61: for (i = 0; i < limb_cnt && round == 0; i++)
62: round = up[i];
63:
64: wp = PTR(w);
65: cnt %= GMP_NUMB_BITS;
66: if (cnt != 0)
67: {
68: round |= rmask & mpn_rshift (wp, up + limb_cnt, wsize, cnt);
69: wsize -= (wp[wsize - 1] == 0);
70: }
71: else
72: MPN_COPY_INCR (wp, up + limb_cnt, wsize);
73:
74: if (round != 0)
75: {
76: if (wsize != 0)
77: {
78: mp_limb_t cy;
79: cy = mpn_add_1 (wp, wp, wsize, CNST_LIMB(1));
80: wp[wsize] = cy;
81: wsize += cy;
82: }
83: else
84: {
85: /* We shifted something to zero. */
86: wp[0] = 1;
87: wsize = 1;
88: }
89: }
90: SIZ(w) = (usize >= 0 ? wsize : -wsize);
91: }
92:
93:
94: void
95: mpz_cdiv_q_2exp (mpz_ptr w, mpz_srcptr u, unsigned long cnt)
96: {
97: cfdiv_q_2exp (w, u, cnt, 1);
98: }
99:
100: void
101: mpz_fdiv_q_2exp (mpz_ptr w, mpz_srcptr u, unsigned long cnt)
102: {
103: cfdiv_q_2exp (w, u, cnt, -1);
104: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>