cfdiv_r_2exp.c
上传用户:qaz666999
上传日期:2022-08-06
资源大小:2570k
文件大小:4k
源码类别:

数学计算

开发平台:

Unix_Linux

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