contrib/gmp/mpn/generic/dc_divappr_q.c

   1 /* mpn_dc_divappr_q -- divide-and-conquer division, returning only approximate
   2    quotient.  The quotient retuened is either correct, or unity too large.
   3
   4    Contributed to the GNU project by Torbjörn Granlund.
   5
   6    THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH A MUTABLE INTERFACE.  IT IS
   7    ONLY SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS
   8    ALMOST GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GMP
   9    RELEASE.
  10
  11 Copyright 2006, 2007 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 #include "gmp.h"
  29 #include "gmp-impl.h"
  30
  31
  32 mp_limb_t
  33 mpn_dc_divappr_q_n (mp_ptr qp, mp_ptr np, mp_srcptr dp, mp_size_t n,
  34                     mp_srcptr dip, mp_ptr tp)
  35 {
  36   mp_size_t lo, hi;
  37   mp_limb_t cy, qh, ql;
  38
  39   lo = n >> 1;                  /* floor(n/2) */
  40   hi = n - lo;                  /* ceil(n/2) */
  41
  42   if (BELOW_THRESHOLD (hi, DC_DIV_QR_THRESHOLD))
  43     qh = mpn_sb_div_qr (qp + lo, np + 2 * lo, 2 * hi, dp + lo, hi, dip);
  44   else
  45     qh = mpn_dc_div_qr_n (qp + lo, np + 2 * lo, dp + lo, hi, dip, tp);
  46
  47   mpn_mul (tp, qp + lo, hi, dp, lo);
  48
  49   cy = mpn_sub_n (np + lo, np + lo, tp, n);
  50   if (qh != 0)
  51     cy += mpn_sub_n (np + n, np + n, dp, lo);
  52
  53   while (cy != 0)
  54     {
  55       qh -= mpn_sub_1 (qp + lo, qp + lo, hi, 1);
  56       cy -= mpn_add_n (np + lo, np + lo, dp, n);
  57     }
  58
  59   if (BELOW_THRESHOLD (lo, DC_DIVAPPR_Q_THRESHOLD))
  60     ql = mpn_sb_divappr_q (qp, np + hi, 2 * lo, dp + hi, lo, dip);
  61   else
  62     ql = mpn_dc_divappr_q_n (qp, np + hi, dp + hi, lo, dip, tp);
  63
  64   if (UNLIKELY (ql != 0))
  65     {
  66       mp_size_t i;
  67       for (i = 0; i < lo; i++)
  68         qp[i] = GMP_NUMB_MASK;
  69     }
  70
  71   return qh;
  72 }
  73
  74 mp_limb_t
  75 mpn_preinv_dc_divappr_q (mp_ptr qp,
  76                          mp_ptr np, mp_size_t nn,
  77                          mp_srcptr dp, mp_size_t dn,
  78                          mp_srcptr dip)
  79 {
  80   mp_size_t qn;
  81   mp_limb_t qh, cy, qsave;
  82   mp_ptr tp;
  83   TMP_DECL;
  84
  85   TMP_MARK;
  86
  87   tp = TMP_SALLOC_LIMBS (dn+1);
  88
  89   qn = nn - dn;
  90   qp += qn;
  91   np += nn;
  92   dp += dn;
  93
  94   if (qn > dn)
  95     {
  96       qn++;                     /* pretend we'll need an extra limb */
  97       /* Reduce qn mod dn without division, optimizing small operations.  */
  98       do
  99         qn -= dn;
 100       while (qn > dn);
 101
 102       qp -= qn;                 /* point at low limb of next quotient block */
 103       np -= qn;                 /* point in the middle of partial remainder */
 104
 105       /* Perform the typically smaller block first.  */
 106       if (BELOW_THRESHOLD (qn, DC_DIV_QR_THRESHOLD))
 107         qh = mpn_sb_div_qr (qp, np - qn, 2 * qn, dp - qn, qn, dip);
 108       else
 109         qh = mpn_dc_div_qr_n (qp, np - qn, dp - qn, qn, dip, tp);
 110
 111       if (qn != dn)
 112         {
 113           if (qn > dn - qn)
 114             mpn_mul (tp, qp, qn, dp - dn, dn - qn);
 115           else
 116             mpn_mul (tp, dp - dn, dn - qn, qp, qn);
 117
 118           cy = mpn_sub_n (np - dn, np - dn, tp, dn);
 119           if (qh != 0)
 120             cy += mpn_sub_n (np - dn + qn, np - dn + qn, dp - dn, dn - qn);
 121
 122           while (cy != 0)
 123             {
 124               qh -= mpn_sub_1 (qp, qp, qn, 1);
 125               cy -= mpn_add_n (np - dn, np - dn, dp - dn, dn);
 126             }
 127         }
 128
 129       qn = nn - dn - qn + 1;
 130       while (qn > dn)
 131         {
 132           qp -= dn;
 133           np -= dn;
 134           mpn_dc_div_qr_n (qp, np - dn, dp - dn, dn, dip, tp);
 135           qn -= dn;
 136         }
 137
 138       /* Since we pretended we'd need an extra quotient limb before, we now
 139          have made sure the code above left just dn-1=qn quotient limbs to
 140          develop.  Develop that plus a guard limb. */
 141       qn--;
 142       qp -= qn;
 143       np -= dn;
 144       qsave = qp[qn];
 145       mpn_dc_divappr_q_n (qp, np - dn, dp - dn, dn, dip, tp);
 146       MPN_COPY_INCR (qp, qp + 1, qn);
 147       qp[qn] = qsave;
 148     }
 149   else
 150     {
 151       if (qn == 0)
 152         {
 153           qh = mpn_cmp (np - dn, dp - dn, dn) >= 0;
 154           if (qh)
 155             mpn_sub_n (np - dn, np - dn, dp - dn, dn);
 156           TMP_FREE;
 157           return qh;
 158         }
 159
 160       qp -= qn;                 /* point at low limb of next quotient block */
 161       np -= qn;                 /* point in the middle of partial remainder */
 162
 163       if (BELOW_THRESHOLD (qn, DC_DIVAPPR_Q_THRESHOLD))
 164          /* Full precision.  Optimal?  */
 165         qh = mpn_sb_divappr_q (qp, np - dn, nn, dp - dn, dn, dip);
 166       else
 167         {
 168           /* Put quotient in tp, use qp as temporary, since qp lacks a limb.  */
 169           qh = mpn_dc_divappr_q_n (tp, np - qn - 2, dp - (qn + 1), qn + 1, dip, qp);
 170           MPN_COPY (qp, tp + 1, qn);
 171         }
 172     }
 173
 174   TMP_FREE;
 175   return qh;
 176 }
 177
 178 mp_limb_t
 179 mpn_dc_divappr_q (mp_ptr qp, mp_ptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn)
 180 {
 181   mp_limb_t cy;
 182   mp_limb_t xp[2], dip[2];
 183
 184   ASSERT (dn >= 2);
 185
 186   cy = mpn_add_1 (xp, dp + dn - 2, 2, 1);
 187   if (cy != 0)
 188     dip[0] = dip[1] = 0;
 189   else
 190     {
 191       mp_limb_t scratch[10];    /* FIXME */
 192       mpn_invert (dip, xp, 2, scratch);
 193     }
 194
 195   return mpn_preinv_dc_divappr_q (qp, np, nn, dp, dn, dip);
 196 }