changed the dir layout
[libgcrypt.git] / mpi / mpih-mul.c
1 /* mpihelp-mul.c  -  MPI helper functions
2  *      Copyright (c) 1997 by Werner Koch (dd9jn)
3  *      Copyright (C) 1994, 1996 Free Software Foundation, Inc.
4  *
5  * This file is part of G10.
6  *
7  * G10 is free software; you can redistribute it and/or modify
8  * it under the terms of the GNU General Public License as published by
9  * the Free Software Foundation; either version 2 of the License, or
10  * (at your option) any later version.
11  *
12  * G10 is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15  * GNU General Public License for more details.
16  *
17  * You should have received a copy of the GNU General Public License
18  * along with this program; if not, write to the Free Software
19  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA
20  *
21  * Note: This code is heavily based on the GNU MP Library.
22  *       Actually it's the same code with only minor changes in the
23  *       way the data is stored; this is to support the abstraction
24  *       of an optional secure memory allocation which may be used
25  *       to avoid revealing of sensitive data due to paging etc.
26  *       The GNU MP Library itself is published under the LGPL;
27  *       however I decided to publish this code under the plain GPL.
28  */
29
30 #include <config.h>
31 #include <stdio.h>
32 #include <stdlib.h>
33 #include "mpi-internal.h"
34 #include "longlong.h"
35
36 /* If KARATSUBA_THRESHOLD is not already defined, define it to a
37  * value which is good on most machines.  */
38 #ifndef KARATSUBA_THRESHOLD
39     #define KARATSUBA_THRESHOLD 32
40 #endif
41
42 /* The code can't handle KARATSUBA_THRESHOLD smaller than 2.  */
43 #if KARATSUBA_THRESHOLD < 2
44     #undef KARATSUBA_THRESHOLD
45     #define KARATSUBA_THRESHOLD 2
46 #endif
47
48
49 #define MPN_MUL_N_RECURSE(prodp, up, vp, size, tspace) \
50     do {                                                \
51         if( (size) < KARATSUBA_THRESHOLD )              \
52             mul_n_basecase (prodp, up, vp, size);       \
53         else                                            \
54             mul_n (prodp, up, vp, size, tspace);        \
55     } while (0);
56
57 #define MPN_SQR_N_RECURSE(prodp, up, size, tspace) \
58     do {                                            \
59         if ((size) < KARATSUBA_THRESHOLD)           \
60             sqr_n_basecase (prodp, up, size);       \
61         else                                        \
62             sqr_n (prodp, up, size, tspace);        \
63     } while (0);
64
65
66
67
68 /* Multiply the natural numbers u (pointed to by UP) and v (pointed to by VP),
69  * both with SIZE limbs, and store the result at PRODP.  2 * SIZE limbs are
70  * always stored.  Return the most significant limb.
71  *
72  * Argument constraints:
73  * 1. PRODP != UP and PRODP != VP, i.e. the destination
74  *    must be distinct from the multiplier and the multiplicand.
75  *
76  *
77  * Handle simple cases with traditional multiplication.
78  *
79  * This is the most critical code of multiplication.  All multiplies rely
80  * on this, both small and huge.  Small ones arrive here immediately.  Huge
81  * ones arrive here as this is the base case for Karatsuba's recursive
82  * algorithm below.
83  */
84
85 static mpi_limb_t
86 mul_n_basecase( mpi_ptr_t prodp, mpi_ptr_t up,
87                                  mpi_ptr_t vp, mpi_size_t size)
88 {
89     mpi_size_t i;
90     mpi_limb_t cy;
91     mpi_limb_t v_limb;
92
93     /* Multiply by the first limb in V separately, as the result can be
94      * stored (not added) to PROD.  We also avoid a loop for zeroing.  */
95     v_limb = vp[0];
96     if( v_limb <= 1 ) {
97         if( v_limb == 1 )
98             MPN_COPY( prodp, up, size );
99         else
100             MPN_ZERO( prodp, size );
101         cy = 0;
102     }
103     else
104         cy = mpihelp_mul_1( prodp, up, size, v_limb );
105
106     prodp[size] = cy;
107     prodp++;
108
109     /* For each iteration in the outer loop, multiply one limb from
110      * U with one limb from V, and add it to PROD.  */
111     for( i = 1; i < size; i++ ) {
112         v_limb = vp[i];
113         if( v_limb <= 1 ) {
114             cy = 0;
115             if( v_limb == 1 )
116                cy = mpihelp_add_n(prodp, prodp, up, size);
117         }
118         else
119             cy = mpihelp_addmul_1(prodp, up, size, v_limb);
120
121         prodp[size] = cy;
122         prodp++;
123     }
124
125     return cy;
126 }
127
128
129 static void
130 mul_n( mpi_ptr_t prodp, mpi_ptr_t up, mpi_ptr_t vp,
131                         mpi_size_t size, mpi_ptr_t tspace )
132 {
133     if( size & 1 ) {
134       /* The size is odd, the code code below doesn't handle that.
135        * Multiply the least significant (size - 1) limbs with a recursive
136        * call, and handle the most significant limb of S1 and S2
137        * separately.
138        * A slightly faster way to do this would be to make the Karatsuba
139        * code below behave as if the size were even, and let it check for
140        * odd size in the end.  I.e., in essence move this code to the end.
141        * Doing so would save us a recursive call, and potentially make the
142        * stack grow a lot less.
143        */
144       mpi_size_t esize = size - 1;       /* even size */
145       mpi_limb_t cy_limb;
146
147       MPN_MUL_N_RECURSE( prodp, up, vp, esize, tspace );
148       cy_limb = mpihelp_addmul_1( prodp + esize, up, esize, vp[esize] );
149       prodp[esize + esize] = cy_limb;
150       cy_limb = mpihelp_addmul_1( prodp + esize, vp, size, up[esize] );
151       prodp[esize + size] = cy_limb;
152     }
153     else {
154         /* Anatolij Alekseevich Karatsuba's divide-and-conquer algorithm.
155          *
156          * Split U in two pieces, U1 and U0, such that
157          * U = U0 + U1*(B**n),
158          * and V in V1 and V0, such that
159          * V = V0 + V1*(B**n).
160          *
161          * UV is then computed recursively using the identity
162          *
163          *        2n   n          n                     n
164          * UV = (B  + B )U V  +  B (U -U )(V -V )  +  (B + 1)U V
165          *                1 1        1  0   0  1              0 0
166          *
167          * Where B = 2**BITS_PER_MP_LIMB.
168          */
169         mpi_size_t hsize = size >> 1;
170         mpi_limb_t cy;
171         int negflg;
172
173         /* Product H.      ________________  ________________
174          *                |_____U1 x V1____||____U0 x V0_____|
175          * Put result in upper part of PROD and pass low part of TSPACE
176          * as new TSPACE.
177          */
178         MPN_MUL_N_RECURSE(prodp + size, up + hsize, vp + hsize, hsize, tspace);
179
180         /* Product M.      ________________
181          *                |_(U1-U0)(V0-V1)_|
182          */
183         if( mpihelp_cmp(up + hsize, up, hsize) >= 0 ) {
184             mpihelp_sub_n(prodp, up + hsize, up, hsize);
185             negflg = 0;
186         }
187         else {
188             mpihelp_sub_n(prodp, up, up + hsize, hsize);
189             negflg = 1;
190         }
191         if( mpihelp_cmp(vp + hsize, vp, hsize) >= 0 ) {
192             mpihelp_sub_n(prodp + hsize, vp + hsize, vp, hsize);
193             negflg ^= 1;
194         }
195         else {
196             mpihelp_sub_n(prodp + hsize, vp, vp + hsize, hsize);
197             /* No change of NEGFLG.  */
198         }
199         /* Read temporary operands from low part of PROD.
200          * Put result in low part of TSPACE using upper part of TSPACE
201          * as new TSPACE.
202          */
203         MPN_MUL_N_RECURSE(tspace, prodp, prodp + hsize, hsize, tspace + size);
204
205         /* Add/copy product H. */
206         MPN_COPY (prodp + hsize, prodp + size, hsize);
207         cy = mpihelp_add_n( prodp + size, prodp + size,
208                             prodp + size + hsize, hsize);
209
210         /* Add product M (if NEGFLG M is a negative number) */
211         if(negflg)
212             cy -= mpihelp_sub_n(prodp + hsize, prodp + hsize, tspace, size);
213         else
214             cy += mpihelp_add_n(prodp + hsize, prodp + hsize, tspace, size);
215
216         /* Product L.      ________________  ________________
217          *                |________________||____U0 x V0_____|
218          * Read temporary operands from low part of PROD.
219          * Put result in low part of TSPACE using upper part of TSPACE
220          * as new TSPACE.
221          */
222         MPN_MUL_N_RECURSE(tspace, up, vp, hsize, tspace + size);
223
224         /* Add/copy Product L (twice) */
225
226         cy += mpihelp_add_n(prodp + hsize, prodp + hsize, tspace, size);
227         if( cy )
228           mpihelp_add_1(prodp + hsize + size, prodp + hsize + size, hsize, cy);
229
230         MPN_COPY(prodp, tspace, hsize);
231         cy = mpihelp_add_n(prodp + hsize, prodp + hsize, tspace + hsize, hsize);
232         if( cy )
233             mpihelp_add_1(prodp + size, prodp + size, size, 1);
234     }
235 }
236
237
238 static void
239 sqr_n_basecase( mpi_ptr_t prodp, mpi_ptr_t up, mpi_size_t size )
240 {
241     mpi_size_t i;
242     mpi_limb_t cy_limb;
243     mpi_limb_t v_limb;
244
245     /* Multiply by the first limb in V separately, as the result can be
246      * stored (not added) to PROD.  We also avoid a loop for zeroing.  */
247     v_limb = up[0];
248     if( v_limb <= 1 ) {
249         if( v_limb == 1 )
250             MPN_COPY( prodp, up, size );
251         else
252             MPN_ZERO(prodp, size);
253         cy_limb = 0;
254     }
255     else
256         cy_limb = mpihelp_mul_1( prodp, up, size, v_limb );
257
258     prodp[size] = cy_limb;
259     prodp++;
260
261     /* For each iteration in the outer loop, multiply one limb from
262      * U with one limb from V, and add it to PROD.  */
263     for( i=1; i < size; i++) {
264         v_limb = up[i];
265         if( v_limb <= 1 ) {
266             cy_limb = 0;
267             if( v_limb == 1 )
268                 cy_limb = mpihelp_add_n(prodp, prodp, up, size);
269         }
270         else
271             cy_limb = mpihelp_addmul_1(prodp, up, size, v_limb);
272
273         prodp[size] = cy_limb;
274         prodp++;
275     }
276 }
277
278
279 static void
280 sqr_n( mpi_ptr_t prodp, mpi_ptr_t up, mpi_size_t size, mpi_ptr_t tspace)
281 {
282     if( size & 1 ) {
283         /* The size is odd, the code code below doesn't handle that.
284          * Multiply the least significant (size - 1) limbs with a recursive
285          * call, and handle the most significant limb of S1 and S2
286          * separately.
287          * A slightly faster way to do this would be to make the Karatsuba
288          * code below behave as if the size were even, and let it check for
289          * odd size in the end.  I.e., in essence move this code to the end.
290          * Doing so would save us a recursive call, and potentially make the
291          * stack grow a lot less.
292          */
293         mpi_size_t esize = size - 1;       /* even size */
294         mpi_limb_t cy_limb;
295
296         MPN_SQR_N_RECURSE( prodp, up, esize, tspace );
297         cy_limb = mpihelp_addmul_1( prodp + esize, up, esize, up[esize] );
298         prodp[esize + esize] = cy_limb;
299         cy_limb = mpihelp_addmul_1( prodp + esize, up, size, up[esize] );
300
301         prodp[esize + size] = cy_limb;
302     }
303     else {
304         mpi_size_t hsize = size >> 1;
305         mpi_limb_t cy;
306
307         /* Product H.      ________________  ________________
308          *                |_____U1 x U1____||____U0 x U0_____|
309          * Put result in upper part of PROD and pass low part of TSPACE
310          * as new TSPACE.
311          */
312         MPN_SQR_N_RECURSE(prodp + size, up + hsize, hsize, tspace);
313
314         /* Product M.      ________________
315          *                |_(U1-U0)(U0-U1)_|
316          */
317         if( mpihelp_cmp( up + hsize, up, hsize) >= 0 )
318             mpihelp_sub_n( prodp, up + hsize, up, hsize);
319         else
320             mpihelp_sub_n (prodp, up, up + hsize, hsize);
321
322         /* Read temporary operands from low part of PROD.
323          * Put result in low part of TSPACE using upper part of TSPACE
324          * as new TSPACE.  */
325         MPN_SQR_N_RECURSE(tspace, prodp, hsize, tspace + size);
326
327         /* Add/copy product H  */
328         MPN_COPY(prodp + hsize, prodp + size, hsize);
329         cy = mpihelp_add_n(prodp + size, prodp + size,
330                            prodp + size + hsize, hsize);
331
332         /* Add product M (if NEGFLG M is a negative number).  */
333         cy -= mpihelp_sub_n (prodp + hsize, prodp + hsize, tspace, size);
334
335         /* Product L.      ________________  ________________
336          *                |________________||____U0 x U0_____|
337          * Read temporary operands from low part of PROD.
338          * Put result in low part of TSPACE using upper part of TSPACE
339          * as new TSPACE.  */
340         MPN_SQR_N_RECURSE (tspace, up, hsize, tspace + size);
341
342         /* Add/copy Product L (twice).  */
343         cy += mpihelp_add_n (prodp + hsize, prodp + hsize, tspace, size);
344         if( cy )
345             mpihelp_add_1(prodp + hsize + size, prodp + hsize + size,
346                                                             hsize, cy);
347
348         MPN_COPY(prodp, tspace, hsize);
349         cy = mpihelp_add_n (prodp + hsize, prodp + hsize, tspace + hsize, hsize);
350         if( cy )
351             mpihelp_add_1 (prodp + size, prodp + size, size, 1);
352     }
353 }
354
355
356 /* This should be made into an inline function in gmp.h.  */
357 void
358 mpihelp_mul_n( mpi_ptr_t prodp, mpi_ptr_t up, mpi_ptr_t vp, mpi_size_t size)
359 {
360     /* FIXME: mpi_alloc_limb_space, secure arg is wrong! */
361
362     if( up == vp ) {
363         if( size < KARATSUBA_THRESHOLD )
364             sqr_n_basecase( prodp, up, size );
365         else {
366             mpi_ptr_t tspace;
367             tspace = mpi_alloc_limb_space( 2 * size, 0 );
368             sqr_n( prodp, up, size, tspace );
369             mpi_free_limb_space( tspace );
370         }
371     }
372     else {
373         if( size < KARATSUBA_THRESHOLD )
374             mul_n_basecase( prodp, up, vp, size );
375         else {
376             mpi_ptr_t tspace;
377             tspace = mpi_alloc_limb_space( 2 * size, 0 );
378             mul_n (prodp, up, vp, size, tspace);
379             mpi_free_limb_space( tspace );
380         }
381     }
382 }
383
384
385 /* Multiply the natural numbers u (pointed to by UP, with USIZE limbs)
386  * and v (pointed to by VP, with VSIZE limbs), and store the result at
387  * PRODP.  USIZE + VSIZE limbs are always stored, but if the input
388  * operands are normalized.  Return the most significant limb of the
389  * result.
390  *
391  * NOTE: The space pointed to by PRODP is overwritten before finished
392  * with U and V, so overlap is an error.
393  *
394  * Argument constraints:
395  * 1. USIZE >= VSIZE.
396  * 2. PRODP != UP and PRODP != VP, i.e. the destination
397  *    must be distinct from the multiplier and the multiplicand.
398  */
399
400 mpi_limb_t
401 mpihelp_mul( mpi_ptr_t prodp, mpi_ptr_t up, mpi_size_t usize,
402                               mpi_ptr_t vp, mpi_size_t vsize)
403 {
404     mpi_ptr_t prod_endp = prodp + usize + vsize - 1;
405     mpi_limb_t cy;
406     mpi_ptr_t tspace;
407
408     if( vsize < KARATSUBA_THRESHOLD ) {
409         mpi_size_t i;
410         mpi_limb_t v_limb;
411
412         if( !vsize )
413             return 0;
414
415         /* Multiply by the first limb in V separately, as the result can be
416          * stored (not added) to PROD.  We also avoid a loop for zeroing.  */
417         v_limb = vp[0];
418         if( v_limb <= 1 ) {
419             if( v_limb == 1 )
420                 MPN_COPY( prodp, up, usize );
421             else
422                 MPN_ZERO( prodp, usize );
423             cy = 0;
424         }
425         else
426             cy = mpihelp_mul_1( prodp, up, usize, v_limb );
427
428         prodp[usize] = cy;
429         prodp++;
430
431         /* For each iteration in the outer loop, multiply one limb from
432          * U with one limb from V, and add it to PROD.  */
433         for( i = 1; i < vsize; i++ ) {
434             v_limb = vp[i];
435             if( v_limb <= 1 ) {
436                 cy = 0;
437                 if( v_limb == 1 )
438                    cy = mpihelp_add_n(prodp, prodp, up, usize);
439             }
440             else
441                 cy = mpihelp_addmul_1(prodp, up, usize, v_limb);
442
443             prodp[usize] = cy;
444             prodp++;
445         }
446
447         return cy;
448     }
449
450     /* FIXME: mpi_alloc_limb_space, secure arg is wrong! */
451     tspace = mpi_alloc_limb_space( 2 * vsize, 0 );
452     MPN_MUL_N_RECURSE( prodp, up, vp, vsize, tspace );
453
454     prodp += vsize;
455     up += vsize;
456     usize -= vsize;
457     if( usize >= vsize ) {
458         /* FIXME: mpi_alloc_limb_space, secure arg is wrong! */
459         mpi_ptr_t tp = mpi_alloc_limb_space( 2 * vsize, 0 );
460         do {
461             MPN_MUL_N_RECURSE( tp, up, vp, vsize, tspace );
462             cy = mpihelp_add_n( prodp, prodp, tp, vsize );
463             mpihelp_add_1( prodp + vsize, tp + vsize, vsize, cy );
464             prodp += vsize;
465             up += vsize;
466             usize -= vsize;
467         } while( usize >= vsize );
468         mpi_free_limb_space( tp );
469     }
470
471     if( usize ) {
472         mpihelp_mul( tspace, vp, vsize, up, usize );
473         cy = mpihelp_add_n( prodp, prodp, tspace, vsize);
474         mpihelp_add_1( prodp + vsize, tspace + vsize, usize, cy );
475     }
476
477     mpi_free_limb_space( tspace );
478     return *prod_endp;
479 }
480
481