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

1.1       maekawa     1: /* mpn/bdivmod.c: mpn_bdivmod for computing U/V mod 2^d.
                      2:
1.1.1.3 ! ohara       3: Copyright 1991, 1993, 1994, 1995, 1996, 1999, 2000, 2001, 2002 Free Software
1.1.1.2   maekawa     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:
1.1.1.3 ! ohara      55:
1.1       maekawa    56: mp_limb_t
                     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: {
1.1.1.2   maekawa    60:   mp_limb_t v_inv;
1.1       maekawa    61:
1.1.1.3 ! ohara      62:   ASSERT (usize >= 1);
        !            63:   ASSERT (vsize >= 1);
        !            64:   ASSERT (usize * GMP_NUMB_BITS >= d);
        !            65:   ASSERT (! MPN_OVERLAP_P (up, usize, vp, vsize));
        !            66:   ASSERT (! MPN_OVERLAP_P (qp, d/GMP_NUMB_BITS, vp, vsize));
        !            67:   ASSERT (MPN_SAME_OR_INCR2_P (qp, d/GMP_NUMB_BITS, up, usize));
        !            68:   ASSERT_MPN (up, usize);
        !            69:   ASSERT_MPN (vp, vsize);
        !            70:
        !            71:   /* 1/V mod 2^GMP_NUMB_BITS. */
1.1.1.2   maekawa    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 &&
1.1.1.3 ! ohara      77:       (d == GMP_NUMB_BITS || d == 2*GMP_NUMB_BITS))
1.1       maekawa    78:     {
                     79:       mp_limb_t hi, lo;
1.1.1.3 ! ohara      80:       mp_limb_t q = (up[0] * v_inv) & GMP_NUMB_MASK;
        !            81:       umul_ppmm (hi, lo, q, vp[0] << GMP_NAIL_BITS);
        !            82:       up[0] = 0;
        !            83:       up[1] -= hi + q*vp[1];
        !            84:       qp[0] = q;
        !            85:       if (d == 2*GMP_NUMB_BITS)
        !            86:         {
        !            87:           q = (up[1] * v_inv) & GMP_NUMB_MASK;
        !            88:           up[1] = 0;
        !            89:           qp[1] = q;
        !            90:         }
1.1       maekawa    91:       return 0;
                     92:     }
                     93:
                     94:   /* Main loop.  */
1.1.1.3 ! ohara      95:   while (d >= GMP_NUMB_BITS)
1.1       maekawa    96:     {
1.1.1.3 ! ohara      97:       mp_limb_t q = (up[0] * v_inv) & GMP_NUMB_MASK;
1.1       maekawa    98:       mp_limb_t b = mpn_submul_1 (up, vp, MIN (usize, vsize), q);
                     99:       if (usize > vsize)
                    100:        mpn_sub_1 (up + vsize, up + vsize, usize - vsize, b);
1.1.1.3 ! ohara     101:       d -= GMP_NUMB_BITS;
1.1       maekawa   102:       up += 1, usize -= 1;
                    103:       *qp++ = q;
                    104:     }
                    105:
                    106:   if (d)
                    107:     {
                    108:       mp_limb_t b;
                    109:       mp_limb_t q = (up[0] * v_inv) & (((mp_limb_t)1<<d) - 1);
1.1.1.2   maekawa   110:       if (q <= 1)
1.1       maekawa   111:        {
1.1.1.2   maekawa   112:          if (q == 0)
                    113:            return 0;
                    114:          else
                    115:            b = mpn_sub_n (up, up, vp, MIN (usize, vsize));
1.1       maekawa   116:        }
1.1.1.2   maekawa   117:       else
                    118:        b = mpn_submul_1 (up, vp, MIN (usize, vsize), q);
                    119:
1.1       maekawa   120:       if (usize > vsize)
                    121:        mpn_sub_1 (up + vsize, up + vsize, usize - vsize, b);
                    122:       return q;
                    123:     }
                    124:
                    125:   return 0;
                    126: }