ghc/runtime/gmp/mpz_div.c

   1 /* mpz_div -- divide two integers and produce a quotient.
   2
   3 Copyright (C) 1991 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 General Public License as published by
   9 the Free Software Foundation; either version 2, or (at your option)
  10 any later version.
  11
  12 The GNU MP Library is distributed in the hope that it will be useful,
  13 but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 GNU General Public License for more details.
  16
  17 You should have received a copy of the GNU General Public License
  18 along with the GNU MP Library; see the file COPYING.  If not, write to
  19 the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.  */
  20
  21 #include "gmp.h"
  22 #include "gmp-impl.h"
  23 #include "longlong.h"
  24
  25 void
  26 #ifdef __STDC__
  27 mpz_div (MP_INT *quot, const MP_INT *num, const MP_INT *den)
  28 #else
  29 mpz_div (quot, num, den)
  30      MP_INT *quot;
  31      const MP_INT *num;
  32      const MP_INT *den;
  33 #endif
  34 {
  35   mp_srcptr np, dp;
  36   mp_ptr qp, rp;
  37   mp_size nsize = num->size;
  38   mp_size dsize = den->size;
  39   mp_size qsize, rsize;
  40   mp_size sign_quotient = nsize ^ dsize;
  41   unsigned normalization_steps;
  42
  43   nsize = ABS (nsize);
  44   dsize = ABS (dsize);
  45
  46   /* Ensure space is enough for quotient. */
  47
  48   qsize = nsize - dsize + 1;    /* qsize cannot be bigger than this.  */
  49   if (qsize <= 0)
  50     {
  51       quot->size = 0;
  52       return;
  53     }
  54
  55   if (quot->alloc < qsize)
  56     _mpz_realloc (quot, qsize);
  57
  58   qp = quot->d;
  59   np = num->d;
  60   dp = den->d;
  61   rp = (mp_ptr) alloca ((nsize + 1) * BYTES_PER_MP_LIMB);
  62
  63   count_leading_zeros (normalization_steps, dp[dsize - 1]);
  64
  65   /* Normalize the denominator and the numerator.  */
  66   if (normalization_steps != 0)
  67     {
  68       mp_ptr tp;
  69       mp_limb ndigit;
  70
  71       /* Shift up the denominator setting the most significant bit of
  72          the most significant word.  Use temporary storage not to clobber
  73          the original contents of the denominator.  */
  74       tp = (mp_ptr) alloca (dsize * BYTES_PER_MP_LIMB);
  75       (void) mpn_lshift (tp, dp, dsize, normalization_steps);
  76       dp = tp;
  77
  78       /* Shift up the numerator, possibly introducing a new most
  79          significant word.  Move the shifted numerator in the remainder
  80          meanwhile.  */
  81       ndigit = mpn_lshift (rp, np, nsize, normalization_steps);
  82       if (ndigit != 0)
  83         {
  84           rp[nsize] = ndigit;
  85           rsize = nsize + 1;
  86         }
  87       else
  88         rsize = nsize;
  89     }
  90   else
  91     {
  92       /* The denominator is already normalized, as required.
  93          Copy it to temporary space if it overlaps with the quotient.  */
  94       if (dp == qp)
  95         {
  96           dp = (mp_ptr) alloca (dsize * BYTES_PER_MP_LIMB);
  97           MPN_COPY ((mp_ptr) dp, qp, dsize);
  98         }
  99
 100       /* Move the numerator to the remainder.  */
 101       MPN_COPY (rp, np, nsize);
 102       rsize = nsize;
 103     }
 104
 105   qsize = rsize - dsize + mpn_div (qp, rp, rsize, dp, dsize);
 106
 107   /* Normalize the quotient.  We may have at most one leading
 108      zero-word, so no loop is needed.  */
 109   if (qsize > 0)
 110     qsize -= (qp[qsize - 1] == 0);
 111
 112   if (sign_quotient < 0)
 113     qsize = -qsize;
 114   quot->size = qsize;
 115
 116   alloca (0);
 117 }