contrib/gmp/mpz/cfdiv_r_2exp.c

   1 /* mpz_cdiv_r_2exp, mpz_fdiv_r_2exp -- remainder from mpz divided by 2^n.
   2
   3 Copyright 2001, 2002, 2004 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 3 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.  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
  30 static void __gmpz_cfdiv_r_2exp __GMP_PROTO ((REGPARM_3_1 (mpz_ptr w, mpz_srcptr u, unsigned long cnt, int dir))) REGPARM_ATTR (1);
  31 #define cfdiv_r_2exp(w,u,cnt,dir)  __gmpz_cfdiv_r_2exp (REGPARM_3_1 (w, u, cnt, dir))
  32
  33 REGPARM_ATTR (1) static void
  34 cfdiv_r_2exp (mpz_ptr w, mpz_srcptr u, unsigned long 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_n (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
 145 void
 146 mpz_cdiv_r_2exp (mpz_ptr w, mpz_srcptr u, unsigned long cnt)
 147 {
 148   cfdiv_r_2exp (w, u, cnt, 1);
 149 }
 150
 151 void
 152 mpz_fdiv_r_2exp (mpz_ptr w, mpz_srcptr u, unsigned long cnt)
 153 {
 154   cfdiv_r_2exp (w, u, cnt, -1);
 155 }