Return to cfdiv_r_2exp.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpz

Annotation of OpenXM_contrib/gmp/mpz/cfdiv_r_2exp.c, Revision 1.1

1.1     ! ohara       1: /* mpz_cdiv_r_2exp, mpz_fdiv_r_2exp -- remainder from mpz divided by 2^n.
        !             2:
        !             3: Copyright 2001, 2002 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 2.1 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; see the file COPYING.LIB.  If not, write to
        !            19: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
        !            20: MA 02111-1307, USA. */
        !            21:
        !            22: #include "gmp.h"
        !            23: #include "gmp-impl.h"
        !            24:
        !            25:
        !            26: /* Bit mask of "n" least significant bits of a limb. */
        !            27: #define LOW_MASK(n)   ((CNST_LIMB(1) << (n)) - 1)
        !            28:
        !            29:
        !            30: /* dir==1 for ceil, dir==-1 for floor */
        !            31:
        !            32: static void __gmpz_cfdiv_r_2exp _PROTO ((REGPARM_3_1 (mpz_ptr w, mpz_srcptr u, unsigned long cnt, int dir))) REGPARM_ATTR (1);
        !            33: #define cfdiv_r_2exp(w,u,cnt,dir)  __gmpz_cfdiv_r_2exp (REGPARM_3_1 (w, u, cnt, dir))
        !            34:
        !            35: static void
        !            36: cfdiv_r_2exp (mpz_ptr w, mpz_srcptr u, unsigned long cnt, int dir)
        !            37: {
        !            38:   mp_size_t  usize, abs_usize, limb_cnt, i;
        !            39:   mp_srcptr  up;
        !            40:   mp_ptr     wp;
        !            41:   mp_limb_t  high;
        !            42:
        !            43:   usize = SIZ(u);
        !            44:   if (usize == 0)
        !            45:     {
        !            46:       SIZ(w) = 0;
        !            47:       return;
        !            48:     }
        !            49:
        !            50:   limb_cnt = cnt / GMP_NUMB_BITS;
        !            51:   cnt %= GMP_NUMB_BITS;
        !            52:   abs_usize = ABS (usize);
        !            53:
        !            54:   /* MPZ_REALLOC(w) below is only when w!=u, so we can fetch PTR(u) here
        !            55:      nice and early */
        !            56:   up = PTR(u);
        !            57:
        !            58:   if ((usize ^ dir) < 0)
        !            59:     {
        !            60:       /* Round towards zero, means just truncate */
        !            61:
        !            62:       if (w == u)
        !            63:         {
        !            64:           /* if already smaller than limb_cnt then do nothing */
        !            65:           if (abs_usize <= limb_cnt)
        !            66:             return;
        !            67:           wp = PTR(w);
        !            68:         }
        !            69:       else
        !            70:         {
        !            71:           i = MIN (abs_usize, limb_cnt+1);
        !            72:           MPZ_REALLOC (w, i);
        !            73:           wp = PTR(w);
        !            74:           MPN_COPY (wp, up, i);
        !            75:
        !            76:           /* if smaller than limb_cnt then only the copy is needed */
        !            77:           if (abs_usize <= limb_cnt)
        !            78:             {
        !            79:               SIZ(w) = usize;
        !            80:               return;
        !            81:             }
        !            82:         }
        !            83:     }
        !            84:   else
        !            85:     {
        !            86:       /* Round away from zero, means twos complement if non-zero */
        !            87:
        !            88:       /* if u!=0 and smaller than divisor, then must negate */
        !            89:       if (abs_usize <= limb_cnt)
        !            90:         goto negate;
        !            91:
        !            92:       /* if non-zero low limb, then must negate */
        !            93:       for (i = 0; i < limb_cnt; i++)
        !            94:         if (up[i] != 0)
        !            95:           goto negate;
        !            96:
        !            97:       /* if non-zero partial limb, then must negate */
        !            98:       if ((up[limb_cnt] & LOW_MASK (cnt)) != 0)
        !            99:         goto negate;
        !           100:
        !           101:       /* otherwise low bits of u are zero, so that's the result */
        !           102:       SIZ(w) = 0;
        !           103:       return;
        !           104:
        !           105:     negate:
        !           106:       /* twos complement negation to get 2**cnt-u */
        !           107:
        !           108:       MPZ_REALLOC (w, limb_cnt+1);
        !           109:       up = PTR(u);
        !           110:       wp = PTR(w);
        !           111:
        !           112:       /* Ones complement */
        !           113:       i = MIN (abs_usize, limb_cnt+1);
        !           114:       mpn_com_n (wp, up, i);
        !           115:       for ( ; i <= limb_cnt; i++)
        !           116:         wp[i] = GMP_NUMB_MAX;
        !           117:
        !           118:       /* Twos complement.  Since u!=0 in the relevant part, the twos
        !           119:          complement never gives 0 and a carry, so can use MPN_INCR_U. */
        !           120:       MPN_INCR_U (wp, limb_cnt+1, CNST_LIMB(1));
        !           121:
        !           122:       usize = -usize;
        !           123:     }
        !           124:
        !           125:   /* Mask the high limb */
        !           126:   high = wp[limb_cnt];
        !           127:   high &= LOW_MASK (cnt);
        !           128:   wp[limb_cnt] = high;
        !           129:
        !           130:   /* Strip any consequent high zeros */
        !           131:   while (high == 0)
        !           132:     {
        !           133:       limb_cnt--;
        !           134:       if (limb_cnt < 0)
        !           135:         {
        !           136:           SIZ(w) = 0;
        !           137:           return;
        !           138:         }
        !           139:       high = wp[limb_cnt];
        !           140:     }
        !           141:
        !           142:   limb_cnt++;
        !           143:   SIZ(w) = (usize >= 0 ? limb_cnt : -limb_cnt);
        !           144: }
        !           145:
        !           146:
        !           147: void
        !           148: mpz_cdiv_r_2exp (mpz_ptr w, mpz_srcptr u, unsigned long cnt)
        !           149: {
        !           150:   cfdiv_r_2exp (w, u, cnt, 1);
        !           151: }
        !           152:
        !           153: void
        !           154: mpz_fdiv_r_2exp (mpz_ptr w, mpz_srcptr u, unsigned long cnt)
        !           155: {
        !           156:   cfdiv_r_2exp (w, u, cnt, -1);
        !           157: }