stdlib/divrem.c

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