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

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: }