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.2

1.1       maekawa     1: /* mpn/bdivmod.c: mpn_bdivmod for computing U/V mod 2^d.
                      2:
1.1.1.2 ! maekawa     3: Copyright (C) 1991, 1993, 1994, 1995, 1996, 1999, 2000 Free Software
        !             4: Foundation, Inc.
1.1       maekawa     5:
                      6: This file is part of the GNU MP Library.
                      7:
                      8: The GNU MP Library is free software; you can redistribute it and/or modify
1.1.1.2 ! maekawa     9: it under the terms of the GNU Lesser General Public License as published by
        !            10: the Free Software Foundation; either version 2.1 of the License, or (at your
1.1       maekawa    11: option) any later version.
                     12:
                     13: The GNU MP Library is distributed in the hope that it will be useful, but
                     14: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.1.1.2 ! maekawa    15: or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
1.1       maekawa    16: License for more details.
                     17:
1.1.1.2 ! maekawa    18: You should have received a copy of the GNU Lesser General Public License
1.1       maekawa    19: along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
                     20: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
                     21: MA 02111-1307, USA. */
                     22:
                     23: /* q_high = mpn_bdivmod (qp, up, usize, vp, vsize, d).
                     24:
                     25:    Puts the low d/BITS_PER_MP_LIMB limbs of Q = U / V mod 2^d at qp, and
                     26:    returns the high d%BITS_PER_MP_LIMB bits of Q as the result.
                     27:
                     28:    Also, U - Q * V mod 2^(usize*BITS_PER_MP_LIMB) is placed at up.  Since the
                     29:    low d/BITS_PER_MP_LIMB limbs of this difference are zero, the code allows
                     30:    the limb vectors at qp to overwrite the low limbs at up, provided qp <= up.
                     31:
                     32:    Preconditions:
                     33:    1.  V is odd.
                     34:    2.  usize * BITS_PER_MP_LIMB >= d.
                     35:    3.  If Q and U overlap, qp <= up.
                     36:
                     37:    Ken Weber (kweber@mat.ufrgs.br, kweber@mcs.kent.edu)
                     38:
                     39:    Funding for this work has been partially provided by Conselho Nacional
                     40:    de Desenvolvimento Cienti'fico e Tecnolo'gico (CNPq) do Brazil, Grant
                     41:    301314194-2, and was done while I was a visiting reseacher in the Instituto
                     42:    de Matema'tica at Universidade Federal do Rio Grande do Sul (UFRGS).
                     43:
                     44:    References:
                     45:        T. Jebelean, An algorithm for exact division, Journal of Symbolic
                     46:        Computation, v. 15, 1993, pp. 169-180.
                     47:
                     48:        K. Weber, The accelerated integer GCD algorithm, ACM Transactions on
                     49:        Mathematical Software, v. 21 (March), 1995, pp. 111-122.  */
                     50:
                     51: #include "gmp.h"
                     52: #include "gmp-impl.h"
                     53: #include "longlong.h"
                     54:
                     55: mp_limb_t
                     56: #if __STDC__
                     57: mpn_bdivmod (mp_ptr qp, mp_ptr up, mp_size_t usize,
                     58:             mp_srcptr vp, mp_size_t vsize, unsigned long int d)
                     59: #else
                     60: mpn_bdivmod (qp, up, usize, vp, vsize, d)
                     61:      mp_ptr qp;
                     62:      mp_ptr up;
                     63:      mp_size_t usize;
                     64:      mp_srcptr vp;
                     65:      mp_size_t vsize;
                     66:      unsigned long int d;
                     67: #endif
                     68: {
1.1.1.2 ! maekawa    69:   mp_limb_t v_inv;
1.1       maekawa    70:
1.1.1.2 ! maekawa    71:   /* 1/V mod 2^BITS_PER_MP_LIMB. */
        !            72:   modlimb_invert (v_inv, vp[0]);
1.1       maekawa    73:
1.1.1.2 ! maekawa    74:   /* Fast code for two cases previously used by the accel part of mpn_gcd.
        !            75:      (Could probably remove this now it's inlined there.) */
1.1       maekawa    76:   if (usize == 2 && vsize == 2 &&
                     77:       (d == BITS_PER_MP_LIMB || d == 2*BITS_PER_MP_LIMB))
                     78:     {
                     79:       mp_limb_t hi, lo;
                     80:       mp_limb_t q = up[0] * v_inv;
                     81:       umul_ppmm (hi, lo, q, vp[0]);
                     82:       up[0] = 0, up[1] -= hi + q*vp[1], qp[0] = q;
                     83:       if (d == 2*BITS_PER_MP_LIMB)
                     84:        q = up[1] * v_inv, up[1] = 0, qp[1] = q;
                     85:       return 0;
                     86:     }
                     87:
                     88:   /* Main loop.  */
                     89:   while (d >= BITS_PER_MP_LIMB)
                     90:     {
                     91:       mp_limb_t q = up[0] * v_inv;
                     92:       mp_limb_t b = mpn_submul_1 (up, vp, MIN (usize, vsize), q);
                     93:       if (usize > vsize)
                     94:        mpn_sub_1 (up + vsize, up + vsize, usize - vsize, b);
                     95:       d -= BITS_PER_MP_LIMB;
                     96:       up += 1, usize -= 1;
                     97:       *qp++ = q;
                     98:     }
                     99:
                    100:   if (d)
                    101:     {
                    102:       mp_limb_t b;
                    103:       mp_limb_t q = (up[0] * v_inv) & (((mp_limb_t)1<<d) - 1);
1.1.1.2 ! maekawa   104:       if (q <= 1)
1.1       maekawa   105:        {
1.1.1.2 ! maekawa   106:          if (q == 0)
        !           107:            return 0;
        !           108:          else
        !           109:            b = mpn_sub_n (up, up, vp, MIN (usize, vsize));
1.1       maekawa   110:        }
1.1.1.2 ! maekawa   111:       else
        !           112:        b = mpn_submul_1 (up, vp, MIN (usize, vsize), q);
        !           113:
1.1       maekawa   114:       if (usize > vsize)
                    115:        mpn_sub_1 (up + vsize, up + vsize, usize - vsize, b);
                    116:       return q;
                    117:     }
                    118:
                    119:   return 0;
                    120: }