lib/gmp/mpn/divrem.c

   1 /* mpn_divrem -- Divide natural numbers, producing both remainder and
   2    quotient.
   3
   4 Copyright (C) 1993, 1994, 1995, 1996 Free Software Foundation, Inc.
   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
   9 it under the terms of the GNU Library General Public License as published by
  10 the Free Software Foundation; either version 2 of the License, or (at your
  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
  15 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
  16 License for more details.
  17
  18 You should have received a copy of the GNU Library General Public License
  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 #include <config.h>
  24 #include "gmp.h"
  25 #include "gmp-impl.h"
  26 #include "longlong.h"
  27
  28 /* Divide num (NP/NSIZE) by den (DP/DSIZE) and write
  29    the NSIZE-DSIZE least significant quotient limbs at QP
  30    and the DSIZE long remainder at NP.  If QEXTRA_LIMBS is
  31    non-zero, generate that many fraction bits and append them after the
  32    other quotient limbs.
  33    Return the most significant limb of the quotient, this is always 0 or 1.
  34
  35    Preconditions:
  36    0. NSIZE >= DSIZE.
  37    1. The most significant bit of the divisor must be set.
  38    2. QP must either not overlap with the input operands at all, or
  39       QP + DSIZE >= NP must hold true.  (This means that it's
  40       possible to put the quotient in the high part of NUM, right after the
  41       remainder in NUM.
  42    3. NSIZE >= DSIZE, even if QEXTRA_LIMBS is non-zero.  */
  43
  44 mp_limb_t
  45 #if __STDC__
  46 mpn_divrem (mp_ptr qp, mp_size_t qextra_limbs,
  47             mp_ptr np, mp_size_t nsize,
  48             mp_srcptr dp, mp_size_t dsize)
  49 #else
  50 mpn_divrem (qp, qextra_limbs, np, nsize, dp, dsize)
  51      mp_ptr qp;
  52      mp_size_t qextra_limbs;
  53      mp_ptr np;
  54      mp_size_t nsize;
  55      mp_srcptr dp;
  56      mp_size_t dsize;
  57 #endif
  58 {
  59   mp_limb_t most_significant_q_limb = 0;
  60
  61   switch (dsize)
  62     {
  63     case 0:
  64       /* We are asked to divide by zero, so go ahead and do it!  (To make
  65          the compiler not remove this statement, return the value.)  */
  66       return 1 / dsize;
  67
  68     case 1:
  69       {
  70         mp_size_t i;
  71         mp_limb_t n1;
  72         mp_limb_t d;
  73
  74         d = dp[0];
  75         n1 = np[nsize - 1];
  76
  77         if (n1 >= d)
  78           {
  79             n1 -= d;
  80             most_significant_q_limb = 1;
  81           }
  82
  83         qp += qextra_limbs;
  84         for (i = nsize - 2; i >= 0; i--)
  85           udiv_qrnnd (qp[i], n1, n1, np[i], d);
  86         qp -= qextra_limbs;
  87
  88         for (i = qextra_limbs - 1; i >= 0; i--)
  89           udiv_qrnnd (qp[i], n1, n1, 0, d);
  90
  91         np[0] = n1;
  92       }
  93       break;
  94
  95     case 2:
  96       {
  97         mp_size_t i;
  98         mp_limb_t n1, n0, n2;
  99         mp_limb_t d1, d0;
 100
 101         np += nsize - 2;
 102         d1 = dp[1];
 103         d0 = dp[0];
 104         n1 = np[1];
 105         n0 = np[0];
 106
 107         if (n1 >= d1 && (n1 > d1 || n0 >= d0))
 108           {
 109             sub_ddmmss (n1, n0, n1, n0, d1, d0);
 110             most_significant_q_limb = 1;
 111           }
 112
 113         for (i = qextra_limbs + nsize - 2 - 1; i >= 0; i--)
 114           {
 115             mp_limb_t q;
 116             mp_limb_t r;
 117
 118             if (i >= qextra_limbs)
 119               np--;
 120             else
 121               np[0] = 0;
 122
 123             if (n1 == d1)
 124               {
 125                 /* Q should be either 111..111 or 111..110.  Need special
 126                    treatment of this rare case as normal division would
 127                    give overflow.  */
 128                 q = ~(mp_limb_t) 0;
 129
 130                 r = n0 + d1;
 131                 if (r < d1)     /* Carry in the addition? */
 132                   {
 133                     add_ssaaaa (n1, n0, r - d0, np[0], 0, d0);
 134                     qp[i] = q;
 135                     continue;
 136                   }
 137                 n1 = d0 - (d0 != 0);
 138                 n0 = -d0;
 139               }
 140             else
 141               {
 142                 udiv_qrnnd (q, r, n1, n0, d1);
 143                 umul_ppmm (n1, n0, d0, q);
 144               }
 145
 146             n2 = np[0];
 147           q_test:
 148             if (n1 > r || (n1 == r && n0 > n2))
 149               {
 150                 /* The estimated Q was too large.  */
 151                 q--;
 152
 153                 sub_ddmmss (n1, n0, n1, n0, 0, d0);
 154                 r += d1;
 155                 if (r >= d1)    /* If not carry, test Q again.  */
 156                   goto q_test;
 157               }
 158
 159             qp[i] = q;
 160             sub_ddmmss (n1, n0, r, n2, n1, n0);
 161           }
 162         np[1] = n1;
 163         np[0] = n0;
 164       }
 165       break;
 166
 167     default:
 168       {
 169         mp_size_t i;
 170         mp_limb_t dX, d1, n0;
 171
 172         np += nsize - dsize;
 173         dX = dp[dsize - 1];
 174         d1 = dp[dsize - 2];
 175         n0 = np[dsize - 1];
 176
 177         if (n0 >= dX)
 178           {
 179             if (n0 > dX || mpn_cmp (np, dp, dsize - 1) >= 0)
 180               {
 181                 mpn_sub_n (np, np, dp, dsize);
 182                 n0 = np[dsize - 1];
 183                 most_significant_q_limb = 1;
 184               }
 185           }
 186
 187         for (i = qextra_limbs + nsize - dsize - 1; i >= 0; i--)
 188           {
 189             mp_limb_t q;
 190             mp_limb_t n1, n2;
 191             mp_limb_t cy_limb;
 192
 193             if (i >= qextra_limbs)
 194               {
 195                 np--;
 196                 n2 = np[dsize];
 197               }
 198             else
 199               {
 200                 n2 = np[dsize - 1];
 201                 MPN_COPY_DECR (np + 1, np, dsize);
 202                 np[0] = 0;
 203               }
 204
 205             if (n0 == dX)
 206               /* This might over-estimate q, but it's probably not worth
 207                  the extra code here to find out.  */
 208               q = ~(mp_limb_t) 0;
 209             else
 210               {
 211                 mp_limb_t r;
 212
 213                 udiv_qrnnd (q, r, n0, np[dsize - 1], dX);
 214                 umul_ppmm (n1, n0, d1, q);
 215
 216                 while (n1 > r || (n1 == r && n0 > np[dsize - 2]))
 217                   {
 218                     q--;
 219                     r += dX;
 220                     if (r < dX) /* I.e. "carry in previous addition?"  */
 221                       break;
 222                     n1 -= n0 < d1;
 223                     n0 -= d1;
 224                   }
 225               }
 226
 227             /* Possible optimization: We already have (q * n0) and (1 * n1)
 228                after the calculation of q.  Taking advantage of that, we
 229                could make this loop make two iterations less.  */
 230
 231             cy_limb = mpn_submul_1 (np, dp, dsize, q);
 232
 233             if (n2 != cy_limb)
 234               {
 235                 mpn_add_n (np, np, dp, dsize);
 236                 q--;
 237               }
 238
 239             qp[i] = q;
 240             n0 = np[dsize - 1];
 241           }
 242       }
 243     }
 244
 245   return most_significant_q_limb;
 246 }