Return to divexact.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpz

Annotation of OpenXM_contrib/gmp/mpz/divexact.c, Revision 1.1.1.1

1.1       maekawa     1: /* mpz_divexact -- finds quotient when known that quot * den == num && den != 0.
                      2:
                      3: Copyright (C) 1991, 1993, 1994, 1995 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: /*  Ken Weber (kweber@mat.ufrgs.br, kweber@mcs.kent.edu)
                     23:
                     24:     Funding for this work has been partially provided by Conselho Nacional
                     25:     de Desenvolvimento Cienti'fico e Tecnolo'gico (CNPq) do Brazil, Grant
                     26:     301314194-2, and was done while I was a visiting reseacher in the Instituto
                     27:     de Matema'tica at Universidade Federal do Rio Grande do Sul (UFRGS).
                     28:
                     29:     References:
                     30:         T. Jebelean, An algorithm for exact division, Journal of Symbolic
                     31:         Computation, v. 15, 1993, pp. 169-180.  */
                     32:
                     33: #include "gmp.h"
                     34: #include "gmp-impl.h"
                     35: #include "longlong.h"
                     36:
                     37: void
                     38: #if __STDC__
                     39: mpz_divexact (mpz_ptr quot, mpz_srcptr num, mpz_srcptr den)
                     40: #else
                     41: mpz_divexact (quot, num, den)
                     42:      mpz_ptr quot;
                     43:      mpz_srcptr num;
                     44:      mpz_srcptr den;
                     45: #endif
                     46: {
                     47:   mp_ptr qp, tp;
                     48:   mp_size_t qsize, tsize;
                     49:
                     50:   mp_srcptr np = num->_mp_d;
                     51:   mp_srcptr dp = den->_mp_d;
                     52:   mp_size_t nsize = ABS (num->_mp_size);
                     53:   mp_size_t dsize = ABS (den->_mp_size);
                     54:   TMP_DECL (marker);
                     55:
                     56:   /*  Generate divide-by-zero error if dsize == 0.  */
                     57:   if (dsize == 0)
                     58:     {
                     59:       quot->_mp_size = 1 / dsize;
                     60:       return;
                     61:     }
                     62:
                     63:   if (nsize == 0)
                     64:     {
                     65:       quot->_mp_size = 0;
                     66:       return;
                     67:     }
                     68:
                     69:   qsize = nsize - dsize + 1;
                     70:   if (quot->_mp_alloc < qsize)
                     71:     _mpz_realloc (quot, qsize);
                     72:   qp = quot->_mp_d;
                     73:
                     74:   TMP_MARK (marker);
                     75:
                     76:   /*  QUOT <-- NUM/2^r, T <-- DEN/2^r where = r number of twos in DEN.  */
                     77:   while (dp[0] == 0)
                     78:     np += 1, nsize -= 1, dp += 1, dsize -= 1;
                     79:   tsize = MIN (qsize, dsize);
                     80:   if (dp[0] & 1)
                     81:     {
                     82:       if (qp != dp)
                     83:        MPN_COPY (qp, np, qsize);
                     84:       if (qp == dp)            /*  QUOT and DEN overlap.  */
                     85:        {
                     86:          tp = (mp_ptr) TMP_ALLOC (sizeof (mp_limb_t) * tsize);
                     87:          MPN_COPY (tp, dp, tsize);
                     88:        }
                     89:       else
                     90:        tp = (mp_ptr) dp;
                     91:     }
                     92:   else
                     93:     {
                     94:       unsigned long int r;
                     95:       tp = (mp_ptr) TMP_ALLOC (sizeof (mp_limb_t) * tsize);
                     96:       count_trailing_zeros (r, dp[0]);
                     97:       mpn_rshift (tp, dp, tsize, r);
                     98:       if (dsize > tsize)
                     99:        tp[tsize-1] |= dp[tsize] << (BITS_PER_MP_LIMB - r);
                    100:       mpn_rshift (qp, np, qsize, r);
                    101:       if (nsize > qsize)
                    102:        qp[qsize-1] |= np[qsize] << (BITS_PER_MP_LIMB - r);
                    103:     }
                    104:
                    105:   /*  Now QUOT <-- QUOT/T.  */
                    106:   mpn_bdivmod (qp, qp, qsize, tp, tsize, qsize * BITS_PER_MP_LIMB);
                    107:   MPN_NORMALIZE (qp, qsize);
                    108:
                    109:   quot->_mp_size = (num->_mp_size < 0) == (den->_mp_size < 0) ? qsize : -qsize;
                    110:
                    111:   TMP_FREE (marker);
                    112: }