Return to pre_divrem_1.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpn / generic

Annotation of OpenXM_contrib/gmp/mpn/generic/pre_divrem_1.c, Revision 1.1.1.1

1.1       ohara       1: /* mpn_preinv_divrem_1 -- mpn by limb division with pre-inverted divisor.
                      2:
                      3: Copyright 2000, 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: #include "longlong.h"
                     25:
                     26:
                     27: /* Don't bloat a shared library with unused code. */
                     28: #if USE_PREINV_DIVREM_1
                     29:
                     30: /* Same test here for skipping one divide step as in mpn_divrem_1.
                     31:
                     32:    The main reason for a separate shift==0 case is that not all CPUs give
                     33:    zero for "n0 >> BITS_PER_MP_LIMB" which would arise in the general case
                     34:    code used on shift==0.  shift==0 is also reasonably common in __mp_bases
                     35:    big_base, for instance base==10 on a 64-bit limb.
                     36:
                     37:    Under shift!=0 it would be possible to call mpn_lshift to adjust the
                     38:    dividend all in one go (into the quotient space say), rather than
                     39:    limb-by-limb in the loop.  This might help if mpn_lshift is a lot faster
                     40:    than what the compiler can generate for EXTRACT.  But this is left to CPU
                     41:    specific implementations to consider, especially since EXTRACT isn't on
                     42:    the dependent chain.
                     43:
                     44:    If size==0 then the result is simply xsize limbs of zeros, but nothing
                     45:    special is done for that, since it wouldn't be a usual call, and
                     46:    certainly never arises from mpn_get_str which is our main caller.  */
                     47:
                     48: mp_limb_t
                     49: mpn_preinv_divrem_1 (mp_ptr qp, mp_size_t xsize,
                     50:                     mp_srcptr ap, mp_size_t size, mp_limb_t d_unnorm,
                     51:                     mp_limb_t dinv, int shift)
                     52: {
                     53:   mp_limb_t  ahigh, qhigh, r;
                     54:   mp_size_t  i;
                     55:   mp_limb_t  n1, n0;
                     56:   mp_limb_t  d;
                     57:
                     58:   ASSERT (xsize >= 0);
                     59:   ASSERT (size >= 1);
                     60:   ASSERT (d_unnorm != 0);
                     61: #if WANT_ASSERT
                     62:   {
                     63:     int        want_shift;
                     64:     mp_limb_t  want_dinv;
                     65:     count_leading_zeros (want_shift, d_unnorm);
                     66:     ASSERT (shift == want_shift);
                     67:     invert_limb (want_dinv, d_unnorm << shift);
                     68:     ASSERT (dinv == want_dinv);
                     69:   }
                     70: #endif
                     71:   /* FIXME: What's the correct overlap rule when xsize!=0? */
                     72:   ASSERT (MPN_SAME_OR_SEPARATE_P (qp+xsize, ap, size));
                     73:
                     74:   ahigh = ap[size-1];
                     75:   d = d_unnorm << shift;
                     76:   qp += (size + xsize - 1);   /* dest high limb */
                     77:
                     78:   if (shift == 0)
                     79:     {
                     80:       /* High quotient limb is 0 or 1, and skip a divide step. */
                     81:       r = ahigh;
                     82:       qhigh = (r >= d);
                     83:       r = (qhigh ? r-d : r);
                     84:       *qp-- = qhigh;
                     85:       size--;
                     86:
                     87:       for (i = size-1; i >= 0; i--)
                     88:        {
                     89:          n0 = ap[i];
                     90:          udiv_qrnnd_preinv (*qp, r, r, n0, d, dinv);
                     91:          qp--;
                     92:        }
                     93:     }
                     94:   else
                     95:     {
                     96:       r = 0;
                     97:       if (ahigh < d_unnorm)
                     98:        {
                     99:          r = ahigh << shift;
                    100:          *qp-- = 0;
                    101:          size--;
                    102:          if (size == 0)
                    103:            goto done_integer;
                    104:        }
                    105:
                    106:       n1 = ap[size-1];
                    107:       r |= n1 >> (BITS_PER_MP_LIMB - shift);
                    108:
                    109:       for (i = size-2; i >= 0; i--)
                    110:        {
                    111:          ASSERT (r < d);
                    112:          n0 = ap[i];
                    113:          udiv_qrnnd_preinv (*qp, r, r,
                    114:                             ((n1 << shift) | (n0 >> (BITS_PER_MP_LIMB - shift))),
                    115:                             d, dinv);
                    116:          qp--;
                    117:          n1 = n0;
                    118:        }
                    119:       udiv_qrnnd_preinv (*qp, r, r, n1 << shift, d, dinv);
                    120:       qp--;
                    121:     }
                    122:
                    123:  done_integer:
                    124:   for (i = 0; i < xsize; i++)
                    125:     {
                    126:       udiv_qrnnd_preinv (*qp, r, r, 0, d, dinv);
                    127:       qp--;
                    128:     }
                    129:
                    130:   return r >> shift;
                    131: }
                    132:
                    133: #endif /* USE_PREINV_DIVREM_1 */