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