contrib/gmp/mpn/generic/sb_div_q.c

   1 /* mpn_sb_div_q -- schoolbook division with 2-limb sloppy non-greater
   2    precomputed inverse, returning an accurate quotient.
   3
   4    Contributed to the GNU project by Torbjörn Granlund.
   5
   6    THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH A MUTABLE INTERFACE.  IT IS
   7    ONLY SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS
   8    ALMOST GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GMP
   9    RELEASE.
  10
  11 Copyright 2006, 2007 Free Software Foundation, Inc.
  12
  13 This file is part of the GNU MP Library.
  14
  15 The GNU MP Library is free software; you can redistribute it and/or modify
  16 it under the terms of the GNU Lesser General Public License as published by
  17 the Free Software Foundation; either version 3 of the License, or (at your
  18 option) any later version.
  19
  20 The GNU MP Library is distributed in the hope that it will be useful, but
  21 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  22 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  23 License for more details.
  24
  25 You should have received a copy of the GNU Lesser General Public License
  26 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  27
  28 #include "gmp.h"
  29 #include "gmp-impl.h"
  30 #include "longlong.h"
  31
  32 /*
  33   CAVEATS:
  34   1. Should it demand normalized operands like now, or normalize on-the-fly?
  35   2. Overwrites {np,nn}.
  36   3. Uses mpn_submul_1.  It would be nice to somehow make it use mpn_addmul_1
  37      instead.  (That would open for mpn_addmul_2 straightforwardly.)
  38 */
  39
  40 mp_limb_t
  41 mpn_sb_div_q (mp_ptr qp,
  42               mp_ptr np, mp_size_t nn,
  43               mp_srcptr dp, mp_size_t dn,
  44               mp_srcptr dip)
  45 {
  46   mp_limb_t q, q10, q01a, q00a, q01b, q00b;
  47   mp_limb_t cy;
  48   mp_size_t i;
  49   mp_limb_t qh;
  50   mp_limb_t di1, di0;
  51   mp_size_t qn;
  52
  53   mp_size_t dn_orig = dn;
  54   mp_srcptr dp_orig = dp;
  55   mp_ptr np_orig = np;
  56
  57   ASSERT (dn > 0);
  58   ASSERT (nn >= dn);
  59   ASSERT ((dp[dn-1] & GMP_NUMB_HIGHBIT) != 0);
  60   ASSERT (! MPN_OVERLAP_P (np, nn, dp, dn));
  61   ASSERT (! MPN_OVERLAP_P (qp, nn-dn, dp, dn));
  62   ASSERT (! MPN_OVERLAP_P (qp, nn-dn, np, nn) || qp+dn >= np);
  63   ASSERT_MPN (np, nn);
  64   ASSERT_MPN (dp, dn);
  65
  66   np += nn;
  67   qn = nn - dn;
  68   if (qn + 1 < dn)
  69     {
  70       dp += dn - (qn + 1);
  71       dn = qn + 1;
  72     }
  73
  74   qh = mpn_cmp (np - dn, dp, dn) >= 0;
  75   if (qh != 0)
  76     mpn_sub_n (np - dn, np - dn, dp, dn);
  77
  78   qp += qn;
  79   di1 = dip[1]; di0 = dip[0];
  80   for (i = qn; i >= dn; i--)
  81     {
  82       np--;
  83       umul_ppmm (q, q10, np[0], di1);
  84       umul_ppmm (q01a, q00a, np[-1], di1);
  85       add_ssaaaa (q, q10, q, q10, np[0], q01a);
  86       umul_ppmm (q01b, q00b, np[0], di0);
  87       add_ssaaaa (q, q10, q, q10, 0, q01b);
  88       add_ssaaaa (q, q10, q, q10, 0, np[-1]);
  89
  90       cy = mpn_submul_1 (np - dn, dp, dn, q);
  91
  92       if (UNLIKELY (np[0] > cy || mpn_cmp (np - dn, dp, dn) >= 0))
  93         {
  94           q = q + 1;
  95           mpn_sub_n (np - dn, np - dn, dp, dn);
  96         }
  97
  98       *--qp = q;
  99     }
 100
 101   for (i = dn - 1; i > 0; i--)
 102     {
 103       np--;
 104       umul_ppmm (q, q10, np[0], di1);
 105       umul_ppmm (q01a, q00a, np[-1], di1);
 106       add_ssaaaa (q, q10, q, q10, np[0], q01a);
 107       umul_ppmm (q01b, q00b, np[0], di0);
 108       add_ssaaaa (q, q10, q, q10, 0, q01b);
 109       add_ssaaaa (q, q10, q, q10, 0, np[-1]);
 110
 111       cy = mpn_submul_1 (np - dn, dp, dn, q);
 112
 113       if (UNLIKELY (np[0] > cy || mpn_cmp (np - dn, dp, dn) >= 0))
 114         {
 115           q = q + 1;
 116           if (q == 0)
 117             q = GMP_NUMB_MAX;
 118           else
 119             mpn_sub_n (np - dn, np - dn, dp, dn);
 120         }
 121
 122       *--qp = q;
 123
 124       /* Truncate operands.  */
 125       dn--;
 126       dp++;
 127
 128       /* The partial remainder might be equal to the truncated divisor,
 129          thus non-canonical.  When that happens, the rest of the quotient
 130          should be all ones.  */
 131       if (UNLIKELY (mpn_cmp (np - dn, dp, dn) == 0))
 132         {
 133           while (--i)
 134             *--qp = GMP_NUMB_MAX;
 135           break;
 136         }
 137     }
 138
 139   dn = dn_orig;
 140   if (UNLIKELY (np[-1] < dn))
 141     {
 142       mp_limb_t q, x;
 143
 144       /* The quotient may be too large if the remainder is small.  Recompute
 145          for above ignored operand parts, until the remainder spills.
 146
 147          FIXME: The quality of this code isn't the same as the code above.
 148          1. We don't compute things in an optimal order, high-to-low, in order
 149             to terminate as quickly as possible.
 150          2. We mess with pointers and sizes, adding and subtracting and
 151             adjusting to get things right.  It surely could be streamlined.
 152          3. The only termination criteria are that we determine that the
 153             quotient needs to be adjusted, or that we have recomputed
 154             everything.  We should stop when the remainder is so large
 155             that no additional subtracting could make it spill.
 156          4. If nothing else, we should not do two loops of submul_1 over the
 157             data, instead handle both the triangularization and chopping at
 158             once.  */
 159
 160       x = np[-1];
 161
 162       if (dn > 2)
 163         {
 164           /* Compensate for triangularization.  */
 165           mp_limb_t y;
 166
 167           dp = dp_orig;
 168           if (qn + 1 < dn)
 169             {
 170               dp += dn - (qn + 1);
 171               dn = qn + 1;
 172             }
 173
 174           y = np[-2];
 175
 176           for (i = dn - 3; i >= 0; i--)
 177             {
 178               q = qp[i];
 179               cy = mpn_submul_1 (np - (dn - i), dp, dn - i - 2, q);
 180
 181               if (y < cy)
 182                 {
 183                   if (x == 0)
 184                     {
 185                       cy = mpn_sub_1 (qp, qp, qn, 1);
 186                       ASSERT_ALWAYS (cy == 0);
 187                       return qh - cy;
 188                     }
 189                   x--;
 190                 }
 191               y -= cy;
 192             }
 193           np[-2] = y;
 194         }
 195
 196       dn = dn_orig;
 197       if (qn + 1 < dn)
 198         {
 199           /* Compensate for ignored dividend and divisor tails.  */
 200
 201           if (qn == 0)
 202             return qh;
 203
 204           dp = dp_orig;
 205           np = np_orig;
 206
 207           if (qh != 0)
 208             {
 209               cy = mpn_sub_n (np + qn, np + qn, dp, dn - (qn + 1));
 210               if (cy != 0)
 211                 {
 212                   if (x == 0)
 213                     {
 214                       cy = mpn_sub_1 (qp, qp, qn, 1);
 215                       return qh - cy;
 216                     }
 217                   x--;
 218                 }
 219             }
 220
 221           for (i = dn - qn - 2; i >= 0; i--)
 222             {
 223               cy = mpn_submul_1 (np + i, qp, qn, dp[i]);
 224               cy = mpn_sub_1 (np + qn + i, np + qn + i, dn - qn - i - 1, cy);
 225               if (cy != 0)
 226                 {
 227                   if (x == 0)
 228                     {
 229                       cy = mpn_sub_1 (qp, qp, qn, 1);
 230                       ASSERT_ALWAYS (cy == 0);
 231                       return qh - cy;
 232                     }
 233                   x--;
 234                 }
 235             }
 236         }
 237     }
 238
 239   return qh;
 240 }