LCOV - code coverage report
Current view: top level - source4/heimdal/lib/hcrypto/libtommath - bn_mp_exptmod_fast.c (source / functions) Hit Total Coverage
Test: coverage report for abartlet/fix-coverage dd10fb34 Lines: 78 117 66.7 %
Date: 2021-09-23 10:06:22 Functions: 1 1 100.0 %

          Line data    Source code
       1             : #include <tommath.h>
       2             : #ifdef BN_MP_EXPTMOD_FAST_C
       3             : /* LibTomMath, multiple-precision integer library -- Tom St Denis
       4             :  *
       5             :  * LibTomMath is a library that provides multiple-precision
       6             :  * integer arithmetic as well as number theoretic functionality.
       7             :  *
       8             :  * The library was designed directly after the MPI library by
       9             :  * Michael Fromberger but has been written from scratch with
      10             :  * additional optimizations in place.
      11             :  *
      12             :  * The library is free for all purposes without any express
      13             :  * guarantee it works.
      14             :  *
      15             :  * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
      16             :  */
      17             : 
      18             : /* computes Y == G**X mod P, HAC pp.616, Algorithm 14.85
      19             :  *
      20             :  * Uses a left-to-right k-ary sliding window to compute the modular exponentiation.
      21             :  * The value of k changes based on the size of the exponent.
      22             :  *
      23             :  * Uses Montgomery or Diminished Radix reduction [whichever appropriate]
      24             :  */
      25             : 
      26             : #ifdef MP_LOW_MEM
      27             :    #define TAB_SIZE 32
      28             : #else
      29             :    #define TAB_SIZE 256
      30             : #endif
      31             : 
      32         872 : int mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode)
      33             : {
      34             :   mp_int  M[TAB_SIZE], res;
      35             :   mp_digit buf, mp;
      36             :   int     err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
      37             : 
      38             :   /* use a pointer to the reduction algorithm.  This allows us to use
      39             :    * one of many reduction algorithms without modding the guts of
      40             :    * the code with if statements everywhere.
      41             :    */
      42             :   int     (*redux)(mp_int*,mp_int*,mp_digit);
      43             : 
      44             :   /* find window size */
      45         872 :   x = mp_count_bits (X);
      46         872 :   if (x <= 7) {
      47           0 :     winsize = 2;
      48         872 :   } else if (x <= 36) {
      49         418 :     winsize = 3;
      50         438 :   } else if (x <= 140) {
      51           0 :     winsize = 4;
      52         438 :   } else if (x <= 450) {
      53           0 :     winsize = 5;
      54         438 :   } else if (x <= 1303) {
      55         278 :     winsize = 6;
      56         160 :   } else if (x <= 3529) {
      57         144 :     winsize = 7;
      58             :   } else {
      59           0 :     winsize = 8;
      60             :   }
      61             : 
      62             : #ifdef MP_LOW_MEM
      63             :   if (winsize > 5) {
      64             :      winsize = 5;
      65             :   }
      66             : #endif
      67             : 
      68             :   /* init M array */
      69             :   /* init first cell */
      70         872 :   if ((err = mp_init(&M[1])) != MP_OKAY) {
      71           0 :      return err;
      72             :   }
      73             : 
      74             :   /* now init the second half of the array */
      75       21744 :   for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
      76       20872 :     if ((err = mp_init(&M[x])) != MP_OKAY) {
      77           0 :       for (y = 1<<(winsize-1); y < x; y++) {
      78           0 :         mp_clear (&M[y]);
      79             :       }
      80           0 :       mp_clear(&M[1]);
      81           0 :       return err;
      82             :     }
      83             :   }
      84             : 
      85             :   /* determine and setup reduction code */
      86         872 :   if (redmode == 0) {
      87             : #ifdef BN_MP_MONTGOMERY_SETUP_C
      88             :      /* now setup montgomery  */
      89         872 :      if ((err = mp_montgomery_setup (P, &mp)) != MP_OKAY) {
      90           0 :         goto LBL_M;
      91             :      }
      92             : #else
      93             :      err = MP_VAL;
      94             :      goto LBL_M;
      95             : #endif
      96             : 
      97             :      /* automatically pick the comba one if available (saves quite a few calls/ifs) */
      98             : #ifdef BN_FAST_MP_MONTGOMERY_REDUCE_C
      99        1712 :      if (((P->used * 2 + 1) < MP_WARRAY) &&
     100         840 :           P->used < (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
     101         840 :         redux = fast_mp_montgomery_reduce;
     102             :      } else
     103             : #endif
     104             :      {
     105             : #ifdef BN_MP_MONTGOMERY_REDUCE_C
     106             :         /* use slower baseline Montgomery method */
     107           0 :         redux = mp_montgomery_reduce;
     108             : #else
     109             :         err = MP_VAL;
     110             :         goto LBL_M;
     111             : #endif
     112             :      }
     113           0 :   } else if (redmode == 1) {
     114             : #if defined(BN_MP_DR_SETUP_C) && defined(BN_MP_DR_REDUCE_C)
     115             :      /* setup DR reduction for moduli of the form B**k - b */
     116           0 :      mp_dr_setup(P, &mp);
     117           0 :      redux = mp_dr_reduce;
     118             : #else
     119             :      err = MP_VAL;
     120             :      goto LBL_M;
     121             : #endif
     122             :   } else {
     123             : #if defined(BN_MP_REDUCE_2K_SETUP_C) && defined(BN_MP_REDUCE_2K_C)
     124             :      /* setup DR reduction for moduli of the form 2**k - b */
     125           0 :      if ((err = mp_reduce_2k_setup(P, &mp)) != MP_OKAY) {
     126           0 :         goto LBL_M;
     127             :      }
     128           0 :      redux = mp_reduce_2k;
     129             : #else
     130             :      err = MP_VAL;
     131             :      goto LBL_M;
     132             : #endif
     133             :   }
     134             : 
     135             :   /* setup result */
     136         872 :   if ((err = mp_init (&res)) != MP_OKAY) {
     137           0 :     goto LBL_M;
     138             :   }
     139             : 
     140             :   /* create M table
     141             :    *
     142             : 
     143             :    *
     144             :    * The first half of the table is not computed though accept for M[0] and M[1]
     145             :    */
     146             : 
     147         872 :   if (redmode == 0) {
     148             : #ifdef BN_MP_MONTGOMERY_CALC_NORMALIZATION_C
     149             :      /* now we need R mod m */
     150         872 :      if ((err = mp_montgomery_calc_normalization (&res, P)) != MP_OKAY) {
     151           0 :        goto LBL_RES;
     152             :      }
     153             : #else
     154             :      err = MP_VAL;
     155             :      goto LBL_RES;
     156             : #endif
     157             : 
     158             :      /* now set M[1] to G * R mod m */
     159         872 :      if ((err = mp_mulmod (G, &res, P, &M[1])) != MP_OKAY) {
     160           0 :        goto LBL_RES;
     161             :      }
     162             :   } else {
     163           0 :      mp_set(&res, 1);
     164           0 :      if ((err = mp_mod(G, P, &M[1])) != MP_OKAY) {
     165           0 :         goto LBL_RES;
     166             :      }
     167             :   }
     168             : 
     169             :   /* compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times */
     170         872 :   if ((err = mp_copy (&M[1], &M[1 << (winsize - 1)])) != MP_OKAY) {
     171           0 :     goto LBL_RES;
     172             :   }
     173             : 
     174        4058 :   for (x = 0; x < (winsize - 1); x++) {
     175        3218 :     if ((err = mp_sqr (&M[1 << (winsize - 1)], &M[1 << (winsize - 1)])) != MP_OKAY) {
     176           0 :       goto LBL_RES;
     177             :     }
     178        3218 :     if ((err = redux (&M[1 << (winsize - 1)], P, mp)) != MP_OKAY) {
     179           0 :       goto LBL_RES;
     180             :     }
     181             :   }
     182             : 
     183             :   /* create upper table */
     184       20872 :   for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
     185       20000 :     if ((err = mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) {
     186           0 :       goto LBL_RES;
     187             :     }
     188       20000 :     if ((err = redux (&M[x], P, mp)) != MP_OKAY) {
     189           0 :       goto LBL_RES;
     190             :     }
     191             :   }
     192             : 
     193             :   /* set initial mode and bit cnt */
     194         872 :   mode   = 0;
     195         872 :   bitcnt = 1;
     196         872 :   buf    = 0;
     197         872 :   digidx = X->used - 1;
     198         872 :   bitcpy = 0;
     199         872 :   bitbuf = 0;
     200             : 
     201             :   for (;;) {
     202             :     /* grab next digit as required */
     203     1290872 :     if (--bitcnt == 0) {
     204             :       /* if digidx == -1 we are out of digits so break */
     205       11910 :       if (digidx == -1) {
     206         840 :         break;
     207             :       }
     208             :       /* read next digit and reset bitcnt */
     209       11038 :       buf    = X->dp[digidx--];
     210       11038 :       bitcnt = (int)DIGIT_BIT;
     211             :     }
     212             : 
     213             :     /* grab the next msb from the exponent */
     214      662280 :     y     = (mp_digit)(buf >> (DIGIT_BIT - 1)) & 1;
     215      662280 :     buf <<= (mp_digit)1;
     216             : 
     217             :     /* if the bit is zero and mode == 0 then we ignore it
     218             :      * These represent the leading zero bits before the first 1 bit
     219             :      * in the exponent.  Technically this opt is not required but it
     220             :      * does lower the # of trivial squaring/reductions used
     221             :      */
     222      662280 :     if (mode == 0 && y == 0) {
     223       43124 :       continue;
     224             :     }
     225             : 
     226             :     /* if the bit is zero and mode == 1 then we square */
     227      619156 :     if (mode == 1 && y == 0) {
     228       86459 :       if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
     229           0 :         goto LBL_RES;
     230             :       }
     231       86459 :       if ((err = redux (&res, P, mp)) != MP_OKAY) {
     232           0 :         goto LBL_RES;
     233             :       }
     234       86459 :       continue;
     235             :     }
     236             : 
     237             :     /* else we add it to the window */
     238      532697 :     bitbuf |= (y << (winsize - ++bitcpy));
     239      532697 :     mode    = 2;
     240             : 
     241      532697 :     if (bitcpy == winsize) {
     242             :       /* ok window is filled so square as required and multiply  */
     243             :       /* square first */
     244      608682 :       for (x = 0; x < winsize; x++) {
     245      530901 :         if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
     246           0 :           goto LBL_RES;
     247             :         }
     248      530901 :         if ((err = redux (&res, P, mp)) != MP_OKAY) {
     249           0 :           goto LBL_RES;
     250             :         }
     251             :       }
     252             : 
     253             :       /* then multiply */
     254       81897 :       if ((err = mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) {
     255           0 :         goto LBL_RES;
     256             :       }
     257       81897 :       if ((err = redux (&res, P, mp)) != MP_OKAY) {
     258           0 :         goto LBL_RES;
     259             :       }
     260             : 
     261             :       /* empty window and reset */
     262       77781 :       bitcpy = 0;
     263       77781 :       bitbuf = 0;
     264       77781 :       mode   = 1;
     265             :     }
     266             :   }
     267             : 
     268             :   /* if bits remain then square/multiply */
     269         872 :   if (mode == 2 && bitcpy > 0) {
     270             :     /* square then multiply if the bit is set */
     271        2604 :     for (x = 0; x < bitcpy; x++) {
     272        1796 :       if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
     273           0 :         goto LBL_RES;
     274             :       }
     275        1796 :       if ((err = redux (&res, P, mp)) != MP_OKAY) {
     276           0 :         goto LBL_RES;
     277             :       }
     278             : 
     279             :       /* get next bit of the window */
     280        1796 :       bitbuf <<= 1;
     281        1796 :       if ((bitbuf & (1 << winsize)) != 0) {
     282             :         /* then multiply */
     283        1278 :         if ((err = mp_mul (&res, &M[1], &res)) != MP_OKAY) {
     284           0 :           goto LBL_RES;
     285             :         }
     286        1278 :         if ((err = redux (&res, P, mp)) != MP_OKAY) {
     287           0 :           goto LBL_RES;
     288             :         }
     289             :       }
     290             :     }
     291             :   }
     292             : 
     293         872 :   if (redmode == 0) {
     294             :      /* fixup result if Montgomery reduction is used
     295             :       * recall that any value in a Montgomery system is
     296             :       * actually multiplied by R mod n.  So we have
     297             :       * to reduce one more time to cancel out the factor
     298             :       * of R.
     299             :       */
     300         872 :      if ((err = redux(&res, P, mp)) != MP_OKAY) {
     301           0 :        goto LBL_RES;
     302             :      }
     303             :   }
     304             : 
     305             :   /* swap res with Y */
     306         872 :   mp_exch (&res, Y);
     307         872 :   err = MP_OKAY;
     308         872 : LBL_RES:mp_clear (&res);
     309         872 : LBL_M:
     310         872 :   mp_clear(&M[1]);
     311       21744 :   for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
     312       20872 :     mp_clear (&M[x]);
     313             :   }
     314         840 :   return err;
     315             : }
     316             : #endif
     317             : 
     318             : 
     319             : /* $Source: /cvs/libtom/libtommath/bn_mp_exptmod_fast.c,v $ */
     320             : /* $Revision: 1.4 $ */
     321             : /* $Date: 2006/12/28 01:25:13 $ */

Generated by: LCOV version 1.13