fcf500c3660212582687a69dc0db062cc4eddc21
[gnupg.git] / mpi / mpi-pow.c
1 /* mpi-pow.c  -  MPI 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 #include <assert.h>
36
37
38 /****************
39  * RES = BASE ^ EXP mod MOD
40  */
41 void
42 mpi_powm( MPI res, MPI base, MPI exp, MPI mod)
43 {
44     mpi_ptr_t  rp, ep, mp, bp;
45     mpi_size_t esize, msize, bsize, rsize;
46     int        esign, msign, bsign, rsign;
47     int        esec,  msec,  bsec,  rsec;
48     mpi_size_t size;
49     int mod_shift_cnt;
50     int negative_result;
51     mpi_ptr_t mp_marker=NULL, bp_marker=NULL, ep_marker=NULL;
52     mpi_ptr_t xp_marker=NULL;
53     int assign_rp=0;
54     mpi_ptr_t tspace = NULL;
55     mpi_size_t tsize;
56
57     esize = exp->nlimbs;
58     msize = mod->nlimbs;
59     size = 2 * msize;
60     esign = exp->sign;
61     msign = mod->sign;
62
63     esec = exp->secure;
64     msec = mod->secure;
65     bsec = base->secure;
66     rsec = res->secure;
67
68     rp = res->d;
69     ep = exp->d;
70
71     if( !msize )
72         msize = 1 / msize;          /* provoke a signal */
73
74     if( !esize ) {
75         /* Exponent is zero, result is 1 mod MOD, i.e., 1 or 0
76          * depending on if MOD equals 1.  */
77         rp[0] = 1;
78         res->nlimbs = (msize == 1 && mod->d[0] == 1) ? 0 : 1;
79         res->sign = 0;
80         goto leave;
81     }
82
83     /* Normalize MOD (i.e. make its most significant bit set) as required by
84      * mpn_divrem.  This will make the intermediate values in the calculation
85      * slightly larger, but the correct result is obtained after a final
86      * reduction using the original MOD value.  */
87     mp = mp_marker = mpi_alloc_limb_space(msize, msec);
88     count_leading_zeros( mod_shift_cnt, mod->d[msize-1] );
89     if( mod_shift_cnt )
90         mpihelp_lshift( mp, mod->d, msize, mod_shift_cnt );
91     else
92         MPN_COPY( mp, mod->d, msize );
93
94     bsize = base->nlimbs;
95     bsign = base->sign;
96     if( bsize > msize ) { /* The base is larger than the module. Reduce it. */
97         /* Allocate (BSIZE + 1) with space for remainder and quotient.
98          * (The quotient is (bsize - msize + 1) limbs.)  */
99         bp = bp_marker = mpi_alloc_limb_space( bsize + 1, bsec );
100         MPN_COPY( bp, base->d, bsize );
101         /* We don't care about the quotient, store it above the remainder,
102          * at BP + MSIZE.  */
103         mpihelp_divrem( bp + msize, 0, bp, bsize, mp, msize );
104         bsize = msize;
105         /* Canonicalize the base, since we are going to multiply with it
106          * quite a few times.  */
107         MPN_NORMALIZE( bp, bsize );
108     }
109     else
110         bp = base->d;
111
112     if( !bsize ) {
113         res->nlimbs = 0;
114         res->sign = 0;
115         goto leave;
116     }
117
118     if( res->alloced < size ) {
119         /* We have to allocate more space for RES.  If any of the input
120          * parameters are identical to RES, defer deallocation of the old
121          * space.  */
122         if( rp == ep || rp == mp || rp == bp ) {
123             rp = mpi_alloc_limb_space( size, rsec );
124             assign_rp = 1;
125         }
126         else {
127             mpi_resize( res, size );
128             rp = res->d;
129         }
130     }
131     else { /* Make BASE, EXP and MOD not overlap with RES.  */
132         if( rp == bp ) {
133             /* RES and BASE are identical.  Allocate temp. space for BASE.  */
134             assert( !bp_marker );
135             bp = bp_marker = mpi_alloc_limb_space( bsize, bsec );
136             MPN_COPY(bp, rp, bsize);
137         }
138         if( rp == ep ) {
139             /* RES and EXP are identical.  Allocate temp. space for EXP.  */
140             ep = ep_marker = mpi_alloc_limb_space( esize, esec );
141             MPN_COPY(ep, rp, esize);
142         }
143         if( rp == mp ) {
144             /* RES and MOD are identical.  Allocate temporary space for MOD.*/
145             assert( !mp_marker );
146             mp = mp_marker = mpi_alloc_limb_space( msize, msec );
147             MPN_COPY(mp, rp, msize);
148         }
149     }
150
151     MPN_COPY( rp, bp, bsize );
152     rsize = bsize;
153     rsign = bsign;
154
155     {
156         mpi_size_t i;
157         mpi_ptr_t xp = xp_marker = mpi_alloc_limb_space( 2 * (msize + 1), msec );
158         int c;
159         mpi_limb_t e;
160         mpi_limb_t carry_limb;
161
162         negative_result = (ep[0] & 1) && base->sign;
163
164         i = esize - 1;
165         e = ep[i];
166         count_leading_zeros (c, e);
167         e = (e << c) << 1;     /* shift the exp bits to the left, lose msb */
168         c = BITS_PER_MPI_LIMB - 1 - c;
169
170         /* Main loop.
171          *
172          * Make the result be pointed to alternately by XP and RP.  This
173          * helps us avoid block copying, which would otherwise be necessary
174          * with the overlap restrictions of mpihelp_divmod. With 50% probability
175          * the result after this loop will be in the area originally pointed
176          * by RP (==RES->d), and with 50% probability in the area originally
177          * pointed to by XP.
178          */
179         for(;;) {
180             while( c ) {
181                 mpi_ptr_t tp;
182                 mpi_size_t xsize;
183
184                 /*mpihelp_mul_n(xp, rp, rp, rsize);*/
185                 if( rsize < KARATSUBA_THRESHOLD )
186                     mpih_sqr_n_basecase( xp, rp, rsize );
187                 else {
188                     if( !tspace ) {
189                         tsize = 2 * rsize;
190                         tspace = mpi_alloc_limb_space( tsize, 0 );
191                     }
192                     else if( tsize < (2*rsize) ) {
193                         mpi_free_limb_space( tspace );
194                         tsize = 2 * rsize;
195                         tspace = mpi_alloc_limb_space( tsize, 0 );
196
197                     }
198                     mpih_sqr_n( xp, rp, rsize, tspace );
199                 }
200
201                 xsize = 2 * rsize;
202                 if( xsize > msize ) {
203                     mpihelp_divrem(xp + msize, 0, xp, xsize, mp, msize);
204                     xsize = msize;
205                 }
206
207                 tp = rp; rp = xp; xp = tp;
208                 rsize = xsize;
209
210                 if( (mpi_limb_signed_t)e < 0 ) {
211                     mpihelp_mul( xp, rp, rsize, bp, bsize );
212                     xsize = rsize + bsize;
213                     if( xsize > msize ) {
214                         mpihelp_divrem(xp + msize, 0, xp, xsize, mp, msize);
215                         xsize = msize;
216                     }
217
218                     tp = rp; rp = xp; xp = tp;
219                     rsize = xsize;
220                 }
221                 e <<= 1;
222                 c--;
223             }
224
225             i--;
226             if( i < 0 )
227                 break;
228             e = ep[i];
229             c = BITS_PER_MPI_LIMB;
230         }
231
232         /* We shifted MOD, the modulo reduction argument, left MOD_SHIFT_CNT
233          * steps.  Adjust the result by reducing it with the original MOD.
234          *
235          * Also make sure the result is put in RES->d (where it already
236          * might be, see above).
237          */
238         if( mod_shift_cnt ) {
239             carry_limb = mpihelp_lshift( res->d, rp, rsize, mod_shift_cnt);
240             rp = res->d;
241             if( carry_limb ) {
242                 rp[rsize] = carry_limb;
243                 rsize++;
244             }
245         }
246         else {
247             MPN_COPY( res->d, rp, rsize);
248             rp = res->d;
249         }
250
251         if( rsize >= msize ) {
252             mpihelp_divrem(rp + msize, 0, rp, rsize, mp, msize);
253             rsize = msize;
254         }
255
256         /* Remove any leading zero words from the result.  */
257         if( mod_shift_cnt )
258             mpihelp_rshift( rp, rp, rsize, mod_shift_cnt);
259         MPN_NORMALIZE (rp, rsize);
260     }
261
262     if( negative_result && rsize ) {
263         if( mod_shift_cnt )
264             mpihelp_rshift( mp, mp, msize, mod_shift_cnt);
265         mpihelp_sub( rp, mp, msize, rp, rsize);
266         rsize = msize;
267         rsign = msign;
268         MPN_NORMALIZE(rp, rsize);
269     }
270     res->nlimbs = rsize;
271     res->sign = rsign;
272
273   leave:
274     if( assign_rp ) mpi_assign_limb_space( res, rp, size );
275     if( mp_marker ) mpi_free_limb_space( mp_marker );
276     if( bp_marker ) mpi_free_limb_space( bp_marker );
277     if( ep_marker ) mpi_free_limb_space( ep_marker );
278     if( xp_marker ) mpi_free_limb_space( xp_marker );
279     if( tspace )    mpi_free_limb_space( tspace );
280 }
281