contrib/gmp/mpn/generic/mu_div_q.c

   1 /* mpn_mu_div_q, mpn_preinv_mu_div_q.
   2
   3    Contributed to the GNU project by Torbjörn Granlund.
   4
   5    THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH A MUTABLE INTERFACE.  IT IS
   6    ONLY SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS
   7    ALMOST GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GMP
   8    RELEASE.
   9
  10 Copyright 2005, 2006, 2007 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 the GNU Lesser General Public License as published by
  16 the Free Software Foundation; either version 3 of the License, or (at your
  17 option) any later version.
  18
  19 The GNU MP Library is distributed in the hope that it will be useful, but
  20 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  21 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  22 License for more details.
  23
  24 You should have received a copy of the GNU Lesser General Public License
  25 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  26
  27
  28 /*
  29   Things to work on:
  30
  31   1. This is a rudimentary implementation of mpn_mu_div_q.  The algorithm is
  32      probably close to optimal, except when mpn_mu_divappr_q fails.
  33
  34      An alternative which could be considered for much simpler code for the
  35      complex qn>=dn arm would be to allocate a temporary nn+1 limb buffer, then
  36      simply call mpn_mu_divappr_q.  Such a temporary allocation is
  37      unfortunately very large.
  38
  39   2. Instead of falling back to mpn_mu_div_qr when we detect a possible
  40      mpn_mu_divappr_q rounding problem, we could multiply and compare.
  41      Unfortunately, since mpn_mu_divappr_q does not return the partial
  42      remainder, this also doesn't become optimal.  A mpn_mu_divappr_qr
  43      could solve that.
  44
  45   3. The allocations done here should be made from the scratch area.
  46 */
  47
  48 #include <stdlib.h>             /* for NULL */
  49 #include "gmp.h"
  50 #include "gmp-impl.h"
  51
  52
  53 mp_limb_t
  54 mpn_mu_div_q (mp_ptr qp,
  55               mp_ptr np, mp_size_t nn,
  56               mp_srcptr dp, mp_size_t dn,
  57               mp_ptr scratch)
  58 {
  59   mp_ptr tp, rp, ip, this_ip;
  60   mp_size_t qn, in, this_in;
  61   mp_limb_t cy;
  62   TMP_DECL;
  63
  64   TMP_MARK;
  65
  66   qn = nn - dn;
  67
  68   tp = TMP_BALLOC_LIMBS (qn + 1);
  69
  70   if (qn >= dn)                 /* nn >= 2*dn + 1 */
  71     {
  72       /* Find max inverse size needed by the two preinv calls.  */
  73       if (dn != qn)
  74         {
  75           mp_size_t in1, in2;
  76
  77           in1 = mpn_mu_div_qr_choose_in (qn - dn, dn, 0);
  78           in2 = mpn_mu_divappr_q_choose_in (dn + 1, dn, 0);
  79           in = MAX (in1, in2);
  80         }
  81       else
  82         {
  83           in = mpn_mu_divappr_q_choose_in (dn + 1, dn, 0);
  84         }
  85
  86       ip = TMP_BALLOC_LIMBS (in + 1);
  87
  88       if (dn == in)
  89         {
  90           MPN_COPY (scratch + 1, dp, in);
  91           scratch[0] = 1;
  92           mpn_invert (ip, scratch, in + 1, NULL);
  93           MPN_COPY_INCR (ip, ip + 1, in);
  94         }
  95       else
  96         {
  97           cy = mpn_add_1 (scratch, dp + dn - (in + 1), in + 1, 1);
  98           if (UNLIKELY (cy != 0))
  99             MPN_ZERO (ip, in);
 100           else
 101             {
 102               mpn_invert (ip, scratch, in + 1, NULL);
 103               MPN_COPY_INCR (ip, ip + 1, in);
 104             }
 105         }
 106
 107        /* |_______________________|   dividend
 108                          |________|   divisor  */
 109       rp = TMP_BALLOC_LIMBS (2 * dn + 1);
 110       if (dn != qn)             /* FIXME: perhaps mpn_mu_div_qr should DTRT */
 111         {
 112           this_in = mpn_mu_div_qr_choose_in (qn - dn, dn, 0);
 113           this_ip = ip + in - this_in;
 114           mpn_preinv_mu_div_qr (tp + dn + 1, rp + dn + 1, np + dn, qn, dp, dn,
 115                                 this_ip, this_in, scratch);
 116         }
 117       else
 118         MPN_COPY (rp + dn + 1, np + dn, dn);
 119
 120       MPN_COPY (rp + 1, np, dn);
 121       rp[0] = 0;
 122       this_in = mpn_mu_divappr_q_choose_in (dn + 1, dn, 0);
 123       this_ip = ip + in - this_in;
 124       mpn_preinv_mu_divappr_q (tp, rp, 2*dn + 1, dp, dn, this_ip, this_in, scratch);
 125
 126       /* The max error of mpn_mu_divappr_q is +4.  If the low quotient limb is
 127          greater than the max error, we cannot trust the quotient.  */
 128       if (tp[0] > 4)
 129         {
 130           MPN_COPY (qp, tp + 1, qn);
 131         }
 132       else
 133         {
 134           /* Fall back to plain mpn_mu_div_qr.  */
 135           mpn_mu_div_qr (qp, rp, np, nn, dp, dn, scratch);
 136         }
 137     }
 138   else
 139     {
 140        /* |_______________________|   dividend
 141                  |________________|   divisor  */
 142       mpn_mu_divappr_q (tp, np + nn - (2*qn + 2), 2*qn + 2, dp + dn - (qn + 1), qn + 1, scratch);
 143
 144       if (tp[0] > 4)
 145         {
 146           MPN_COPY (qp, tp + 1, qn);
 147         }
 148       else
 149         {
 150           rp = TMP_BALLOC_LIMBS (dn);
 151           mpn_mu_div_qr (qp, rp, np, nn, dp, dn, scratch);
 152         }
 153     }
 154
 155   TMP_FREE;
 156   return 0;
 157 }