sysdeps/generic/divmod.c

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