Annotation of OpenXM_contrib/gmp/mpz/cfdiv_r_2exp.c, Revision 1.1.1.1
1.1 ohara 1: /* mpz_cdiv_r_2exp, mpz_fdiv_r_2exp -- remainder from mpz divided by 2^n.
2:
3: Copyright 2001, 2002 Free Software Foundation, Inc.
4:
5: This file is part of the GNU MP Library.
6:
7: The GNU MP 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 GNU MP 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 GNU MP 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:
25:
26: /* Bit mask of "n" least significant bits of a limb. */
27: #define LOW_MASK(n) ((CNST_LIMB(1) << (n)) - 1)
28:
29:
30: /* dir==1 for ceil, dir==-1 for floor */
31:
32: static void __gmpz_cfdiv_r_2exp _PROTO ((REGPARM_3_1 (mpz_ptr w, mpz_srcptr u, unsigned long cnt, int dir))) REGPARM_ATTR (1);
33: #define cfdiv_r_2exp(w,u,cnt,dir) __gmpz_cfdiv_r_2exp (REGPARM_3_1 (w, u, cnt, dir))
34:
35: static void
36: cfdiv_r_2exp (mpz_ptr w, mpz_srcptr u, unsigned long cnt, int dir)
37: {
38: mp_size_t usize, abs_usize, limb_cnt, i;
39: mp_srcptr up;
40: mp_ptr wp;
41: mp_limb_t high;
42:
43: usize = SIZ(u);
44: if (usize == 0)
45: {
46: SIZ(w) = 0;
47: return;
48: }
49:
50: limb_cnt = cnt / GMP_NUMB_BITS;
51: cnt %= GMP_NUMB_BITS;
52: abs_usize = ABS (usize);
53:
54: /* MPZ_REALLOC(w) below is only when w!=u, so we can fetch PTR(u) here
55: nice and early */
56: up = PTR(u);
57:
58: if ((usize ^ dir) < 0)
59: {
60: /* Round towards zero, means just truncate */
61:
62: if (w == u)
63: {
64: /* if already smaller than limb_cnt then do nothing */
65: if (abs_usize <= limb_cnt)
66: return;
67: wp = PTR(w);
68: }
69: else
70: {
71: i = MIN (abs_usize, limb_cnt+1);
72: MPZ_REALLOC (w, i);
73: wp = PTR(w);
74: MPN_COPY (wp, up, i);
75:
76: /* if smaller than limb_cnt then only the copy is needed */
77: if (abs_usize <= limb_cnt)
78: {
79: SIZ(w) = usize;
80: return;
81: }
82: }
83: }
84: else
85: {
86: /* Round away from zero, means twos complement if non-zero */
87:
88: /* if u!=0 and smaller than divisor, then must negate */
89: if (abs_usize <= limb_cnt)
90: goto negate;
91:
92: /* if non-zero low limb, then must negate */
93: for (i = 0; i < limb_cnt; i++)
94: if (up[i] != 0)
95: goto negate;
96:
97: /* if non-zero partial limb, then must negate */
98: if ((up[limb_cnt] & LOW_MASK (cnt)) != 0)
99: goto negate;
100:
101: /* otherwise low bits of u are zero, so that's the result */
102: SIZ(w) = 0;
103: return;
104:
105: negate:
106: /* twos complement negation to get 2**cnt-u */
107:
108: MPZ_REALLOC (w, limb_cnt+1);
109: up = PTR(u);
110: wp = PTR(w);
111:
112: /* Ones complement */
113: i = MIN (abs_usize, limb_cnt+1);
114: mpn_com_n (wp, up, i);
115: for ( ; i <= limb_cnt; i++)
116: wp[i] = GMP_NUMB_MAX;
117:
118: /* Twos complement. Since u!=0 in the relevant part, the twos
119: complement never gives 0 and a carry, so can use MPN_INCR_U. */
120: MPN_INCR_U (wp, limb_cnt+1, CNST_LIMB(1));
121:
122: usize = -usize;
123: }
124:
125: /* Mask the high limb */
126: high = wp[limb_cnt];
127: high &= LOW_MASK (cnt);
128: wp[limb_cnt] = high;
129:
130: /* Strip any consequent high zeros */
131: while (high == 0)
132: {
133: limb_cnt--;
134: if (limb_cnt < 0)
135: {
136: SIZ(w) = 0;
137: return;
138: }
139: high = wp[limb_cnt];
140: }
141:
142: limb_cnt++;
143: SIZ(w) = (usize >= 0 ? limb_cnt : -limb_cnt);
144: }
145:
146:
147: void
148: mpz_cdiv_r_2exp (mpz_ptr w, mpz_srcptr u, unsigned long cnt)
149: {
150: cfdiv_r_2exp (w, u, cnt, 1);
151: }
152:
153: void
154: mpz_fdiv_r_2exp (mpz_ptr w, mpz_srcptr u, unsigned long cnt)
155: {
156: cfdiv_r_2exp (w, u, cnt, -1);
157: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>