contrib/gmp/mpn/generic/mu_div_qr.c

   1 /* mpn_mu_div_qr, mpn_preinv_mu_div_qr.
   2
   3    Compute Q = floor(N / D) and R = N-QD.  N is nn limbs and D is dn limbs and
   4    must be normalized, and Q must be nn-dn limbs.  The requirement that Q is
   5    nn-dn limbs (and not nn-dn+1 limbs) was put in place in order to allow us to
   6    let N be unmodified during the operation.
   7
   8    Contributed to the GNU project by Torbjorn Granlund.
   9
  10    THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES.  IT IS ONLY
  11    SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
  12    GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GMP RELEASE.
  13
  14 Copyright 2005, 2006, 2007, 2009, 2010 Free Software Foundation, Inc.
  15
  16 This file is part of the GNU MP Library.
  17
  18 The GNU MP Library is free software; you can redistribute it and/or modify
  19 it under the terms of the GNU Lesser General Public License as published by
  20 the Free Software Foundation; either version 3 of the License, or (at your
  21 option) any later version.
  22
  23 The GNU MP Library is distributed in the hope that it will be useful, but
  24 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  25 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  26 License for more details.
  27
  28 You should have received a copy of the GNU Lesser General Public License
  29 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  30
  31
  32 /*
  33    The idea of the algorithm used herein is to compute a smaller inverted value
  34    than used in the standard Barrett algorithm, and thus save time in the
  35    Newton iterations, and pay just a small price when using the inverted value
  36    for developing quotient bits.  This algorithm was presented at ICMS 2006.
  37 */
  38
  39 /* CAUTION: This code and the code in mu_divappr_q.c should be edited in sync.
  40
  41  Things to work on:
  42
  43   * This isn't optimal when the quotient isn't needed, as it might take a lot
  44     of space.  The computation is always needed, though, so there is no time to
  45     save with special code.
  46
  47   * The itch/scratch scheme isn't perhaps such a good idea as it once seemed,
  48     demonstrated by the fact that the mpn_invertappr function's scratch needs
  49     mean that we need to keep a large allocation long after it is needed.
  50     Things are worse as mpn_mul_fft does not accept any scratch parameter,
  51     which means we'll have a large memory hole while in mpn_mul_fft.  In
  52     general, a peak scratch need in the beginning of a function isn't
  53     well-handled by the itch/scratch scheme.
  54 */
  55
  56 #ifdef STAT
  57 #undef STAT
  58 #define STAT(x) x
  59 #else
  60 #define STAT(x)
  61 #endif
  62
  63 #include <stdlib.h>             /* for NULL */
  64 #include "gmp.h"
  65 #include "gmp-impl.h"
  66
  67
  68 /* FIXME: The MU_DIV_QR_SKEW_THRESHOLD was not analysed properly.  It gives a
  69    speedup according to old measurements, but does the decision mechanism
  70    really make sense?  It seem like the quotient between dn and qn might be
  71    what we really should be checking.  */
  72 #ifndef MU_DIV_QR_SKEW_THRESHOLD
  73 #define MU_DIV_QR_SKEW_THRESHOLD 100
  74 #endif
  75
  76 #ifdef CHECK                            /* FIXME: Enable in minithres */
  77 #undef  MU_DIV_QR_SKEW_THRESHOLD
  78 #define MU_DIV_QR_SKEW_THRESHOLD 1
  79 #endif
  80
  81
  82 static mp_limb_t mpn_mu_div_qr2 (mp_ptr, mp_ptr, mp_srcptr, mp_size_t, mp_srcptr, mp_size_t, mp_ptr);
  83
  84
  85 mp_limb_t
  86 mpn_mu_div_qr (mp_ptr qp,
  87                mp_ptr rp,
  88                mp_srcptr np,
  89                mp_size_t nn,
  90                mp_srcptr dp,
  91                mp_size_t dn,
  92                mp_ptr scratch)
  93 {
  94   mp_size_t qn;
  95   mp_limb_t cy, qh;
  96
  97   qn = nn - dn;
  98   if (qn + MU_DIV_QR_SKEW_THRESHOLD < dn)
  99     {
 100       /* |______________|_ign_first__|   dividend                         nn
 101                 |_______|_ign_first__|   divisor                          dn
 102
 103                 |______|             quotient (prel)                      qn
 104
 105                  |___________________|   quotient * ignored-divisor-part  dn-1
 106       */
 107
 108       /* Compute a preliminary quotient and a partial remainder by dividing the
 109          most significant limbs of each operand.  */
 110       qh = mpn_mu_div_qr2 (qp, rp + nn - (2 * qn + 1),
 111                            np + nn - (2 * qn + 1), 2 * qn + 1,
 112                            dp + dn - (qn + 1), qn + 1,
 113                            scratch);
 114
 115       /* Multiply the quotient by the divisor limbs ignored above.  */
 116       if (dn - (qn + 1) > qn)
 117         mpn_mul (scratch, dp, dn - (qn + 1), qp, qn);  /* prod is dn-1 limbs */
 118       else
 119         mpn_mul (scratch, qp, qn, dp, dn - (qn + 1));  /* prod is dn-1 limbs */
 120
 121       if (qh)
 122         cy = mpn_add_n (scratch + qn, scratch + qn, dp, dn - (qn + 1));
 123       else
 124         cy = 0;
 125       scratch[dn - 1] = cy;
 126
 127       cy = mpn_sub_n (rp, np, scratch, nn - (2 * qn + 1));
 128       cy = mpn_sub_nc (rp + nn - (2 * qn + 1),
 129                        rp + nn - (2 * qn + 1),
 130                        scratch + nn - (2 * qn + 1),
 131                        qn + 1, cy);
 132       if (cy)
 133         {
 134           qh -= mpn_sub_1 (qp, qp, qn, 1);
 135           mpn_add_n (rp, rp, dp, dn);
 136         }
 137     }
 138   else
 139     {
 140       qh = mpn_mu_div_qr2 (qp, rp, np, nn, dp, dn, scratch);
 141     }
 142
 143   return qh;
 144 }
 145
 146 static mp_limb_t
 147 mpn_mu_div_qr2 (mp_ptr qp,
 148                 mp_ptr rp,
 149                 mp_srcptr np,
 150                 mp_size_t nn,
 151                 mp_srcptr dp,
 152                 mp_size_t dn,
 153                 mp_ptr scratch)
 154 {
 155   mp_size_t qn, in;
 156   mp_limb_t cy, qh;
 157   mp_ptr ip, tp;
 158
 159   ASSERT (dn > 1);
 160
 161   qn = nn - dn;
 162
 163   /* Compute the inverse size.  */
 164   in = mpn_mu_div_qr_choose_in (qn, dn, 0);
 165   ASSERT (in <= dn);
 166
 167 #if 1
 168   /* This alternative inverse computation method gets slightly more accurate
 169      results.  FIXMEs: (1) Temp allocation needs not analysed (2) itch function
 170      not adapted (3) mpn_invertappr scratch needs not met.  */
 171   ip = scratch;
 172   tp = scratch + in + 1;
 173
 174   /* compute an approximate inverse on (in+1) limbs */
 175   if (dn == in)
 176     {
 177       MPN_COPY (tp + 1, dp, in);
 178       tp[0] = 1;
 179       mpn_invertappr (ip, tp, in + 1, NULL);
 180       MPN_COPY_INCR (ip, ip + 1, in);
 181     }
 182   else
 183     {
 184       cy = mpn_add_1 (tp, dp + dn - (in + 1), in + 1, 1);
 185       if (UNLIKELY (cy != 0))
 186         MPN_ZERO (ip, in);
 187       else
 188         {
 189           mpn_invertappr (ip, tp, in + 1, NULL);
 190           MPN_COPY_INCR (ip, ip + 1, in);
 191         }
 192     }
 193 #else
 194   /* This older inverse computation method gets slightly worse results than the
 195      one above.  */
 196   ip = scratch;
 197   tp = scratch + in;
 198
 199   /* Compute inverse of D to in+1 limbs, then round to 'in' limbs.  Ideally the
 200      inversion function should do this automatically.  */
 201   if (dn == in)
 202     {
 203       tp[in + 1] = 0;
 204       MPN_COPY (tp + in + 2, dp, in);
 205       mpn_invertappr (tp, tp + in + 1, in + 1, NULL);
 206     }
 207   else
 208     {
 209       mpn_invertappr (tp, dp + dn - (in + 1), in + 1, NULL);
 210     }
 211   cy = mpn_sub_1 (tp, tp, in + 1, GMP_NUMB_HIGHBIT);
 212   if (UNLIKELY (cy != 0))
 213     MPN_ZERO (tp + 1, in);
 214   MPN_COPY (ip, tp + 1, in);
 215 #endif
 216
 217   qh = mpn_preinv_mu_div_qr (qp, rp, np, nn, dp, dn, ip, in, scratch + in);
 218
 219   return qh;
 220 }
 221
 222 mp_limb_t
 223 mpn_preinv_mu_div_qr (mp_ptr qp,
 224                       mp_ptr rp,
 225                       mp_srcptr np,
 226                       mp_size_t nn,
 227                       mp_srcptr dp,
 228                       mp_size_t dn,
 229                       mp_srcptr ip,
 230                       mp_size_t in,
 231                       mp_ptr scratch)
 232 {
 233   mp_size_t qn;
 234   mp_limb_t cy, cx, qh;
 235   mp_limb_t r;
 236   mp_size_t tn, wn;
 237
 238 #define tp           scratch
 239 #define scratch_out  (scratch + tn)
 240
 241   qn = nn - dn;
 242
 243   np += qn;
 244   qp += qn;
 245
 246   qh = mpn_cmp (np, dp, dn) >= 0;
 247   if (qh != 0)
 248     mpn_sub_n (rp, np, dp, dn);
 249   else
 250     MPN_COPY (rp, np, dn);
 251
 252   if (qn == 0)
 253     return qh;                  /* Degenerate use.  Should we allow this? */
 254
 255   while (qn > 0)
 256     {
 257       if (qn < in)
 258         {
 259           ip += in - qn;
 260           in = qn;
 261         }
 262       np -= in;
 263       qp -= in;
 264
 265       /* Compute the next block of quotient limbs by multiplying the inverse I
 266          by the upper part of the partial remainder R.  */
 267       mpn_mul_n (tp, rp + dn - in, ip, in);             /* mulhi  */
 268       cy = mpn_add_n (qp, tp + in, rp + dn - in, in);   /* I's msb implicit */
 269       ASSERT_ALWAYS (cy == 0);
 270
 271       qn -= in;
 272
 273       /* Compute the product of the quotient block and the divisor D, to be
 274          subtracted from the partial remainder combined with new limbs from the
 275          dividend N.  We only really need the low dn+1 limbs.  */
 276
 277       if (BELOW_THRESHOLD (in, MUL_TO_MULMOD_BNM1_FOR_2NXN_THRESHOLD))
 278         mpn_mul (tp, dp, dn, qp, in);           /* dn+in limbs, high 'in' cancels */
 279       else
 280         {
 281           tn = mpn_mulmod_bnm1_next_size (dn + 1);
 282           mpn_mulmod_bnm1 (tp, tn, dp, dn, qp, in, scratch_out);
 283           wn = dn + in - tn;                    /* number of wrapped limbs */
 284           if (wn > 0)
 285             {
 286               cy = mpn_sub_n (tp, tp, rp + dn - wn, wn);
 287               cy = mpn_sub_1 (tp + wn, tp + wn, tn - wn, cy);
 288               cx = mpn_cmp (rp + dn - in, tp + dn, tn - dn) < 0;
 289               ASSERT_ALWAYS (cx >= cy);
 290               mpn_incr_u (tp, cx - cy);
 291             }
 292         }
 293
 294       r = rp[dn - in] - tp[dn];
 295
 296       /* Subtract the product from the partial remainder combined with new
 297          limbs from the dividend N, generating a new partial remainder R.  */
 298       if (dn != in)
 299         {
 300           cy = mpn_sub_n (tp, np, tp, in);      /* get next 'in' limbs from N */
 301           cy = mpn_sub_nc (tp + in, rp, tp + in, dn - in, cy);
 302           MPN_COPY (rp, tp, dn);                /* FIXME: try to avoid this */
 303         }
 304       else
 305         {
 306           cy = mpn_sub_n (rp, np, tp, in);      /* get next 'in' limbs from N */
 307         }
 308
 309       STAT (int i; int err = 0;
 310             static int errarr[5]; static int err_rec; static int tot);
 311
 312       /* Check the remainder R and adjust the quotient as needed.  */
 313       r -= cy;
 314       while (r != 0)
 315         {
 316           /* We loop 0 times with about 69% probability, 1 time with about 31%
 317              probability, 2 times with about 0.6% probability, if inverse is
 318              computed as recommended.  */
 319           mpn_incr_u (qp, 1);
 320           cy = mpn_sub_n (rp, rp, dp, dn);
 321           r -= cy;
 322           STAT (err++);
 323         }
 324       if (mpn_cmp (rp, dp, dn) >= 0)
 325         {
 326           /* This is executed with about 76% probability.  */
 327           mpn_incr_u (qp, 1);
 328           cy = mpn_sub_n (rp, rp, dp, dn);
 329           STAT (err++);
 330         }
 331
 332       STAT (
 333             tot++;
 334             errarr[err]++;
 335             if (err > err_rec)
 336               err_rec = err;
 337             if (tot % 0x10000 == 0)
 338               {
 339                 for (i = 0; i <= err_rec; i++)
 340                   printf ("  %d(%.1f%%)", errarr[i], 100.0*errarr[i]/tot);
 341                 printf ("\n");
 342               }
 343             );
 344     }
 345
 346   return qh;
 347 }
 348
 349 /* In case k=0 (automatic choice), we distinguish 3 cases:
 350    (a) dn < qn:         in = ceil(qn / ceil(qn/dn))
 351    (b) dn/3 < qn <= dn: in = ceil(qn / 2)
 352    (c) qn < dn/3:       in = qn
 353    In all cases we have in <= dn.
 354  */
 355 mp_size_t
 356 mpn_mu_div_qr_choose_in (mp_size_t qn, mp_size_t dn, int k)
 357 {
 358   mp_size_t in;
 359
 360   if (k == 0)
 361     {
 362       mp_size_t b;
 363       if (qn > dn)
 364         {
 365           /* Compute an inverse size that is a nice partition of the quotient.  */
 366           b = (qn - 1) / dn + 1;        /* ceil(qn/dn), number of blocks */
 367           in = (qn - 1) / b + 1;        /* ceil(qn/b) = ceil(qn / ceil(qn/dn)) */
 368         }
 369       else if (3 * qn > dn)
 370         {
 371           in = (qn - 1) / 2 + 1;        /* b = 2 */
 372         }
 373       else
 374         {
 375           in = (qn - 1) / 1 + 1;        /* b = 1 */
 376         }
 377     }
 378   else
 379     {
 380       mp_size_t xn;
 381       xn = MIN (dn, qn);
 382       in = (xn - 1) / k + 1;
 383     }
 384
 385   return in;
 386 }
 387
 388 mp_size_t
 389 mpn_mu_div_qr_itch (mp_size_t nn, mp_size_t dn, int mua_k)
 390 {
 391   mp_size_t itch_local = mpn_mulmod_bnm1_next_size (dn + 1);
 392   mp_size_t in = mpn_mu_div_qr_choose_in (nn - dn, dn, mua_k);
 393   mp_size_t itch_out = mpn_mulmod_bnm1_itch (itch_local, dn, in);
 394
 395   return in + itch_local + itch_out;
 396 }
 397
 398 mp_size_t
 399 mpn_preinv_mu_div_qr_itch (mp_size_t nn, mp_size_t dn, mp_size_t in)
 400 {
 401   mp_size_t itch_local = mpn_mulmod_bnm1_next_size (dn + 1);
 402   mp_size_t itch_out = mpn_mulmod_bnm1_itch (itch_local, dn, in);
 403
 404   return itch_local + itch_out;
 405 }