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

   1 /* mpn_sbpi1_divappr_q -- Schoolbook division using the Möller-Granlund 3/2
   2    division algorithm, returning approximate quotient.  The quotient returned
   3    is either correct, or one too large.
   4
   5    Contributed to the GNU project by Torbjorn Granlund.
   6
   7    THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE.  IT IS ONLY
   8    SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
   9    GUARANTEED THAT IT WILL CHANGE OR DISAPPEAR IN A FUTURE GMP RELEASE.
  10
  11 Copyright 2007, 2009 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 either:
  17
  18   * the GNU Lesser General Public License as published by the Free
  19     Software Foundation; either version 3 of the License, or (at your
  20     option) any later version.
  21
  22 or
  23
  24   * the GNU General Public License as published by the Free Software
  25     Foundation; either version 2 of the License, or (at your option) any
  26     later version.
  27
  28 or both in parallel, as here.
  29
  30 The GNU MP Library is distributed in the hope that it will be useful, but
  31 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  32 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  33 for more details.
  34
  35 You should have received copies of the GNU General Public License and the
  36 GNU Lesser General Public License along with the GNU MP Library.  If not,
  37 see https://www.gnu.org/licenses/.  */
  38
  39
  40 #include "gmp.h"
  41 #include "gmp-impl.h"
  42 #include "longlong.h"
  43
  44 mp_limb_t
  45 mpn_sbpi1_divappr_q (mp_ptr qp,
  46                      mp_ptr np, mp_size_t nn,
  47                      mp_srcptr dp, mp_size_t dn,
  48                      mp_limb_t dinv)
  49 {
  50   mp_limb_t qh;
  51   mp_size_t qn, i;
  52   mp_limb_t n1, n0;
  53   mp_limb_t d1, d0;
  54   mp_limb_t cy, cy1;
  55   mp_limb_t q;
  56   mp_limb_t flag;
  57
  58   ASSERT (dn > 2);
  59   ASSERT (nn >= dn);
  60   ASSERT ((dp[dn-1] & GMP_NUMB_HIGHBIT) != 0);
  61
  62   np += nn;
  63
  64   qn = nn - dn;
  65   if (qn + 1 < dn)
  66     {
  67       dp += dn - (qn + 1);
  68       dn = qn + 1;
  69     }
  70
  71   qh = mpn_cmp (np - dn, dp, dn) >= 0;
  72   if (qh != 0)
  73     mpn_sub_n (np - dn, np - dn, dp, dn);
  74
  75   qp += qn;
  76
  77   dn -= 2;                      /* offset dn by 2 for main division loops,
  78                                    saving two iterations in mpn_submul_1.  */
  79   d1 = dp[dn + 1];
  80   d0 = dp[dn + 0];
  81
  82   np -= 2;
  83
  84   n1 = np[1];
  85
  86   for (i = qn - (dn + 2); i >= 0; i--)
  87     {
  88       np--;
  89       if (UNLIKELY (n1 == d1) && np[1] == d0)
  90         {
  91           q = GMP_NUMB_MASK;
  92           mpn_submul_1 (np - dn, dp, dn + 2, q);
  93           n1 = np[1];           /* update n1, last loop's value will now be invalid */
  94         }
  95       else
  96         {
  97           udiv_qr_3by2 (q, n1, n0, n1, np[1], np[0], d1, d0, dinv);
  98
  99           cy = mpn_submul_1 (np - dn, dp, dn, q);
 100
 101           cy1 = n0 < cy;
 102           n0 = (n0 - cy) & GMP_NUMB_MASK;
 103           cy = n1 < cy1;
 104           n1 -= cy1;
 105           np[0] = n0;
 106
 107           if (UNLIKELY (cy != 0))
 108             {
 109               n1 += d1 + mpn_add_n (np - dn, np - dn, dp, dn + 1);
 110               q--;
 111             }
 112         }
 113
 114       *--qp = q;
 115     }
 116
 117   flag = ~CNST_LIMB(0);
 118
 119   if (dn >= 0)
 120     {
 121       for (i = dn; i > 0; i--)
 122         {
 123           np--;
 124           if (UNLIKELY (n1 >= (d1 & flag)))
 125             {
 126               q = GMP_NUMB_MASK;
 127               cy = mpn_submul_1 (np - dn, dp, dn + 2, q);
 128
 129               if (UNLIKELY (n1 != cy))
 130                 {
 131                   if (n1 < (cy & flag))
 132                     {
 133                       q--;
 134                       mpn_add_n (np - dn, np - dn, dp, dn + 2);
 135                     }
 136                   else
 137                     flag = 0;
 138                 }
 139               n1 = np[1];
 140             }
 141           else
 142             {
 143               udiv_qr_3by2 (q, n1, n0, n1, np[1], np[0], d1, d0, dinv);
 144
 145               cy = mpn_submul_1 (np - dn, dp, dn, q);
 146
 147               cy1 = n0 < cy;
 148               n0 = (n0 - cy) & GMP_NUMB_MASK;
 149               cy = n1 < cy1;
 150               n1 -= cy1;
 151               np[0] = n0;
 152
 153               if (UNLIKELY (cy != 0))
 154                 {
 155                   n1 += d1 + mpn_add_n (np - dn, np - dn, dp, dn + 1);
 156                   q--;
 157                 }
 158             }
 159
 160           *--qp = q;
 161
 162           /* Truncate operands.  */
 163           dn--;
 164           dp++;
 165         }
 166
 167       np--;
 168       if (UNLIKELY (n1 >= (d1 & flag)))
 169         {
 170           q = GMP_NUMB_MASK;
 171           cy = mpn_submul_1 (np, dp, 2, q);
 172
 173           if (UNLIKELY (n1 != cy))
 174             {
 175               if (n1 < (cy & flag))
 176                 {
 177                   q--;
 178                   add_ssaaaa (np[1], np[0], np[1], np[0], dp[1], dp[0]);
 179                 }
 180               else
 181                 flag = 0;
 182             }
 183           n1 = np[1];
 184         }
 185       else
 186         {
 187           udiv_qr_3by2 (q, n1, n0, n1, np[1], np[0], d1, d0, dinv);
 188
 189           np[1] = n1;
 190           np[0] = n0;
 191         }
 192
 193       *--qp = q;
 194     }
 195
 196   ASSERT_ALWAYS (np[1] == n1);
 197
 198   return qh;
 199 }