mpi: mpi-pow improvement.
[libgcrypt.git] / mpi / mpi-pow.c
1 /* mpi-pow.c  -  MPI functions for exponentiation
2  * Copyright (C) 1994, 1996, 1998, 2000, 2002
3  *               2003  Free Software Foundation, Inc.
4  *               2013  g10 Code GmbH
5  *
6  * This file is part of Libgcrypt.
7  *
8  * Libgcrypt is free software; you can redistribute it and/or modify
9  * it under the terms of the GNU Lesser General Public License as
10  * published by the Free Software Foundation; either version 2.1 of
11  * the License, or (at your option) any later version.
12  *
13  * Libgcrypt is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16  * GNU Lesser General Public License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public
19  * License along with this program; if not, see <http://www.gnu.org/licenses/>.
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  */
27
28 #include <config.h>
29 #include <stdio.h>
30 #include <stdlib.h>
31 #include <string.h>
32
33 #include "mpi-internal.h"
34 #include "longlong.h"
35
36
37 /*
38  * When you need old implementation, please add compilation option
39  * -DUSE_ALGORITHM_SIMPLE_EXPONENTIATION
40  * or expose this line:
41 #define USE_ALGORITHM_SIMPLE_EXPONENTIATION 1
42  */
43
44 #if defined(USE_ALGORITHM_SIMPLE_EXPONENTIATION)
45 /****************
46  * RES = BASE ^ EXPO mod MOD
47  */
48 void
49 gcry_mpi_powm (gcry_mpi_t res,
50                gcry_mpi_t base, gcry_mpi_t expo, gcry_mpi_t mod)
51 {
52   /* Pointer to the limbs of the arguments, their size and signs. */
53   mpi_ptr_t  rp, ep, mp, bp;
54   mpi_size_t esize, msize, bsize, rsize;
55   int               msign, bsign, rsign;
56   /* Flags telling the secure allocation status of the arguments.  */
57   int        esec,  msec,  bsec;
58   /* Size of the result including space for temporary values.  */
59   mpi_size_t size;
60   /* Helper.  */
61   int mod_shift_cnt;
62   int negative_result;
63   mpi_ptr_t mp_marker = NULL;
64   mpi_ptr_t bp_marker = NULL;
65   mpi_ptr_t ep_marker = NULL;
66   mpi_ptr_t xp_marker = NULL;
67   unsigned int mp_nlimbs = 0;
68   unsigned int bp_nlimbs = 0;
69   unsigned int ep_nlimbs = 0;
70   unsigned int xp_nlimbs = 0;
71   mpi_ptr_t tspace = NULL;
72   mpi_size_t tsize = 0;
73
74
75   esize = expo->nlimbs;
76   msize = mod->nlimbs;
77   size = 2 * msize;
78   msign = mod->sign;
79
80   esec = mpi_is_secure(expo);
81   msec = mpi_is_secure(mod);
82   bsec = mpi_is_secure(base);
83
84   rp = res->d;
85   ep = expo->d;
86
87   if (!msize)
88     _gcry_divide_by_zero();
89
90   if (!esize)
91     {
92       /* Exponent is zero, result is 1 mod MOD, i.e., 1 or 0 depending
93          on if MOD equals 1.  */
94       res->nlimbs = (msize == 1 && mod->d[0] == 1) ? 0 : 1;
95       if (res->nlimbs)
96         {
97           RESIZE_IF_NEEDED (res, 1);
98           rp = res->d;
99           rp[0] = 1;
100         }
101       res->sign = 0;
102       goto leave;
103     }
104
105   /* Normalize MOD (i.e. make its most significant bit set) as
106      required by mpn_divrem.  This will make the intermediate values
107      in the calculation slightly larger, but the correct result is
108      obtained after a final reduction using the original MOD value. */
109   mp_nlimbs = msec? msize:0;
110   mp = mp_marker = mpi_alloc_limb_space(msize, msec);
111   count_leading_zeros (mod_shift_cnt, mod->d[msize-1]);
112   if (mod_shift_cnt)
113     _gcry_mpih_lshift (mp, mod->d, msize, mod_shift_cnt);
114   else
115     MPN_COPY( mp, mod->d, msize );
116
117   bsize = base->nlimbs;
118   bsign = base->sign;
119   if (bsize > msize)
120     {
121       /* The base is larger than the module.  Reduce it.
122
123          Allocate (BSIZE + 1) with space for remainder and quotient.
124          (The quotient is (bsize - msize + 1) limbs.)  */
125       bp_nlimbs = bsec ? (bsize + 1):0;
126       bp = bp_marker = mpi_alloc_limb_space( bsize + 1, bsec );
127       MPN_COPY ( bp, base->d, bsize );
128       /* We don't care about the quotient, store it above the
129        * remainder, at BP + MSIZE.  */
130       _gcry_mpih_divrem( bp + msize, 0, bp, bsize, mp, msize );
131       bsize = msize;
132       /* Canonicalize the base, since we are going to multiply with it
133          quite a few times.  */
134       MPN_NORMALIZE( bp, bsize );
135     }
136   else
137     bp = base->d;
138
139   if (!bsize)
140     {
141       res->nlimbs = 0;
142       res->sign = 0;
143       goto leave;
144     }
145
146
147   /* Make BASE, EXPO and MOD not overlap with RES.  */
148   if ( rp == bp )
149     {
150       /* RES and BASE are identical.  Allocate temp. space for BASE.  */
151       gcry_assert (!bp_marker);
152       bp_nlimbs = bsec? bsize:0;
153       bp = bp_marker = mpi_alloc_limb_space( bsize, bsec );
154       MPN_COPY(bp, rp, bsize);
155     }
156   if ( rp == ep )
157     {
158       /* RES and EXPO are identical.  Allocate temp. space for EXPO.  */
159       ep_nlimbs = esec? esize:0;
160       ep = ep_marker = mpi_alloc_limb_space( esize, esec );
161       MPN_COPY(ep, rp, esize);
162     }
163   if ( rp == mp )
164     {
165       /* RES and MOD are identical.  Allocate temporary space for MOD.*/
166       gcry_assert (!mp_marker);
167       mp_nlimbs = msec?msize:0;
168       mp = mp_marker = mpi_alloc_limb_space( msize, msec );
169       MPN_COPY(mp, rp, msize);
170     }
171
172   /* Copy base to the result.  */
173   if (res->alloced < size)
174     {
175       mpi_resize (res, size);
176       rp = res->d;
177     }
178   MPN_COPY ( rp, bp, bsize );
179   rsize = bsize;
180   rsign = bsign;
181
182   /* Main processing.  */
183   {
184     mpi_size_t i;
185     mpi_ptr_t xp;
186     int c;
187     mpi_limb_t e;
188     mpi_limb_t carry_limb;
189     struct karatsuba_ctx karactx;
190
191     xp_nlimbs = msec? (2 * (msize + 1)):0;
192     xp = xp_marker = mpi_alloc_limb_space( 2 * (msize + 1), msec );
193
194     memset( &karactx, 0, sizeof karactx );
195     negative_result = (ep[0] & 1) && base->sign;
196
197     i = esize - 1;
198     e = ep[i];
199     count_leading_zeros (c, e);
200     e = (e << c) << 1;     /* Shift the expo bits to the left, lose msb.  */
201     c = BITS_PER_MPI_LIMB - 1 - c;
202
203     /* Main loop.
204
205        Make the result be pointed to alternately by XP and RP.  This
206        helps us avoid block copying, which would otherwise be
207        necessary with the overlap restrictions of
208        _gcry_mpih_divmod. With 50% probability the result after this
209        loop will be in the area originally pointed by RP (==RES->d),
210        and with 50% probability in the area originally pointed to by XP. */
211     for (;;)
212       {
213         while (c)
214           {
215             mpi_ptr_t tp;
216             mpi_size_t xsize;
217
218             /*mpih_mul_n(xp, rp, rp, rsize);*/
219             if ( rsize < KARATSUBA_THRESHOLD )
220               _gcry_mpih_sqr_n_basecase( xp, rp, rsize );
221             else
222               {
223                 if ( !tspace )
224                   {
225                     tsize = 2 * rsize;
226                     tspace = mpi_alloc_limb_space( tsize, 0 );
227                   }
228                 else if ( tsize < (2*rsize) )
229                   {
230                     _gcry_mpi_free_limb_space (tspace, 0);
231                     tsize = 2 * rsize;
232                     tspace = mpi_alloc_limb_space (tsize, 0 );
233                   }
234                 _gcry_mpih_sqr_n (xp, rp, rsize, tspace);
235               }
236
237             xsize = 2 * rsize;
238             if ( xsize > msize )
239               {
240                 _gcry_mpih_divrem(xp + msize, 0, xp, xsize, mp, msize);
241                 xsize = msize;
242               }
243
244             tp = rp; rp = xp; xp = tp;
245             rsize = xsize;
246
247             /* To mitigate the Yarom/Falkner flush+reload cache
248              * side-channel attack on the RSA secret exponent, we do
249              * the multiplication regardless of the value of the
250              * high-bit of E.  But to avoid this performance penalty
251              * we do it only if the exponent has been stored in secure
252              * memory and we can thus assume it is a secret exponent.  */
253             if (esec || (mpi_limb_signed_t)e < 0)
254               {
255                 /*mpih_mul( xp, rp, rsize, bp, bsize );*/
256                 if( bsize < KARATSUBA_THRESHOLD )
257                   _gcry_mpih_mul ( xp, rp, rsize, bp, bsize );
258                 else
259                   _gcry_mpih_mul_karatsuba_case (xp, rp, rsize, bp, bsize,
260                                                  &karactx);
261
262                 xsize = rsize + bsize;
263                 if ( xsize > msize )
264                   {
265                     _gcry_mpih_divrem(xp + msize, 0, xp, xsize, mp, msize);
266                     xsize = msize;
267                   }
268               }
269             if ( (mpi_limb_signed_t)e < 0 )
270               {
271                 tp = rp; rp = xp; xp = tp;
272                 rsize = xsize;
273               }
274             e <<= 1;
275             c--;
276           }
277
278         i--;
279         if ( i < 0 )
280           break;
281         e = ep[i];
282         c = BITS_PER_MPI_LIMB;
283       }
284
285     /* We shifted MOD, the modulo reduction argument, left
286        MOD_SHIFT_CNT steps.  Adjust the result by reducing it with the
287        original MOD.
288
289        Also make sure the result is put in RES->d (where it already
290        might be, see above).  */
291     if ( mod_shift_cnt )
292       {
293         carry_limb = _gcry_mpih_lshift( res->d, rp, rsize, mod_shift_cnt);
294         rp = res->d;
295         if ( carry_limb )
296           {
297             rp[rsize] = carry_limb;
298             rsize++;
299           }
300       }
301     else if (res->d != rp)
302       {
303         MPN_COPY (res->d, rp, rsize);
304         rp = res->d;
305       }
306
307     if ( rsize >= msize )
308       {
309         _gcry_mpih_divrem(rp + msize, 0, rp, rsize, mp, msize);
310         rsize = msize;
311       }
312
313     /* Remove any leading zero words from the result.  */
314     if ( mod_shift_cnt )
315       _gcry_mpih_rshift( rp, rp, rsize, mod_shift_cnt);
316     MPN_NORMALIZE (rp, rsize);
317
318     _gcry_mpih_release_karatsuba_ctx (&karactx );
319   }
320
321   /* Fixup for negative results.  */
322   if ( negative_result && rsize )
323     {
324       if ( mod_shift_cnt )
325         _gcry_mpih_rshift( mp, mp, msize, mod_shift_cnt);
326       _gcry_mpih_sub( rp, mp, msize, rp, rsize);
327       rsize = msize;
328       rsign = msign;
329       MPN_NORMALIZE(rp, rsize);
330     }
331   gcry_assert (res->d == rp);
332   res->nlimbs = rsize;
333   res->sign = rsign;
334
335  leave:
336   if (mp_marker)
337     _gcry_mpi_free_limb_space( mp_marker, mp_nlimbs );
338   if (bp_marker)
339     _gcry_mpi_free_limb_space( bp_marker, bp_nlimbs );
340   if (ep_marker)
341     _gcry_mpi_free_limb_space( ep_marker, ep_nlimbs );
342   if (xp_marker)
343     _gcry_mpi_free_limb_space( xp_marker, xp_nlimbs );
344   if (tspace)
345     _gcry_mpi_free_limb_space( tspace, 0 );
346 }
347 #else
348 /**
349  * Internal function to compute
350  *
351  *    X = R * S mod M
352  *
353  * and set the size of X at the pointer XSIZE_P.
354  * Use karatsuba structure at KARACTX_P.
355  *
356  * Condition:
357  *   RSIZE >= SSIZE
358  *   Enough space for X is allocated beforehand.
359  *
360  * For generic cases, we can/should use gcry_mpi_mulm.
361  * This function is use for specific internal case.
362  */
363 static void
364 mul_mod (mpi_ptr_t xp, mpi_size_t *xsize_p,
365          mpi_ptr_t rp, mpi_size_t rsize,
366          mpi_ptr_t sp, mpi_size_t ssize,
367          mpi_ptr_t mp, mpi_size_t msize,
368          struct karatsuba_ctx *karactx_p)
369 {
370   if( ssize < KARATSUBA_THRESHOLD )
371     _gcry_mpih_mul ( xp, rp, rsize, sp, ssize );
372   else
373     _gcry_mpih_mul_karatsuba_case (xp, rp, rsize, sp, ssize, karactx_p);
374
375    if (rsize + ssize > msize)
376     {
377       _gcry_mpih_divrem (xp + msize, 0, xp, rsize + ssize, mp, msize);
378       *xsize_p = msize;
379     }
380    else
381      *xsize_p = rsize + ssize;
382 }
383
384 #define SIZE_B_2I3 ((1 << (5 - 1)) - 1)
385
386 /****************
387  * RES = BASE ^ EXPO mod MOD
388  *
389  * To mitigate the Yarom/Falkner flush+reload cache side-channel
390  * attack on the RSA secret exponent, we don't use the square
391  * routine but multiplication.
392  *
393  * Reference:
394  *   Handbook of Applied Cryptography
395  *       Algorithm 14.83: Modified left-to-right k-ary exponentiation
396  */
397 void
398 gcry_mpi_powm (gcry_mpi_t res,
399                gcry_mpi_t base, gcry_mpi_t expo, gcry_mpi_t mod)
400 {
401   /* Pointer to the limbs of the arguments, their size and signs. */
402   mpi_ptr_t  rp, ep, mp, bp;
403   mpi_size_t esize, msize, bsize, rsize;
404   int               msign, bsign, rsign;
405   /* Flags telling the secure allocation status of the arguments.  */
406   int        esec,  msec,  bsec;
407   /* Size of the result including space for temporary values.  */
408   mpi_size_t size;
409   /* Helper.  */
410   int mod_shift_cnt;
411   int negative_result;
412   mpi_ptr_t mp_marker = NULL;
413   mpi_ptr_t bp_marker = NULL;
414   mpi_ptr_t ep_marker = NULL;
415   mpi_ptr_t xp_marker = NULL;
416   unsigned int mp_nlimbs = 0;
417   unsigned int bp_nlimbs = 0;
418   unsigned int ep_nlimbs = 0;
419   unsigned int xp_nlimbs = 0;
420   mpi_ptr_t b_2i3[SIZE_B_2I3]; /* Pre-computed array: BASE^3, ^5, ^7, ... */
421   mpi_size_t b_2i3size[SIZE_B_2I3];
422   mpi_size_t W;
423   mpi_ptr_t base_u;
424   mpi_size_t base_u_size;
425
426   esize = expo->nlimbs;
427   msize = mod->nlimbs;
428   size = 2 * msize;
429   msign = mod->sign;
430
431   if (esize * BITS_PER_MPI_LIMB > 512)
432     W = 5;
433   else if (esize * BITS_PER_MPI_LIMB > 256)
434     W = 4;
435   else if (esize * BITS_PER_MPI_LIMB > 128)
436     W = 3;
437   else if (esize * BITS_PER_MPI_LIMB > 64)
438     W = 2;
439   else
440     W = 1;
441
442   esec = mpi_is_secure(expo);
443   msec = mpi_is_secure(mod);
444   bsec = mpi_is_secure(base);
445
446   rp = res->d;
447   ep = expo->d;
448
449   if (!msize)
450     _gcry_divide_by_zero();
451
452   if (!esize)
453     {
454       /* Exponent is zero, result is 1 mod MOD, i.e., 1 or 0 depending
455          on if MOD equals 1.  */
456       res->nlimbs = (msize == 1 && mod->d[0] == 1) ? 0 : 1;
457       if (res->nlimbs)
458         {
459           RESIZE_IF_NEEDED (res, 1);
460           rp = res->d;
461           rp[0] = 1;
462         }
463       res->sign = 0;
464       goto leave;
465     }
466
467   /* Normalize MOD (i.e. make its most significant bit set) as
468      required by mpn_divrem.  This will make the intermediate values
469      in the calculation slightly larger, but the correct result is
470      obtained after a final reduction using the original MOD value. */
471   mp_nlimbs = msec? msize:0;
472   mp = mp_marker = mpi_alloc_limb_space(msize, msec);
473   count_leading_zeros (mod_shift_cnt, mod->d[msize-1]);
474   if (mod_shift_cnt)
475     _gcry_mpih_lshift (mp, mod->d, msize, mod_shift_cnt);
476   else
477     MPN_COPY( mp, mod->d, msize );
478
479   bsize = base->nlimbs;
480   bsign = base->sign;
481   if (bsize > msize)
482     {
483       /* The base is larger than the module.  Reduce it.
484
485          Allocate (BSIZE + 1) with space for remainder and quotient.
486          (The quotient is (bsize - msize + 1) limbs.)  */
487       bp_nlimbs = bsec ? (bsize + 1):0;
488       bp = bp_marker = mpi_alloc_limb_space( bsize + 1, bsec );
489       MPN_COPY ( bp, base->d, bsize );
490       /* We don't care about the quotient, store it above the
491        * remainder, at BP + MSIZE.  */
492       _gcry_mpih_divrem( bp + msize, 0, bp, bsize, mp, msize );
493       bsize = msize;
494       /* Canonicalize the base, since we are going to multiply with it
495          quite a few times.  */
496       MPN_NORMALIZE( bp, bsize );
497     }
498   else
499     bp = base->d;
500
501   if (!bsize)
502     {
503       res->nlimbs = 0;
504       res->sign = 0;
505       goto leave;
506     }
507
508
509   /* Make BASE, EXPO and MOD not overlap with RES.  */
510   if ( rp == bp )
511     {
512       /* RES and BASE are identical.  Allocate temp. space for BASE.  */
513       gcry_assert (!bp_marker);
514       bp_nlimbs = bsec? bsize:0;
515       bp = bp_marker = mpi_alloc_limb_space( bsize, bsec );
516       MPN_COPY(bp, rp, bsize);
517     }
518   if ( rp == ep )
519     {
520       /* RES and EXPO are identical.  Allocate temp. space for EXPO.  */
521       ep_nlimbs = esec? esize:0;
522       ep = ep_marker = mpi_alloc_limb_space( esize, esec );
523       MPN_COPY(ep, rp, esize);
524     }
525   if ( rp == mp )
526     {
527       /* RES and MOD are identical.  Allocate temporary space for MOD.*/
528       gcry_assert (!mp_marker);
529       mp_nlimbs = msec?msize:0;
530       mp = mp_marker = mpi_alloc_limb_space( msize, msec );
531       MPN_COPY(mp, rp, msize);
532     }
533
534   /* Copy base to the result.  */
535   if (res->alloced < size)
536     {
537       mpi_resize (res, size);
538       rp = res->d;
539     }
540
541   /* Main processing.  */
542   {
543     mpi_size_t i, j;
544     mpi_ptr_t xp;
545     mpi_size_t xsize;
546     int c;
547     mpi_limb_t e;
548     mpi_limb_t carry_limb;
549     struct karatsuba_ctx karactx;
550     mpi_ptr_t tp;
551
552     xp_nlimbs = msec? (2 * (msize + 1)):0;
553     xp = xp_marker = mpi_alloc_limb_space( 2 * (msize + 1), msec );
554
555     memset( &karactx, 0, sizeof karactx );
556     negative_result = (ep[0] & 1) && bsign;
557
558     /* Precompute B_2I3[], BASE^(2 * i + 3), BASE^3, ^5, ^7, ... */
559     if (W > 1)                  /* X := BASE^2 */
560       mul_mod (xp, &xsize, bp, bsize, bp, bsize, mp, msize, &karactx);
561     for (i = 0; i < (1 << (W - 1)) - 1; i++)
562       {                         /* B_2I3[i] = BASE^(2 * i + 3) */
563         if (i == 0)
564           {
565             base_u = bp;
566             base_u_size = bsize;
567           }
568         else
569           {
570             base_u = b_2i3[i-1];
571             base_u_size = b_2i3size[i-1];
572           }
573
574         if (xsize >= base_u_size)
575           mul_mod (rp, &rsize, xp, xsize, base_u, base_u_size,
576                    mp, msize, &karactx);
577         else
578           mul_mod (rp, &rsize, base_u, base_u_size, xp, xsize,
579                    mp, msize, &karactx);
580         b_2i3[i] = mpi_alloc_limb_space (rsize, esec);
581         b_2i3size[i] = rsize;
582         MPN_COPY (b_2i3[i], rp, rsize);
583       }
584
585     i = esize - 1;
586
587     /* Main loop.
588
589        Make the result be pointed to alternately by XP and RP.  This
590        helps us avoid block copying, which would otherwise be
591        necessary with the overlap restrictions of
592        _gcry_mpih_divmod. With 50% probability the result after this
593        loop will be in the area originally pointed by RP (==RES->d),
594        and with 50% probability in the area originally pointed to by XP. */
595     rsign = 0;
596     if (W == 1)
597       {
598         rsize = bsize;
599       }
600     else
601       {
602         rsize = msize;
603         MPN_ZERO (rp, rsize);
604       }
605     MPN_COPY ( rp, bp, bsize );
606
607     e = ep[i];
608     count_leading_zeros (c, e);
609     e = (e << c) << 1;
610     c = BITS_PER_MPI_LIMB - 1 - c;
611
612     j = 0;
613
614     for (;;)
615       if (e == 0)
616         {
617           j += c;
618           i--;
619           if ( i < 0 )
620             {
621               c = 0;
622               break;
623             }
624
625           e = ep[i];
626           c = BITS_PER_MPI_LIMB;
627         }
628       else
629         {
630           int c0;
631           mpi_limb_t e0;
632
633           count_leading_zeros (c0, e);
634           e = (e << c0);
635           c -= c0;
636           j += c0;
637
638           if (c >= W)
639             {
640               e0 = (e >> (BITS_PER_MPI_LIMB - W));
641               e = (e << W);
642               c -= W;
643             }
644           else
645             {
646               i--;
647               if ( i < 0 )
648                 {
649                   e = (e >> (BITS_PER_MPI_LIMB - c));
650                   break;
651                 }
652
653               c0 = c;
654               e0 = (e >> (BITS_PER_MPI_LIMB - W))
655                 | (ep[i] >> (BITS_PER_MPI_LIMB - W + c0));
656               e = (ep[i] << (W - c0));
657               c = BITS_PER_MPI_LIMB - W + c0;
658             }
659
660           count_trailing_zeros (c0, e0);
661           e0 = (e0 >> c0) >> 1;
662
663           for (j += W - c0; j; j--)
664             {
665               mul_mod (xp, &xsize, rp, rsize, rp, rsize, mp, msize, &karactx);
666               tp = rp; rp = xp; xp = tp;
667               rsize = xsize;
668             }
669
670           if (e0 == 0)
671             {
672               base_u = bp;
673               base_u_size = bsize;
674             }
675           else
676             {
677               base_u = b_2i3[e0 - 1];
678               base_u_size = b_2i3size[e0 -1];
679             }
680
681           mul_mod (xp, &xsize, rp, rsize, base_u, base_u_size,
682                    mp, msize, &karactx);
683           tp = rp; rp = xp; xp = tp;
684           rsize = xsize;
685
686           j = c0;
687         }
688
689     if (c != 0)
690       {
691         j += c;
692         count_trailing_zeros (c, e);
693         e = (e >> c);
694         j -= c;
695       }
696
697     while (j--)
698       {
699         mul_mod (xp, &xsize, rp, rsize, rp, rsize, mp, msize, &karactx);
700         tp = rp; rp = xp; xp = tp;
701         rsize = xsize;
702       }
703
704     if (e != 0)
705       {
706         if ((e>>1) == 0)
707           {
708             base_u = bp;
709             base_u_size = bsize;
710           }
711         else
712           {
713             base_u = b_2i3[(e>>1) - 1];
714             base_u_size = b_2i3size[(e>>1) -1];
715           }
716
717         mul_mod (xp, &xsize, rp, rsize, base_u, base_u_size,
718                  mp, msize, &karactx);
719         tp = rp; rp = xp; xp = tp;
720         rsize = xsize;
721
722         for (; c; c--)
723           {
724             mul_mod (xp, &xsize, rp, rsize, rp, rsize, mp, msize, &karactx);
725             tp = rp; rp = xp; xp = tp;
726             rsize = xsize;
727           }
728       }
729
730     /* We shifted MOD, the modulo reduction argument, left
731        MOD_SHIFT_CNT steps.  Adjust the result by reducing it with the
732        original MOD.
733
734        Also make sure the result is put in RES->d (where it already
735        might be, see above).  */
736     if ( mod_shift_cnt )
737       {
738         carry_limb = _gcry_mpih_lshift( res->d, rp, rsize, mod_shift_cnt);
739         rp = res->d;
740         if ( carry_limb )
741           {
742             rp[rsize] = carry_limb;
743             rsize++;
744           }
745       }
746     else if (res->d != rp)
747       {
748         MPN_COPY (res->d, rp, rsize);
749         rp = res->d;
750       }
751
752     if ( rsize >= msize )
753       {
754         _gcry_mpih_divrem(rp + msize, 0, rp, rsize, mp, msize);
755         rsize = msize;
756       }
757
758     /* Remove any leading zero words from the result.  */
759     if ( mod_shift_cnt )
760       _gcry_mpih_rshift( rp, rp, rsize, mod_shift_cnt);
761     MPN_NORMALIZE (rp, rsize);
762
763     _gcry_mpih_release_karatsuba_ctx (&karactx );
764     for (i = 0; i < (1 << (W - 1)) - 1; i++)
765       _gcry_mpi_free_limb_space( b_2i3[i], esec ? b_2i3size[i] : 0 );
766   }
767
768   /* Fixup for negative results.  */
769   if ( negative_result && rsize )
770     {
771       if ( mod_shift_cnt )
772         _gcry_mpih_rshift( mp, mp, msize, mod_shift_cnt);
773       _gcry_mpih_sub( rp, mp, msize, rp, rsize);
774       rsize = msize;
775       rsign = msign;
776       MPN_NORMALIZE(rp, rsize);
777     }
778   gcry_assert (res->d == rp);
779   res->nlimbs = rsize;
780   res->sign = rsign;
781
782  leave:
783   if (mp_marker)
784     _gcry_mpi_free_limb_space( mp_marker, mp_nlimbs );
785   if (bp_marker)
786     _gcry_mpi_free_limb_space( bp_marker, bp_nlimbs );
787   if (ep_marker)
788     _gcry_mpi_free_limb_space( ep_marker, ep_nlimbs );
789   if (xp_marker)
790     _gcry_mpi_free_limb_space( xp_marker, xp_nlimbs );
791 }
792 #endif