Back to index

glibc  2.9
divrem.c
Go to the documentation of this file.
00001 /* mpn_divrem -- Divide natural numbers, producing both remainder and
00002    quotient.
00003 
00004 Copyright (C) 1993, 1994, 1995, 1996 Free Software Foundation, Inc.
00005 
00006 This file is part of the GNU MP Library.
00007 
00008 The GNU MP Library is free software; you can redistribute it and/or modify
00009 it under the terms of the GNU Lesser General Public License as published by
00010 the Free Software Foundation; either version 2.1 of the License, or (at your
00011 option) any later version.
00012 
00013 The GNU MP Library is distributed in the hope that it will be useful, but
00014 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
00015 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00016 License for more details.
00017 
00018 You should have received a copy of the GNU Lesser General Public License
00019 along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
00020 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
00021 MA 02111-1307, USA. */
00022 
00023 #include <gmp.h>
00024 #include "gmp-impl.h"
00025 #include "longlong.h"
00026 
00027 /* Divide num (NP/NSIZE) by den (DP/DSIZE) and write
00028    the NSIZE-DSIZE least significant quotient limbs at QP
00029    and the DSIZE long remainder at NP.  If QEXTRA_LIMBS is
00030    non-zero, generate that many fraction bits and append them after the
00031    other quotient limbs.
00032    Return the most significant limb of the quotient, this is always 0 or 1.
00033 
00034    Preconditions:
00035    0. NSIZE >= DSIZE.
00036    1. The most significant bit of the divisor must be set.
00037    2. QP must either not overlap with the input operands at all, or
00038       QP + DSIZE >= NP must hold true.  (This means that it's
00039       possible to put the quotient in the high part of NUM, right after the
00040       remainder in NUM.
00041    3. NSIZE >= DSIZE, even if QEXTRA_LIMBS is non-zero.  */
00042 
00043 mp_limb_t
00044 #if __STDC__
00045 mpn_divrem (mp_ptr qp, mp_size_t qextra_limbs,
00046            mp_ptr np, mp_size_t nsize,
00047            mp_srcptr dp, mp_size_t dsize)
00048 #else
00049 mpn_divrem (qp, qextra_limbs, np, nsize, dp, dsize)
00050      mp_ptr qp;
00051      mp_size_t qextra_limbs;
00052      mp_ptr np;
00053      mp_size_t nsize;
00054      mp_srcptr dp;
00055      mp_size_t dsize;
00056 #endif
00057 {
00058   mp_limb_t most_significant_q_limb = 0;
00059 
00060   switch (dsize)
00061     {
00062     case 0:
00063       /* We are asked to divide by zero, so go ahead and do it!  (To make
00064         the compiler not remove this statement, return the value.)  */
00065       return 1 / dsize;
00066 
00067     case 1:
00068       {
00069        mp_size_t i;
00070        mp_limb_t n1;
00071        mp_limb_t d;
00072 
00073        d = dp[0];
00074        n1 = np[nsize - 1];
00075 
00076        if (n1 >= d)
00077          {
00078            n1 -= d;
00079            most_significant_q_limb = 1;
00080          }
00081 
00082        qp += qextra_limbs;
00083        for (i = nsize - 2; i >= 0; i--)
00084          udiv_qrnnd (qp[i], n1, n1, np[i], d);
00085        qp -= qextra_limbs;
00086 
00087        for (i = qextra_limbs - 1; i >= 0; i--)
00088          udiv_qrnnd (qp[i], n1, n1, 0, d);
00089 
00090        np[0] = n1;
00091       }
00092       break;
00093 
00094     case 2:
00095       {
00096        mp_size_t i;
00097        mp_limb_t n1, n0, n2;
00098        mp_limb_t d1, d0;
00099 
00100        np += nsize - 2;
00101        d1 = dp[1];
00102        d0 = dp[0];
00103        n1 = np[1];
00104        n0 = np[0];
00105 
00106        if (n1 >= d1 && (n1 > d1 || n0 >= d0))
00107          {
00108            sub_ddmmss (n1, n0, n1, n0, d1, d0);
00109            most_significant_q_limb = 1;
00110          }
00111 
00112        for (i = qextra_limbs + nsize - 2 - 1; i >= 0; i--)
00113          {
00114            mp_limb_t q;
00115            mp_limb_t r;
00116 
00117            if (i >= qextra_limbs)
00118              np--;
00119            else
00120              np[0] = 0;
00121 
00122            if (n1 == d1)
00123              {
00124               /* Q should be either 111..111 or 111..110.  Need special
00125                  treatment of this rare case as normal division would
00126                  give overflow.  */
00127               q = ~(mp_limb_t) 0;
00128 
00129               r = n0 + d1;
00130               if (r < d1)   /* Carry in the addition? */
00131                 {
00132                   add_ssaaaa (n1, n0, r - d0, np[0], 0, d0);
00133                   qp[i] = q;
00134                   continue;
00135                 }
00136               n1 = d0 - (d0 != 0);
00137               n0 = -d0;
00138              }
00139            else
00140              {
00141               udiv_qrnnd (q, r, n1, n0, d1);
00142               umul_ppmm (n1, n0, d0, q);
00143              }
00144 
00145            n2 = np[0];
00146          q_test:
00147            if (n1 > r || (n1 == r && n0 > n2))
00148              {
00149               /* The estimated Q was too large.  */
00150               q--;
00151 
00152               sub_ddmmss (n1, n0, n1, n0, 0, d0);
00153               r += d1;
00154               if (r >= d1)  /* If not carry, test Q again.  */
00155                 goto q_test;
00156              }
00157 
00158            qp[i] = q;
00159            sub_ddmmss (n1, n0, r, n2, n1, n0);
00160          }
00161        np[1] = n1;
00162        np[0] = n0;
00163       }
00164       break;
00165 
00166     default:
00167       {
00168        mp_size_t i;
00169        mp_limb_t dX, d1, n0;
00170 
00171        np += nsize - dsize;
00172        dX = dp[dsize - 1];
00173        d1 = dp[dsize - 2];
00174        n0 = np[dsize - 1];
00175 
00176        if (n0 >= dX)
00177          {
00178            if (n0 > dX || mpn_cmp (np, dp, dsize - 1) >= 0)
00179              {
00180               mpn_sub_n (np, np, dp, dsize);
00181               n0 = np[dsize - 1];
00182               most_significant_q_limb = 1;
00183              }
00184          }
00185 
00186        for (i = qextra_limbs + nsize - dsize - 1; i >= 0; i--)
00187          {
00188            mp_limb_t q;
00189            mp_limb_t n1, n2;
00190            mp_limb_t cy_limb;
00191 
00192            if (i >= qextra_limbs)
00193              {
00194               np--;
00195               n2 = np[dsize];
00196              }
00197            else
00198              {
00199               n2 = np[dsize - 1];
00200               MPN_COPY_DECR (np + 1, np, dsize);
00201               np[0] = 0;
00202              }
00203 
00204            if (n0 == dX)
00205              /* This might over-estimate q, but it's probably not worth
00206                the extra code here to find out.  */
00207              q = ~(mp_limb_t) 0;
00208            else
00209              {
00210               mp_limb_t r;
00211 
00212               udiv_qrnnd (q, r, n0, np[dsize - 1], dX);
00213               umul_ppmm (n1, n0, d1, q);
00214 
00215               while (n1 > r || (n1 == r && n0 > np[dsize - 2]))
00216                 {
00217                   q--;
00218                   r += dX;
00219                   if (r < dX)      /* I.e. "carry in previous addition?"  */
00220                     break;
00221                   n1 -= n0 < d1;
00222                   n0 -= d1;
00223                 }
00224              }
00225 
00226            /* Possible optimization: We already have (q * n0) and (1 * n1)
00227               after the calculation of q.  Taking advantage of that, we
00228               could make this loop make two iterations less.  */
00229 
00230            cy_limb = mpn_submul_1 (np, dp, dsize, q);
00231 
00232            if (n2 != cy_limb)
00233              {
00234               mpn_add_n (np, np, dp, dsize);
00235               q--;
00236              }
00237 
00238            qp[i] = q;
00239            n0 = np[dsize - 1];
00240          }
00241       }
00242     }
00243 
00244   return most_significant_q_limb;
00245 }