14a1767b4Smrg /* mpz_cdiv_r_2exp, mpz_fdiv_r_2exp -- remainder from mpz divided by 2^n.
24a1767b4Smrg
3d25e02daSmrg Copyright 2001, 2002, 2004, 2012 Free Software Foundation, Inc.
44a1767b4Smrg
54a1767b4Smrg This file is part of the GNU MP Library.
64a1767b4Smrg
74a1767b4Smrg The GNU MP Library is free software; you can redistribute it and/or modify
8*f81b1c5bSmrg it under the terms of either:
9*f81b1c5bSmrg
10*f81b1c5bSmrg * the GNU Lesser General Public License as published by the Free
11*f81b1c5bSmrg Software Foundation; either version 3 of the License, or (at your
124a1767b4Smrg option) any later version.
134a1767b4Smrg
14*f81b1c5bSmrg or
15*f81b1c5bSmrg
16*f81b1c5bSmrg * the GNU General Public License as published by the Free Software
17*f81b1c5bSmrg Foundation; either version 2 of the License, or (at your option) any
18*f81b1c5bSmrg later version.
19*f81b1c5bSmrg
20*f81b1c5bSmrg or both in parallel, as here.
21*f81b1c5bSmrg
224a1767b4Smrg The GNU MP Library is distributed in the hope that it will be useful, but
234a1767b4Smrg WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
24*f81b1c5bSmrg or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
25*f81b1c5bSmrg for more details.
264a1767b4Smrg
27*f81b1c5bSmrg You should have received copies of the GNU General Public License and the
28*f81b1c5bSmrg GNU Lesser General Public License along with the GNU MP Library. If not,
29*f81b1c5bSmrg see https://www.gnu.org/licenses/. */
304a1767b4Smrg
314a1767b4Smrg #include "gmp-impl.h"
324a1767b4Smrg
334a1767b4Smrg
344a1767b4Smrg /* Bit mask of "n" least significant bits of a limb. */
354a1767b4Smrg #define LOW_MASK(n) ((CNST_LIMB(1) << (n)) - 1)
364a1767b4Smrg
374a1767b4Smrg
384a1767b4Smrg /* dir==1 for ceil, dir==-1 for floor */
394a1767b4Smrg
40d25e02daSmrg static void __gmpz_cfdiv_r_2exp (REGPARM_3_1 (mpz_ptr, mpz_srcptr, mp_bitcnt_t, int)) REGPARM_ATTR (1);
414a1767b4Smrg #define cfdiv_r_2exp(w,u,cnt,dir) __gmpz_cfdiv_r_2exp (REGPARM_3_1 (w, u, cnt, dir))
424a1767b4Smrg
434a1767b4Smrg REGPARM_ATTR (1) static void
cfdiv_r_2exp(mpz_ptr w,mpz_srcptr u,mp_bitcnt_t cnt,int dir)444a1767b4Smrg cfdiv_r_2exp (mpz_ptr w, mpz_srcptr u, mp_bitcnt_t cnt, int dir)
454a1767b4Smrg {
464a1767b4Smrg mp_size_t usize, abs_usize, limb_cnt, i;
474a1767b4Smrg mp_srcptr up;
484a1767b4Smrg mp_ptr wp;
494a1767b4Smrg mp_limb_t high;
504a1767b4Smrg
514a1767b4Smrg usize = SIZ(u);
524a1767b4Smrg if (usize == 0)
534a1767b4Smrg {
544a1767b4Smrg SIZ(w) = 0;
554a1767b4Smrg return;
564a1767b4Smrg }
574a1767b4Smrg
584a1767b4Smrg limb_cnt = cnt / GMP_NUMB_BITS;
594a1767b4Smrg cnt %= GMP_NUMB_BITS;
604a1767b4Smrg abs_usize = ABS (usize);
614a1767b4Smrg
624a1767b4Smrg /* MPZ_REALLOC(w) below is only when w!=u, so we can fetch PTR(u) here
634a1767b4Smrg nice and early */
644a1767b4Smrg up = PTR(u);
654a1767b4Smrg
664a1767b4Smrg if ((usize ^ dir) < 0)
674a1767b4Smrg {
684a1767b4Smrg /* Round towards zero, means just truncate */
694a1767b4Smrg
704a1767b4Smrg if (w == u)
714a1767b4Smrg {
724a1767b4Smrg /* if already smaller than limb_cnt then do nothing */
734a1767b4Smrg if (abs_usize <= limb_cnt)
744a1767b4Smrg return;
75*f81b1c5bSmrg wp = (mp_ptr) up;
764a1767b4Smrg }
774a1767b4Smrg else
784a1767b4Smrg {
794a1767b4Smrg i = MIN (abs_usize, limb_cnt+1);
80*f81b1c5bSmrg wp = MPZ_NEWALLOC (w, i);
814a1767b4Smrg MPN_COPY (wp, up, i);
824a1767b4Smrg
834a1767b4Smrg /* if smaller than limb_cnt then only the copy is needed */
844a1767b4Smrg if (abs_usize <= limb_cnt)
854a1767b4Smrg {
864a1767b4Smrg SIZ(w) = usize;
874a1767b4Smrg return;
884a1767b4Smrg }
894a1767b4Smrg }
904a1767b4Smrg }
914a1767b4Smrg else
924a1767b4Smrg {
934a1767b4Smrg /* Round away from zero, means twos complement if non-zero */
944a1767b4Smrg
954a1767b4Smrg /* if u!=0 and smaller than divisor, then must negate */
964a1767b4Smrg if (abs_usize <= limb_cnt)
974a1767b4Smrg goto negate;
984a1767b4Smrg
994a1767b4Smrg /* if non-zero low limb, then must negate */
1004a1767b4Smrg for (i = 0; i < limb_cnt; i++)
1014a1767b4Smrg if (up[i] != 0)
1024a1767b4Smrg goto negate;
1034a1767b4Smrg
1044a1767b4Smrg /* if non-zero partial limb, then must negate */
1054a1767b4Smrg if ((up[limb_cnt] & LOW_MASK (cnt)) != 0)
1064a1767b4Smrg goto negate;
1074a1767b4Smrg
1084a1767b4Smrg /* otherwise low bits of u are zero, so that's the result */
1094a1767b4Smrg SIZ(w) = 0;
1104a1767b4Smrg return;
1114a1767b4Smrg
1124a1767b4Smrg negate:
1134a1767b4Smrg /* twos complement negation to get 2**cnt-u */
1144a1767b4Smrg
115d25e02daSmrg wp = MPZ_REALLOC (w, limb_cnt+1);
1164a1767b4Smrg up = PTR(u);
1174a1767b4Smrg
1184a1767b4Smrg /* Ones complement */
1194a1767b4Smrg i = MIN (abs_usize, limb_cnt+1);
120*f81b1c5bSmrg ASSERT_CARRY (mpn_neg (wp, up, i));
1214a1767b4Smrg for ( ; i <= limb_cnt; i++)
1224a1767b4Smrg wp[i] = GMP_NUMB_MAX;
1234a1767b4Smrg
1244a1767b4Smrg usize = -usize;
1254a1767b4Smrg }
1264a1767b4Smrg
1274a1767b4Smrg /* Mask the high limb */
1284a1767b4Smrg high = wp[limb_cnt];
1294a1767b4Smrg high &= LOW_MASK (cnt);
1304a1767b4Smrg wp[limb_cnt] = high;
1314a1767b4Smrg
1324a1767b4Smrg /* Strip any consequent high zeros */
1334a1767b4Smrg while (high == 0)
1344a1767b4Smrg {
1354a1767b4Smrg limb_cnt--;
1364a1767b4Smrg if (limb_cnt < 0)
1374a1767b4Smrg {
1384a1767b4Smrg SIZ(w) = 0;
1394a1767b4Smrg return;
1404a1767b4Smrg }
1414a1767b4Smrg high = wp[limb_cnt];
1424a1767b4Smrg }
1434a1767b4Smrg
1444a1767b4Smrg limb_cnt++;
1454a1767b4Smrg SIZ(w) = (usize >= 0 ? limb_cnt : -limb_cnt);
1464a1767b4Smrg }
1474a1767b4Smrg
1484a1767b4Smrg
1494a1767b4Smrg void
mpz_cdiv_r_2exp(mpz_ptr w,mpz_srcptr u,mp_bitcnt_t cnt)1504a1767b4Smrg mpz_cdiv_r_2exp (mpz_ptr w, mpz_srcptr u, mp_bitcnt_t cnt)
1514a1767b4Smrg {
1524a1767b4Smrg cfdiv_r_2exp (w, u, cnt, 1);
1534a1767b4Smrg }
1544a1767b4Smrg
1554a1767b4Smrg void
mpz_fdiv_r_2exp(mpz_ptr w,mpz_srcptr u,mp_bitcnt_t cnt)1564a1767b4Smrg mpz_fdiv_r_2exp (mpz_ptr w, mpz_srcptr u, mp_bitcnt_t cnt)
1574a1767b4Smrg {
1584a1767b4Smrg cfdiv_r_2exp (w, u, cnt, -1);
1594a1767b4Smrg }
160