update from main archive 961220
[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, 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 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 }