1 /* mpz_dmincl.c -- include file for mpz_dm.c, mpz_mod.c, mdiv.c.
3 Copyright (C) 1991 Free Software Foundation, Inc.
5 This file is part of the GNU MP Library.
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)
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.
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. */
21 /* THIS CODE IS OBSOLETE. IT WILL SOON BE REPLACED BY CLEANER CODE WITH
22 LESS MEMORY ALLOCATION OVERHEAD. */
24 /* If den == quot, den needs temporary storage.
25 If den == rem, den needs temporary storage.
26 If num == quot, num needs temporary storage.
27 If den has temporary storage, it can be normalized while being copied,
28 i.e no extra storage should be allocated. */
30 /* This is the function body of mdiv, mpz_divmod, and mpz_mod.
32 If COMPUTE_QUOTIENT is defined, the quotient is put in the MP_INT
33 object quot, otherwise that variable is not referenced at all.
35 The remainder is always computed, and the result is put in the MP_INT
41 mp_size nsize = num->size;
42 mp_size dsize = den->size;
44 mp_size sign_remainder = nsize;
45 #ifdef COMPUTE_QUOTIENT
46 mp_size sign_quotient = nsize ^ dsize;
48 unsigned normalization_steps;
53 /* Ensure space is enough for quotient and remainder. */
55 /* We need space for an extra limb in the remainder, because it's
56 up-shifted (normalized) below. */
58 if (rem->alloc < rsize)
59 _mpz_realloc (rem, rsize);
61 qsize = nsize - dsize + 1; /* qsize cannot be bigger than this. */
64 #ifdef COMPUTE_QUOTIENT
69 rem->size = num->size;
70 MPN_COPY (rem->d, num->d, nsize);
75 #ifdef COMPUTE_QUOTIENT
76 if (quot->alloc < qsize)
77 _mpz_realloc (quot, qsize);
80 qp = (mp_ptr) alloca (qsize * BYTES_PER_MP_LIMB);
86 /* Make sure quot and num are different. Otherwise the numerator
87 would be successively overwritten by the quotient digits. */
90 np = (mp_ptr) alloca (nsize * BYTES_PER_MP_LIMB);
91 MPN_COPY (np, qp, nsize);
94 count_leading_zeros (normalization_steps, dp[dsize - 1]);
96 /* Normalize the denominator, i.e. make its most significant bit set by
97 shifting it NORMALIZATION_STEPS bits to the left. Also shift the
98 numerator the same number of steps (to keep the quotient the same!). */
99 if (normalization_steps != 0)
104 /* Shift up the denominator setting the most significant bit of
105 the most significant word. Use temporary storage not to clobber
106 the original contents of the denominator. */
107 tp = (mp_ptr) alloca (dsize * BYTES_PER_MP_LIMB);
108 (void) mpn_lshift (tp, dp, dsize, normalization_steps);
111 /* Shift up the numerator, possibly introducing a new most
112 significant word. Move the shifted numerator in the remainder
114 ndigit = mpn_lshift (rp, np, nsize, normalization_steps);
125 #ifdef COMPUTE_QUOTIENT
126 if (rem == den || quot == den)
133 tp = (mp_ptr) alloca (dsize * BYTES_PER_MP_LIMB);
134 MPN_COPY (tp, dp, dsize);
138 /* Move the numerator to the remainder. */
140 MPN_COPY (rp, np, nsize);
145 qsize = rsize - dsize + mpn_div (qp, rp, rsize, dp, dsize);
149 /* Normalize the remainder. */
152 if (rp[rsize - 1] != 0)
157 if (normalization_steps != 0)
158 rsize = mpn_rshift (rp, rp, rsize, normalization_steps);
160 rem->size = (sign_remainder >= 0) ? rsize : -rsize;
162 #ifdef COMPUTE_QUOTIENT
163 /* Normalize the quotient. We may have at most one leading
164 zero-word, so no loop is needed. */
166 qsize -= (qp[qsize - 1] == 0);
168 quot->size = (sign_quotient >= 0) ? qsize : -qsize;