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