contrib/gmp/mpn/generic/mu_bdiv_qr.c

   1 /* mpn_mu_bdiv_qr(qp,rp,np,nn,dp,dn,tp) -- Compute {np,nn} / {dp,dn} mod B^qn,
   2    where qn = nn-dn, storing the result in {qp,qn}.  Overlap allowed between Q
   3    and N; all other overlap disallowed.
   4
   5    Contributed to the GNU project by Torbjorn Granlund.
   6
   7    THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES.  IT IS ONLY
   8    SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
   9    GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GMP RELEASE.
  10
  11 Copyright 2005, 2006, 2007, 2009, 2010 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 /*
  30    The idea of the algorithm used herein is to compute a smaller inverted value
  31    than used in the standard Barrett algorithm, and thus save time in the
  32    Newton iterations, and pay just a small price when using the inverted value
  33    for developing quotient bits.  This algorithm was presented at ICMS 2006.
  34 */
  35
  36 #include "gmp.h"
  37 #include "gmp-impl.h"
  38
  39
  40 /* N = {np,nn}
  41    D = {dp,dn}
  42
  43    Requirements: N >= D
  44                  D >= 1
  45                  D odd
  46                  dn >= 2
  47                  nn >= 2
  48                  scratch space as determined by mpn_mu_bdiv_qr_itch(nn,dn).
  49
  50    Write quotient to Q = {qp,nn-dn}.
  51
  52    FIXME: When iterating, perhaps do the small step before loop, not after.
  53    FIXME: Try to avoid the scalar divisions when computing inverse size.
  54    FIXME: Trim allocation for (qn > dn) case, 3*dn might be possible.  In
  55           particular, when dn==in, tp and rp could use the same space.
  56 */
  57 mp_limb_t
  58 mpn_mu_bdiv_qr (mp_ptr qp,
  59                 mp_ptr rp,
  60                 mp_srcptr np, mp_size_t nn,
  61                 mp_srcptr dp, mp_size_t dn,
  62                 mp_ptr scratch)
  63 {
  64   mp_size_t qn;
  65   mp_size_t in;
  66   mp_limb_t cy, c0;
  67   int k;
  68   mp_size_t tn, wn;
  69   mp_size_t i;
  70
  71   qn = nn - dn;
  72
  73   ASSERT (dn >= 2);
  74   ASSERT (qn >= 2);
  75
  76   if (qn > dn)
  77     {
  78       mp_size_t b;
  79
  80       /* |_______________________|   dividend
  81                         |________|   divisor  */
  82
  83 #define ip           scratch            /* in */
  84 #define tp           (scratch + in)     /* dn+in or next_size(dn) or rest >= binvert_itch(in) */
  85 #define scratch_out  (scratch + in + tn)/* mulmod_bnm1_itch(next_size(dn)) */
  86
  87       /* Compute an inverse size that is a nice partition of the quotient.  */
  88       b = (qn - 1) / dn + 1;    /* ceil(qn/dn), number of blocks */
  89       in = (qn - 1) / b + 1;    /* ceil(qn/b) = ceil(qn / ceil(qn/dn)) */
  90
  91       /* Some notes on allocation:
  92
  93          When in = dn, R dies when mpn_mullo returns, if in < dn the low in
  94          limbs of R dies at that point.  We could save memory by letting T live
  95          just under R, and let the upper part of T expand into R. These changes
  96          should reduce itch to perhaps 3dn.
  97        */
  98
  99       mpn_binvert (ip, dp, in, tp);
 100
 101       MPN_COPY (rp, np, dn);
 102       np += dn;
 103       cy = 0;
 104
 105       while (qn > in)
 106         {
 107           mpn_mullo_n (qp, rp, ip, in);
 108
 109           if (BELOW_THRESHOLD (in, MUL_TO_MULMOD_BNM1_FOR_2NXN_THRESHOLD))
 110             mpn_mul (tp, dp, dn, qp, in);       /* mulhi, need tp[dn+in-1...in] */
 111           else
 112             {
 113               tn = mpn_mulmod_bnm1_next_size (dn);
 114               mpn_mulmod_bnm1 (tp, tn, dp, dn, qp, in, scratch_out);
 115               wn = dn + in - tn;                /* number of wrapped limbs */
 116               if (wn > 0)
 117                 {
 118                   c0 = mpn_sub_n (tp + tn, tp, rp, wn);
 119                   mpn_decr_u (tp + wn, c0);
 120                 }
 121             }
 122
 123           qp += in;
 124           qn -= in;
 125
 126           if (dn != in)
 127             {
 128               /* Subtract tp[dn-1...in] from partial remainder.  */
 129               cy += mpn_sub_n (rp, rp + in, tp + in, dn - in);
 130               if (cy == 2)
 131                 {
 132                   mpn_incr_u (tp + dn, 1);
 133                   cy = 1;
 134                 }
 135             }
 136           /* Subtract tp[dn+in-1...dn] from dividend.  */
 137           cy = mpn_sub_nc (rp + dn - in, np, tp + dn, in, cy);
 138           np += in;
 139         }
 140
 141       /* Generate last qn limbs.  */
 142       mpn_mullo_n (qp, rp, ip, qn);
 143
 144       if (BELOW_THRESHOLD (qn, MUL_TO_MULMOD_BNM1_FOR_2NXN_THRESHOLD))
 145         mpn_mul (tp, dp, dn, qp, qn);           /* mulhi, need tp[qn+in-1...in] */
 146       else
 147         {
 148           tn = mpn_mulmod_bnm1_next_size (dn);
 149           mpn_mulmod_bnm1 (tp, tn, dp, dn, qp, qn, scratch_out);
 150           wn = dn + qn - tn;                    /* number of wrapped limbs */
 151           if (wn > 0)
 152             {
 153               c0 = mpn_sub_n (tp + tn, tp, rp, wn);
 154               mpn_decr_u (tp + wn, c0);
 155             }
 156         }
 157
 158       if (dn != qn)
 159         {
 160           cy += mpn_sub_n (rp, rp + qn, tp + qn, dn - qn);
 161           if (cy == 2)
 162             {
 163               mpn_incr_u (tp + dn, 1);
 164               cy = 1;
 165             }
 166         }
 167       return mpn_sub_nc (rp + dn - qn, np, tp + dn, qn, cy);
 168
 169 #undef ip
 170 #undef tp
 171 #undef scratch_out
 172     }
 173   else
 174     {
 175       /* |_______________________|   dividend
 176                 |________________|   divisor  */
 177
 178 #define ip           scratch            /* in */
 179 #define tp           (scratch + in)     /* dn+in or next_size(dn) or rest >= binvert_itch(in) */
 180 #define scratch_out  (scratch + in + tn)/* mulmod_bnm1_itch(next_size(dn)) */
 181
 182       /* Compute half-sized inverse.  */
 183       in = qn - (qn >> 1);
 184
 185       mpn_binvert (ip, dp, in, tp);
 186
 187       mpn_mullo_n (qp, np, ip, in);             /* low `in' quotient limbs */
 188
 189       if (BELOW_THRESHOLD (in, MUL_TO_MULMOD_BNM1_FOR_2NXN_THRESHOLD))
 190         mpn_mul (tp, dp, dn, qp, in);           /* mulhigh */
 191       else
 192         {
 193           tn = mpn_mulmod_bnm1_next_size (dn);
 194           mpn_mulmod_bnm1 (tp, tn, dp, dn, qp, in, scratch_out);
 195           wn = dn + in - tn;                    /* number of wrapped limbs */
 196           if (wn > 0)
 197             {
 198               c0 = mpn_sub_n (tp + tn, tp, np, wn);
 199               mpn_decr_u (tp + wn, c0);
 200             }
 201         }
 202
 203       qp += in;
 204       qn -= in;
 205
 206       cy = mpn_sub_n (rp, np + in, tp + in, dn);
 207       mpn_mullo_n (qp, rp, ip, qn);             /* high qn quotient limbs */
 208
 209       if (BELOW_THRESHOLD (qn, MUL_TO_MULMOD_BNM1_FOR_2NXN_THRESHOLD))
 210         mpn_mul (tp, dp, dn, qp, qn);           /* mulhigh */
 211       else
 212         {
 213           tn = mpn_mulmod_bnm1_next_size (dn);
 214           mpn_mulmod_bnm1 (tp, tn, dp, dn, qp, qn, scratch_out);
 215           wn = dn + qn - tn;                    /* number of wrapped limbs */
 216           if (wn > 0)
 217             {
 218               c0 = mpn_sub_n (tp + tn, tp, rp, wn);
 219               mpn_decr_u (tp + wn, c0);
 220             }
 221         }
 222
 223       cy += mpn_sub_n (rp, rp + qn, tp + qn, dn - qn);
 224       if (cy == 2)
 225         {
 226           mpn_incr_u (tp + dn, 1);
 227           cy = 1;
 228         }
 229       return mpn_sub_nc (rp + dn - qn, np + dn + in, tp + dn, qn, cy);
 230
 231 #undef ip
 232 #undef tp
 233 #undef scratch_out
 234     }
 235 }
 236
 237 mp_size_t
 238 mpn_mu_bdiv_qr_itch (mp_size_t nn, mp_size_t dn)
 239 {
 240   mp_size_t qn, in, tn, itch_binvert, itch_out, itches;
 241   mp_size_t b;
 242
 243   qn = nn - dn;
 244
 245   if (qn > dn)
 246     {
 247       b = (qn - 1) / dn + 1;    /* ceil(qn/dn), number of blocks */
 248       in = (qn - 1) / b + 1;    /* ceil(qn/b) = ceil(qn / ceil(qn/dn)) */
 249       if (BELOW_THRESHOLD (in, MUL_TO_MULMOD_BNM1_FOR_2NXN_THRESHOLD))
 250         {
 251           tn = dn + in;
 252           itch_out = 0;
 253         }
 254       else
 255         {
 256           tn = mpn_mulmod_bnm1_next_size (dn);
 257           itch_out = mpn_mulmod_bnm1_itch (tn, dn, in);
 258         }
 259       itch_binvert = mpn_binvert_itch (in);
 260       itches = tn + itch_out;
 261       return in + MAX (itches, itch_binvert);
 262     }
 263   else
 264     {
 265       in = qn - (qn >> 1);
 266       if (BELOW_THRESHOLD (in, MUL_TO_MULMOD_BNM1_FOR_2NXN_THRESHOLD))
 267         {
 268           tn = dn + in;
 269           itch_out = 0;
 270         }
 271       else
 272         {
 273           tn = mpn_mulmod_bnm1_next_size (dn);
 274           itch_out = mpn_mulmod_bnm1_itch (tn, dn, in);
 275         }
 276     }
 277   itch_binvert = mpn_binvert_itch (in);
 278   itches = tn + itch_out;
 279   return in + MAX (itches, itch_binvert);
 280 }