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

          Line data    Source code
       1             : #include <tommath.h>
       2             : #ifdef BN_FAST_MP_MONTGOMERY_REDUCE_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 xR**-1 == x (mod N) via Montgomery Reduction
      19             :  *
      20             :  * This is an optimized implementation of montgomery_reduce
      21             :  * which uses the comba method to quickly calculate the columns of the
      22             :  * reduction.
      23             :  *
      24             :  * Based on Algorithm 14.32 on pp.601 of HAC.
      25             : */
      26      726421 : int fast_mp_montgomery_reduce (mp_int * x, mp_int * n, mp_digit rho)
      27             : {
      28             :   int     ix, res, olduse;
      29             :   mp_word W[MP_WARRAY];
      30             : 
      31             :   /* get old used count */
      32      726421 :   olduse = x->used;
      33             : 
      34             :   /* grow a as required */
      35      726421 :   if (x->alloc < n->used + 1) {
      36           0 :     if ((res = mp_grow (x, n->used + 1)) != MP_OKAY) {
      37           0 :       return res;
      38             :     }
      39             :   }
      40             : 
      41             :   /* first we have to get the digits of the input into
      42             :    * an array of double precision words W[...]
      43             :    */
      44             :   {
      45             :     register mp_word *_W;
      46             :     register mp_digit *tmpx;
      47             : 
      48             :     /* alias for the W[] array */
      49      726421 :     _W   = W;
      50             : 
      51             :     /* alias for the digits of  x*/
      52      726421 :     tmpx = x->dp;
      53             : 
      54             :     /* copy the digits of a into W[0..a->used-1] */
      55    40251773 :     for (ix = 0; ix < x->used; ix++) {
      56    39525352 :       *_W++ = *tmpx++;
      57             :     }
      58             : 
      59             :     /* zero the high words of W[a->used..m->used*2] */
      60     2206026 :     for (; ix < n->used * 2 + 1; ix++) {
      61     1517985 :       *_W++ = 0;
      62             :     }
      63             :   }
      64             : 
      65             :   /* now we proceed to zero successive digits
      66             :    * from the least significant upwards
      67             :    */
      68    20846499 :   for (ix = 0; ix < n->used; ix++) {
      69             :     /* mu = ai * m' mod b
      70             :      *
      71             :      * We avoid a double precision multiplication (which isn't required)
      72             :      * by casting the value down to a mp_digit.  Note this requires
      73             :      * that W[ix-1] have  the carry cleared (see after the inner loop)
      74             :      */
      75             :     register mp_digit mu;
      76    20158458 :     mu = (mp_digit) (((W[ix] & MP_MASK) * rho) & MP_MASK);
      77             : 
      78             :     /* a = a + mu * m * b**i
      79             :      *
      80             :      * This is computed in place and on the fly.  The multiplication
      81             :      * by b**i is handled by offseting which columns the results
      82             :      * are added to.
      83             :      *
      84             :      * Note the comba method normally doesn't handle carries in the
      85             :      * inner loop In this case we fix the carry from the previous
      86             :      * column since the Montgomery reduction requires digits of the
      87             :      * result (so far) [see above] to work.  This is
      88             :      * handled by fixing up one carry after the inner loop.  The
      89             :      * carry fixups are done in order so after these loops the
      90             :      * first m->used words of W[] have the carries fixed
      91             :      */
      92             :     {
      93             :       register int iy;
      94             :       register mp_digit *tmpn;
      95             :       register mp_word *_W;
      96             : 
      97             :       /* alias for the digits of the modulus */
      98    20158458 :       tmpn = n->dp;
      99             : 
     100             :       /* Alias for the columns set by an offset of ix */
     101    20158458 :       _W = W + ix;
     102             : 
     103             :       /* inner loop */
     104   674261502 :       for (iy = 0; iy < n->used; iy++) {
     105   654103044 :           *_W++ += ((mp_word)mu) * ((mp_word)*tmpn++);
     106             :       }
     107             :     }
     108             : 
     109             :     /* now fix carry for next digit, W[ix+1] */
     110    20158458 :     W[ix + 1] += W[ix] >> ((mp_word) DIGIT_BIT);
     111             :   }
     112             : 
     113             :   /* now we have to propagate the carries and
     114             :    * shift the words downward [all those least
     115             :    * significant digits we zeroed].
     116             :    */
     117             :   {
     118             :     register mp_digit *tmpx;
     119             :     register mp_word *_W, *_W1;
     120             : 
     121             :     /* nox fix rest of carries */
     122             : 
     123             :     /* alias for current word */
     124      726421 :     _W1 = W + ix;
     125             : 
     126             :     /* alias for next word, where the carry goes */
     127      726421 :     _W = W + ++ix;
     128             : 
     129    21611300 :     for (; ix <= n->used * 2 + 1; ix++) {
     130    20884879 :       *_W++ += *_W1++ >> ((mp_word) DIGIT_BIT);
     131             :     }
     132             : 
     133             :     /* copy out, A = A/b**n
     134             :      *
     135             :      * The result is A/b**n but instead of converting from an
     136             :      * array of mp_word to mp_digit than calling mp_rshd
     137             :      * we just copy them in the right order
     138             :      */
     139             : 
     140             :     /* alias for destination word */
     141      726421 :     tmpx = x->dp;
     142             : 
     143             :     /* alias for shifted double precision result */
     144      726421 :     _W = W + n->used;
     145             : 
     146    21611300 :     for (ix = 0; ix < n->used + 1; ix++) {
     147    20884879 :       *tmpx++ = (mp_digit)(*_W++ & ((mp_word) MP_MASK));
     148             :     }
     149             : 
     150             :     /* zero oldused digits, if the input a was larger than
     151             :      * m->used+1 we'll have to clear the digits
     152             :      */
     153    19329398 :     for (; ix < olduse; ix++) {
     154    18641357 :       *tmpx++ = 0;
     155             :     }
     156             :   }
     157             : 
     158             :   /* set the max used and clamp */
     159      726421 :   x->used = n->used + 1;
     160      726421 :   mp_clamp (x);
     161             : 
     162             :   /* if A >= m then A = A - m */
     163      726421 :   if (mp_cmp_mag (x, n) != MP_LT) {
     164           0 :     return s_mp_sub (x, n, x);
     165             :   }
     166      688041 :   return MP_OKAY;
     167             : }
     168             : #endif
     169             : 
     170             : /* $Source: /cvs/libtom/libtommath/bn_fast_mp_montgomery_reduce.c,v $ */
     171             : /* $Revision: 1.4 $ */
     172             : /* $Date: 2006/12/28 01:25:13 $ */

Generated by: LCOV version 1.13