source/libs/gmp/gmp-src/mpn/generic/sbpi1_div_q.c

   1 /* mpn_sbpi1_div_q -- Schoolbook division using the Möller-Granlund 3/2
   2    division algorithm.
   3
   4    Contributed to the GNU project by Torbjorn Granlund.
   5
   6    THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE.  IT IS ONLY
   7    SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
   8    GUARANTEED THAT IT WILL CHANGE OR DISAPPEAR IN A FUTURE GMP RELEASE.
   9
  10 Copyright 2007, 2009 Free Software Foundation, Inc.
  11
  12 This file is part of the GNU MP Library.
  13
  14 The GNU MP Library is free software; you can redistribute it and/or modify
  15 it under the terms of either:
  16
  17   * the GNU Lesser General Public License as published by the Free
  18     Software Foundation; either version 3 of the License, or (at your
  19     option) any later version.
  20
  21 or
  22
  23   * the GNU General Public License as published by the Free Software
  24     Foundation; either version 2 of the License, or (at your option) any
  25     later version.
  26
  27 or both in parallel, as here.
  28
  29 The GNU MP Library is distributed in the hope that it will be useful, but
  30 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  31 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  32 for more details.
  33
  34 You should have received copies of the GNU General Public License and the
  35 GNU Lesser General Public License along with the GNU MP Library.  If not,
  36 see https://www.gnu.org/licenses/.  */
  37
  38
  39 #include "gmp.h"
  40 #include "gmp-impl.h"
  41 #include "longlong.h"
  42
  43 mp_limb_t
  44 mpn_sbpi1_div_q (mp_ptr qp,
  45                  mp_ptr np, mp_size_t nn,
  46                  mp_srcptr dp, mp_size_t dn,
  47                  mp_limb_t dinv)
  48 {
  49   mp_limb_t qh;
  50   mp_size_t qn, i;
  51   mp_limb_t n1, n0;
  52   mp_limb_t d1, d0;
  53   mp_limb_t cy, cy1;
  54   mp_limb_t q;
  55   mp_limb_t flag;
  56
  57   mp_size_t dn_orig = dn;
  58   mp_srcptr dp_orig = dp;
  59   mp_ptr np_orig = np;
  60
  61   ASSERT (dn > 2);
  62   ASSERT (nn >= dn);
  63   ASSERT ((dp[dn-1] & GMP_NUMB_HIGHBIT) != 0);
  64
  65   np += nn;
  66
  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
  80   dn -= 2;                      /* offset dn by 2 for main division loops,
  81                                    saving two iterations in mpn_submul_1.  */
  82   d1 = dp[dn + 1];
  83   d0 = dp[dn + 0];
  84
  85   np -= 2;
  86
  87   n1 = np[1];
  88
  89   for (i = qn - (dn + 2); i >= 0; i--)
  90     {
  91       np--;
  92       if (UNLIKELY (n1 == d1) && np[1] == d0)
  93         {
  94           q = GMP_NUMB_MASK;
  95           mpn_submul_1 (np - dn, dp, dn + 2, q);
  96           n1 = np[1];           /* update n1, last loop's value will now be invalid */
  97         }
  98       else
  99         {
 100           udiv_qr_3by2 (q, n1, n0, n1, np[1], np[0], d1, d0, dinv);
 101
 102           cy = mpn_submul_1 (np - dn, dp, dn, q);
 103
 104           cy1 = n0 < cy;
 105           n0 = (n0 - cy) & GMP_NUMB_MASK;
 106           cy = n1 < cy1;
 107           n1 -= cy1;
 108           np[0] = n0;
 109
 110           if (UNLIKELY (cy != 0))
 111             {
 112               n1 += d1 + mpn_add_n (np - dn, np - dn, dp, dn + 1);
 113               q--;
 114             }
 115         }
 116
 117       *--qp = q;
 118     }
 119
 120   flag = ~CNST_LIMB(0);
 121
 122   if (dn >= 0)
 123     {
 124       for (i = dn; i > 0; i--)
 125         {
 126           np--;
 127           if (UNLIKELY (n1 >= (d1 & flag)))
 128             {
 129               q = GMP_NUMB_MASK;
 130               cy = mpn_submul_1 (np - dn, dp, dn + 2, q);
 131
 132               if (UNLIKELY (n1 != cy))
 133                 {
 134                   if (n1 < (cy & flag))
 135                     {
 136                       q--;
 137                       mpn_add_n (np - dn, np - dn, dp, dn + 2);
 138                     }
 139                   else
 140                     flag = 0;
 141                 }
 142               n1 = np[1];
 143             }
 144           else
 145             {
 146               udiv_qr_3by2 (q, n1, n0, n1, np[1], np[0], d1, d0, dinv);
 147
 148               cy = mpn_submul_1 (np - dn, dp, dn, q);
 149
 150               cy1 = n0 < cy;
 151               n0 = (n0 - cy) & GMP_NUMB_MASK;
 152               cy = n1 < cy1;
 153               n1 -= cy1;
 154               np[0] = n0;
 155
 156               if (UNLIKELY (cy != 0))
 157                 {
 158                   n1 += d1 + mpn_add_n (np - dn, np - dn, dp, dn + 1);
 159                   q--;
 160                 }
 161             }
 162
 163           *--qp = q;
 164
 165           /* Truncate operands.  */
 166           dn--;
 167           dp++;
 168         }
 169
 170       np--;
 171       if (UNLIKELY (n1 >= (d1 & flag)))
 172         {
 173           q = GMP_NUMB_MASK;
 174           cy = mpn_submul_1 (np, dp, 2, q);
 175
 176           if (UNLIKELY (n1 != cy))
 177             {
 178               if (n1 < (cy & flag))
 179                 {
 180                   q--;
 181                   add_ssaaaa (np[1], np[0], np[1], np[0], dp[1], dp[0]);
 182                 }
 183               else
 184                 flag = 0;
 185             }
 186           n1 = np[1];
 187         }
 188       else
 189         {
 190           udiv_qr_3by2 (q, n1, n0, n1, np[1], np[0], d1, d0, dinv);
 191
 192           np[0] = n0;
 193           np[1] = n1;
 194         }
 195
 196       *--qp = q;
 197     }
 198   ASSERT_ALWAYS (np[1] == n1);
 199   np += 2;
 200
 201
 202   dn = dn_orig;
 203   if (UNLIKELY (n1 < (dn & flag)))
 204     {
 205       mp_limb_t q, x;
 206
 207       /* The quotient may be too large if the remainder is small.  Recompute
 208          for above ignored operand parts, until the remainder spills.
 209
 210          FIXME: The quality of this code isn't the same as the code above.
 211          1. We don't compute things in an optimal order, high-to-low, in order
 212             to terminate as quickly as possible.
 213          2. We mess with pointers and sizes, adding and subtracting and
 214             adjusting to get things right.  It surely could be streamlined.
 215          3. The only termination criteria are that we determine that the
 216             quotient needs to be adjusted, or that we have recomputed
 217             everything.  We should stop when the remainder is so large
 218             that no additional subtracting could make it spill.
 219          4. If nothing else, we should not do two loops of submul_1 over the
 220             data, instead handle both the triangularization and chopping at
 221             once.  */
 222
 223       x = n1;
 224
 225       if (dn > 2)
 226         {
 227           /* Compensate for triangularization.  */
 228           mp_limb_t y;
 229
 230           dp = dp_orig;
 231           if (qn + 1 < dn)
 232             {
 233               dp += dn - (qn + 1);
 234               dn = qn + 1;
 235             }
 236
 237           y = np[-2];
 238
 239           for (i = dn - 3; i >= 0; i--)
 240             {
 241               q = qp[i];
 242               cy = mpn_submul_1 (np - (dn - i), dp, dn - i - 2, q);
 243
 244               if (y < cy)
 245                 {
 246                   if (x == 0)
 247                     {
 248                       cy = mpn_sub_1 (qp, qp, qn, 1);
 249                       ASSERT_ALWAYS (cy == 0);
 250                       return qh - cy;
 251                     }
 252                   x--;
 253                 }
 254               y -= cy;
 255             }
 256           np[-2] = y;
 257         }
 258
 259       dn = dn_orig;
 260       if (qn + 1 < dn)
 261         {
 262           /* Compensate for ignored dividend and divisor tails.  */
 263
 264           dp = dp_orig;
 265           np = np_orig;
 266
 267           if (qh != 0)
 268             {
 269               cy = mpn_sub_n (np + qn, np + qn, dp, dn - (qn + 1));
 270               if (cy != 0)
 271                 {
 272                   if (x == 0)
 273                     {
 274                       if (qn != 0)
 275                         cy = mpn_sub_1 (qp, qp, qn, 1);
 276                       return qh - cy;
 277                     }
 278                   x--;
 279                 }
 280             }
 281
 282           if (qn == 0)
 283             return qh;
 284
 285           for (i = dn - qn - 2; i >= 0; i--)
 286             {
 287               cy = mpn_submul_1 (np + i, qp, qn, dp[i]);
 288               cy = mpn_sub_1 (np + qn + i, np + qn + i, dn - qn - i - 1, cy);
 289               if (cy != 0)
 290                 {
 291                   if (x == 0)
 292                     {
 293                       cy = mpn_sub_1 (qp, qp, qn, 1);
 294                       return qh;
 295                     }
 296                   x--;
 297                 }
 298             }
 299         }
 300     }
 301
 302   return qh;
 303 }