Annotation of OpenXM_contrib/gmp/mpn/generic/bz_divrem_n.c, Revision 1.1.1.1
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:
38: static mp_limb_t mpn_bz_div_3_halves_by_2 _PROTO ((mp_ptr, mp_ptr, mp_srcptr,
39: mp_size_t, mp_ptr));
40:
41: static mp_limb_t mpn_bz_divrem_aux _PROTO ((mp_ptr, mp_ptr, mp_srcptr,
42: mp_size_t, mp_ptr));
43:
44: /* mpn_bz_divrem_n(n) calls 2*mul(n/2)+2*div(n/2), thus to be faster than
45: div(n) = 4*div(n/2), we need mul(n/2) to be faster than the classic way,
46: i.e. n/2 >= KARATSUBA_MUL_THRESHOLD */
47: #ifndef BZ_THRESHOLD
48: #define BZ_THRESHOLD (7 * KARATSUBA_MUL_THRESHOLD)
49: #endif
50:
51: #if 0
52: static
53: unused_mpn_divrem (qp, qxn, np, nn, dp, dn)
54: mp_ptr qp;
55: mp_size_t qxn;
56: mp_ptr np;
57: mp_size_t nn;
58: mp_srcptr dp;
59: mp_size_t dn;
60: {
61: /* This might be useful: */
62: if (qxn != 0)
63: {
64: mp_limb_t c;
65: mp_ptr tp = alloca ((nn + qxn) * BYTES_PER_MP_LIMB);
66: MPN_COPY (tp + qxn - nn, np, nn);
67: MPN_ZERO (tp, qxn);
68: c = mpn_divrem (qp, 0L, tp, nn + qxn, dp, dn);
69: /* Maybe copy proper part of tp to np? Documentation is unclear about
70: the returned np value when qxn != 0 */
71: return c;
72: }
73: }
74: #endif
75:
76: /* mpn_bz_divrem_n - Implements algorithm of page 8 in [1]: divides (np,2n)
77: by (dp,n) and puts the quotient in (qp,n), the remainder in (np,n).
78: Returns most significant limb of the quotient, which is 0 or 1.
79: Requires that the most significant bit of the divisor is set. */
80:
81: mp_limb_t
82: #if __STDC__
83: mpn_bz_divrem_n (mp_ptr qp, mp_ptr np, mp_srcptr dp, mp_size_t n)
84: #else
85: mpn_bz_divrem_n (qp, np, dp, n)
86: mp_ptr qp;
87: mp_ptr np;
88: mp_srcptr dp;
89: mp_size_t n;
90: #endif
91: {
92: mp_limb_t qhl = 0;
93: if (mpn_cmp (np + n, dp, n) >= 0)
94: {
95: qhl = 1;
96: mpn_sub_n (np + n, np + n, dp, n);
97: abort ();
98: }
99: if (n % 2 != 0)
100: {
101: /* divide (2n - 2) most significant limbs from np by those (n - 1) from dp */
102: if (n < BZ_THRESHOLD)
103: qhl += mpn_sb_divrem_mn (qp + 1, np + 2, 2 * (n - 1), dp + 1, n - 1);
104: else
105: qhl += mpn_bz_divrem_n (qp + 1, np + 2, dp + 1, n - 1);
106: /* now (qp + 1, n - 1) contains the quotient of (np + 2, 2n - 2) by
107: (dp + 1, n - 1) and (np + 2, n - 1) contains the remainder */
108: if (mpn_sub_1 (np + n, np + n, 1,
109: mpn_submul_1 (np + 1, qp + 1, n - 1, dp[0])))
110: {
111: /* quotient too large */
112: qhl -= mpn_sub_1 (qp + 1, qp + 1, n - 1, 1L);
113: if (mpn_add_n (np + 1, np + 1, dp, n) == 0)
114: { /* still too large */
115: qhl -= mpn_sub_1 (qp + 1, qp + 1, n - 1, 1L);
116: mpn_add_n (np + 1, np + 1, dp, n); /* always carry here */
117: }
118: }
119: /* now divide (np, n + 1) by (dp, n) */
120: qhl += mpn_add_1 (qp + 1, qp + 1, n - 1,
121: mpn_sb_divrem_mn (qp, np, n + 1, dp, n));
122: }
123: else
124: {
125: mp_ptr tmp;
126: mp_size_t n2 = n/2;
127: TMP_DECL (marker);
128: TMP_MARK (marker);
129: tmp = (mp_ptr) TMP_ALLOC (n * BYTES_PER_MP_LIMB);
130: qhl = mpn_bz_div_3_halves_by_2 (qp + n2, np + n2, dp, n2, tmp);
131: qhl += mpn_add_1 (qp + n2, qp + n2, n2,
132: mpn_bz_div_3_halves_by_2 (qp, np, dp, n2, tmp));
133: TMP_FREE (marker);
134: }
135: return qhl;
136: }
137:
138: /* Like mpn_bz_divrem_n, but without memory allocation. Also
139: assumes mpn_cmp (np + n, dp, n) < 0 */
140:
141: static mp_limb_t
142: #if __STDC__
143: mpn_bz_divrem_aux (mp_ptr qp, mp_ptr np, mp_srcptr dp, mp_size_t n, mp_ptr tmp)
144: #else
145: mpn_bz_divrem_aux (qp, np, dp, n, tmp)
146: mp_ptr qp;
147: mp_ptr np;
148: mp_srcptr dp;
149: mp_size_t n;
150: mp_ptr tmp;
151: #endif
152: {
153: mp_limb_t qhl;
154:
155: if (n % 2 != 0)
156: {
157: /* divide (2n - 2) most significant limbs from np by those (n - 1) from dp */
158: qhl = mpn_bz_divrem_aux (qp + 1, np + 2, dp + 1, n - 1, tmp);
159: /* now (qp + 1, n - 1) contains the quotient of (np + 2, 2n - 2) by
160: (dp + 1, n - 1) and (np + 2, n - 1) contains the remainder */
161: if (mpn_sub_1 (np + n, np + n, 1,
162: mpn_submul_1 (np + 1, qp + 1, n - 1, dp[0])))
163: {
164: /* quotient too large */
165: qhl -= mpn_sub_1 (qp + 1, qp + 1, n - 1, 1L);
166: if (mpn_add_n (np + 1, np + 1, dp, n) == 0)
167: { /* still too large */
168: qhl -= mpn_sub_1 (qp + 1, qp + 1, n - 1, 1L);
169: mpn_add_n (np + 1, np + 1, dp, n); /* always carry here */
170: }
171: }
172: /* now divide (np, n + 1) by (dp, n) */
173: qhl += mpn_add_1 (qp + 1, qp + 1, n - 1,
174: mpn_sb_divrem_mn (qp, np, n + 1, dp, n));
175: }
176: else
177: {
178: mp_size_t n2 = n/2;
179: qhl = mpn_bz_div_3_halves_by_2 (qp + n2, np + n2, dp, n2, tmp);
180: qhl += mpn_add_1 (qp + n2, qp + n2, n2,
181: mpn_bz_div_3_halves_by_2 (qp, np, dp, n2, tmp));
182: }
183: return qhl;
184: }
185:
186: /* divides (np, 3n) by (dp, 2n) and puts the quotient in (qp, n),
187: the remainder in (np, 2n) */
188:
189: static mp_limb_t
190: #if __STDC__
191: mpn_bz_div_3_halves_by_2 (mp_ptr qp, mp_ptr np, mp_srcptr dp, mp_size_t n,
192: mp_ptr tmp)
193: #else
194: mpn_bz_div_3_halves_by_2 (qp, np, dp, n, tmp)
195: mp_ptr qp;
196: mp_ptr np;
197: mp_srcptr dp;
198: mp_size_t n;
199: mp_ptr tmp;
200: #endif
201: {
202: mp_size_t twon = n + n;
203: mp_limb_t qhl;
204:
205: if (n < BZ_THRESHOLD)
206: qhl = mpn_sb_divrem_mn (qp, np + n, twon, dp + n, n);
207: else
208: qhl = mpn_bz_divrem_aux (qp, np + n, dp + n, n, tmp);
209: /* q = (qp, n), c = (np + n, n) with the notations of [1] */
210: mpn_mul_n (tmp, qp, dp, n);
211: if (qhl != 0)
212: mpn_add_n (tmp + n, tmp + n, dp, n);
213: if (mpn_sub_n (np, np, tmp, twon)) /* R = (np, 2n) */
214: {
215: qhl -= mpn_sub_1 (qp, qp, n, 1L);
216: if (mpn_add_n (np, np, dp, twon) == 0)
217: { /* qp still too large */
218: qhl -= mpn_sub_1 (qp, qp, n, 1L);
219: mpn_add_n (np, np, dp, twon); /* always carry here */
220: }
221: }
222: return qhl;
223: }
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