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 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
  29 #include "gmp.h"
  30 #include "gmp-impl.h"
  31 #include "longlong.h"
  32
  33 mp_limb_t
  34 mpn_sbpi1_divappr_q (mp_ptr qp,
  35                      mp_ptr np, mp_size_t nn,
  36                      mp_srcptr dp, mp_size_t dn,
  37                      mp_limb_t dinv)
  38 {
  39   mp_limb_t qh;
  40   mp_size_t qn, i;
  41   mp_limb_t n1, n0;
  42   mp_limb_t d1, d0;
  43   mp_limb_t cy, cy1;
  44   mp_limb_t q;
  45   mp_limb_t flag;
  46
  47   ASSERT (dn > 2);
  48   ASSERT (nn >= dn);
  49   ASSERT ((dp[dn-1] & GMP_NUMB_HIGHBIT) != 0);
  50
  51   np += nn;
  52
  53   qn = nn - dn;
  54   if (qn + 1 < dn)
  55     {
  56       dp += dn - (qn + 1);
  57       dn = qn + 1;
  58     }
  59
  60   qh = mpn_cmp (np - dn, dp, dn) >= 0;
  61   if (qh != 0)
  62     mpn_sub_n (np - dn, np - dn, dp, dn);
  63
  64   qp += qn;
  65
  66   dn -= 2;                      /* offset dn by 2 for main division loops,
  67                                    saving two iterations in mpn_submul_1.  */
  68   d1 = dp[dn + 1];
  69   d0 = dp[dn + 0];
  70
  71   np -= 2;
  72
  73   n1 = np[1];
  74
  75   for (i = qn - (dn + 2); i >= 0; i--)
  76     {
  77       np--;
  78       if (UNLIKELY (n1 == d1) && np[1] == d0)
  79         {
  80           q = GMP_NUMB_MASK;
  81           mpn_submul_1 (np - dn, dp, dn + 2, q);
  82           n1 = np[1];           /* update n1, last loop's value will now be invalid */
  83         }
  84       else
  85         {
  86           udiv_qr_3by2 (q, n1, n0, n1, np[1], np[0], d1, d0, dinv);
  87
  88           cy = mpn_submul_1 (np - dn, dp, dn, q);
  89
  90           cy1 = n0 < cy;
  91           n0 = (n0 - cy) & GMP_NUMB_MASK;
  92           cy = n1 < cy1;
  93           n1 -= cy1;
  94           np[0] = n0;
  95
  96           if (UNLIKELY (cy != 0))
  97             {
  98               n1 += d1 + mpn_add_n (np - dn, np - dn, dp, dn + 1);
  99               q--;
 100             }
 101         }
 102
 103       *--qp = q;
 104     }
 105
 106   flag = ~CNST_LIMB(0);
 107
 108   if (dn >= 0)
 109     {
 110       for (i = dn; i > 0; i--)
 111         {
 112           np--;
 113           if (UNLIKELY (n1 >= (d1 & flag)))
 114             {
 115               q = GMP_NUMB_MASK;
 116               cy = mpn_submul_1 (np - dn, dp, dn + 2, q);
 117
 118               if (UNLIKELY (n1 != cy))
 119                 {
 120                   if (n1 < (cy & flag))
 121                     {
 122                       q--;
 123                       mpn_add_n (np - dn, np - dn, dp, dn + 2);
 124                     }
 125                   else
 126                     flag = 0;
 127                 }
 128               n1 = np[1];
 129             }
 130           else
 131             {
 132               udiv_qr_3by2 (q, n1, n0, n1, np[1], np[0], d1, d0, dinv);
 133
 134               cy = mpn_submul_1 (np - dn, dp, dn, q);
 135
 136               cy1 = n0 < cy;
 137               n0 = (n0 - cy) & GMP_NUMB_MASK;
 138               cy = n1 < cy1;
 139               n1 -= cy1;
 140               np[0] = n0;
 141
 142               if (UNLIKELY (cy != 0))
 143                 {
 144                   n1 += d1 + mpn_add_n (np - dn, np - dn, dp, dn + 1);
 145                   q--;
 146                 }
 147             }
 148
 149           *--qp = q;
 150
 151           /* Truncate operands.  */
 152           dn--;
 153           dp++;
 154         }
 155
 156       np--;
 157       if (UNLIKELY (n1 >= (d1 & flag)))
 158         {
 159           q = GMP_NUMB_MASK;
 160           cy = mpn_submul_1 (np, dp, 2, q);
 161
 162           if (UNLIKELY (n1 != cy))
 163             {
 164               if (n1 < (cy & flag))
 165                 {
 166                   q--;
 167                   add_ssaaaa (np[1], np[0], np[1], np[0], dp[1], dp[0]);
 168                 }
 169               else
 170                 flag = 0;
 171             }
 172           n1 = np[1];
 173         }
 174       else
 175         {
 176           udiv_qr_3by2 (q, n1, n0, n1, np[1], np[0], d1, d0, dinv);
 177
 178           np[1] = n1;
 179           np[0] = n0;
 180         }
 181
 182       *--qp = q;
 183     }
 184
 185   ASSERT_ALWAYS (np[1] == n1);
 186
 187   return qh;
 188 }