Annotation of OpenXM_contrib/gmp/mpz/fdiv_r_2exp.c, Revision 1.1.1.2
1.1 maekawa 1: /* mpz_fdiv_r_2exp -- Divide a integer by 2**CNT and produce a remainder.
2:
1.1.1.2 ! maekawa 3: Copyright (C) 1991, 1993, 1994, 1995, 1998, 1999 Free Software Foundation,
! 4: Inc.
1.1 maekawa 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
1.1.1.2 ! maekawa 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
1.1 maekawa 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
1.1.1.2 ! maekawa 15: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
1.1 maekawa 16: License for more details.
17:
1.1.1.2 ! maekawa 18: You should have received a copy of the GNU Lesser General Public License
1.1 maekawa 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_r_2exp (mpz_ptr res, mpz_srcptr in, unsigned long int cnt)
29: #else
30: mpz_fdiv_r_2exp (res, in, cnt)
31: mpz_ptr res;
32: mpz_srcptr in;
33: unsigned long int cnt;
34: #endif
35: {
36: mp_size_t in_size = ABS (in->_mp_size);
37: mp_size_t res_size;
38: mp_size_t limb_cnt = cnt / BITS_PER_MP_LIMB;
39: mp_srcptr in_ptr = in->_mp_d;
40:
41: if (in_size > limb_cnt)
42: {
43: /* The input operand is (probably) greater than 2**CNT. */
44: mp_limb_t x;
45:
46: x = in_ptr[limb_cnt] & (((mp_limb_t) 1 << cnt % BITS_PER_MP_LIMB) - 1);
47: if (x != 0)
48: {
49: res_size = limb_cnt + 1;
50: if (res->_mp_alloc < res_size)
51: _mpz_realloc (res, res_size);
52:
53: res->_mp_d[limb_cnt] = x;
54: }
55: else
56: {
57: res_size = limb_cnt;
58: MPN_NORMALIZE (in_ptr, res_size);
59:
60: if (res->_mp_alloc < res_size)
61: _mpz_realloc (res, res_size);
62:
63: limb_cnt = res_size;
64: }
65: }
66: else
67: {
68: /* The input operand is smaller than 2**CNT. We perform a no-op,
1.1.1.2 ! maekawa 69: apart from that we might need to copy IN to RES, and may need
! 70: to round the result. */
1.1 maekawa 71: res_size = in_size;
72: if (res->_mp_alloc < res_size)
73: _mpz_realloc (res, res_size);
74:
75: limb_cnt = res_size;
76: }
77:
78: if (res != in)
79: MPN_COPY (res->_mp_d, in->_mp_d, limb_cnt);
1.1.1.2 ! maekawa 80: in_size = in->_mp_size;
1.1 maekawa 81: res->_mp_size = res_size;
1.1.1.2 ! maekawa 82: if (in_size < 0 && res_size != 0)
1.1 maekawa 83: {
84: /* Result should be 2^CNT - RES */
85: mpz_t tmp;
1.1.1.2 ! maekawa 86: MPZ_TMP_INIT (tmp, cnt/BITS_PER_MP_LIMB + 2);
1.1 maekawa 87: mpz_set_ui (tmp, 1L);
88: mpz_mul_2exp (tmp, tmp, cnt);
89: mpz_sub (res, tmp, res);
90: }
91: }
1.1.1.2 ! maekawa 92:
! 93: /* This is an alternative ending of the above function using just low-level
! 94: functions. Tested, but perhaps excessive? */
! 95: #if 0
! 96: if (in->_mp_size < 0 && res_size != 0)
! 97: {
! 98: /* Result should be 2^CNT - RES */
! 99:
! 100: mp_ptr rp;
! 101:
! 102: limb_cnt = cnt / BITS_PER_MP_LIMB;
! 103:
! 104: if (res->_mp_alloc <= limb_cnt)
! 105: _mpz_realloc (res, limb_cnt + 1);
! 106: rp = PTR(res);
! 107: if (res_size > limb_cnt)
! 108: {
! 109: mpn_nz_neg (rp, rp, res_size);
! 110: rp[limb_cnt] &= ~(~(mp_limb_t) 0 << cnt % BITS_PER_MP_LIMB);
! 111: MPN_NORMALIZE_NOT_ZERO (rp, res_size);
! 112: }
! 113: else
! 114: {
! 115: mp_size_t i;
! 116: mpn_nz_neg (rp, rp, res_size);
! 117: for (i = res_size; i < limb_cnt; i++)
! 118: rp[i] = ~ (mp_limb_t) 0;
! 119: res_size = limb_cnt;
! 120: if (cnt % BITS_PER_MP_LIMB != 0)
! 121: {
! 122: rp[res_size] = ((mp_limb_t) 1 << (cnt % BITS_PER_MP_LIMB)) - 1;
! 123: res_size++;
! 124: }
! 125: else
! 126: MPN_NORMALIZE_NOT_ZERO (rp, res_size);
! 127: }
! 128: }
! 129: SIZ(res) = res_size;
! 130: }
! 131:
! 132: static void
! 133: mpn_nz_neg (rp, sp, n)
! 134: mp_ptr rp, sp;
! 135: mp_size_t n;
! 136: {
! 137: mp_size_t i;
! 138: mp_limb_t x;
! 139:
! 140: x = sp[0];
! 141: rp[0] = -x;
! 142: for (i = 1; x == 0; i++)
! 143: {
! 144: x = sp[i];
! 145: rp[i] = -x;
! 146: }
! 147:
! 148: for (; i < n; i++)
! 149: {
! 150: rp[i] = ~sp[i];
! 151: }
! 152: }
! 153: #endif
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>