Return to dc_divrem_n.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpn / generic

Annotation of OpenXM_contrib/gmp/mpn/generic/dc_divrem_n.c, Revision 1.1.1.1

1.1       ohara       1: /* mpn_dc_divrem_n and auxilliary routines.
                      2:
                      3:    THE FUNCTIONS IN THIS FILE ARE INTERNAL FUNCTIONS WITH MUTABLE
                      4:    INTERFACES.  IT IS ONLY SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.
                      5:    IN FACT, IT IS ALMOST GUARANTEED THAT THEY'LL CHANGE OR DISAPPEAR IN A
                      6:    FUTURE GNU MP RELEASE.
                      7:
                      8:
                      9: Copyright 2000, 2001, 2002 Free Software Foundation, Inc.
                     10: Contributed by Paul Zimmermann.
                     11:
                     12: This file is part of the GNU MP Library.
                     13:
                     14: The GNU MP Library is free software; you can redistribute it and/or modify
                     15: it under the terms of the GNU Lesser General Public License as published by
                     16: the Free Software Foundation; either version 2.1 of the License, or (at your
                     17: option) any later version.
                     18:
                     19: The GNU MP Library is distributed in the hope that it will be useful, but
                     20: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
                     21: or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
                     22: License for more details.
                     23:
                     24: You should have received a copy of the GNU Lesser General Public License
                     25: along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
                     26: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
                     27: MA 02111-1307, USA. */
                     28:
                     29: #include "gmp.h"
                     30: #include "gmp-impl.h"
                     31:
                     32: /*
                     33: [1] Fast Recursive Division, by Christoph Burnikel and Joachim Ziegler,
                     34:     Technical report MPI-I-98-1-022, october 1998.
                     35:     http://www.mpi-sb.mpg.de/~ziegler/TechRep.ps.gz
                     36: */
                     37:
                     38: static mp_limb_t mpn_dc_div_3_halves_by_2
                     39:   _PROTO ((mp_ptr qp, mp_ptr np, mp_srcptr dp, mp_size_t n));
                     40:
                     41:
                     42: #if 0
                     43: static
                     44: unused_mpn_divrem (qp, qxn, np, nn, dp, dn)
                     45:      mp_ptr qp;
                     46:      mp_size_t qxn;
                     47:      mp_ptr np;
                     48:      mp_size_t nn;
                     49:      mp_srcptr dp;
                     50:      mp_size_t dn;
                     51: {
                     52:   /* This might be useful: */
                     53:   if (qxn != 0)
                     54:     {
                     55:       mp_limb_t c;
                     56:       mp_ptr tp = alloca ((nn + qxn) * BYTES_PER_MP_LIMB);
                     57:       MPN_COPY (tp + qxn - nn, np, nn);
                     58:       MPN_ZERO (tp, qxn);
                     59:       c = mpn_divrem (qp, 0L, tp, nn + qxn, dp, dn);
                     60:       /* Maybe copy proper part of tp to np?  Documentation is unclear about
                     61:         the returned np value when qxn != 0 */
                     62:       return c;
                     63:     }
                     64: }
                     65: #endif
                     66:
                     67:
                     68: /* mpn_dc_divrem_n - Implements algorithm of page 8 in [1]: divides (np,2n)
                     69:    by (dp,n) and puts the quotient in (qp,n), the remainder in (np,n).
                     70:    Returns most significant limb of the quotient, which is 0 or 1.
                     71:    Requires that the most significant bit of the divisor is set.  */
                     72:
                     73: mp_limb_t
                     74: mpn_dc_divrem_n (mp_ptr qp, mp_ptr np, mp_srcptr dp, mp_size_t n)
                     75: {
                     76:   mp_limb_t qhl, cc;
                     77:
                     78:   if (n % 2 != 0)
                     79:     {
                     80:       qhl = mpn_dc_divrem_n (qp + 1, np + 2, dp + 1, n - 1);
                     81:       cc = mpn_submul_1 (np + 1, qp + 1, n - 1, dp[0]);
                     82:       cc = mpn_sub_1 (np + n, np + n, 1, cc);
                     83:       if (qhl) cc += mpn_sub_1 (np + n, np + n, 1, dp[0]);
                     84:       while (cc)
                     85:        {
                     86:          qhl -= mpn_sub_1 (qp + 1, qp + 1, n - 1, (mp_limb_t) 1);
                     87:          cc -= mpn_add_n (np + 1, np + 1, dp, n);
                     88:        }
                     89:       qhl += mpn_add_1 (qp + 1, qp + 1, n - 1,
                     90:                        mpn_sb_divrem_mn (qp, np, n + 1, dp, n));
                     91:     }
                     92:   else
                     93:     {
                     94:       mp_size_t n2 = n/2;
                     95:       qhl = mpn_dc_div_3_halves_by_2 (qp + n2, np + n2, dp, n2);
                     96:       qhl += mpn_add_1 (qp + n2, qp + n2, n2,
                     97:                        mpn_dc_div_3_halves_by_2 (qp, np, dp, n2));
                     98:     }
                     99:   return qhl;
                    100: }
                    101:
                    102:
                    103: /* divides (np, 3n) by (dp, 2n) and puts the quotient in (qp, n),
                    104:    the remainder in (np, 2n) */
                    105:
                    106: static mp_limb_t
                    107: mpn_dc_div_3_halves_by_2 (mp_ptr qp, mp_ptr np, mp_srcptr dp, mp_size_t n)
                    108: {
                    109:   mp_size_t twon = n + n;
                    110:   mp_limb_t qhl, cc;
                    111:   mp_ptr tmp;
                    112:   TMP_DECL (marker);
                    113:
                    114:   TMP_MARK (marker);
                    115:   if (n < DIV_DC_THRESHOLD)
                    116:     qhl = mpn_sb_divrem_mn (qp, np + n, twon, dp + n, n);
                    117:   else
                    118:     qhl = mpn_dc_divrem_n (qp, np + n, dp + n, n);
                    119:   tmp = (mp_ptr) TMP_ALLOC (twon * BYTES_PER_MP_LIMB);
                    120:   mpn_mul_n (tmp, qp, dp, n);
                    121:   cc = mpn_sub_n (np, np, tmp, twon);
                    122:   TMP_FREE (marker);
                    123:   if (qhl) cc += mpn_sub_n (np + n, np + n, dp, n);
                    124:   while (cc)
                    125:     {
                    126:       qhl -= mpn_sub_1 (qp, qp, n, (mp_limb_t) 1);
                    127:       cc -= mpn_add_n (np, np, dp, twon);
                    128:     }
                    129:   return qhl;
                    130: }