Return to bdivmod.c CVS log

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

Annotation of OpenXM_contrib/gmp/mpn/generic/bdivmod.c, Revision 1.1

1.1     ! maekawa     1: /* mpn/bdivmod.c: mpn_bdivmod for computing U/V mod 2^d.
        !             2:
        !             3: Copyright (C) 1991, 1993, 1994, 1995, 1996 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 Library General Public License as published by
        !             9: the Free Software Foundation; either version 2 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 Library General Public
        !            15: License for more details.
        !            16:
        !            17: You should have received a copy of the GNU Library 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: /* q_high = mpn_bdivmod (qp, up, usize, vp, vsize, d).
        !            23:
        !            24:    Puts the low d/BITS_PER_MP_LIMB limbs of Q = U / V mod 2^d at qp, and
        !            25:    returns the high d%BITS_PER_MP_LIMB bits of Q as the result.
        !            26:
        !            27:    Also, U - Q * V mod 2^(usize*BITS_PER_MP_LIMB) is placed at up.  Since the
        !            28:    low d/BITS_PER_MP_LIMB limbs of this difference are zero, the code allows
        !            29:    the limb vectors at qp to overwrite the low limbs at up, provided qp <= up.
        !            30:
        !            31:    Preconditions:
        !            32:    1.  V is odd.
        !            33:    2.  usize * BITS_PER_MP_LIMB >= d.
        !            34:    3.  If Q and U overlap, qp <= up.
        !            35:
        !            36:    Ken Weber (kweber@mat.ufrgs.br, kweber@mcs.kent.edu)
        !            37:
        !            38:    Funding for this work has been partially provided by Conselho Nacional
        !            39:    de Desenvolvimento Cienti'fico e Tecnolo'gico (CNPq) do Brazil, Grant
        !            40:    301314194-2, and was done while I was a visiting reseacher in the Instituto
        !            41:    de Matema'tica at Universidade Federal do Rio Grande do Sul (UFRGS).
        !            42:
        !            43:    References:
        !            44:        T. Jebelean, An algorithm for exact division, Journal of Symbolic
        !            45:        Computation, v. 15, 1993, pp. 169-180.
        !            46:
        !            47:        K. Weber, The accelerated integer GCD algorithm, ACM Transactions on
        !            48:        Mathematical Software, v. 21 (March), 1995, pp. 111-122.  */
        !            49:
        !            50: #include "gmp.h"
        !            51: #include "gmp-impl.h"
        !            52: #include "longlong.h"
        !            53:
        !            54: mp_limb_t
        !            55: #if __STDC__
        !            56: mpn_bdivmod (mp_ptr qp, mp_ptr up, mp_size_t usize,
        !            57:             mp_srcptr vp, mp_size_t vsize, unsigned long int d)
        !            58: #else
        !            59: mpn_bdivmod (qp, up, usize, vp, vsize, d)
        !            60:      mp_ptr qp;
        !            61:      mp_ptr up;
        !            62:      mp_size_t usize;
        !            63:      mp_srcptr vp;
        !            64:      mp_size_t vsize;
        !            65:      unsigned long int d;
        !            66: #endif
        !            67: {
        !            68:   /* Cache for v_inv is used to make mpn_accelgcd faster.  */
        !            69:   static mp_limb_t previous_low_vlimb = 0;
        !            70:   static mp_limb_t v_inv;              /* 1/V mod 2^BITS_PER_MP_LIMB.  */
        !            71:
        !            72:   if (vp[0] != previous_low_vlimb)     /* Cache miss; compute v_inv.  */
        !            73:     {
        !            74:       mp_limb_t v = previous_low_vlimb = vp[0];
        !            75:       mp_limb_t make_zero = 1;
        !            76:       mp_limb_t two_i = 1;
        !            77:       v_inv = 0;
        !            78:       do
        !            79:        {
        !            80:          while ((two_i & make_zero) == 0)
        !            81:            two_i <<= 1, v <<= 1;
        !            82:          v_inv += two_i;
        !            83:          make_zero -= v;
        !            84:        }
        !            85:       while (make_zero);
        !            86:     }
        !            87:
        !            88:   /* Need faster computation for some common cases in mpn_accelgcd.  */
        !            89:   if (usize == 2 && vsize == 2 &&
        !            90:       (d == BITS_PER_MP_LIMB || d == 2*BITS_PER_MP_LIMB))
        !            91:     {
        !            92:       mp_limb_t hi, lo;
        !            93:       mp_limb_t q = up[0] * v_inv;
        !            94:       umul_ppmm (hi, lo, q, vp[0]);
        !            95:       up[0] = 0, up[1] -= hi + q*vp[1], qp[0] = q;
        !            96:       if (d == 2*BITS_PER_MP_LIMB)
        !            97:        q = up[1] * v_inv, up[1] = 0, qp[1] = q;
        !            98:       return 0;
        !            99:     }
        !           100:
        !           101:   /* Main loop.  */
        !           102:   while (d >= BITS_PER_MP_LIMB)
        !           103:     {
        !           104:       mp_limb_t q = up[0] * v_inv;
        !           105:       mp_limb_t b = mpn_submul_1 (up, vp, MIN (usize, vsize), q);
        !           106:       if (usize > vsize)
        !           107:        mpn_sub_1 (up + vsize, up + vsize, usize - vsize, b);
        !           108:       d -= BITS_PER_MP_LIMB;
        !           109:       up += 1, usize -= 1;
        !           110:       *qp++ = q;
        !           111:     }
        !           112:
        !           113:   if (d)
        !           114:     {
        !           115:       mp_limb_t b;
        !           116:       mp_limb_t q = (up[0] * v_inv) & (((mp_limb_t)1<<d) - 1);
        !           117:       switch (q)
        !           118:        {
        !           119:          case 0:  return 0;
        !           120:          case 1:  b = mpn_sub_n (up, up, vp, MIN (usize, vsize));   break;
        !           121:          default: b = mpn_submul_1 (up, vp, MIN (usize, vsize), q); break;
        !           122:        }
        !           123:       if (usize > vsize)
        !           124:        mpn_sub_1 (up + vsize, up + vsize, usize - vsize, b);
        !           125:       return q;
        !           126:     }
        !           127:
        !           128:   return 0;
        !           129: }