Annotation of OpenXM_contrib/gmp/mpn/generic/mul_n.c, Revision 1.1.1.1
1.1 maekawa 1: /* mpn_mul_n -- Multiply two natural numbers of length n.
2:
3: Copyright (C) 1991, 1992, 1993, 1994, 1996 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: #include "gmp.h"
23: #include "gmp-impl.h"
24:
25: /* Multiply the natural numbers u (pointed to by UP) and v (pointed to by VP),
26: both with SIZE limbs, and store the result at PRODP. 2 * SIZE limbs are
27: always stored. Return the most significant limb.
28:
29: Argument constraints:
30: 1. PRODP != UP and PRODP != VP, i.e. the destination
31: must be distinct from the multiplier and the multiplicand. */
32:
33: /* If KARATSUBA_THRESHOLD is not already defined, define it to a
34: value which is good on most machines. */
35: #ifndef KARATSUBA_THRESHOLD
36: #define KARATSUBA_THRESHOLD 32
37: #endif
38:
39: /* The code can't handle KARATSUBA_THRESHOLD smaller than 2. */
40: #if KARATSUBA_THRESHOLD < 2
41: #undef KARATSUBA_THRESHOLD
42: #define KARATSUBA_THRESHOLD 2
43: #endif
44:
45: /* Handle simple cases with traditional multiplication.
46:
47: This is the most critical code of multiplication. All multiplies rely
48: on this, both small and huge. Small ones arrive here immediately. Huge
49: ones arrive here as this is the base case for Karatsuba's recursive
50: algorithm below. */
51:
52: void
53: #if __STDC__
54: impn_mul_n_basecase (mp_ptr prodp, mp_srcptr up, mp_srcptr vp, mp_size_t size)
55: #else
56: impn_mul_n_basecase (prodp, up, vp, size)
57: mp_ptr prodp;
58: mp_srcptr up;
59: mp_srcptr vp;
60: mp_size_t size;
61: #endif
62: {
63: mp_size_t i;
64: mp_limb_t cy_limb;
65: mp_limb_t v_limb;
66:
67: /* Multiply by the first limb in V separately, as the result can be
68: stored (not added) to PROD. We also avoid a loop for zeroing. */
69: v_limb = vp[0];
70: if (v_limb <= 1)
71: {
72: if (v_limb == 1)
73: MPN_COPY (prodp, up, size);
74: else
75: MPN_ZERO (prodp, size);
76: cy_limb = 0;
77: }
78: else
79: cy_limb = mpn_mul_1 (prodp, up, size, v_limb);
80:
81: prodp[size] = cy_limb;
82: prodp++;
83:
84: /* For each iteration in the outer loop, multiply one limb from
85: U with one limb from V, and add it to PROD. */
86: for (i = 1; i < size; i++)
87: {
88: v_limb = vp[i];
89: if (v_limb <= 1)
90: {
91: cy_limb = 0;
92: if (v_limb == 1)
93: cy_limb = mpn_add_n (prodp, prodp, up, size);
94: }
95: else
96: cy_limb = mpn_addmul_1 (prodp, up, size, v_limb);
97:
98: prodp[size] = cy_limb;
99: prodp++;
100: }
101: }
102:
103: void
104: #if __STDC__
105: impn_mul_n (mp_ptr prodp,
106: mp_srcptr up, mp_srcptr vp, mp_size_t size, mp_ptr tspace)
107: #else
108: impn_mul_n (prodp, up, vp, size, tspace)
109: mp_ptr prodp;
110: mp_srcptr up;
111: mp_srcptr vp;
112: mp_size_t size;
113: mp_ptr tspace;
114: #endif
115: {
116: if ((size & 1) != 0)
117: {
118: /* The size is odd, the code code below doesn't handle that.
119: Multiply the least significant (size - 1) limbs with a recursive
120: call, and handle the most significant limb of S1 and S2
121: separately. */
122: /* A slightly faster way to do this would be to make the Karatsuba
123: code below behave as if the size were even, and let it check for
124: odd size in the end. I.e., in essence move this code to the end.
125: Doing so would save us a recursive call, and potentially make the
126: stack grow a lot less. */
127:
128: mp_size_t esize = size - 1; /* even size */
129: mp_limb_t cy_limb;
130:
131: MPN_MUL_N_RECURSE (prodp, up, vp, esize, tspace);
132: cy_limb = mpn_addmul_1 (prodp + esize, up, esize, vp[esize]);
133: prodp[esize + esize] = cy_limb;
134: cy_limb = mpn_addmul_1 (prodp + esize, vp, size, up[esize]);
135:
136: prodp[esize + size] = cy_limb;
137: }
138: else
139: {
140: /* Anatolij Alekseevich Karatsuba's divide-and-conquer algorithm.
141:
142: Split U in two pieces, U1 and U0, such that
143: U = U0 + U1*(B**n),
144: and V in V1 and V0, such that
145: V = V0 + V1*(B**n).
146:
147: UV is then computed recursively using the identity
148:
149: 2n n n n
150: UV = (B + B )U V + B (U -U )(V -V ) + (B + 1)U V
151: 1 1 1 0 0 1 0 0
152:
153: Where B = 2**BITS_PER_MP_LIMB. */
154:
155: mp_size_t hsize = size >> 1;
156: mp_limb_t cy;
157: int negflg;
158:
159: /*** Product H. ________________ ________________
160: |_____U1 x V1____||____U0 x V0_____| */
161: /* Put result in upper part of PROD and pass low part of TSPACE
162: as new TSPACE. */
163: MPN_MUL_N_RECURSE (prodp + size, up + hsize, vp + hsize, hsize, tspace);
164:
165: /*** Product M. ________________
166: |_(U1-U0)(V0-V1)_| */
167: if (mpn_cmp (up + hsize, up, hsize) >= 0)
168: {
169: mpn_sub_n (prodp, up + hsize, up, hsize);
170: negflg = 0;
171: }
172: else
173: {
174: mpn_sub_n (prodp, up, up + hsize, hsize);
175: negflg = 1;
176: }
177: if (mpn_cmp (vp + hsize, vp, hsize) >= 0)
178: {
179: mpn_sub_n (prodp + hsize, vp + hsize, vp, hsize);
180: negflg ^= 1;
181: }
182: else
183: {
184: mpn_sub_n (prodp + hsize, vp, vp + hsize, hsize);
185: /* No change of NEGFLG. */
186: }
187: /* Read temporary operands from low part of PROD.
188: Put result in low part of TSPACE using upper part of TSPACE
189: as new TSPACE. */
190: MPN_MUL_N_RECURSE (tspace, prodp, prodp + hsize, hsize, tspace + size);
191:
192: /*** Add/copy product H. */
193: MPN_COPY (prodp + hsize, prodp + size, hsize);
194: cy = mpn_add_n (prodp + size, prodp + size, prodp + size + hsize, hsize);
195:
196: /*** Add product M (if NEGFLG M is a negative number). */
197: if (negflg)
198: cy -= mpn_sub_n (prodp + hsize, prodp + hsize, tspace, size);
199: else
200: cy += mpn_add_n (prodp + hsize, prodp + hsize, tspace, size);
201:
202: /*** Product L. ________________ ________________
203: |________________||____U0 x V0_____| */
204: /* Read temporary operands from low part of PROD.
205: Put result in low part of TSPACE using upper part of TSPACE
206: as new TSPACE. */
207: MPN_MUL_N_RECURSE (tspace, up, vp, hsize, tspace + size);
208:
209: /*** Add/copy Product L (twice). */
210:
211: cy += mpn_add_n (prodp + hsize, prodp + hsize, tspace, size);
212: if (cy)
213: mpn_add_1 (prodp + hsize + size, prodp + hsize + size, hsize, cy);
214:
215: MPN_COPY (prodp, tspace, hsize);
216: cy = mpn_add_n (prodp + hsize, prodp + hsize, tspace + hsize, hsize);
217: if (cy)
218: mpn_add_1 (prodp + size, prodp + size, size, 1);
219: }
220: }
221:
222: void
223: #if __STDC__
224: impn_sqr_n_basecase (mp_ptr prodp, mp_srcptr up, mp_size_t size)
225: #else
226: impn_sqr_n_basecase (prodp, up, size)
227: mp_ptr prodp;
228: mp_srcptr up;
229: mp_size_t size;
230: #endif
231: {
232: mp_size_t i;
233: mp_limb_t cy_limb;
234: mp_limb_t v_limb;
235:
236: /* Multiply by the first limb in V separately, as the result can be
237: stored (not added) to PROD. We also avoid a loop for zeroing. */
238: v_limb = up[0];
239: if (v_limb <= 1)
240: {
241: if (v_limb == 1)
242: MPN_COPY (prodp, up, size);
243: else
244: MPN_ZERO (prodp, size);
245: cy_limb = 0;
246: }
247: else
248: cy_limb = mpn_mul_1 (prodp, up, size, v_limb);
249:
250: prodp[size] = cy_limb;
251: prodp++;
252:
253: /* For each iteration in the outer loop, multiply one limb from
254: U with one limb from V, and add it to PROD. */
255: for (i = 1; i < size; i++)
256: {
257: v_limb = up[i];
258: if (v_limb <= 1)
259: {
260: cy_limb = 0;
261: if (v_limb == 1)
262: cy_limb = mpn_add_n (prodp, prodp, up, size);
263: }
264: else
265: cy_limb = mpn_addmul_1 (prodp, up, size, v_limb);
266:
267: prodp[size] = cy_limb;
268: prodp++;
269: }
270: }
271:
272: void
273: #if __STDC__
274: impn_sqr_n (mp_ptr prodp,
275: mp_srcptr up, mp_size_t size, mp_ptr tspace)
276: #else
277: impn_sqr_n (prodp, up, size, tspace)
278: mp_ptr prodp;
279: mp_srcptr up;
280: mp_size_t size;
281: mp_ptr tspace;
282: #endif
283: {
284: if ((size & 1) != 0)
285: {
286: /* The size is odd, the code code below doesn't handle that.
287: Multiply the least significant (size - 1) limbs with a recursive
288: call, and handle the most significant limb of S1 and S2
289: separately. */
290: /* A slightly faster way to do this would be to make the Karatsuba
291: code below behave as if the size were even, and let it check for
292: odd size in the end. I.e., in essence move this code to the end.
293: Doing so would save us a recursive call, and potentially make the
294: stack grow a lot less. */
295:
296: mp_size_t esize = size - 1; /* even size */
297: mp_limb_t cy_limb;
298:
299: MPN_SQR_N_RECURSE (prodp, up, esize, tspace);
300: cy_limb = mpn_addmul_1 (prodp + esize, up, esize, up[esize]);
301: prodp[esize + esize] = cy_limb;
302: cy_limb = mpn_addmul_1 (prodp + esize, up, size, up[esize]);
303:
304: prodp[esize + size] = cy_limb;
305: }
306: else
307: {
308: mp_size_t hsize = size >> 1;
309: mp_limb_t cy;
310:
311: /*** Product H. ________________ ________________
312: |_____U1 x U1____||____U0 x U0_____| */
313: /* Put result in upper part of PROD and pass low part of TSPACE
314: as new TSPACE. */
315: MPN_SQR_N_RECURSE (prodp + size, up + hsize, hsize, tspace);
316:
317: /*** Product M. ________________
318: |_(U1-U0)(U0-U1)_| */
319: if (mpn_cmp (up + hsize, up, hsize) >= 0)
320: {
321: mpn_sub_n (prodp, up + hsize, up, hsize);
322: }
323: else
324: {
325: mpn_sub_n (prodp, up, up + hsize, hsize);
326: }
327:
328: /* Read temporary operands from low part of PROD.
329: Put result in low part of TSPACE using upper part of TSPACE
330: as new TSPACE. */
331: MPN_SQR_N_RECURSE (tspace, prodp, hsize, tspace + size);
332:
333: /*** Add/copy product H. */
334: MPN_COPY (prodp + hsize, prodp + size, hsize);
335: cy = mpn_add_n (prodp + size, prodp + size, prodp + size + hsize, hsize);
336:
337: /*** Add product M (if NEGFLG M is a negative number). */
338: cy -= mpn_sub_n (prodp + hsize, prodp + hsize, tspace, size);
339:
340: /*** Product L. ________________ ________________
341: |________________||____U0 x U0_____| */
342: /* Read temporary operands from low part of PROD.
343: Put result in low part of TSPACE using upper part of TSPACE
344: as new TSPACE. */
345: MPN_SQR_N_RECURSE (tspace, up, hsize, tspace + size);
346:
347: /*** Add/copy Product L (twice). */
348:
349: cy += mpn_add_n (prodp + hsize, prodp + hsize, tspace, size);
350: if (cy)
351: mpn_add_1 (prodp + hsize + size, prodp + hsize + size, hsize, cy);
352:
353: MPN_COPY (prodp, tspace, hsize);
354: cy = mpn_add_n (prodp + hsize, prodp + hsize, tspace + hsize, hsize);
355: if (cy)
356: mpn_add_1 (prodp + size, prodp + size, size, 1);
357: }
358: }
359:
360: /* This should be made into an inline function in gmp.h. */
361: inline void
362: #if __STDC__
363: mpn_mul_n (mp_ptr prodp, mp_srcptr up, mp_srcptr vp, mp_size_t size)
364: #else
365: mpn_mul_n (prodp, up, vp, size)
366: mp_ptr prodp;
367: mp_srcptr up;
368: mp_srcptr vp;
369: mp_size_t size;
370: #endif
371: {
372: TMP_DECL (marker);
373: TMP_MARK (marker);
374: if (up == vp)
375: {
376: if (size < KARATSUBA_THRESHOLD)
377: {
378: impn_sqr_n_basecase (prodp, up, size);
379: }
380: else
381: {
382: mp_ptr tspace;
383: tspace = (mp_ptr) TMP_ALLOC (2 * size * BYTES_PER_MP_LIMB);
384: impn_sqr_n (prodp, up, size, tspace);
385: }
386: }
387: else
388: {
389: if (size < KARATSUBA_THRESHOLD)
390: {
391: impn_mul_n_basecase (prodp, up, vp, size);
392: }
393: else
394: {
395: mp_ptr tspace;
396: tspace = (mp_ptr) TMP_ALLOC (2 * size * BYTES_PER_MP_LIMB);
397: impn_mul_n (prodp, up, vp, size, tspace);
398: }
399: }
400: TMP_FREE (marker);
401: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to mul_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