contrib/libgmp/mpz/tdiv_q.c

   1 /* mpz_tdiv_q -- divide two integers and produce a quotient.
   2
   3 Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc.
   4
   5 This file is part of the GNU MP Library.
   6
   7 The GNU MP Library is free software; you can redistribute it and/or modify
   8 it under the terms of the GNU Library General Public License as published by
   9 the Free Software Foundation; either version 2 of the License, or (at your
  10 option) any later version.
  11
  12 The GNU MP Library is distributed in the hope that it will be useful, but
  13 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  14 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
  15 License for more details.
  16
  17 You should have received a copy of the GNU Library General Public License
  18 along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
  19 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  20 MA 02111-1307, USA. */
  21
  22 #include "gmp.h"
  23 #include "gmp-impl.h"
  24 #include "longlong.h"
  25
  26 void
  27 #if __STDC__
  28 mpz_tdiv_q (mpz_ptr quot, mpz_srcptr num, mpz_srcptr den)
  29 #else
  30 mpz_tdiv_q (quot, num, den)
  31      mpz_ptr quot;
  32      mpz_srcptr num;
  33      mpz_srcptr den;
  34 #endif
  35 {
  36   mp_srcptr np, dp;
  37   mp_ptr qp, rp;
  38   mp_size_t nsize = num->_mp_size;
  39   mp_size_t dsize = den->_mp_size;
  40   mp_size_t qsize, rsize;
  41   mp_size_t sign_quotient = nsize ^ dsize;
  42   unsigned normalization_steps;
  43   mp_limb_t q_limb;
  44   TMP_DECL (marker);
  45
  46   nsize = ABS (nsize);
  47   dsize = ABS (dsize);
  48
  49   /* Ensure space is enough for quotient. */
  50
  51   qsize = nsize - dsize + 1;    /* qsize cannot be bigger than this.  */
  52   if (qsize <= 0)
  53     {
  54       quot->_mp_size = 0;
  55       return;
  56     }
  57
  58   if (quot->_mp_alloc < qsize)
  59     _mpz_realloc (quot, qsize);
  60
  61   qp = quot->_mp_d;
  62   np = num->_mp_d;
  63   dp = den->_mp_d;
  64
  65   /* Optimize division by a single-limb divisor.  */
  66   if (dsize == 1)
  67     {
  68       mpn_divmod_1 (qp, np, nsize, dp[0]);
  69       qsize -= qp[qsize - 1] == 0;
  70       quot->_mp_size = sign_quotient >= 0 ? qsize : -qsize;
  71       return;
  72     }
  73
  74   TMP_MARK (marker);
  75
  76   rp = (mp_ptr) TMP_ALLOC ((nsize + 1) * BYTES_PER_MP_LIMB);
  77
  78   count_leading_zeros (normalization_steps, dp[dsize - 1]);
  79
  80   /* Normalize the denominator, i.e. make its most significant bit set by
  81      shifting it NORMALIZATION_STEPS bits to the left.  Also shift the
  82      numerator the same number of steps (to keep the quotient the same!).  */
  83   if (normalization_steps != 0)
  84     {
  85       mp_ptr tp;
  86       mp_limb_t nlimb;
  87
  88       /* Shift up the denominator setting the most significant bit of
  89          the most significant word.  Use temporary storage not to clobber
  90          the original contents of the denominator.  */
  91       tp = (mp_ptr) TMP_ALLOC (dsize * BYTES_PER_MP_LIMB);
  92       mpn_lshift (tp, dp, dsize, normalization_steps);
  93       dp = tp;
  94
  95       /* Shift up the numerator, possibly introducing a new most
  96          significant word.  Move the shifted numerator in the remainder
  97          meanwhile.  */
  98       nlimb = mpn_lshift (rp, np, nsize, normalization_steps);
  99       if (nlimb != 0)
 100         {
 101           rp[nsize] = nlimb;
 102           rsize = nsize + 1;
 103         }
 104       else
 105         rsize = nsize;
 106     }
 107   else
 108     {
 109       /* The denominator is already normalized, as required.  Copy it to
 110          temporary space if it overlaps with the quotient.  */
 111       if (dp == qp)
 112         {
 113           dp = (mp_ptr) TMP_ALLOC (dsize * BYTES_PER_MP_LIMB);
 114           MPN_COPY ((mp_ptr) dp, qp, dsize);
 115         }
 116
 117       /* Move the numerator to the remainder.  */
 118       MPN_COPY (rp, np, nsize);
 119       rsize = nsize;
 120     }
 121
 122   q_limb = mpn_divmod (qp, rp, rsize, dp, dsize);
 123
 124   qsize = rsize - dsize;
 125   if (q_limb)
 126     {
 127       qp[qsize] = q_limb;
 128       qsize += 1;
 129     }
 130
 131   quot->_mp_size = sign_quotient >= 0 ? qsize : -qsize;
 132   TMP_FREE (marker);
 133 }