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

   1 /* mpn_dcpi1_div_qr_n -- recursive divide-and-conquer division for arbitrary
   2    size operands.
   3
   4    Contributed to the GNU project by Torbjorn Granlund.
   5
   6    THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES.  IT IS ONLY
   7    SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
   8    GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GMP RELEASE.
   9
  10 Copyright 2006, 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 #include "gmp.h"
  39 #include "gmp-impl.h"
  40 #include "longlong.h"
  41
  42
  43 mp_limb_t
  44 mpn_dcpi1_div_qr_n (mp_ptr qp, mp_ptr np, mp_srcptr dp, mp_size_t n,
  45                     gmp_pi1_t *dinv, mp_ptr tp)
  46 {
  47   mp_size_t lo, hi;
  48   mp_limb_t cy, qh, ql;
  49
  50   lo = n >> 1;                  /* floor(n/2) */
  51   hi = n - lo;                  /* ceil(n/2) */
  52
  53   if (BELOW_THRESHOLD (hi, DC_DIV_QR_THRESHOLD))
  54     qh = mpn_sbpi1_div_qr (qp + lo, np + 2 * lo, 2 * hi, dp + lo, hi, dinv->inv32);
  55   else
  56     qh = mpn_dcpi1_div_qr_n (qp + lo, np + 2 * lo, dp + lo, hi, dinv, tp);
  57
  58   mpn_mul (tp, qp + lo, hi, dp, lo);
  59
  60   cy = mpn_sub_n (np + lo, np + lo, tp, n);
  61   if (qh != 0)
  62     cy += mpn_sub_n (np + n, np + n, dp, lo);
  63
  64   while (cy != 0)
  65     {
  66       qh -= mpn_sub_1 (qp + lo, qp + lo, hi, 1);
  67       cy -= mpn_add_n (np + lo, np + lo, dp, n);
  68     }
  69
  70   if (BELOW_THRESHOLD (lo, DC_DIV_QR_THRESHOLD))
  71     ql = mpn_sbpi1_div_qr (qp, np + hi, 2 * lo, dp + hi, lo, dinv->inv32);
  72   else
  73     ql = mpn_dcpi1_div_qr_n (qp, np + hi, dp + hi, lo, dinv, tp);
  74
  75   mpn_mul (tp, dp, hi, qp, lo);
  76
  77   cy = mpn_sub_n (np, np, tp, n);
  78   if (ql != 0)
  79     cy += mpn_sub_n (np + lo, np + lo, dp, hi);
  80
  81   while (cy != 0)
  82     {
  83       mpn_sub_1 (qp, qp, lo, 1);
  84       cy -= mpn_add_n (np, np, dp, n);
  85     }
  86
  87   return qh;
  88 }
  89
  90 mp_limb_t
  91 mpn_dcpi1_div_qr (mp_ptr qp,
  92                   mp_ptr np, mp_size_t nn,
  93                   mp_srcptr dp, mp_size_t dn,
  94                   gmp_pi1_t *dinv)
  95 {
  96   mp_size_t qn;
  97   mp_limb_t qh, cy;
  98   mp_ptr tp;
  99   TMP_DECL;
 100
 101   TMP_MARK;
 102
 103   ASSERT (dn >= 6);             /* to adhere to mpn_sbpi1_div_qr's limits */
 104   ASSERT (nn - dn >= 3);        /* to adhere to mpn_sbpi1_div_qr's limits */
 105   ASSERT (dp[dn-1] & GMP_NUMB_HIGHBIT);
 106
 107   tp = TMP_ALLOC_LIMBS (dn);
 108
 109   qn = nn - dn;
 110   qp += qn;
 111   np += nn;
 112   dp += dn;
 113
 114   if (qn > dn)
 115     {
 116       /* Reduce qn mod dn without division, optimizing small operations.  */
 117       do
 118         qn -= dn;
 119       while (qn > dn);
 120
 121       qp -= qn;                 /* point at low limb of next quotient block */
 122       np -= qn;                 /* point in the middle of partial remainder */
 123
 124       /* Perform the typically smaller block first.  */
 125       if (qn == 1)
 126         {
 127           mp_limb_t q, n2, n1, n0, d1, d0;
 128
 129           /* Handle qh up front, for simplicity. */
 130           qh = mpn_cmp (np - dn + 1, dp - dn, dn) >= 0;
 131           if (qh)
 132             ASSERT_NOCARRY (mpn_sub_n (np - dn + 1, np - dn + 1, dp - dn, dn));
 133
 134           /* A single iteration of schoolbook: One 3/2 division,
 135              followed by the bignum update and adjustment. */
 136           n2 = np[0];
 137           n1 = np[-1];
 138           n0 = np[-2];
 139           d1 = dp[-1];
 140           d0 = dp[-2];
 141
 142           ASSERT (n2 < d1 || (n2 == d1 && n1 <= d0));
 143
 144           if (UNLIKELY (n2 == d1) && n1 == d0)
 145             {
 146               q = GMP_NUMB_MASK;
 147               cy = mpn_submul_1 (np - dn, dp - dn, dn, q);
 148               ASSERT (cy == n2);
 149             }
 150           else
 151             {
 152               udiv_qr_3by2 (q, n1, n0, n2, n1, n0, d1, d0, dinv->inv32);
 153
 154               if (dn > 2)
 155                 {
 156                   mp_limb_t cy, cy1;
 157                   cy = mpn_submul_1 (np - dn, dp - dn, dn - 2, q);
 158
 159                   cy1 = n0 < cy;
 160                   n0 = (n0 - cy) & GMP_NUMB_MASK;
 161                   cy = n1 < cy1;
 162                   n1 = (n1 - cy1) & GMP_NUMB_MASK;
 163                   np[-2] = n0;
 164
 165                   if (UNLIKELY (cy != 0))
 166                     {
 167                       n1 += d1 + mpn_add_n (np - dn, np - dn, dp - dn, dn - 1);
 168                       qh -= (q == 0);
 169                       q = (q - 1) & GMP_NUMB_MASK;
 170                     }
 171                 }
 172               else
 173                 np[-2] = n0;
 174
 175               np[-1] = n1;
 176             }
 177           qp[0] = q;
 178         }
 179       else
 180         {
 181           /* Do a 2qn / qn division */
 182           if (qn == 2)
 183             qh = mpn_divrem_2 (qp, 0L, np - 2, 4, dp - 2); /* FIXME: obsolete function. Use 5/3 division? */
 184           else if (BELOW_THRESHOLD (qn, DC_DIV_QR_THRESHOLD))
 185             qh = mpn_sbpi1_div_qr (qp, np - qn, 2 * qn, dp - qn, qn, dinv->inv32);
 186           else
 187             qh = mpn_dcpi1_div_qr_n (qp, np - qn, dp - qn, qn, dinv, tp);
 188
 189           if (qn != dn)
 190             {
 191               if (qn > dn - qn)
 192                 mpn_mul (tp, qp, qn, dp - dn, dn - qn);
 193               else
 194                 mpn_mul (tp, dp - dn, dn - qn, qp, qn);
 195
 196               cy = mpn_sub_n (np - dn, np - dn, tp, dn);
 197               if (qh != 0)
 198                 cy += mpn_sub_n (np - dn + qn, np - dn + qn, dp - dn, dn - qn);
 199
 200               while (cy != 0)
 201                 {
 202                   qh -= mpn_sub_1 (qp, qp, qn, 1);
 203                   cy -= mpn_add_n (np - dn, np - dn, dp - dn, dn);
 204                 }
 205             }
 206         }
 207
 208       qn = nn - dn - qn;
 209       do
 210         {
 211           qp -= dn;
 212           np -= dn;
 213           mpn_dcpi1_div_qr_n (qp, np - dn, dp - dn, dn, dinv, tp);
 214           qn -= dn;
 215         }
 216       while (qn > 0);
 217     }
 218   else
 219     {
 220       qp -= qn;                 /* point at low limb of next quotient block */
 221       np -= qn;                 /* point in the middle of partial remainder */
 222
 223       if (BELOW_THRESHOLD (qn, DC_DIV_QR_THRESHOLD))
 224         qh = mpn_sbpi1_div_qr (qp, np - qn, 2 * qn, dp - qn, qn, dinv->inv32);
 225       else
 226         qh = mpn_dcpi1_div_qr_n (qp, np - qn, dp - qn, qn, dinv, tp);
 227
 228       if (qn != dn)
 229         {
 230           if (qn > dn - qn)
 231             mpn_mul (tp, qp, qn, dp - dn, dn - qn);
 232           else
 233             mpn_mul (tp, dp - dn, dn - qn, qp, qn);
 234
 235           cy = mpn_sub_n (np - dn, np - dn, tp, dn);
 236           if (qh != 0)
 237             cy += mpn_sub_n (np - dn + qn, np - dn + qn, dp - dn, dn - qn);
 238
 239           while (cy != 0)
 240             {
 241               qh -= mpn_sub_1 (qp, qp, qn, 1);
 242               cy -= mpn_add_n (np - dn, np - dn, dp - dn, dn);
 243             }
 244         }
 245     }
 246
 247   TMP_FREE;
 248   return qh;
 249 }