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