Annotation of OpenXM_contrib/gmp/mpn/generic/bz_divrem_n.c, Revision 1.1.1.2
1.1 maekawa 1: /* mpn_bz_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 (C) 2000 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:
1.1.1.2 ! maekawa 38: static mp_limb_t mpn_bz_div_3_halves_by_2
! 39: _PROTO ((mp_ptr qp, mp_ptr np, mp_srcptr dp, mp_size_t n));
1.1 maekawa 40:
41:
42: /* mpn_bz_divrem_n(n) calls 2*mul(n/2)+2*div(n/2), thus to be faster than
43: div(n) = 4*div(n/2), we need mul(n/2) to be faster than the classic way,
44: i.e. n/2 >= KARATSUBA_MUL_THRESHOLD */
45: #ifndef BZ_THRESHOLD
46: #define BZ_THRESHOLD (7 * KARATSUBA_MUL_THRESHOLD)
47: #endif
48:
49: #if 0
50: static
51: unused_mpn_divrem (qp, qxn, np, nn, dp, dn)
52: mp_ptr qp;
53: mp_size_t qxn;
54: mp_ptr np;
55: mp_size_t nn;
56: mp_srcptr dp;
57: mp_size_t dn;
58: {
59: /* This might be useful: */
60: if (qxn != 0)
61: {
62: mp_limb_t c;
63: mp_ptr tp = alloca ((nn + qxn) * BYTES_PER_MP_LIMB);
64: MPN_COPY (tp + qxn - nn, np, nn);
65: MPN_ZERO (tp, qxn);
66: c = mpn_divrem (qp, 0L, tp, nn + qxn, dp, dn);
67: /* Maybe copy proper part of tp to np? Documentation is unclear about
68: the returned np value when qxn != 0 */
69: return c;
70: }
71: }
72: #endif
73:
1.1.1.2 ! maekawa 74:
1.1 maekawa 75: /* mpn_bz_divrem_n - Implements algorithm of page 8 in [1]: divides (np,2n)
76: by (dp,n) and puts the quotient in (qp,n), the remainder in (np,n).
77: Returns most significant limb of the quotient, which is 0 or 1.
78: Requires that the most significant bit of the divisor is set. */
79:
80: mp_limb_t
81: #if __STDC__
82: mpn_bz_divrem_n (mp_ptr qp, mp_ptr np, mp_srcptr dp, mp_size_t n)
83: #else
84: mpn_bz_divrem_n (qp, np, dp, n)
85: mp_ptr qp;
86: mp_ptr np;
87: mp_srcptr dp;
88: mp_size_t n;
89: #endif
90: {
1.1.1.2 ! maekawa 91: mp_limb_t qhl, cc;
1.1 maekawa 92:
93: if (n % 2 != 0)
94: {
1.1.1.2 ! maekawa 95: qhl = mpn_bz_divrem_n (qp + 1, np + 2, dp + 1, n - 1);
! 96: cc = mpn_submul_1 (np + 1, qp + 1, n - 1, dp[0]);
! 97: cc = mpn_sub_1 (np + n, np + n, 1, cc);
! 98: if (qhl) cc += mpn_sub_1 (np + n, np + n, 1, dp[0]);
! 99: while (cc)
! 100: {
! 101: qhl -= mpn_sub_1 (qp + 1, qp + 1, n - 1, (mp_limb_t) 1);
! 102: cc -= mpn_add_n (np + 1, np + 1, dp, n);
! 103: }
1.1 maekawa 104: qhl += mpn_add_1 (qp + 1, qp + 1, n - 1,
1.1.1.2 ! maekawa 105: mpn_sb_divrem_mn (qp, np, n + 1, dp, n));
1.1 maekawa 106: }
107: else
108: {
109: mp_size_t n2 = n/2;
1.1.1.2 ! maekawa 110: qhl = mpn_bz_div_3_halves_by_2 (qp + n2, np + n2, dp, n2);
1.1 maekawa 111: qhl += mpn_add_1 (qp + n2, qp + n2, n2,
1.1.1.2 ! maekawa 112: mpn_bz_div_3_halves_by_2 (qp, np, dp, n2));
1.1 maekawa 113: }
114: return qhl;
115: }
116:
1.1.1.2 ! maekawa 117:
1.1 maekawa 118: /* divides (np, 3n) by (dp, 2n) and puts the quotient in (qp, n),
119: the remainder in (np, 2n) */
120:
121: static mp_limb_t
122: #if __STDC__
1.1.1.2 ! maekawa 123: mpn_bz_div_3_halves_by_2 (mp_ptr qp, mp_ptr np, mp_srcptr dp, mp_size_t n)
1.1 maekawa 124: #else
1.1.1.2 ! maekawa 125: mpn_bz_div_3_halves_by_2 (qp, np, dp, n)
1.1 maekawa 126: mp_ptr qp;
127: mp_ptr np;
128: mp_srcptr dp;
129: mp_size_t n;
130: #endif
131: {
1.1.1.2 ! maekawa 132: mp_size_t twon = n + n;
! 133: mp_limb_t qhl, cc;
! 134: mp_ptr tmp;
! 135: TMP_DECL (marker);
1.1 maekawa 136:
1.1.1.2 ! maekawa 137: TMP_MARK (marker);
1.1 maekawa 138: if (n < BZ_THRESHOLD)
139: qhl = mpn_sb_divrem_mn (qp, np + n, twon, dp + n, n);
140: else
1.1.1.2 ! maekawa 141: qhl = mpn_bz_divrem_n (qp, np + n, dp + n, n);
! 142: tmp = (mp_ptr) TMP_ALLOC (twon * BYTES_PER_MP_LIMB);
1.1 maekawa 143: mpn_mul_n (tmp, qp, dp, n);
1.1.1.2 ! maekawa 144: cc = mpn_sub_n (np, np, tmp, twon);
! 145: TMP_FREE (marker);
! 146: if (qhl) cc += mpn_sub_n (np + n, np + n, dp, n);
! 147: while (cc)
1.1 maekawa 148: {
1.1.1.2 ! maekawa 149: qhl -= mpn_sub_1 (qp, qp, n, (mp_limb_t) 1);
! 150: cc -= mpn_add_n (np, np, dp, twon);
1.1 maekawa 151: }
152: return qhl;
153: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to bz_divrem_n.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpn / generic