1/* mpz_cdiv_r_2exp, mpz_fdiv_r_2exp -- remainder from mpz divided by 2^n. 2 3Copyright 2001, 2002, 2004 Free Software Foundation, Inc. 4 5This file is part of the GNU MP Library. 6 7The GNU MP Library is free software; you can redistribute it and/or modify 8it under the terms of the GNU Lesser General Public License as published by 9the Free Software Foundation; either version 3 of the License, or (at your 10option) any later version. 11 12The GNU MP Library is distributed in the hope that it will be useful, but 13WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 14or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 15License for more details. 16 17You should have received a copy of the GNU Lesser General Public License 18along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ 19 20#include "gmp.h" 21#include "gmp-impl.h" 22 23 24/* Bit mask of "n" least significant bits of a limb. */ 25#define LOW_MASK(n) ((CNST_LIMB(1) << (n)) - 1) 26 27 28/* dir==1 for ceil, dir==-1 for floor */ 29 30static void __gmpz_cfdiv_r_2exp __GMP_PROTO ((REGPARM_3_1 (mpz_ptr, mpz_srcptr, mp_bitcnt_t, int))) REGPARM_ATTR (1); 31#define cfdiv_r_2exp(w,u,cnt,dir) __gmpz_cfdiv_r_2exp (REGPARM_3_1 (w, u, cnt, dir)) 32 33REGPARM_ATTR (1) static void 34cfdiv_r_2exp (mpz_ptr w, mpz_srcptr u, mp_bitcnt_t cnt, int dir) 35{ 36 mp_size_t usize, abs_usize, limb_cnt, i; 37 mp_srcptr up; 38 mp_ptr wp; 39 mp_limb_t high; 40 41 usize = SIZ(u); 42 if (usize == 0) 43 { 44 SIZ(w) = 0; 45 return; 46 } 47 48 limb_cnt = cnt / GMP_NUMB_BITS; 49 cnt %= GMP_NUMB_BITS; 50 abs_usize = ABS (usize); 51 52 /* MPZ_REALLOC(w) below is only when w!=u, so we can fetch PTR(u) here 53 nice and early */ 54 up = PTR(u); 55 56 if ((usize ^ dir) < 0) 57 { 58 /* Round towards zero, means just truncate */ 59 60 if (w == u) 61 { 62 /* if already smaller than limb_cnt then do nothing */ 63 if (abs_usize <= limb_cnt) 64 return; 65 wp = PTR(w); 66 } 67 else 68 { 69 i = MIN (abs_usize, limb_cnt+1); 70 MPZ_REALLOC (w, i); 71 wp = PTR(w); 72 MPN_COPY (wp, up, i); 73 74 /* if smaller than limb_cnt then only the copy is needed */ 75 if (abs_usize <= limb_cnt) 76 { 77 SIZ(w) = usize; 78 return; 79 } 80 } 81 } 82 else 83 { 84 /* Round away from zero, means twos complement if non-zero */ 85 86 /* if u!=0 and smaller than divisor, then must negate */ 87 if (abs_usize <= limb_cnt) 88 goto negate; 89 90 /* if non-zero low limb, then must negate */ 91 for (i = 0; i < limb_cnt; i++) 92 if (up[i] != 0) 93 goto negate; 94 95 /* if non-zero partial limb, then must negate */ 96 if ((up[limb_cnt] & LOW_MASK (cnt)) != 0) 97 goto negate; 98 99 /* otherwise low bits of u are zero, so that's the result */ 100 SIZ(w) = 0; 101 return; 102 103 negate: 104 /* twos complement negation to get 2**cnt-u */ 105 106 MPZ_REALLOC (w, limb_cnt+1); 107 up = PTR(u); 108 wp = PTR(w); 109 110 /* Ones complement */ 111 i = MIN (abs_usize, limb_cnt+1); 112 mpn_com (wp, up, i); 113 for ( ; i <= limb_cnt; i++) 114 wp[i] = GMP_NUMB_MAX; 115 116 /* Twos complement. Since u!=0 in the relevant part, the twos 117 complement never gives 0 and a carry, so can use MPN_INCR_U. */ 118 MPN_INCR_U (wp, limb_cnt+1, CNST_LIMB(1)); 119 120 usize = -usize; 121 } 122 123 /* Mask the high limb */ 124 high = wp[limb_cnt]; 125 high &= LOW_MASK (cnt); 126 wp[limb_cnt] = high; 127 128 /* Strip any consequent high zeros */ 129 while (high == 0) 130 { 131 limb_cnt--; 132 if (limb_cnt < 0) 133 { 134 SIZ(w) = 0; 135 return; 136 } 137 high = wp[limb_cnt]; 138 } 139 140 limb_cnt++; 141 SIZ(w) = (usize >= 0 ? limb_cnt : -limb_cnt); 142} 143 144 145void 146mpz_cdiv_r_2exp (mpz_ptr w, mpz_srcptr u, mp_bitcnt_t cnt) 147{ 148 cfdiv_r_2exp (w, u, cnt, 1); 149} 150 151void 152mpz_fdiv_r_2exp (mpz_ptr w, mpz_srcptr u, mp_bitcnt_t cnt) 153{ 154 cfdiv_r_2exp (w, u, cnt, -1); 155} 156