rts/gmp/mpn/generic/divrem_2.c

   1 /* mpn_divrem_2 -- Divide natural numbers, producing both remainder and
   2    quotient.  The divisor is two limbs.
   3
   4    THIS FILE CONTAINS INTERNAL FUNCTIONS WITH MUTABLE INTERFACES.  IT IS
   5    ONLY SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS
   6    ALMOST GUARANTEED THAT THEY'LL CHANGE OR DISAPPEAR IN A FUTURE GNU MP
   7    RELEASE.
   8
   9
  10 Copyright (C) 1993, 1994, 1995, 1996, 1999, 2000 Free Software Foundation,
  11 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 2.1 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; see the file COPYING.LIB.  If not, write to
  27 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  28 MA 02111-1307, USA. */
  29
  30 #include "gmp.h"
  31 #include "gmp-impl.h"
  32 #include "longlong.h"
  33
  34 /* Divide num (NP/NSIZE) by den (DP/2) and write
  35    the NSIZE-2 least significant quotient limbs at QP
  36    and the 2 long remainder at NP.  If QEXTRA_LIMBS is
  37    non-zero, generate that many fraction bits and append them after the
  38    other quotient limbs.
  39    Return the most significant limb of the quotient, this is always 0 or 1.
  40
  41    Preconditions:
  42    0. NSIZE >= 2.
  43    1. The most significant bit of the divisor must be set.
  44    2. QP must either not overlap with the input operands at all, or
  45       QP + 2 >= NP must hold true.  (This means that it's
  46       possible to put the quotient in the high part of NUM, right after the
  47       remainder in NUM.
  48    3. NSIZE >= 2, even if QEXTRA_LIMBS is non-zero.  */
  49
  50 mp_limb_t
  51 #if __STDC__
  52 mpn_divrem_2 (mp_ptr qp, mp_size_t qxn,
  53               mp_ptr np, mp_size_t nsize,
  54               mp_srcptr dp)
  55 #else
  56 mpn_divrem_2 (qp, qxn, np, nsize, dp)
  57      mp_ptr qp;
  58      mp_size_t qxn;
  59      mp_ptr np;
  60      mp_size_t nsize;
  61      mp_srcptr dp;
  62 #endif
  63 {
  64   mp_limb_t most_significant_q_limb = 0;
  65   mp_size_t i;
  66   mp_limb_t n1, n0, n2;
  67   mp_limb_t d1, d0;
  68   mp_limb_t d1inv;
  69   int have_preinv;
  70
  71   np += nsize - 2;
  72   d1 = dp[1];
  73   d0 = dp[0];
  74   n1 = np[1];
  75   n0 = np[0];
  76
  77   if (n1 >= d1 && (n1 > d1 || n0 >= d0))
  78     {
  79       sub_ddmmss (n1, n0, n1, n0, d1, d0);
  80       most_significant_q_limb = 1;
  81     }
  82
  83   /* If multiplication is much faster than division, preinvert the most
  84      significant divisor limb before entering the loop.  */
  85   if (UDIV_TIME > 2 * UMUL_TIME + 6)
  86     {
  87       have_preinv = 0;
  88       if ((UDIV_TIME - (2 * UMUL_TIME + 6)) * (nsize - 2) > UDIV_TIME)
  89         {
  90           invert_limb (d1inv, d1);
  91           have_preinv = 1;
  92         }
  93     }
  94
  95   for (i = qxn + nsize - 2 - 1; i >= 0; i--)
  96     {
  97       mp_limb_t q;
  98       mp_limb_t r;
  99
 100       if (i >= qxn)
 101         np--;
 102       else
 103         np[0] = 0;
 104
 105       if (n1 == d1)
 106         {
 107           /* Q should be either 111..111 or 111..110.  Need special treatment
 108              of this rare case as normal division would give overflow.  */
 109           q = ~(mp_limb_t) 0;
 110
 111           r = n0 + d1;
 112           if (r < d1)   /* Carry in the addition? */
 113             {
 114               add_ssaaaa (n1, n0, r - d0, np[0], 0, d0);
 115               qp[i] = q;
 116               continue;
 117             }
 118           n1 = d0 - (d0 != 0);
 119           n0 = -d0;
 120         }
 121       else
 122         {
 123           if (UDIV_TIME > 2 * UMUL_TIME + 6 && have_preinv)
 124             udiv_qrnnd_preinv (q, r, n1, n0, d1, d1inv);
 125           else
 126             udiv_qrnnd (q, r, n1, n0, d1);
 127           umul_ppmm (n1, n0, d0, q);
 128         }
 129
 130       n2 = np[0];
 131
 132     q_test:
 133       if (n1 > r || (n1 == r && n0 > n2))
 134         {
 135           /* The estimated Q was too large.  */
 136           q--;
 137
 138           sub_ddmmss (n1, n0, n1, n0, 0, d0);
 139           r += d1;
 140           if (r >= d1)  /* If not carry, test Q again.  */
 141             goto q_test;
 142         }
 143
 144       qp[i] = q;
 145       sub_ddmmss (n1, n0, r, n2, n1, n0);
 146     }
 147   np[1] = n1;
 148   np[0] = n0;
 149
 150   return most_significant_q_limb;
 151 }