1/* mpi-mpow.c - MPI functions 2 * Copyright (C) 1998, 1999, 2001, 2002, 2003 Free Software Foundation, Inc. 3 * 4 * This file is part of Libgcrypt. 5 * 6 * Libgcrypt is free software; you can redistribute it and/or modify 7 * it under the terms of the GNU Lesser General Public License as 8 * published by the Free Software Foundation; either version 2.1 of 9 * the License, or (at your option) any later version. 10 * 11 * Libgcrypt 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 Lesser General Public License for more details. 15 * 16 * You should have received a copy of the GNU Lesser General Public 17 * License 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 25#include "mpi-internal.h" 26#include "longlong.h" 27#include "g10lib.h" 28 29 30/* Barrett is slower than the classical way. It can be tweaked by 31 * using partial multiplications 32 */ 33/*#define USE_BARRETT*/ 34 35 36 37#ifdef USE_BARRETT 38static void barrett_mulm( gcry_mpi_t w, gcry_mpi_t u, gcry_mpi_t v, gcry_mpi_t m, gcry_mpi_t y, int k, gcry_mpi_t r1, gcry_mpi_t r2 ); 39static gcry_mpi_t init_barrett( gcry_mpi_t m, int *k, gcry_mpi_t *r1, gcry_mpi_t *r2 ); 40static int calc_barrett( gcry_mpi_t r, gcry_mpi_t x, gcry_mpi_t m, gcry_mpi_t y, int k, gcry_mpi_t r1, gcry_mpi_t r2 ); 41#else 42#define barrett_mulm( w, u, v, m, y, k, r1, r2 ) gcry_mpi_mulm( (w), (u), (v), (m) ) 43#endif 44 45 46static int 47build_index( gcry_mpi_t *exparray, int k, int i, int t ) 48{ 49 int j, bitno; 50 int idx = 0; 51 52 bitno = t-i; 53 for(j=k-1; j >= 0; j-- ) { 54 idx <<= 1; 55 if( mpi_test_bit( exparray[j], bitno ) ) 56 idx |= 1; 57 } 58 /*log_debug("t=%d i=%d idx=%d\n", t, i, idx );*/ 59 return idx; 60} 61 62/**************** 63 * RES = (BASE[0] ^ EXP[0]) * (BASE[1] ^ EXP[1]) * ... * mod M 64 */ 65void 66_gcry_mpi_mulpowm( gcry_mpi_t res, gcry_mpi_t *basearray, gcry_mpi_t *exparray, gcry_mpi_t m) 67{ 68 int k; /* number of elements */ 69 int t; /* bit size of largest exponent */ 70 int i, j, idx; 71 gcry_mpi_t *G; /* table with precomputed values of size 2^k */ 72 gcry_mpi_t tmp; 73#ifdef USE_BARRETT 74 gcry_mpi_t barrett_y, barrett_r1, barrett_r2; 75 int barrett_k; 76#endif 77 78 for(k=0; basearray[k]; k++ ) 79 ; 80 gcry_assert(k); 81 for(t=0, i=0; (tmp=exparray[i]); i++ ) { 82 /*log_mpidump("exp: ", tmp );*/ 83 j = mpi_get_nbits(tmp); 84 if( j > t ) 85 t = j; 86 } 87 /*log_mpidump("mod: ", m );*/ 88 gcry_assert (i==k); 89 gcry_assert (t); 90 gcry_assert (k < 10); 91 92 G = gcry_xcalloc( (1<<k) , sizeof *G ); 93#ifdef USE_BARRETT 94 barrett_y = init_barrett( m, &barrett_k, &barrett_r1, &barrett_r2 ); 95#endif 96 /* and calculate */ 97 tmp = mpi_alloc( mpi_get_nlimbs(m)+1 ); 98 mpi_set_ui( res, 1 ); 99 for(i = 1; i <= t; i++ ) { 100 barrett_mulm(tmp, res, res, m, barrett_y, barrett_k, 101 barrett_r1, barrett_r2 ); 102 idx = build_index( exparray, k, i, t ); 103 gcry_assert (idx >= 0 && idx < (1<<k)); 104 if( !G[idx] ) { 105 if( !idx ) 106 G[0] = mpi_alloc_set_ui( 1 ); 107 else { 108 for(j=0; j < k; j++ ) { 109 if( (idx & (1<<j) ) ) { 110 if( !G[idx] ) 111 G[idx] = mpi_copy( basearray[j] ); 112 else 113 barrett_mulm( G[idx], G[idx], basearray[j], 114 m, barrett_y, barrett_k, barrett_r1, barrett_r2 ); 115 } 116 } 117 if( !G[idx] ) 118 G[idx] = mpi_alloc(0); 119 } 120 } 121 barrett_mulm(res, tmp, G[idx], m, barrett_y, barrett_k, barrett_r1, barrett_r2 ); 122 } 123 124 /* cleanup */ 125 mpi_free(tmp); 126#ifdef USE_BARRETT 127 mpi_free(barrett_y); 128 mpi_free(barrett_r1); 129 mpi_free(barrett_r2); 130#endif 131 for(i=0; i < (1<<k); i++ ) 132 mpi_free(G[i]); 133 gcry_free(G); 134} 135 136 137 138#ifdef USE_BARRETT 139static void 140barrett_mulm( gcry_mpi_t w, gcry_mpi_t u, gcry_mpi_t v, gcry_mpi_t m, gcry_mpi_t y, int k, gcry_mpi_t r1, gcry_mpi_t r2 ) 141{ 142 mpi_mul(w, u, v); 143 if( calc_barrett( w, w, m, y, k, r1, r2 ) ) 144 mpi_fdiv_r( w, w, m ); 145} 146 147/**************** 148 * Barrett precalculation: y = floor(b^(2k) / m) 149 */ 150static gcry_mpi_t 151init_barrett( gcry_mpi_t m, int *k, gcry_mpi_t *r1, gcry_mpi_t *r2 ) 152{ 153 gcry_mpi_t tmp; 154 155 mpi_normalize( m ); 156 *k = mpi_get_nlimbs( m ); 157 tmp = mpi_alloc( *k + 1 ); 158 mpi_set_ui( tmp, 1 ); 159 mpi_lshift_limbs( tmp, 2 * *k ); 160 mpi_fdiv_q( tmp, tmp, m ); 161 *r1 = mpi_alloc( 2* *k + 1 ); 162 *r2 = mpi_alloc( 2* *k + 1 ); 163 return tmp; 164} 165 166/**************** 167 * Barrett reduction: We assume that these conditions are met: 168 * Given x =(x_2k-1 ...x_0)_b 169 * m =(m_k-1 ....m_0)_b with m_k-1 != 0 170 * Output r = x mod m 171 * Before using this function init_barret must be used to calucalte y and k. 172 * Returns: false = no error 173 * true = can't perform barret reduction 174 */ 175static int 176calc_barrett( gcry_mpi_t r, gcry_mpi_t x, gcry_mpi_t m, gcry_mpi_t y, int k, gcry_mpi_t r1, gcry_mpi_t r2 ) 177{ 178 int xx = k > 3 ? k-3:0; 179 180 mpi_normalize( x ); 181 if( mpi_get_nlimbs(x) > 2*k ) 182 return 1; /* can't do it */ 183 184 /* 1. q1 = floor( x / b^k-1) 185 * q2 = q1 * y 186 * q3 = floor( q2 / b^k+1 ) 187 * Actually, we don't need qx, we can work direct on r2 188 */ 189 mpi_set( r2, x ); 190 mpi_rshift_limbs( r2, k-1 ); 191 mpi_mul( r2, r2, y ); 192 mpi_rshift_limbs( r2, k+1 ); 193 194 /* 2. r1 = x mod b^k+1 195 * r2 = q3 * m mod b^k+1 196 * r = r1 - r2 197 * 3. if r < 0 then r = r + b^k+1 198 */ 199 mpi_set( r1, x ); 200 if( r1->nlimbs > k+1 ) /* quick modulo operation */ 201 r1->nlimbs = k+1; 202 mpi_mul( r2, r2, m ); 203 if( r2->nlimbs > k+1 ) /* quick modulo operation */ 204 r2->nlimbs = k+1; 205 mpi_sub( r, r1, r2 ); 206 207 if( mpi_is_neg( r ) ) { 208 gcry_mpi_t tmp; 209 210 tmp = mpi_alloc( k + 2 ); 211 mpi_set_ui( tmp, 1 ); 212 mpi_lshift_limbs( tmp, k+1 ); 213 mpi_add( r, r, tmp ); 214 mpi_free(tmp); 215 } 216 217 /* 4. while r >= m do r = r - m */ 218 while( mpi_cmp( r, m ) >= 0 ) 219 mpi_sub( r, r, m ); 220 221 return 0; 222} 223#endif /* USE_BARRETT */ 224