PolarSSL v1.2.9
bignum.c
Go to the documentation of this file.
1 /*
2  * Multi-precision integer library
3  *
4  * Copyright (C) 2006-2010, Brainspark B.V.
5  *
6  * This file is part of PolarSSL (http://www.polarssl.org)
7  * Lead Maintainer: Paul Bakker <polarssl_maintainer at polarssl.org>
8  *
9  * All rights reserved.
10  *
11  * This program is free software; you can redistribute it and/or modify
12  * it under the terms of the GNU General Public License as published by
13  * the Free Software Foundation; either version 2 of the License, or
14  * (at your option) any later version.
15  *
16  * This program is distributed in the hope that it will be useful,
17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19  * GNU General Public License for more details.
20  *
21  * You should have received a copy of the GNU General Public License along
22  * with this program; if not, write to the Free Software Foundation, Inc.,
23  * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24  */
25 /*
26  * This MPI implementation is based on:
27  *
28  * http://www.cacr.math.uwaterloo.ca/hac/about/chap14.pdf
29  * http://www.stillhq.com/extracted/gnupg-api/mpi/
30  * http://math.libtomcrypt.com/files/tommath.pdf
31  */
32 
33 #include "polarssl/config.h"
34 
35 #if defined(POLARSSL_BIGNUM_C)
36 
37 #include "polarssl/bignum.h"
38 #include "polarssl/bn_mul.h"
39 
40 #include <stdlib.h>
41 
42 #define ciL (sizeof(t_uint)) /* chars in limb */
43 #define biL (ciL << 3) /* bits in limb */
44 #define biH (ciL << 2) /* half limb size */
45 
46 /*
47  * Convert between bits/chars and number of limbs
48  */
49 #define BITS_TO_LIMBS(i) (((i) + biL - 1) / biL)
50 #define CHARS_TO_LIMBS(i) (((i) + ciL - 1) / ciL)
51 
52 /*
53  * Initialize one MPI
54  */
55 void mpi_init( mpi *X )
56 {
57  if( X == NULL )
58  return;
59 
60  X->s = 1;
61  X->n = 0;
62  X->p = NULL;
63 }
64 
65 /*
66  * Unallocate one MPI
67  */
68 void mpi_free( mpi *X )
69 {
70  if( X == NULL )
71  return;
72 
73  if( X->p != NULL )
74  {
75  memset( X->p, 0, X->n * ciL );
76  free( X->p );
77  }
78 
79  X->s = 1;
80  X->n = 0;
81  X->p = NULL;
82 }
83 
84 /*
85  * Enlarge to the specified number of limbs
86  */
87 int mpi_grow( mpi *X, size_t nblimbs )
88 {
89  t_uint *p;
90 
91  if( nblimbs > POLARSSL_MPI_MAX_LIMBS )
93 
94  if( X->n < nblimbs )
95  {
96  if( ( p = (t_uint *) malloc( nblimbs * ciL ) ) == NULL )
98 
99  memset( p, 0, nblimbs * ciL );
100 
101  if( X->p != NULL )
102  {
103  memcpy( p, X->p, X->n * ciL );
104  memset( X->p, 0, X->n * ciL );
105  free( X->p );
106  }
107 
108  X->n = nblimbs;
109  X->p = p;
110  }
111 
112  return( 0 );
113 }
114 
115 /*
116  * Copy the contents of Y into X
117  */
118 int mpi_copy( mpi *X, const mpi *Y )
119 {
120  int ret;
121  size_t i;
122 
123  if( X == Y )
124  return( 0 );
125 
126  for( i = Y->n - 1; i > 0; i-- )
127  if( Y->p[i] != 0 )
128  break;
129  i++;
130 
131  X->s = Y->s;
132 
133  MPI_CHK( mpi_grow( X, i ) );
134 
135  memset( X->p, 0, X->n * ciL );
136  memcpy( X->p, Y->p, i * ciL );
137 
138 cleanup:
139 
140  return( ret );
141 }
142 
143 /*
144  * Swap the contents of X and Y
145  */
146 void mpi_swap( mpi *X, mpi *Y )
147 {
148  mpi T;
149 
150  memcpy( &T, X, sizeof( mpi ) );
151  memcpy( X, Y, sizeof( mpi ) );
152  memcpy( Y, &T, sizeof( mpi ) );
153 }
154 
155 /*
156  * Set value from integer
157  */
158 int mpi_lset( mpi *X, t_sint z )
159 {
160  int ret;
161 
162  MPI_CHK( mpi_grow( X, 1 ) );
163  memset( X->p, 0, X->n * ciL );
164 
165  X->p[0] = ( z < 0 ) ? -z : z;
166  X->s = ( z < 0 ) ? -1 : 1;
167 
168 cleanup:
169 
170  return( ret );
171 }
172 
173 /*
174  * Get a specific bit
175  */
176 int mpi_get_bit( const mpi *X, size_t pos )
177 {
178  if( X->n * biL <= pos )
179  return( 0 );
180 
181  return ( X->p[pos / biL] >> ( pos % biL ) ) & 0x01;
182 }
183 
184 /*
185  * Set a bit to a specific value of 0 or 1
186  */
187 int mpi_set_bit( mpi *X, size_t pos, unsigned char val )
188 {
189  int ret = 0;
190  size_t off = pos / biL;
191  size_t idx = pos % biL;
192 
193  if( val != 0 && val != 1 )
195 
196  if( X->n * biL <= pos )
197  {
198  if( val == 0 )
199  return ( 0 );
200 
201  MPI_CHK( mpi_grow( X, off + 1 ) );
202  }
203 
204  X->p[off] = ( X->p[off] & ~( 0x01 << idx ) ) | ( val << idx );
205 
206 cleanup:
207 
208  return( ret );
209 }
210 
211 /*
212  * Return the number of least significant bits
213  */
214 size_t mpi_lsb( const mpi *X )
215 {
216  size_t i, j, count = 0;
217 
218  for( i = 0; i < X->n; i++ )
219  for( j = 0; j < biL; j++, count++ )
220  if( ( ( X->p[i] >> j ) & 1 ) != 0 )
221  return( count );
222 
223  return( 0 );
224 }
225 
226 /*
227  * Return the number of most significant bits
228  */
229 size_t mpi_msb( const mpi *X )
230 {
231  size_t i, j;
232 
233  for( i = X->n - 1; i > 0; i-- )
234  if( X->p[i] != 0 )
235  break;
236 
237  for( j = biL; j > 0; j-- )
238  if( ( ( X->p[i] >> ( j - 1 ) ) & 1 ) != 0 )
239  break;
240 
241  return( ( i * biL ) + j );
242 }
243 
244 /*
245  * Return the total size in bytes
246  */
247 size_t mpi_size( const mpi *X )
248 {
249  return( ( mpi_msb( X ) + 7 ) >> 3 );
250 }
251 
252 /*
253  * Convert an ASCII character to digit value
254  */
255 static int mpi_get_digit( t_uint *d, int radix, char c )
256 {
257  *d = 255;
258 
259  if( c >= 0x30 && c <= 0x39 ) *d = c - 0x30;
260  if( c >= 0x41 && c <= 0x46 ) *d = c - 0x37;
261  if( c >= 0x61 && c <= 0x66 ) *d = c - 0x57;
262 
263  if( *d >= (t_uint) radix )
265 
266  return( 0 );
267 }
268 
269 /*
270  * Import from an ASCII string
271  */
272 int mpi_read_string( mpi *X, int radix, const char *s )
273 {
274  int ret;
275  size_t i, j, slen, n;
276  t_uint d;
277  mpi T;
278 
279  if( radix < 2 || radix > 16 )
281 
282  mpi_init( &T );
283 
284  slen = strlen( s );
285 
286  if( radix == 16 )
287  {
288  n = BITS_TO_LIMBS( slen << 2 );
289 
290  MPI_CHK( mpi_grow( X, n ) );
291  MPI_CHK( mpi_lset( X, 0 ) );
292 
293  for( i = slen, j = 0; i > 0; i--, j++ )
294  {
295  if( i == 1 && s[i - 1] == '-' )
296  {
297  X->s = -1;
298  break;
299  }
300 
301  MPI_CHK( mpi_get_digit( &d, radix, s[i - 1] ) );
302  X->p[j / (2 * ciL)] |= d << ( (j % (2 * ciL)) << 2 );
303  }
304  }
305  else
306  {
307  MPI_CHK( mpi_lset( X, 0 ) );
308 
309  for( i = 0; i < slen; i++ )
310  {
311  if( i == 0 && s[i] == '-' )
312  {
313  X->s = -1;
314  continue;
315  }
316 
317  MPI_CHK( mpi_get_digit( &d, radix, s[i] ) );
318  MPI_CHK( mpi_mul_int( &T, X, radix ) );
319 
320  if( X->s == 1 )
321  {
322  MPI_CHK( mpi_add_int( X, &T, d ) );
323  }
324  else
325  {
326  MPI_CHK( mpi_sub_int( X, &T, d ) );
327  }
328  }
329  }
330 
331 cleanup:
332 
333  mpi_free( &T );
334 
335  return( ret );
336 }
337 
338 /*
339  * Helper to write the digits high-order first
340  */
341 static int mpi_write_hlp( mpi *X, int radix, char **p )
342 {
343  int ret;
344  t_uint r;
345 
346  if( radix < 2 || radix > 16 )
348 
349  MPI_CHK( mpi_mod_int( &r, X, radix ) );
350  MPI_CHK( mpi_div_int( X, NULL, X, radix ) );
351 
352  if( mpi_cmp_int( X, 0 ) != 0 )
353  MPI_CHK( mpi_write_hlp( X, radix, p ) );
354 
355  if( r < 10 )
356  *(*p)++ = (char)( r + 0x30 );
357  else
358  *(*p)++ = (char)( r + 0x37 );
359 
360 cleanup:
361 
362  return( ret );
363 }
364 
365 /*
366  * Export into an ASCII string
367  */
368 int mpi_write_string( const mpi *X, int radix, char *s, size_t *slen )
369 {
370  int ret = 0;
371  size_t n;
372  char *p;
373  mpi T;
374 
375  if( radix < 2 || radix > 16 )
377 
378  n = mpi_msb( X );
379  if( radix >= 4 ) n >>= 1;
380  if( radix >= 16 ) n >>= 1;
381  n += 3;
382 
383  if( *slen < n )
384  {
385  *slen = n;
387  }
388 
389  p = s;
390  mpi_init( &T );
391 
392  if( X->s == -1 )
393  *p++ = '-';
394 
395  if( radix == 16 )
396  {
397  int c;
398  size_t i, j, k;
399 
400  for( i = X->n, k = 0; i > 0; i-- )
401  {
402  for( j = ciL; j > 0; j-- )
403  {
404  c = ( X->p[i - 1] >> ( ( j - 1 ) << 3) ) & 0xFF;
405 
406  if( c == 0 && k == 0 && ( i + j + 3 ) != 0 )
407  continue;
408 
409  *(p++) = "0123456789ABCDEF" [c / 16];
410  *(p++) = "0123456789ABCDEF" [c % 16];
411  k = 1;
412  }
413  }
414  }
415  else
416  {
417  MPI_CHK( mpi_copy( &T, X ) );
418 
419  if( T.s == -1 )
420  T.s = 1;
421 
422  MPI_CHK( mpi_write_hlp( &T, radix, &p ) );
423  }
424 
425  *p++ = '\0';
426  *slen = p - s;
427 
428 cleanup:
429 
430  mpi_free( &T );
431 
432  return( ret );
433 }
434 
435 #if defined(POLARSSL_FS_IO)
436 /*
437  * Read X from an opened file
438  */
439 int mpi_read_file( mpi *X, int radix, FILE *fin )
440 {
441  t_uint d;
442  size_t slen;
443  char *p;
444  /*
445  * Buffer should have space for (short) label and decimal formatted MPI,
446  * newline characters and '\0'
447  */
448  char s[ POLARSSL_MPI_RW_BUFFER_SIZE ];
449 
450  memset( s, 0, sizeof( s ) );
451  if( fgets( s, sizeof( s ) - 1, fin ) == NULL )
453 
454  slen = strlen( s );
455  if( slen == sizeof( s ) - 2 )
457 
458  if( s[slen - 1] == '\n' ) { slen--; s[slen] = '\0'; }
459  if( s[slen - 1] == '\r' ) { slen--; s[slen] = '\0'; }
460 
461  p = s + slen;
462  while( --p >= s )
463  if( mpi_get_digit( &d, radix, *p ) != 0 )
464  break;
465 
466  return( mpi_read_string( X, radix, p + 1 ) );
467 }
468 
469 /*
470  * Write X into an opened file (or stdout if fout == NULL)
471  */
472 int mpi_write_file( const char *p, const mpi *X, int radix, FILE *fout )
473 {
474  int ret;
475  size_t n, slen, plen;
476  /*
477  * Buffer should have space for (short) label and decimal formatted MPI,
478  * newline characters and '\0'
479  */
480  char s[ POLARSSL_MPI_RW_BUFFER_SIZE ];
481 
482  n = sizeof( s );
483  memset( s, 0, n );
484  n -= 2;
485 
486  MPI_CHK( mpi_write_string( X, radix, s, (size_t *) &n ) );
487 
488  if( p == NULL ) p = "";
489 
490  plen = strlen( p );
491  slen = strlen( s );
492  s[slen++] = '\r';
493  s[slen++] = '\n';
494 
495  if( fout != NULL )
496  {
497  if( fwrite( p, 1, plen, fout ) != plen ||
498  fwrite( s, 1, slen, fout ) != slen )
500  }
501  else
502  printf( "%s%s", p, s );
503 
504 cleanup:
505 
506  return( ret );
507 }
508 #endif /* POLARSSL_FS_IO */
509 
510 /*
511  * Import X from unsigned binary data, big endian
512  */
513 int mpi_read_binary( mpi *X, const unsigned char *buf, size_t buflen )
514 {
515  int ret;
516  size_t i, j, n;
517 
518  for( n = 0; n < buflen; n++ )
519  if( buf[n] != 0 )
520  break;
521 
522  MPI_CHK( mpi_grow( X, CHARS_TO_LIMBS( buflen - n ) ) );
523  MPI_CHK( mpi_lset( X, 0 ) );
524 
525  for( i = buflen, j = 0; i > n; i--, j++ )
526  X->p[j / ciL] |= ((t_uint) buf[i - 1]) << ((j % ciL) << 3);
527 
528 cleanup:
529 
530  return( ret );
531 }
532 
533 /*
534  * Export X into unsigned binary data, big endian
535  */
536 int mpi_write_binary( const mpi *X, unsigned char *buf, size_t buflen )
537 {
538  size_t i, j, n;
539 
540  n = mpi_size( X );
541 
542  if( buflen < n )
544 
545  memset( buf, 0, buflen );
546 
547  for( i = buflen - 1, j = 0; n > 0; i--, j++, n-- )
548  buf[i] = (unsigned char)( X->p[j / ciL] >> ((j % ciL) << 3) );
549 
550  return( 0 );
551 }
552 
553 /*
554  * Left-shift: X <<= count
555  */
556 int mpi_shift_l( mpi *X, size_t count )
557 {
558  int ret;
559  size_t i, v0, t1;
560  t_uint r0 = 0, r1;
561 
562  v0 = count / (biL );
563  t1 = count & (biL - 1);
564 
565  i = mpi_msb( X ) + count;
566 
567  if( X->n * biL < i )
568  MPI_CHK( mpi_grow( X, BITS_TO_LIMBS( i ) ) );
569 
570  ret = 0;
571 
572  /*
573  * shift by count / limb_size
574  */
575  if( v0 > 0 )
576  {
577  for( i = X->n; i > v0; i-- )
578  X->p[i - 1] = X->p[i - v0 - 1];
579 
580  for( ; i > 0; i-- )
581  X->p[i - 1] = 0;
582  }
583 
584  /*
585  * shift by count % limb_size
586  */
587  if( t1 > 0 )
588  {
589  for( i = v0; i < X->n; i++ )
590  {
591  r1 = X->p[i] >> (biL - t1);
592  X->p[i] <<= t1;
593  X->p[i] |= r0;
594  r0 = r1;
595  }
596  }
597 
598 cleanup:
599 
600  return( ret );
601 }
602 
603 /*
604  * Right-shift: X >>= count
605  */
606 int mpi_shift_r( mpi *X, size_t count )
607 {
608  size_t i, v0, v1;
609  t_uint r0 = 0, r1;
610 
611  v0 = count / biL;
612  v1 = count & (biL - 1);
613 
614  if( v0 > X->n || ( v0 == X->n && v1 > 0 ) )
615  return mpi_lset( X, 0 );
616 
617  /*
618  * shift by count / limb_size
619  */
620  if( v0 > 0 )
621  {
622  for( i = 0; i < X->n - v0; i++ )
623  X->p[i] = X->p[i + v0];
624 
625  for( ; i < X->n; i++ )
626  X->p[i] = 0;
627  }
628 
629  /*
630  * shift by count % limb_size
631  */
632  if( v1 > 0 )
633  {
634  for( i = X->n; i > 0; i-- )
635  {
636  r1 = X->p[i - 1] << (biL - v1);
637  X->p[i - 1] >>= v1;
638  X->p[i - 1] |= r0;
639  r0 = r1;
640  }
641  }
642 
643  return( 0 );
644 }
645 
646 /*
647  * Compare unsigned values
648  */
649 int mpi_cmp_abs( const mpi *X, const mpi *Y )
650 {
651  size_t i, j;
652 
653  for( i = X->n; i > 0; i-- )
654  if( X->p[i - 1] != 0 )
655  break;
656 
657  for( j = Y->n; j > 0; j-- )
658  if( Y->p[j - 1] != 0 )
659  break;
660 
661  if( i == 0 && j == 0 )
662  return( 0 );
663 
664  if( i > j ) return( 1 );
665  if( j > i ) return( -1 );
666 
667  for( ; i > 0; i-- )
668  {
669  if( X->p[i - 1] > Y->p[i - 1] ) return( 1 );
670  if( X->p[i - 1] < Y->p[i - 1] ) return( -1 );
671  }
672 
673  return( 0 );
674 }
675 
676 /*
677  * Compare signed values
678  */
679 int mpi_cmp_mpi( const mpi *X, const mpi *Y )
680 {
681  size_t i, j;
682 
683  for( i = X->n; i > 0; i-- )
684  if( X->p[i - 1] != 0 )
685  break;
686 
687  for( j = Y->n; j > 0; j-- )
688  if( Y->p[j - 1] != 0 )
689  break;
690 
691  if( i == 0 && j == 0 )
692  return( 0 );
693 
694  if( i > j ) return( X->s );
695  if( j > i ) return( -Y->s );
696 
697  if( X->s > 0 && Y->s < 0 ) return( 1 );
698  if( Y->s > 0 && X->s < 0 ) return( -1 );
699 
700  for( ; i > 0; i-- )
701  {
702  if( X->p[i - 1] > Y->p[i - 1] ) return( X->s );
703  if( X->p[i - 1] < Y->p[i - 1] ) return( -X->s );
704  }
705 
706  return( 0 );
707 }
708 
709 /*
710  * Compare signed values
711  */
712 int mpi_cmp_int( const mpi *X, t_sint z )
713 {
714  mpi Y;
715  t_uint p[1];
716 
717  *p = ( z < 0 ) ? -z : z;
718  Y.s = ( z < 0 ) ? -1 : 1;
719  Y.n = 1;
720  Y.p = p;
721 
722  return( mpi_cmp_mpi( X, &Y ) );
723 }
724 
725 /*
726  * Unsigned addition: X = |A| + |B| (HAC 14.7)
727  */
728 int mpi_add_abs( mpi *X, const mpi *A, const mpi *B )
729 {
730  int ret;
731  size_t i, j;
732  t_uint *o, *p, c;
733 
734  if( X == B )
735  {
736  const mpi *T = A; A = X; B = T;
737  }
738 
739  if( X != A )
740  MPI_CHK( mpi_copy( X, A ) );
741 
742  /*
743  * X should always be positive as a result of unsigned additions.
744  */
745  X->s = 1;
746 
747  for( j = B->n; j > 0; j-- )
748  if( B->p[j - 1] != 0 )
749  break;
750 
751  MPI_CHK( mpi_grow( X, j ) );
752 
753  o = B->p; p = X->p; c = 0;
754 
755  for( i = 0; i < j; i++, o++, p++ )
756  {
757  *p += c; c = ( *p < c );
758  *p += *o; c += ( *p < *o );
759  }
760 
761  while( c != 0 )
762  {
763  if( i >= X->n )
764  {
765  MPI_CHK( mpi_grow( X, i + 1 ) );
766  p = X->p + i;
767  }
768 
769  *p += c; c = ( *p < c ); i++; p++;
770  }
771 
772 cleanup:
773 
774  return( ret );
775 }
776 
777 /*
778  * Helper for mpi substraction
779  */
780 static void mpi_sub_hlp( size_t n, t_uint *s, t_uint *d )
781 {
782  size_t i;
783  t_uint c, z;
784 
785  for( i = c = 0; i < n; i++, s++, d++ )
786  {
787  z = ( *d < c ); *d -= c;
788  c = ( *d < *s ) + z; *d -= *s;
789  }
790 
791  while( c != 0 )
792  {
793  z = ( *d < c ); *d -= c;
794  c = z; i++; d++;
795  }
796 }
797 
798 /*
799  * Unsigned substraction: X = |A| - |B| (HAC 14.9)
800  */
801 int mpi_sub_abs( mpi *X, const mpi *A, const mpi *B )
802 {
803  mpi TB;
804  int ret;
805  size_t n;
806 
807  if( mpi_cmp_abs( A, B ) < 0 )
809 
810  mpi_init( &TB );
811 
812  if( X == B )
813  {
814  MPI_CHK( mpi_copy( &TB, B ) );
815  B = &TB;
816  }
817 
818  if( X != A )
819  MPI_CHK( mpi_copy( X, A ) );
820 
821  /*
822  * X should always be positive as a result of unsigned substractions.
823  */
824  X->s = 1;
825 
826  ret = 0;
827 
828  for( n = B->n; n > 0; n-- )
829  if( B->p[n - 1] != 0 )
830  break;
831 
832  mpi_sub_hlp( n, B->p, X->p );
833 
834 cleanup:
835 
836  mpi_free( &TB );
837 
838  return( ret );
839 }
840 
841 /*
842  * Signed addition: X = A + B
843  */
844 int mpi_add_mpi( mpi *X, const mpi *A, const mpi *B )
845 {
846  int ret, s = A->s;
847 
848  if( A->s * B->s < 0 )
849  {
850  if( mpi_cmp_abs( A, B ) >= 0 )
851  {
852  MPI_CHK( mpi_sub_abs( X, A, B ) );
853  X->s = s;
854  }
855  else
856  {
857  MPI_CHK( mpi_sub_abs( X, B, A ) );
858  X->s = -s;
859  }
860  }
861  else
862  {
863  MPI_CHK( mpi_add_abs( X, A, B ) );
864  X->s = s;
865  }
866 
867 cleanup:
868 
869  return( ret );
870 }
871 
872 /*
873  * Signed substraction: X = A - B
874  */
875 int mpi_sub_mpi( mpi *X, const mpi *A, const mpi *B )
876 {
877  int ret, s = A->s;
878 
879  if( A->s * B->s > 0 )
880  {
881  if( mpi_cmp_abs( A, B ) >= 0 )
882  {
883  MPI_CHK( mpi_sub_abs( X, A, B ) );
884  X->s = s;
885  }
886  else
887  {
888  MPI_CHK( mpi_sub_abs( X, B, A ) );
889  X->s = -s;
890  }
891  }
892  else
893  {
894  MPI_CHK( mpi_add_abs( X, A, B ) );
895  X->s = s;
896  }
897 
898 cleanup:
899 
900  return( ret );
901 }
902 
903 /*
904  * Signed addition: X = A + b
905  */
906 int mpi_add_int( mpi *X, const mpi *A, t_sint b )
907 {
908  mpi _B;
909  t_uint p[1];
910 
911  p[0] = ( b < 0 ) ? -b : b;
912  _B.s = ( b < 0 ) ? -1 : 1;
913  _B.n = 1;
914  _B.p = p;
915 
916  return( mpi_add_mpi( X, A, &_B ) );
917 }
918 
919 /*
920  * Signed substraction: X = A - b
921  */
922 int mpi_sub_int( mpi *X, const mpi *A, t_sint b )
923 {
924  mpi _B;
925  t_uint p[1];
926 
927  p[0] = ( b < 0 ) ? -b : b;
928  _B.s = ( b < 0 ) ? -1 : 1;
929  _B.n = 1;
930  _B.p = p;
931 
932  return( mpi_sub_mpi( X, A, &_B ) );
933 }
934 
935 /*
936  * Helper for mpi multiplication
937  */
938 static
939 #if defined(__APPLE__) && defined(__arm__)
940 /*
941  * Apple LLVM version 4.2 (clang-425.0.24) (based on LLVM 3.2svn)
942  * appears to need this to prevent bad ARM code generation at -O3.
943  */
944 __attribute__ ((noinline))
945 #endif
946 void mpi_mul_hlp( size_t i, t_uint *s, t_uint *d, t_uint b )
947 {
948  t_uint c = 0, t = 0;
949 
950 #if defined(MULADDC_HUIT)
951  for( ; i >= 8; i -= 8 )
952  {
954  MULADDC_HUIT
956  }
957 
958  for( ; i > 0; i-- )
959  {
963  }
964 #else
965  for( ; i >= 16; i -= 16 )
966  {
972 
978  }
979 
980  for( ; i >= 8; i -= 8 )
981  {
985 
989  }
990 
991  for( ; i > 0; i-- )
992  {
996  }
997 #endif
998 
999  t++;
1000 
1001  do {
1002  *d += c; c = ( *d < c ); d++;
1003  }
1004  while( c != 0 );
1005 }
1006 
1007 /*
1008  * Baseline multiplication: X = A * B (HAC 14.12)
1009  */
1010 int mpi_mul_mpi( mpi *X, const mpi *A, const mpi *B )
1011 {
1012  int ret;
1013  size_t i, j;
1014  mpi TA, TB;
1015 
1016  mpi_init( &TA ); mpi_init( &TB );
1017 
1018  if( X == A ) { MPI_CHK( mpi_copy( &TA, A ) ); A = &TA; }
1019  if( X == B ) { MPI_CHK( mpi_copy( &TB, B ) ); B = &TB; }
1020 
1021  for( i = A->n; i > 0; i-- )
1022  if( A->p[i - 1] != 0 )
1023  break;
1024 
1025  for( j = B->n; j > 0; j-- )
1026  if( B->p[j - 1] != 0 )
1027  break;
1028 
1029  MPI_CHK( mpi_grow( X, i + j ) );
1030  MPI_CHK( mpi_lset( X, 0 ) );
1031 
1032  for( i++; j > 0; j-- )
1033  mpi_mul_hlp( i - 1, A->p, X->p + j - 1, B->p[j - 1] );
1034 
1035  X->s = A->s * B->s;
1036 
1037 cleanup:
1038 
1039  mpi_free( &TB ); mpi_free( &TA );
1040 
1041  return( ret );
1042 }
1043 
1044 /*
1045  * Baseline multiplication: X = A * b
1046  */
1047 int mpi_mul_int( mpi *X, const mpi *A, t_sint b )
1048 {
1049  mpi _B;
1050  t_uint p[1];
1051 
1052  _B.s = 1;
1053  _B.n = 1;
1054  _B.p = p;
1055  p[0] = b;
1056 
1057  return( mpi_mul_mpi( X, A, &_B ) );
1058 }
1059 
1060 /*
1061  * Division by mpi: A = Q * B + R (HAC 14.20)
1062  */
1063 int mpi_div_mpi( mpi *Q, mpi *R, const mpi *A, const mpi *B )
1064 {
1065  int ret;
1066  size_t i, n, t, k;
1067  mpi X, Y, Z, T1, T2;
1068 
1069  if( mpi_cmp_int( B, 0 ) == 0 )
1071 
1072  mpi_init( &X ); mpi_init( &Y ); mpi_init( &Z );
1073  mpi_init( &T1 ); mpi_init( &T2 );
1074 
1075  if( mpi_cmp_abs( A, B ) < 0 )
1076  {
1077  if( Q != NULL ) MPI_CHK( mpi_lset( Q, 0 ) );
1078  if( R != NULL ) MPI_CHK( mpi_copy( R, A ) );
1079  return( 0 );
1080  }
1081 
1082  MPI_CHK( mpi_copy( &X, A ) );
1083  MPI_CHK( mpi_copy( &Y, B ) );
1084  X.s = Y.s = 1;
1085 
1086  MPI_CHK( mpi_grow( &Z, A->n + 2 ) );
1087  MPI_CHK( mpi_lset( &Z, 0 ) );
1088  MPI_CHK( mpi_grow( &T1, 2 ) );
1089  MPI_CHK( mpi_grow( &T2, 3 ) );
1090 
1091  k = mpi_msb( &Y ) % biL;
1092  if( k < biL - 1 )
1093  {
1094  k = biL - 1 - k;
1095  MPI_CHK( mpi_shift_l( &X, k ) );
1096  MPI_CHK( mpi_shift_l( &Y, k ) );
1097  }
1098  else k = 0;
1099 
1100  n = X.n - 1;
1101  t = Y.n - 1;
1102  MPI_CHK( mpi_shift_l( &Y, biL * (n - t) ) );
1103 
1104  while( mpi_cmp_mpi( &X, &Y ) >= 0 )
1105  {
1106  Z.p[n - t]++;
1107  mpi_sub_mpi( &X, &X, &Y );
1108  }
1109  mpi_shift_r( &Y, biL * (n - t) );
1110 
1111  for( i = n; i > t ; i-- )
1112  {
1113  if( X.p[i] >= Y.p[t] )
1114  Z.p[i - t - 1] = ~0;
1115  else
1116  {
1117 #if defined(POLARSSL_HAVE_UDBL)
1118  t_udbl r;
1119 
1120  r = (t_udbl) X.p[i] << biL;
1121  r |= (t_udbl) X.p[i - 1];
1122  r /= Y.p[t];
1123  if( r > ((t_udbl) 1 << biL) - 1)
1124  r = ((t_udbl) 1 << biL) - 1;
1125 
1126  Z.p[i - t - 1] = (t_uint) r;
1127 #else
1128  /*
1129  * __udiv_qrnnd_c, from gmp/longlong.h
1130  */
1131  t_uint q0, q1, r0, r1;
1132  t_uint d0, d1, d, m;
1133 
1134  d = Y.p[t];
1135  d0 = ( d << biH ) >> biH;
1136  d1 = ( d >> biH );
1137 
1138  q1 = X.p[i] / d1;
1139  r1 = X.p[i] - d1 * q1;
1140  r1 <<= biH;
1141  r1 |= ( X.p[i - 1] >> biH );
1142 
1143  m = q1 * d0;
1144  if( r1 < m )
1145  {
1146  q1--, r1 += d;
1147  while( r1 >= d && r1 < m )
1148  q1--, r1 += d;
1149  }
1150  r1 -= m;
1151 
1152  q0 = r1 / d1;
1153  r0 = r1 - d1 * q0;
1154  r0 <<= biH;
1155  r0 |= ( X.p[i - 1] << biH ) >> biH;
1156 
1157  m = q0 * d0;
1158  if( r0 < m )
1159  {
1160  q0--, r0 += d;
1161  while( r0 >= d && r0 < m )
1162  q0--, r0 += d;
1163  }
1164  r0 -= m;
1165 
1166  Z.p[i - t - 1] = ( q1 << biH ) | q0;
1167 #endif
1168  }
1169 
1170  Z.p[i - t - 1]++;
1171  do
1172  {
1173  Z.p[i - t - 1]--;
1174 
1175  MPI_CHK( mpi_lset( &T1, 0 ) );
1176  T1.p[0] = (t < 1) ? 0 : Y.p[t - 1];
1177  T1.p[1] = Y.p[t];
1178  MPI_CHK( mpi_mul_int( &T1, &T1, Z.p[i - t - 1] ) );
1179 
1180  MPI_CHK( mpi_lset( &T2, 0 ) );
1181  T2.p[0] = (i < 2) ? 0 : X.p[i - 2];
1182  T2.p[1] = (i < 1) ? 0 : X.p[i - 1];
1183  T2.p[2] = X.p[i];
1184  }
1185  while( mpi_cmp_mpi( &T1, &T2 ) > 0 );
1186 
1187  MPI_CHK( mpi_mul_int( &T1, &Y, Z.p[i - t - 1] ) );
1188  MPI_CHK( mpi_shift_l( &T1, biL * (i - t - 1) ) );
1189  MPI_CHK( mpi_sub_mpi( &X, &X, &T1 ) );
1190 
1191  if( mpi_cmp_int( &X, 0 ) < 0 )
1192  {
1193  MPI_CHK( mpi_copy( &T1, &Y ) );
1194  MPI_CHK( mpi_shift_l( &T1, biL * (i - t - 1) ) );
1195  MPI_CHK( mpi_add_mpi( &X, &X, &T1 ) );
1196  Z.p[i - t - 1]--;
1197  }
1198  }
1199 
1200  if( Q != NULL )
1201  {
1202  mpi_copy( Q, &Z );
1203  Q->s = A->s * B->s;
1204  }
1205 
1206  if( R != NULL )
1207  {
1208  mpi_shift_r( &X, k );
1209  X.s = A->s;
1210  mpi_copy( R, &X );
1211 
1212  if( mpi_cmp_int( R, 0 ) == 0 )
1213  R->s = 1;
1214  }
1215 
1216 cleanup:
1217 
1218  mpi_free( &X ); mpi_free( &Y ); mpi_free( &Z );
1219  mpi_free( &T1 ); mpi_free( &T2 );
1220 
1221  return( ret );
1222 }
1223 
1224 /*
1225  * Division by int: A = Q * b + R
1226  */
1227 int mpi_div_int( mpi *Q, mpi *R, const mpi *A, t_sint b )
1228 {
1229  mpi _B;
1230  t_uint p[1];
1231 
1232  p[0] = ( b < 0 ) ? -b : b;
1233  _B.s = ( b < 0 ) ? -1 : 1;
1234  _B.n = 1;
1235  _B.p = p;
1236 
1237  return( mpi_div_mpi( Q, R, A, &_B ) );
1238 }
1239 
1240 /*
1241  * Modulo: R = A mod B
1242  */
1243 int mpi_mod_mpi( mpi *R, const mpi *A, const mpi *B )
1244 {
1245  int ret;
1246 
1247  if( mpi_cmp_int( B, 0 ) < 0 )
1249 
1250  MPI_CHK( mpi_div_mpi( NULL, R, A, B ) );
1251 
1252  while( mpi_cmp_int( R, 0 ) < 0 )
1253  MPI_CHK( mpi_add_mpi( R, R, B ) );
1254 
1255  while( mpi_cmp_mpi( R, B ) >= 0 )
1256  MPI_CHK( mpi_sub_mpi( R, R, B ) );
1257 
1258 cleanup:
1259 
1260  return( ret );
1261 }
1262 
1263 /*
1264  * Modulo: r = A mod b
1265  */
1266 int mpi_mod_int( t_uint *r, const mpi *A, t_sint b )
1267 {
1268  size_t i;
1269  t_uint x, y, z;
1270 
1271  if( b == 0 )
1273 
1274  if( b < 0 )
1276 
1277  /*
1278  * handle trivial cases
1279  */
1280  if( b == 1 )
1281  {
1282  *r = 0;
1283  return( 0 );
1284  }
1285 
1286  if( b == 2 )
1287  {
1288  *r = A->p[0] & 1;
1289  return( 0 );
1290  }
1291 
1292  /*
1293  * general case
1294  */
1295  for( i = A->n, y = 0; i > 0; i-- )
1296  {
1297  x = A->p[i - 1];
1298  y = ( y << biH ) | ( x >> biH );
1299  z = y / b;
1300  y -= z * b;
1301 
1302  x <<= biH;
1303  y = ( y << biH ) | ( x >> biH );
1304  z = y / b;
1305  y -= z * b;
1306  }
1307 
1308  /*
1309  * If A is negative, then the current y represents a negative value.
1310  * Flipping it to the positive side.
1311  */
1312  if( A->s < 0 && y != 0 )
1313  y = b - y;
1314 
1315  *r = y;
1316 
1317  return( 0 );
1318 }
1319 
1320 /*
1321  * Fast Montgomery initialization (thanks to Tom St Denis)
1322  */
1323 static void mpi_montg_init( t_uint *mm, const mpi *N )
1324 {
1325  t_uint x, m0 = N->p[0];
1326 
1327  x = m0;
1328  x += ( ( m0 + 2 ) & 4 ) << 1;
1329  x *= ( 2 - ( m0 * x ) );
1330 
1331  if( biL >= 16 ) x *= ( 2 - ( m0 * x ) );
1332  if( biL >= 32 ) x *= ( 2 - ( m0 * x ) );
1333  if( biL >= 64 ) x *= ( 2 - ( m0 * x ) );
1334 
1335  *mm = ~x + 1;
1336 }
1337 
1338 /*
1339  * Montgomery multiplication: A = A * B * R^-1 mod N (HAC 14.36)
1340  */
1341 static void mpi_montmul( mpi *A, const mpi *B, const mpi *N, t_uint mm, const mpi *T )
1342 {
1343  size_t i, n, m;
1344  t_uint u0, u1, *d;
1345 
1346  memset( T->p, 0, T->n * ciL );
1347 
1348  d = T->p;
1349  n = N->n;
1350  m = ( B->n < n ) ? B->n : n;
1351 
1352  for( i = 0; i < n; i++ )
1353  {
1354  /*
1355  * T = (T + u0*B + u1*N) / 2^biL
1356  */
1357  u0 = A->p[i];
1358  u1 = ( d[0] + u0 * B->p[0] ) * mm;
1359 
1360  mpi_mul_hlp( m, B->p, d, u0 );
1361  mpi_mul_hlp( n, N->p, d, u1 );
1362 
1363  *d++ = u0; d[n + 1] = 0;
1364  }
1365 
1366  memcpy( A->p, d, (n + 1) * ciL );
1367 
1368  if( mpi_cmp_abs( A, N ) >= 0 )
1369  mpi_sub_hlp( n, N->p, A->p );
1370  else
1371  /* prevent timing attacks */
1372  mpi_sub_hlp( n, A->p, T->p );
1373 }
1374 
1375 /*
1376  * Montgomery reduction: A = A * R^-1 mod N
1377  */
1378 static void mpi_montred( mpi *A, const mpi *N, t_uint mm, const mpi *T )
1379 {
1380  t_uint z = 1;
1381  mpi U;
1382 
1383  U.n = U.s = (int) z;
1384  U.p = &z;
1385 
1386  mpi_montmul( A, &U, N, mm, T );
1387 }
1388 
1389 /*
1390  * Sliding-window exponentiation: X = A^E mod N (HAC 14.85)
1391  */
1392 int mpi_exp_mod( mpi *X, const mpi *A, const mpi *E, const mpi *N, mpi *_RR )
1393 {
1394  int ret;
1395  size_t wbits, wsize, one = 1;
1396  size_t i, j, nblimbs;
1397  size_t bufsize, nbits;
1398  t_uint ei, mm, state;
1399  mpi RR, T, W[ 2 << POLARSSL_MPI_WINDOW_SIZE ], Apos;
1400  int neg;
1401 
1402  if( mpi_cmp_int( N, 0 ) < 0 || ( N->p[0] & 1 ) == 0 )
1404 
1405  if( mpi_cmp_int( E, 0 ) < 0 )
1407 
1408  /*
1409  * Init temps and window size
1410  */
1411  mpi_montg_init( &mm, N );
1412  mpi_init( &RR ); mpi_init( &T );
1413  memset( W, 0, sizeof( W ) );
1414 
1415  i = mpi_msb( E );
1416 
1417  wsize = ( i > 671 ) ? 6 : ( i > 239 ) ? 5 :
1418  ( i > 79 ) ? 4 : ( i > 23 ) ? 3 : 1;
1419 
1420  if( wsize > POLARSSL_MPI_WINDOW_SIZE )
1421  wsize = POLARSSL_MPI_WINDOW_SIZE;
1422 
1423  j = N->n + 1;
1424  MPI_CHK( mpi_grow( X, j ) );
1425  MPI_CHK( mpi_grow( &W[1], j ) );
1426  MPI_CHK( mpi_grow( &T, j * 2 ) );
1427 
1428  /*
1429  * Compensate for negative A (and correct at the end)
1430  */
1431  neg = ( A->s == -1 );
1432 
1433  mpi_init( &Apos );
1434  if( neg )
1435  {
1436  MPI_CHK( mpi_copy( &Apos, A ) );
1437  Apos.s = 1;
1438  A = &Apos;
1439  }
1440 
1441  /*
1442  * If 1st call, pre-compute R^2 mod N
1443  */
1444  if( _RR == NULL || _RR->p == NULL )
1445  {
1446  MPI_CHK( mpi_lset( &RR, 1 ) );
1447  MPI_CHK( mpi_shift_l( &RR, N->n * 2 * biL ) );
1448  MPI_CHK( mpi_mod_mpi( &RR, &RR, N ) );
1449 
1450  if( _RR != NULL )
1451  memcpy( _RR, &RR, sizeof( mpi ) );
1452  }
1453  else
1454  memcpy( &RR, _RR, sizeof( mpi ) );
1455 
1456  /*
1457  * W[1] = A * R^2 * R^-1 mod N = A * R mod N
1458  */
1459  if( mpi_cmp_mpi( A, N ) >= 0 )
1460  mpi_mod_mpi( &W[1], A, N );
1461  else mpi_copy( &W[1], A );
1462 
1463  mpi_montmul( &W[1], &RR, N, mm, &T );
1464 
1465  /*
1466  * X = R^2 * R^-1 mod N = R mod N
1467  */
1468  MPI_CHK( mpi_copy( X, &RR ) );
1469  mpi_montred( X, N, mm, &T );
1470 
1471  if( wsize > 1 )
1472  {
1473  /*
1474  * W[1 << (wsize - 1)] = W[1] ^ (wsize - 1)
1475  */
1476  j = one << (wsize - 1);
1477 
1478  MPI_CHK( mpi_grow( &W[j], N->n + 1 ) );
1479  MPI_CHK( mpi_copy( &W[j], &W[1] ) );
1480 
1481  for( i = 0; i < wsize - 1; i++ )
1482  mpi_montmul( &W[j], &W[j], N, mm, &T );
1483 
1484  /*
1485  * W[i] = W[i - 1] * W[1]
1486  */
1487  for( i = j + 1; i < (one << wsize); i++ )
1488  {
1489  MPI_CHK( mpi_grow( &W[i], N->n + 1 ) );
1490  MPI_CHK( mpi_copy( &W[i], &W[i - 1] ) );
1491 
1492  mpi_montmul( &W[i], &W[1], N, mm, &T );
1493  }
1494  }
1495 
1496  nblimbs = E->n;
1497  bufsize = 0;
1498  nbits = 0;
1499  wbits = 0;
1500  state = 0;
1501 
1502  while( 1 )
1503  {
1504  if( bufsize == 0 )
1505  {
1506  if( nblimbs-- == 0 )
1507  break;
1508 
1509  bufsize = sizeof( t_uint ) << 3;
1510  }
1511 
1512  bufsize--;
1513 
1514  ei = (E->p[nblimbs] >> bufsize) & 1;
1515 
1516  /*
1517  * skip leading 0s
1518  */
1519  if( ei == 0 && state == 0 )
1520  continue;
1521 
1522  if( ei == 0 && state == 1 )
1523  {
1524  /*
1525  * out of window, square X
1526  */
1527  mpi_montmul( X, X, N, mm, &T );
1528  continue;
1529  }
1530 
1531  /*
1532  * add ei to current window
1533  */
1534  state = 2;
1535 
1536  nbits++;
1537  wbits |= (ei << (wsize - nbits));
1538 
1539  if( nbits == wsize )
1540  {
1541  /*
1542  * X = X^wsize R^-1 mod N
1543  */
1544  for( i = 0; i < wsize; i++ )
1545  mpi_montmul( X, X, N, mm, &T );
1546 
1547  /*
1548  * X = X * W[wbits] R^-1 mod N
1549  */
1550  mpi_montmul( X, &W[wbits], N, mm, &T );
1551 
1552  state--;
1553  nbits = 0;
1554  wbits = 0;
1555  }
1556  }
1557 
1558  /*
1559  * process the remaining bits
1560  */
1561  for( i = 0; i < nbits; i++ )
1562  {
1563  mpi_montmul( X, X, N, mm, &T );
1564 
1565  wbits <<= 1;
1566 
1567  if( (wbits & (one << wsize)) != 0 )
1568  mpi_montmul( X, &W[1], N, mm, &T );
1569  }
1570 
1571  /*
1572  * X = A^E * R * R^-1 mod N = A^E mod N
1573  */
1574  mpi_montred( X, N, mm, &T );
1575 
1576  if( neg )
1577  {
1578  X->s = -1;
1579  mpi_add_mpi( X, N, X );
1580  }
1581 
1582 cleanup:
1583 
1584  for( i = (one << (wsize - 1)); i < (one << wsize); i++ )
1585  mpi_free( &W[i] );
1586 
1587  mpi_free( &W[1] ); mpi_free( &T ); mpi_free( &Apos );
1588 
1589  if( _RR == NULL )
1590  mpi_free( &RR );
1591 
1592  return( ret );
1593 }
1594 
1595 /*
1596  * Greatest common divisor: G = gcd(A, B) (HAC 14.54)
1597  */
1598 int mpi_gcd( mpi *G, const mpi *A, const mpi *B )
1599 {
1600  int ret;
1601  size_t lz, lzt;
1602  mpi TG, TA, TB;
1603 
1604  mpi_init( &TG ); mpi_init( &TA ); mpi_init( &TB );
1605 
1606  MPI_CHK( mpi_copy( &TA, A ) );
1607  MPI_CHK( mpi_copy( &TB, B ) );
1608 
1609  lz = mpi_lsb( &TA );
1610  lzt = mpi_lsb( &TB );
1611 
1612  if ( lzt < lz )
1613  lz = lzt;
1614 
1615  MPI_CHK( mpi_shift_r( &TA, lz ) );
1616  MPI_CHK( mpi_shift_r( &TB, lz ) );
1617 
1618  TA.s = TB.s = 1;
1619 
1620  while( mpi_cmp_int( &TA, 0 ) != 0 )
1621  {
1622  MPI_CHK( mpi_shift_r( &TA, mpi_lsb( &TA ) ) );
1623  MPI_CHK( mpi_shift_r( &TB, mpi_lsb( &TB ) ) );
1624 
1625  if( mpi_cmp_mpi( &TA, &TB ) >= 0 )
1626  {
1627  MPI_CHK( mpi_sub_abs( &TA, &TA, &TB ) );
1628  MPI_CHK( mpi_shift_r( &TA, 1 ) );
1629  }
1630  else
1631  {
1632  MPI_CHK( mpi_sub_abs( &TB, &TB, &TA ) );
1633  MPI_CHK( mpi_shift_r( &TB, 1 ) );
1634  }
1635  }
1636 
1637  MPI_CHK( mpi_shift_l( &TB, lz ) );
1638  MPI_CHK( mpi_copy( G, &TB ) );
1639 
1640 cleanup:
1641 
1642  mpi_free( &TG ); mpi_free( &TA ); mpi_free( &TB );
1643 
1644  return( ret );
1645 }
1646 
1647 int mpi_fill_random( mpi *X, size_t size,
1648  int (*f_rng)(void *, unsigned char *, size_t),
1649  void *p_rng )
1650 {
1651  int ret;
1652 
1653  MPI_CHK( mpi_grow( X, CHARS_TO_LIMBS( size ) ) );
1654  MPI_CHK( mpi_lset( X, 0 ) );
1655 
1656  MPI_CHK( f_rng( p_rng, (unsigned char *) X->p, size ) );
1657 
1658 cleanup:
1659  return( ret );
1660 }
1661 
1662 /*
1663  * Modular inverse: X = A^-1 mod N (HAC 14.61 / 14.64)
1664  */
1665 int mpi_inv_mod( mpi *X, const mpi *A, const mpi *N )
1666 {
1667  int ret;
1668  mpi G, TA, TU, U1, U2, TB, TV, V1, V2;
1669 
1670  if( mpi_cmp_int( N, 0 ) <= 0 )
1672 
1673  mpi_init( &TA ); mpi_init( &TU ); mpi_init( &U1 ); mpi_init( &U2 );
1674  mpi_init( &G ); mpi_init( &TB ); mpi_init( &TV );
1675  mpi_init( &V1 ); mpi_init( &V2 );
1676 
1677  MPI_CHK( mpi_gcd( &G, A, N ) );
1678 
1679  if( mpi_cmp_int( &G, 1 ) != 0 )
1680  {
1682  goto cleanup;
1683  }
1684 
1685  MPI_CHK( mpi_mod_mpi( &TA, A, N ) );
1686  MPI_CHK( mpi_copy( &TU, &TA ) );
1687  MPI_CHK( mpi_copy( &TB, N ) );
1688  MPI_CHK( mpi_copy( &TV, N ) );
1689 
1690  MPI_CHK( mpi_lset( &U1, 1 ) );
1691  MPI_CHK( mpi_lset( &U2, 0 ) );
1692  MPI_CHK( mpi_lset( &V1, 0 ) );
1693  MPI_CHK( mpi_lset( &V2, 1 ) );
1694 
1695  do
1696  {
1697  while( ( TU.p[0] & 1 ) == 0 )
1698  {
1699  MPI_CHK( mpi_shift_r( &TU, 1 ) );
1700 
1701  if( ( U1.p[0] & 1 ) != 0 || ( U2.p[0] & 1 ) != 0 )
1702  {
1703  MPI_CHK( mpi_add_mpi( &U1, &U1, &TB ) );
1704  MPI_CHK( mpi_sub_mpi( &U2, &U2, &TA ) );
1705  }
1706 
1707  MPI_CHK( mpi_shift_r( &U1, 1 ) );
1708  MPI_CHK( mpi_shift_r( &U2, 1 ) );
1709  }
1710 
1711  while( ( TV.p[0] & 1 ) == 0 )
1712  {
1713  MPI_CHK( mpi_shift_r( &TV, 1 ) );
1714 
1715  if( ( V1.p[0] & 1 ) != 0 || ( V2.p[0] & 1 ) != 0 )
1716  {
1717  MPI_CHK( mpi_add_mpi( &V1, &V1, &TB ) );
1718  MPI_CHK( mpi_sub_mpi( &V2, &V2, &TA ) );
1719  }
1720 
1721  MPI_CHK( mpi_shift_r( &V1, 1 ) );
1722  MPI_CHK( mpi_shift_r( &V2, 1 ) );
1723  }
1724 
1725  if( mpi_cmp_mpi( &TU, &TV ) >= 0 )
1726  {
1727  MPI_CHK( mpi_sub_mpi( &TU, &TU, &TV ) );
1728  MPI_CHK( mpi_sub_mpi( &U1, &U1, &V1 ) );
1729  MPI_CHK( mpi_sub_mpi( &U2, &U2, &V2 ) );
1730  }
1731  else
1732  {
1733  MPI_CHK( mpi_sub_mpi( &TV, &TV, &TU ) );
1734  MPI_CHK( mpi_sub_mpi( &V1, &V1, &U1 ) );
1735  MPI_CHK( mpi_sub_mpi( &V2, &V2, &U2 ) );
1736  }
1737  }
1738  while( mpi_cmp_int( &TU, 0 ) != 0 );
1739 
1740  while( mpi_cmp_int( &V1, 0 ) < 0 )
1741  MPI_CHK( mpi_add_mpi( &V1, &V1, N ) );
1742 
1743  while( mpi_cmp_mpi( &V1, N ) >= 0 )
1744  MPI_CHK( mpi_sub_mpi( &V1, &V1, N ) );
1745 
1746  MPI_CHK( mpi_copy( X, &V1 ) );
1747 
1748 cleanup:
1749 
1750  mpi_free( &TA ); mpi_free( &TU ); mpi_free( &U1 ); mpi_free( &U2 );
1751  mpi_free( &G ); mpi_free( &TB ); mpi_free( &TV );
1752  mpi_free( &V1 ); mpi_free( &V2 );
1753 
1754  return( ret );
1755 }
1756 
1757 #if defined(POLARSSL_GENPRIME)
1758 
1759 static const int small_prime[] =
1760 {
1761  3, 5, 7, 11, 13, 17, 19, 23,
1762  29, 31, 37, 41, 43, 47, 53, 59,
1763  61, 67, 71, 73, 79, 83, 89, 97,
1764  101, 103, 107, 109, 113, 127, 131, 137,
1765  139, 149, 151, 157, 163, 167, 173, 179,
1766  181, 191, 193, 197, 199, 211, 223, 227,
1767  229, 233, 239, 241, 251, 257, 263, 269,
1768  271, 277, 281, 283, 293, 307, 311, 313,
1769  317, 331, 337, 347, 349, 353, 359, 367,
1770  373, 379, 383, 389, 397, 401, 409, 419,
1771  421, 431, 433, 439, 443, 449, 457, 461,
1772  463, 467, 479, 487, 491, 499, 503, 509,
1773  521, 523, 541, 547, 557, 563, 569, 571,
1774  577, 587, 593, 599, 601, 607, 613, 617,
1775  619, 631, 641, 643, 647, 653, 659, 661,
1776  673, 677, 683, 691, 701, 709, 719, 727,
1777  733, 739, 743, 751, 757, 761, 769, 773,
1778  787, 797, 809, 811, 821, 823, 827, 829,
1779  839, 853, 857, 859, 863, 877, 881, 883,
1780  887, 907, 911, 919, 929, 937, 941, 947,
1781  953, 967, 971, 977, 983, 991, 997, -103
1782 };
1783 
1784 /*
1785  * Miller-Rabin primality test (HAC 4.24)
1786  */
1787 int mpi_is_prime( mpi *X,
1788  int (*f_rng)(void *, unsigned char *, size_t),
1789  void *p_rng )
1790 {
1791  int ret, xs;
1792  size_t i, j, n, s;
1793  mpi W, R, T, A, RR;
1794 
1795  if( mpi_cmp_int( X, 0 ) == 0 ||
1796  mpi_cmp_int( X, 1 ) == 0 )
1798 
1799  if( mpi_cmp_int( X, 2 ) == 0 )
1800  return( 0 );
1801 
1802  mpi_init( &W ); mpi_init( &R ); mpi_init( &T ); mpi_init( &A );
1803  mpi_init( &RR );
1804 
1805  xs = X->s; X->s = 1;
1806 
1807  /*
1808  * test trivial factors first
1809  */
1810  if( ( X->p[0] & 1 ) == 0 )
1812 
1813  for( i = 0; small_prime[i] > 0; i++ )
1814  {
1815  t_uint r;
1816 
1817  if( mpi_cmp_int( X, small_prime[i] ) <= 0 )
1818  return( 0 );
1819 
1820  MPI_CHK( mpi_mod_int( &r, X, small_prime[i] ) );
1821 
1822  if( r == 0 )
1824  }
1825 
1826  /*
1827  * W = |X| - 1
1828  * R = W >> lsb( W )
1829  */
1830  MPI_CHK( mpi_sub_int( &W, X, 1 ) );
1831  s = mpi_lsb( &W );
1832  MPI_CHK( mpi_copy( &R, &W ) );
1833  MPI_CHK( mpi_shift_r( &R, s ) );
1834 
1835  i = mpi_msb( X );
1836  /*
1837  * HAC, table 4.4
1838  */
1839  n = ( ( i >= 1300 ) ? 2 : ( i >= 850 ) ? 3 :
1840  ( i >= 650 ) ? 4 : ( i >= 350 ) ? 8 :
1841  ( i >= 250 ) ? 12 : ( i >= 150 ) ? 18 : 27 );
1842 
1843  for( i = 0; i < n; i++ )
1844  {
1845  /*
1846  * pick a random A, 1 < A < |X| - 1
1847  */
1848  MPI_CHK( mpi_fill_random( &A, X->n * ciL, f_rng, p_rng ) );
1849 
1850  if( mpi_cmp_mpi( &A, &W ) >= 0 )
1851  {
1852  j = mpi_msb( &A ) - mpi_msb( &W );
1853  MPI_CHK( mpi_shift_r( &A, j + 1 ) );
1854  }
1855  A.p[0] |= 3;
1856 
1857  /*
1858  * A = A^R mod |X|
1859  */
1860  MPI_CHK( mpi_exp_mod( &A, &A, &R, X, &RR ) );
1861 
1862  if( mpi_cmp_mpi( &A, &W ) == 0 ||
1863  mpi_cmp_int( &A, 1 ) == 0 )
1864  continue;
1865 
1866  j = 1;
1867  while( j < s && mpi_cmp_mpi( &A, &W ) != 0 )
1868  {
1869  /*
1870  * A = A * A mod |X|
1871  */
1872  MPI_CHK( mpi_mul_mpi( &T, &A, &A ) );
1873  MPI_CHK( mpi_mod_mpi( &A, &T, X ) );
1874 
1875  if( mpi_cmp_int( &A, 1 ) == 0 )
1876  break;
1877 
1878  j++;
1879  }
1880 
1881  /*
1882  * not prime if A != |X| - 1 or A == 1
1883  */
1884  if( mpi_cmp_mpi( &A, &W ) != 0 ||
1885  mpi_cmp_int( &A, 1 ) == 0 )
1886  {
1888  break;
1889  }
1890  }
1891 
1892 cleanup:
1893 
1894  X->s = xs;
1895 
1896  mpi_free( &W ); mpi_free( &R ); mpi_free( &T ); mpi_free( &A );
1897  mpi_free( &RR );
1898 
1899  return( ret );
1900 }
1901 
1902 /*
1903  * Prime number generation
1904  */
1905 int mpi_gen_prime( mpi *X, size_t nbits, int dh_flag,
1906  int (*f_rng)(void *, unsigned char *, size_t),
1907  void *p_rng )
1908 {
1909  int ret;
1910  size_t k, n;
1911  mpi Y;
1912 
1913  if( nbits < 3 || nbits > POLARSSL_MPI_MAX_BITS )
1915 
1916  mpi_init( &Y );
1917 
1918  n = BITS_TO_LIMBS( nbits );
1919 
1920  MPI_CHK( mpi_fill_random( X, n * ciL, f_rng, p_rng ) );
1921 
1922  k = mpi_msb( X );
1923  if( k < nbits ) MPI_CHK( mpi_shift_l( X, nbits - k ) );
1924  if( k > nbits ) MPI_CHK( mpi_shift_r( X, k - nbits ) );
1925 
1926  X->p[0] |= 3;
1927 
1928  if( dh_flag == 0 )
1929  {
1930  while( ( ret = mpi_is_prime( X, f_rng, p_rng ) ) != 0 )
1931  {
1932  if( ret != POLARSSL_ERR_MPI_NOT_ACCEPTABLE )
1933  goto cleanup;
1934 
1935  MPI_CHK( mpi_add_int( X, X, 2 ) );
1936  }
1937  }
1938  else
1939  {
1940  MPI_CHK( mpi_sub_int( &Y, X, 1 ) );
1941  MPI_CHK( mpi_shift_r( &Y, 1 ) );
1942 
1943  while( 1 )
1944  {
1945  if( ( ret = mpi_is_prime( X, f_rng, p_rng ) ) == 0 )
1946  {
1947  if( ( ret = mpi_is_prime( &Y, f_rng, p_rng ) ) == 0 )
1948  break;
1949 
1950  if( ret != POLARSSL_ERR_MPI_NOT_ACCEPTABLE )
1951  goto cleanup;
1952  }
1953 
1954  if( ret != POLARSSL_ERR_MPI_NOT_ACCEPTABLE )
1955  goto cleanup;
1956 
1957  MPI_CHK( mpi_add_int( &Y, X, 1 ) );
1958  MPI_CHK( mpi_add_int( X, X, 2 ) );
1959  MPI_CHK( mpi_shift_r( &Y, 1 ) );
1960  }
1961  }
1962 
1963 cleanup:
1964 
1965  mpi_free( &Y );
1966 
1967  return( ret );
1968 }
1969 
1970 #endif
1971 
1972 #if defined(POLARSSL_SELF_TEST)
1973 
1974 #define GCD_PAIR_COUNT 3
1975 
1976 static const int gcd_pairs[GCD_PAIR_COUNT][3] =
1977 {
1978  { 693, 609, 21 },
1979  { 1764, 868, 28 },
1980  { 768454923, 542167814, 1 }
1981 };
1982 
1983 /*
1984  * Checkup routine
1985  */
1986 int mpi_self_test( int verbose )
1987 {
1988  int ret, i;
1989  mpi A, E, N, X, Y, U, V;
1990 
1991  mpi_init( &A ); mpi_init( &E ); mpi_init( &N ); mpi_init( &X );
1992  mpi_init( &Y ); mpi_init( &U ); mpi_init( &V );
1993 
1994  MPI_CHK( mpi_read_string( &A, 16,
1995  "EFE021C2645FD1DC586E69184AF4A31E" \
1996  "D5F53E93B5F123FA41680867BA110131" \
1997  "944FE7952E2517337780CB0DB80E61AA" \
1998  "E7C8DDC6C5C6AADEB34EB38A2F40D5E6" ) );
1999 
2000  MPI_CHK( mpi_read_string( &E, 16,
2001  "B2E7EFD37075B9F03FF989C7C5051C20" \
2002  "34D2A323810251127E7BF8625A4F49A5" \
2003  "F3E27F4DA8BD59C47D6DAABA4C8127BD" \
2004  "5B5C25763222FEFCCFC38B832366C29E" ) );
2005 
2006  MPI_CHK( mpi_read_string( &N, 16,
2007  "0066A198186C18C10B2F5ED9B522752A" \
2008  "9830B69916E535C8F047518A889A43A5" \
2009  "94B6BED27A168D31D4A52F88925AA8F5" ) );
2010 
2011  MPI_CHK( mpi_mul_mpi( &X, &A, &N ) );
2012 
2013  MPI_CHK( mpi_read_string( &U, 16,
2014  "602AB7ECA597A3D6B56FF9829A5E8B85" \
2015  "9E857EA95A03512E2BAE7391688D264A" \
2016  "A5663B0341DB9CCFD2C4C5F421FEC814" \
2017  "8001B72E848A38CAE1C65F78E56ABDEF" \
2018  "E12D3C039B8A02D6BE593F0BBBDA56F1" \
2019  "ECF677152EF804370C1A305CAF3B5BF1" \
2020  "30879B56C61DE584A0F53A2447A51E" ) );
2021 
2022  if( verbose != 0 )
2023  printf( " MPI test #1 (mul_mpi): " );
2024 
2025  if( mpi_cmp_mpi( &X, &U ) != 0 )
2026  {
2027  if( verbose != 0 )
2028  printf( "failed\n" );
2029 
2030  return( 1 );
2031  }
2032 
2033  if( verbose != 0 )
2034  printf( "passed\n" );
2035 
2036  MPI_CHK( mpi_div_mpi( &X, &Y, &A, &N ) );
2037 
2038  MPI_CHK( mpi_read_string( &U, 16,
2039  "256567336059E52CAE22925474705F39A94" ) );
2040 
2041  MPI_CHK( mpi_read_string( &V, 16,
2042  "6613F26162223DF488E9CD48CC132C7A" \
2043  "0AC93C701B001B092E4E5B9F73BCD27B" \
2044  "9EE50D0657C77F374E903CDFA4C642" ) );
2045 
2046  if( verbose != 0 )
2047  printf( " MPI test #2 (div_mpi): " );
2048 
2049  if( mpi_cmp_mpi( &X, &U ) != 0 ||
2050  mpi_cmp_mpi( &Y, &V ) != 0 )
2051  {
2052  if( verbose != 0 )
2053  printf( "failed\n" );
2054 
2055  return( 1 );
2056  }
2057 
2058  if( verbose != 0 )
2059  printf( "passed\n" );
2060 
2061  MPI_CHK( mpi_exp_mod( &X, &A, &E, &N, NULL ) );
2062 
2063  MPI_CHK( mpi_read_string( &U, 16,
2064  "36E139AEA55215609D2816998ED020BB" \
2065  "BD96C37890F65171D948E9BC7CBAA4D9" \
2066  "325D24D6A3C12710F10A09FA08AB87" ) );
2067 
2068  if( verbose != 0 )
2069  printf( " MPI test #3 (exp_mod): " );
2070 
2071  if( mpi_cmp_mpi( &X, &U ) != 0 )
2072  {
2073  if( verbose != 0 )
2074  printf( "failed\n" );
2075 
2076  return( 1 );
2077  }
2078 
2079  if( verbose != 0 )
2080  printf( "passed\n" );
2081 
2082 #if defined(POLARSSL_GENPRIME)
2083  MPI_CHK( mpi_inv_mod( &X, &A, &N ) );
2084 
2085  MPI_CHK( mpi_read_string( &U, 16,
2086  "003A0AAEDD7E784FC07D8F9EC6E3BFD5" \
2087  "C3DBA76456363A10869622EAC2DD84EC" \
2088  "C5B8A74DAC4D09E03B5E0BE779F2DF61" ) );
2089 
2090  if( verbose != 0 )
2091  printf( " MPI test #4 (inv_mod): " );
2092 
2093  if( mpi_cmp_mpi( &X, &U ) != 0 )
2094  {
2095  if( verbose != 0 )
2096  printf( "failed\n" );
2097 
2098  return( 1 );
2099  }
2100 
2101  if( verbose != 0 )
2102  printf( "passed\n" );
2103 #endif
2104 
2105  if( verbose != 0 )
2106  printf( " MPI test #5 (simple gcd): " );
2107 
2108  for ( i = 0; i < GCD_PAIR_COUNT; i++)
2109  {
2110  MPI_CHK( mpi_lset( &X, gcd_pairs[i][0] ) );
2111  MPI_CHK( mpi_lset( &Y, gcd_pairs[i][1] ) );
2112 
2113  MPI_CHK( mpi_gcd( &A, &X, &Y ) );
2114 
2115  if( mpi_cmp_int( &A, gcd_pairs[i][2] ) != 0 )
2116  {
2117  if( verbose != 0 )
2118  printf( "failed at %d\n", i );
2119 
2120  return( 1 );
2121  }
2122  }
2123 
2124  if( verbose != 0 )
2125  printf( "passed\n" );
2126 
2127 cleanup:
2128 
2129  if( ret != 0 && verbose != 0 )
2130  printf( "Unexpected error, return code = %08X\n", ret );
2131 
2132  mpi_free( &A ); mpi_free( &E ); mpi_free( &N ); mpi_free( &X );
2133  mpi_free( &Y ); mpi_free( &U ); mpi_free( &V );
2134 
2135  if( verbose != 0 )
2136  printf( "\n" );
2137 
2138  return( ret );
2139 }
2140 
2141 #endif
2142 
2143 #endif