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