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, 2004 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 3 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. If not, see http://www.gnu.org/licenses/. */ 20 21 #include "gmp.h" 22 #include "gmp-impl.h" 23 24 25 /* dir==1 for ceil, dir==-1 for floor */ 26 27 static void __gmpz_cfdiv_q_2exp __GMP_PROTO ((REGPARM_3_1 (mpz_ptr, mpz_srcptr, mp_bitcnt_t, int))) REGPARM_ATTR (1); 28 #define cfdiv_q_2exp(w,u,cnt,dir) __gmpz_cfdiv_q_2exp (REGPARM_3_1 (w,u,cnt,dir)) 29 30 REGPARM_ATTR (1) static void 31 cfdiv_q_2exp (mpz_ptr w, mpz_srcptr u, mp_bitcnt_t cnt, int dir) 32 { 33 mp_size_t wsize, usize, abs_usize, limb_cnt, i; 34 mp_srcptr up; 35 mp_ptr wp; 36 mp_limb_t round, rmask; 37 38 usize = SIZ (u); 39 abs_usize = ABS (usize); 40 limb_cnt = cnt / GMP_NUMB_BITS; 41 wsize = abs_usize - limb_cnt; 42 if (wsize <= 0) 43 { 44 /* u < 2**cnt, so result 1, 0 or -1 according to rounding */ 45 PTR(w)[0] = 1; 46 SIZ(w) = (usize == 0 || (usize ^ dir) < 0 ? 0 : dir); 47 return; 48 } 49 50 /* +1 limb to allow for mpn_add_1 below */ 51 MPZ_REALLOC (w, wsize+1); 52 53 /* Check for rounding if direction matches u sign. 54 Set round if we're skipping non-zero limbs. */ 55 up = PTR(u); 56 round = 0; 57 rmask = ((usize ^ dir) >= 0 ? MP_LIMB_T_MAX : 0); 58 if (rmask != 0) 59 for (i = 0; i < limb_cnt && round == 0; i++) 60 round = up[i]; 61 62 wp = PTR(w); 63 cnt %= GMP_NUMB_BITS; 64 if (cnt != 0) 65 { 66 round |= rmask & mpn_rshift (wp, up + limb_cnt, wsize, cnt); 67 wsize -= (wp[wsize - 1] == 0); 68 } 69 else 70 MPN_COPY_INCR (wp, up + limb_cnt, wsize); 71 72 if (round != 0) 73 { 74 if (wsize != 0) 75 { 76 mp_limb_t cy; 77 cy = mpn_add_1 (wp, wp, wsize, CNST_LIMB(1)); 78 wp[wsize] = cy; 79 wsize += cy; 80 } 81 else 82 { 83 /* We shifted something to zero. */ 84 wp[0] = 1; 85 wsize = 1; 86 } 87 } 88 SIZ(w) = (usize >= 0 ? wsize : -wsize); 89 } 90 91 92 void 93 mpz_cdiv_q_2exp (mpz_ptr w, mpz_srcptr u, mp_bitcnt_t cnt) 94 { 95 cfdiv_q_2exp (w, u, cnt, 1); 96 } 97 98 void 99 mpz_fdiv_q_2exp (mpz_ptr w, mpz_srcptr u, mp_bitcnt_t cnt) 100 { 101 cfdiv_q_2exp (w, u, cnt, -1); 102 } 103