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