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

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: }