Return to divrem_2.c CVS log

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

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

1.1     ! maekawa     1: /* mpn_divrem_2 -- Divide natural numbers, producing both remainder and
        !             2:    quotient.  The divisor is two limbs.
        !             3:
        !             4:    THIS FILE CONTAINS INTERNAL FUNCTIONS WITH MUTABLE INTERFACES.  IT IS
        !             5:    ONLY SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS
        !             6:    ALMOST GUARANTEED THAT THEY'LL CHANGE OR DISAPPEAR IN A FUTURE GNU MP
        !             7:    RELEASE.
        !             8:
        !             9:
        !            10: Copyright (C) 1993, 1994, 1995, 1996, 1999, 2000 Free Software Foundation,
        !            11: Inc.
        !            12:
        !            13: This file is part of the GNU MP Library.
        !            14:
        !            15: The GNU MP Library is free software; you can redistribute it and/or modify
        !            16: it under the terms of the GNU Lesser General Public License as published by
        !            17: the Free Software Foundation; either version 2.1 of the License, or (at your
        !            18: option) any later version.
        !            19:
        !            20: The GNU MP Library is distributed in the hope that it will be useful, but
        !            21: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
        !            22: or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
        !            23: License for more details.
        !            24:
        !            25: You should have received a copy of the GNU Lesser General Public License
        !            26: along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
        !            27: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
        !            28: MA 02111-1307, USA. */
        !            29:
        !            30: #include "gmp.h"
        !            31: #include "gmp-impl.h"
        !            32: #include "longlong.h"
        !            33:
        !            34: /* Divide num (NP/NSIZE) by den (DP/2) and write
        !            35:    the NSIZE-2 least significant quotient limbs at QP
        !            36:    and the 2 long remainder at NP.  If QEXTRA_LIMBS is
        !            37:    non-zero, generate that many fraction bits and append them after the
        !            38:    other quotient limbs.
        !            39:    Return the most significant limb of the quotient, this is always 0 or 1.
        !            40:
        !            41:    Preconditions:
        !            42:    0. NSIZE >= 2.
        !            43:    1. The most significant bit of the divisor must be set.
        !            44:    2. QP must either not overlap with the input operands at all, or
        !            45:       QP + 2 >= NP must hold true.  (This means that it's
        !            46:       possible to put the quotient in the high part of NUM, right after the
        !            47:       remainder in NUM.
        !            48:    3. NSIZE >= 2, even if QEXTRA_LIMBS is non-zero.  */
        !            49:
        !            50: mp_limb_t
        !            51: #if __STDC__
        !            52: mpn_divrem_2 (mp_ptr qp, mp_size_t qxn,
        !            53:              mp_ptr np, mp_size_t nsize,
        !            54:              mp_srcptr dp)
        !            55: #else
        !            56: mpn_divrem_2 (qp, qxn, np, nsize, dp)
        !            57:      mp_ptr qp;
        !            58:      mp_size_t qxn;
        !            59:      mp_ptr np;
        !            60:      mp_size_t nsize;
        !            61:      mp_srcptr dp;
        !            62: #endif
        !            63: {
        !            64:   mp_limb_t most_significant_q_limb = 0;
        !            65:   mp_size_t i;
        !            66:   mp_limb_t n1, n0, n2;
        !            67:   mp_limb_t d1, d0;
        !            68:   mp_limb_t d1inv;
        !            69:   int have_preinv;
        !            70:
        !            71:   np += nsize - 2;
        !            72:   d1 = dp[1];
        !            73:   d0 = dp[0];
        !            74:   n1 = np[1];
        !            75:   n0 = np[0];
        !            76:
        !            77:   if (n1 >= d1 && (n1 > d1 || n0 >= d0))
        !            78:     {
        !            79:       sub_ddmmss (n1, n0, n1, n0, d1, d0);
        !            80:       most_significant_q_limb = 1;
        !            81:     }
        !            82:
        !            83:   /* If multiplication is much faster than division, preinvert the most
        !            84:      significant divisor limb before entering the loop.  */
        !            85:   if (UDIV_TIME > 2 * UMUL_TIME + 6)
        !            86:     {
        !            87:       have_preinv = 0;
        !            88:       if ((UDIV_TIME - (2 * UMUL_TIME + 6)) * (nsize - 2) > UDIV_TIME)
        !            89:        {
        !            90:          invert_limb (d1inv, d1);
        !            91:          have_preinv = 1;
        !            92:        }
        !            93:     }
        !            94:
        !            95:   for (i = qxn + nsize - 2 - 1; i >= 0; i--)
        !            96:     {
        !            97:       mp_limb_t q;
        !            98:       mp_limb_t r;
        !            99:
        !           100:       if (i >= qxn)
        !           101:        np--;
        !           102:       else
        !           103:        np[0] = 0;
        !           104:
        !           105:       if (n1 == d1)
        !           106:        {
        !           107:          /* Q should be either 111..111 or 111..110.  Need special treatment
        !           108:             of this rare case as normal division would give overflow.  */
        !           109:          q = ~(mp_limb_t) 0;
        !           110:
        !           111:          r = n0 + d1;
        !           112:          if (r < d1)   /* Carry in the addition? */
        !           113:            {
        !           114:              add_ssaaaa (n1, n0, r - d0, np[0], 0, d0);
        !           115:              qp[i] = q;
        !           116:              continue;
        !           117:            }
        !           118:          n1 = d0 - (d0 != 0);
        !           119:          n0 = -d0;
        !           120:        }
        !           121:       else
        !           122:        {
        !           123:          if (UDIV_TIME > 2 * UMUL_TIME + 6 && have_preinv)
        !           124:            udiv_qrnnd_preinv (q, r, n1, n0, d1, d1inv);
        !           125:          else
        !           126:            udiv_qrnnd (q, r, n1, n0, d1);
        !           127:          umul_ppmm (n1, n0, d0, q);
        !           128:        }
        !           129:
        !           130:       n2 = np[0];
        !           131:
        !           132:     q_test:
        !           133:       if (n1 > r || (n1 == r && n0 > n2))
        !           134:        {
        !           135:          /* The estimated Q was too large.  */
        !           136:          q--;
        !           137:
        !           138:          sub_ddmmss (n1, n0, n1, n0, 0, d0);
        !           139:          r += d1;
        !           140:          if (r >= d1)  /* If not carry, test Q again.  */
        !           141:            goto q_test;
        !           142:        }
        !           143:
        !           144:       qp[i] = q;
        !           145:       sub_ddmmss (n1, n0, r, n2, n1, n0);
        !           146:     }
        !           147:   np[1] = n1;
        !           148:   np[0] = n0;
        !           149:
        !           150:   return most_significant_q_limb;
        !           151: }