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