1 /*

     2 * Portions Copyright (c) 2009 Nokia Corporation and/or its subsidiary(-ies).

     3 * All rights reserved.

     4 * This component and the accompanying materials are made available

     5 * under the terms of "Eclipse Public License v1.0"

     6 * which accompanies this distribution, and is available

     7 * at the URL "http://www.eclipse.org/legal/epl-v10.html".

     8 *

     9 * Initial Contributors:

    10 * Nokia Corporation - initial contribution.

    11 *

    12 * Contributors:

    13 *

    14 * Description: 

    15 * The original NIST Statistical Test Suite code is placed in public domain.

    16 * (http://csrc.nist.gov/groups/ST/toolkit/rng/documentation_software.html) 

    17 * 

    18 * This software was developed at the National Institute of Standards and Technology by 

    19 * employees of the Federal Government in the course of their official duties. Pursuant

    20 * to title 17 Section 105 of the United States Code this software is not subject to 

    21 * copyright protection and is in the public domain. The NIST Statistical Test Suite is

    22 * an experimental system. NIST assumes no responsibility whatsoever for its use by other 

    23 * parties, and makes no guarantees, expressed or implied, about its quality, reliability, 

    24 * or any other characteristic. We would appreciate acknowledgment if the software is used.

    25 */

    26 

    27 /*

    28  * file: mp.c

    29  *

    30  * DESCRIPTION

    31  *

    32  * These functions comprise a multi-precision integer arithmetic

    33  * and discrete function package.

    34  */

    35 

    36 #include	"../include/genutils.h"

    37 

    38 #define	MAXPLEN		384

    39 

    40 

    41 /*****************************************

    42 ** greater - Test if x > y               *

    43 **                                       *

    44 ** Returns TRUE (1) if x greater than y, *

    45 ** otherwise FALSE (0).                  *

    46 **                                       *

    47 ** Parameters:                           *

    48 **                                       *

    49 **  x      Address of array x            *

    50 **  y      Address of array y            *

    51 **  l      Length both x and y in bytes  *

    52 **                                       *

    53 ******************************************/ 

    54 int greater(BYTE *x, BYTE *y, int l)

    55 {

    56 	int		i;

    57 

    58 	for ( i=0; i<l; i++ )

    59 		if ( x[i] != y[i] )

    60 			break;

    61 

    62 	if ( i == l )

    63 		return 0;

    64 

    65 	if ( x[i] > y[i] )

    66 		return 1;

    67 

    68 	return 0;

    69 }

    70 

    71 

    72 /*****************************************

    73 ** less - Test if x < y                  *

    74 **                                       *

    75 ** Returns TRUE (1) if x less than y,    *

    76 ** otherwise FALSE (0).                  *

    77 **                                       *

    78 ** Parameters:                           *

    79 **                                       *

    80 **  x      Address of array x            *

    81 **  y      Address of array y            *

    82 **  l      Length both x and y in bytes  *

    83 **                                       *

    84 ******************************************/ 

    85 int less(BYTE *x, BYTE *y, int l)

    86 {

    87 	int		i;

    88 

    89 	for ( i=0; i<l; i++ )

    90 		if ( x[i] != y[i] )

    91 			break;

    92 

    93 	if ( i == l ) {

    94 		return 0;

    95 	}

    96 

    97 	if ( x[i] < y[i] ) {

    98 		return 1;

    99 	}

   100 

   101 	return 0;

   102 }

   103 

   104 

   105 /*****************************************

   106 ** bshl - shifts array left              *

   107 **                  by one bit.          *

   108 **                                       *	

   109 ** x = x * 2                             *

   110 **                                       *

   111 ** Parameters:                           *	

   112 **                                       *

   113 **  x      Address of array x            *

   114 **  l      Length array x in bytes       *

   115 **                                       *

   116 ******************************************/ 

   117 BYTE bshl(BYTE *x, int l)

   118 {

   119 	BYTE	*p;

   120 	int		c1, c2;

   121 

   122 	p = x + l - 1;

   123 	c1 = 0;

   124 	c2 = 0;

   125 	while ( p != x ) {

   126 		if ( *p & 0x80 )

   127 			c2 = 1;

   128 		*p <<= 1;  /* shift the word left once (ls bit = 0) */

   129 		if ( c1 )

   130 			*p |= 1;

   131 		c1 = c2;

   132 		c2 = 0;

   133 		p--;

   134 	}

   135 

   136 	if ( *p & 0x80 )

   137 		c2 = 1;

   138 	*p <<= 1;  /* shift the word left once (ls bit = 0) */

   139 	if ( c1 )

   140 		*p |= (DIGIT)1;

   141 

   142 	return (BYTE)c2;

   143 }

   144 

   145 

   146 /*****************************************

   147 ** bshr - shifts array right             *

   148 **                   by one bit.         *

   149 **                                       *	

   150 ** x = x / 2                             *

   151 **                                       *

   152 ** Parameters:                           *	

   153 **                                       *

   154 **  x      Address of array x            *

   155 **  l      Length array x in bytes       *	

   156 **                                       *

   157 ******************************************/

   158 void bshr(BYTE *x, int l)	

   159 {

   160 	BYTE	*p;

   161 	int		c1,c2;

   162 

   163 	p = x;

   164 	c1 = 0;

   165 	c2 = 0;

   166 	while ( p != x+l-1 ) {

   167 		if ( *p & 0x01 )

   168 			c2 = 1;

   169 		*p >>= 1;  /* shift the word right once (ms bit = 0) */

   170 		if ( c1 )

   171 			*p |= 0x80;

   172 		c1 = c2;

   173 		c2 = 0;

   174 		p++;

   175 	}

   176 

   177 	*p >>= 1;  /* shift the word right once (ms bit = 0) */

   178 	if ( c1 )

   179 		*p |= 0x80;

   180 }

   181 

   182 

   183 /*****************************************

   184 ** Mult - Multiply two integers          *

   185 **                                       *

   186 ** A = B * C                             *

   187 **                                       *

   188 ** Parameters:                           *	

   189 **                                       *

   190 **  A      Address of the result         *

   191 **  B      Address of the multiplier     *

   192 **  C      Address of the multiplicand   *

   193 **  LB      Length of B in bytes         *

   194 **  LC      Length of C in bytes         *

   195 **                                       *

   196 **  NOTE:  A MUST be LB+LC in length     *

   197 **                                       *

   198 ******************************************/

   199 int Mult(BYTE *A, BYTE *B, int LB, BYTE *C, int LC)

   200 {

   201 	int    i, j;

   202 	int    k = 0;

   203 	DIGIT	result;

   204 

   205 

   206 	for ( i=LB-1; i>=0; i-- ) {

   207 		result = 0;

   208 		for ( j=LC-1; j>=0; j-- ) {

   209 			k = i+j+1;

   210 			result = (DIGIT)((DIGIT)A[k] + ((DIGIT)(B[i] * C[j])) + (result >> 8));

   211 			A[k] = (BYTE)result;

   212 			}

   213 		A[--k] = (BYTE)(result >> 8);

   214 	}

   215 

   216 	return 0;

   217 }

   218 

   219 

   220 void ModSqr(BYTE *A, BYTE *B, int LB, BYTE *M, int LM)

   221 {

   222 

   223 	Square(A, B, LB);

   224 	Mod(A, 2*LB, M, LM);

   225 }

   226 

   227 void ModMult(BYTE *A, BYTE *B, int LB, BYTE *C, int LC, BYTE *M, int LM)

   228 {

   229 	Mult(A, B, LB, C, LC);

   230 	Mod(A, (LB+LC), M, LM);

   231 }

   232 

   233 

   234 /*****************************************

   235 ** smult - Multiply array by a scalar.   *

   236 **                                       *

   237 ** A = b * C                             *

   238 **                                       *

   239 ** Parameters:                           *	

   240 **                                       *

   241 **  A      Address of the result         *

   242 **  b      Scalar (1 BYTE)               *

   243 **  C      Address of the multiplicand   *

   244 **  L      Length of C in bytes          *

   245 **                                       *

   246 **  NOTE:  A MUST be L+1 in length       *

   247 **                                       *

   248 ******************************************/

   249 void smult(BYTE *A, BYTE b, BYTE *C, int L)

   250 {

   251 	int		i;

   252 	DIGIT	result;

   253 

   254 	result = 0;

   255 	for ( i=L-1; i>0; i-- ) {

   256 		result = (DIGIT)(A[i] + ((DIGIT)b * C[i]) + (result >> 8));

   257 		A[i] = (BYTE)(result & 0xff);

   258 		A[i-1] = (BYTE)(result >> 8);

   259 	}

   260 }

   261 

   262 /*****************************************

   263 ** Square() - Square an integer          *

   264 **                                       *

   265 ** A = B^2                               *

   266 **                                       *

   267 ** Parameters:                           *

   268 **                                       *

   269 **  A      Address of the result         *

   270 **  B      Address of the operand        *

   271 **  L      Length of B in bytes          *

   272 **                                       *

   273 **  NOTE:  A MUST be 2*L in length       *

   274 **                                       *

   275 ******************************************/

   276 void Square(BYTE *A, BYTE *B, int L)

   277 {

   278 	Mult(A, B, L, B, L);

   279 }

   280 

   281 /*****************************************

   282 ** ModExp - Modular Exponentiation       *

   283 **                                       *

   284 ** A = B ** C (MOD M)                    *

   285 **                                       *	

   286 ** Parameters:                           *	

   287 **                                       *

   288 **  A      Address of result             *

   289 **  B      Address of mantissa           *

   290 **  C      Address of exponent           *

   291 **  M      Address of modulus            *

   292 **  LB     Length of B in bytes          *

   293 **  LC     Length of C in bytes          *

   294 **  LM     Length of M in bytes          *

   295 **                                       *

   296 **  NOTE: The integer B must be less     *

   297 **        than the modulus M.      	 *

   298 **  NOTE: A must be at least 3*LM        *

   299 **        bytes long.  However, the      *

   300 **        result stored in A will be     *

   301 **        only LM bytes long.            *

   302 ******************************************/

   303 void ModExp(BYTE *A, BYTE *B, int LB, BYTE *C, int LC, BYTE *M, int LM)

   304 {

   305 	BYTE	wmask;

   306 	int		bits;

   307 

   308 	bits = LC*8;

   309 	wmask = 0x80;

   310 

   311 	A[LM-1] = 1;

   312 

   313 	while ( !sniff_bit(C,wmask) ) {

   314 		wmask >>= 1;

   315 		bits--;

   316 		if ( !wmask ) {

   317 			wmask = 0x80;

   318 			C++;

   319 		}

   320 	}

   321 

   322 	while ( bits-- ) {

   323 		memset(A+LM, 0x00, LM*2);

   324 

   325 		/* temp = A*A (MOD M) */

   326 		ModSqr(A+LM, A,LM,  M,LM);

   327 

   328 		/* A = lower L bytes of temp */

   329 		memcpy(A, A+LM*2, LM);

   330 		memset(A+LM, 0x00, 2*LM);

   331 

   332 		if ( sniff_bit(C,wmask) ) {

   333 			memset(A+LM, 0x00, (LM+LB));

   334 			ModMult(A+LM, B,LB, A,LM,  M,LM);       /* temp = B * A (MOD M) */

   335 			memcpy(A, A+LM+(LM+LB)-LM, LM);  /* A = lower LM bytes of temp */

   336 			memset(A+LM, 0x00, 2*LM);

   337 		}

   338  

   339 		wmask >>= 1;

   340 		if ( !wmask ) {

   341 			wmask = 0x80;

   342 			C++;

   343 		}

   344 	}

   345 }

   346 

   347 

   348 /* DivMod:

   349  *

   350  *   computes:

   351  *         quot = x / n

   352  *         rem = x % n

   353  *   returns:

   354  *         length of "quot"

   355  *

   356  *  len of rem is lenx+1

   357  */

   358 int DivMod(BYTE *x, int lenx, BYTE *n, int lenn, BYTE *quot, BYTE *rem)

   359 {

   360 	BYTE	*tx, *tn, *ttx, *ts, bmult[1];

   361 	int		i, shift, lgth_x, lgth_n, t_len, lenq;

   362 	DIGIT	tMSn, mult;

   363 	ULONG	tMSx;

   364 	int		underflow;

   365 

   366 	tx = x;

   367 	tn = n;

   368 	

   369 	/* point to the MSD of n  */

   370 	for ( i=0, lgth_n=lenn; i<lenn; i++, lgth_n-- ) {

   371 		if ( *tn )

   372 			break;

   373 		tn++;

   374 	}

   375 	if ( !lgth_n )

   376 		return 0;

   377 	

   378 	/* point to the MSD of x  */

   379 	for ( i=0, lgth_x=lenx; i<lenx; i++, lgth_x-- ) {

   380 		if ( *tx )

   381 			break;

   382 		tx++;

   383 	}

   384 	if ( !lgth_x )

   385 		return 0;

   386 

   387 	if ( lgth_x < lgth_n )

   388 		lenq = 1;

   389 	else

   390 		lenq = lgth_x - lgth_n + 1;

   391 	memset(quot, 0x00, lenq);

   392 	

   393 	/* Loop while x > n,  WATCH OUT if lgth_x == lgth_n */

   394 	while ( (lgth_x > lgth_n) || ((lgth_x == lgth_n) && !less(tx, tn, lgth_n)) ) {

   395 		shift = 1;

   396 		if ( lgth_n == 1 ) {

   397 			if ( *tx < *tn ) {

   398 				tMSx = (DIGIT) (((*tx) << 8) | *(tx+1));

   399 				tMSn = *tn;

   400 				shift = 0;

   401 			}

   402 			else {

   403 				tMSx = *tx;

   404 				tMSn = *tn;

   405 			}

   406 		}

   407 		else if ( lgth_n > 1 ) {

   408 			tMSx = (DIGIT) (((*tx) << 8) | *(tx+1));

   409 			tMSn = (DIGIT) (((*tn) << 8) | *(tn+1));

   410 			if ( (tMSx < tMSn) || ((tMSx == tMSn) && less(tx, tn, lgth_n)) ) {

   411 				tMSx = (tMSx << 8) | *(tx+2);

   412 				shift = 0;

   413 			}

   414 		}

   415 		else {

   416 			tMSx = (DIGIT) (((*tx) << 8) | *(tx+1));

   417 			tMSn = *tn;

   418 			shift = 0;

   419 		}

   420 

   421 		mult = (DIGIT) (tMSx / tMSn);

   422 		if ( mult > 0xff )

   423 			mult = 0xff;

   424 		bmult[0] = (BYTE)(mult & 0xff);

   425 

   426 		ts = rem;

   427 		do {

   428 			memset(ts, 0x00, lgth_x+1);

   429 			Mult(ts, tn, lgth_n, bmult, 1);

   430 

   431 			underflow = 0;

   432 			if ( shift ) {

   433 				if ( ts[0] != 0 )

   434 					underflow = 1;

   435 				else {

   436 					for ( i=0; i<lgth_x; i++ )

   437 						ts[i] = ts[i+1];

   438 					ts[lgth_x] = 0x00;

   439 				}

   440 			}

   441 			if ( greater(ts, tx, lgth_x) || underflow ) {

   442 				bmult[0]--;

   443 				underflow = 1;

   444 			}

   445 			else

   446 				underflow = 0;

   447 		} while ( underflow );

   448 		sub(tx, lgth_x, ts, lgth_x);

   449 		if ( shift )

   450 			quot[lenq - (lgth_x - lgth_n) - 1] = bmult[0];

   451 		else

   452 			quot[lenq - (lgth_x - lgth_n)] = bmult[0];

   453 		

   454 		ttx = tx;

   455 		t_len = lgth_x;

   456 		for ( i=0, lgth_x=t_len; i<t_len; i++, lgth_x-- ) {

   457 			if ( *ttx )

   458 				break;

   459 			ttx++;

   460 		}

   461 		tx = ttx;

   462 	}

   463 	memset(rem, 0x00, lenn);

   464 	if ( lgth_x )

   465 		memcpy(rem+lenn-lgth_x, tx, lgth_x);

   466 

   467 	return lenq;

   468 }

   469 

   470 

   471 /* 

   472  * Mod - Computes an integer modulo another integer

   473  *

   474  * x = x (mod n)

   475  *

   476  */

   477 void Mod(BYTE *x, int lenx, BYTE *n, int lenn)

   478 {

   479 	BYTE	quot[MAXPLEN+1], rem[2*MAXPLEN+1];

   480 

   481 	memset(quot, 0x00, sizeof(quot));

   482 	memset(rem, 0x00, sizeof(rem));

   483 	if ( DivMod(x, lenx, n, lenn, quot, rem) ) {

   484 		memset(x, 0x00, lenx);

   485 		memcpy(x+lenx-lenn, rem, lenn);

   486 	}

   487 }

   488 

   489 /* 

   490  * Div - Computes the integer division of two numbers

   491  *

   492  * x = x / n

   493  *

   494  */

   495 void Div(BYTE *x, int lenx, BYTE *n, int lenn)

   496 {

   497 	BYTE	quot[MAXPLEN+1], rem[2*MAXPLEN+1];

   498 	int		lenq;

   499 

   500 	memset(quot, 0x00, sizeof(quot));

   501 	memset(rem, 0x00, sizeof(rem));

   502 	if ( (lenq = DivMod(x, lenx, n, lenn, quot, rem)) != 0 ) {

   503 		memset(x, 0x00, lenx);

   504 		memcpy(x+lenx-lenq, quot, lenq);

   505 	}

   506 }

   507 

   508 

   509 /*****************************************

   510 ** sub - Subtract two integers           *

   511 **                                       *

   512 ** A = A - B                             *

   513 **                                       *

   514 **                                       *

   515 ** Parameters:                           *	

   516 **                                       *

   517 **  A      Address of subtrahend integer *

   518 **  B      Address of subtractor integer *

   519 **  L      Length of A and B in bytes    *

   520 **                                       *

   521 **  NOTE: In order to save RAM, B is     *

   522 **        two's complemented twice,      *

   523 **        rather than using a copy of B  *

   524 **                                       *

   525 ******************************************/

   526 void sub(BYTE *A, int LA, BYTE *B, int LB)

   527 {

   528 	BYTE	*tb;

   529 

   530 	tb = (BYTE *)calloc(LA, 1);

   531 	memcpy(tb, B, LB);

   532 	negate(tb, LB);

   533 	add(A, LA, tb, LA);

   534 

   535 	FREE(tb);

   536 }

   537 

   538 

   539 /*****************************************

   540 ** negate - Negate an integer            *

   541 **                                       *

   542 ** A = -A                                *

   543 **                                       *

   544 **                                       *

   545 ** Parameters:                           *	

   546 **                                       *

   547 **  A      Address of integer to negate  *

   548 **  L      Length of A in bytes          *

   549 **                                       *

   550 ******************************************/

   551 int negate(BYTE *A, int L)

   552 {

   553 	int		i, tL;

   554 	DIGIT	accum;

   555 

   556 	/* Take one's complement of A */

   557 	for ( i=0; i<L; i++ )

   558 		A[i] = (BYTE)(~(A[i]));

   559 

   560 	/* Add one to get two's complement of A */

   561 	accum = 1;

   562 	tL = L-1;

   563 	while ( accum && (tL >= 0) ) {

   564 		accum = (DIGIT)(accum + A[tL]);

   565 		A[tL--] = (BYTE)(accum & 0xff);

   566 		accum = (DIGIT)(accum >> 8);

   567 	}

   568 

   569 	return accum;

   570 }

   571 

   572 

   573 /*

   574  * add()

   575  *

   576  * A = A + B

   577  *

   578  * LB must be <= LA

   579  *

   580  */

   581 BYTE add(BYTE *A, int LA, BYTE *B, int LB)

   582 {

   583 	int		i, indexA, indexB;

   584 	DIGIT	accum;

   585 

   586 	indexA = LA - 1; 	/* LSD of result */

   587 	indexB = LB - 1; 	/* LSD of B */

   588 

   589 	accum = 0;

   590 	for ( i = 0; i < LB; i++ ) {

   591 		accum = (DIGIT)(accum + A[indexA]);

   592 		accum = (DIGIT)(accum + B[indexB--]);

   593 		A[indexA--] = (BYTE)(accum & 0xff);

   594 		accum = (DIGIT)(accum >> 8);

   595 	}

   596 

   597 	if ( LA > LB )

   598 		while ( accum  && (indexA >= 0) ) {

   599 			accum = (DIGIT)(accum + A[indexA]);

   600 			A[indexA--] = (BYTE)(accum & 0xff);

   601 			accum = (DIGIT)(accum >> 8);

   602 		}

   603 

   604 	return (BYTE)accum;

   605 }

   606 

   607 

   608 void prettyprintBstr(char *S, BYTE *A, int L)

   609 {

   610 	int		i, extra, ctrb, ctrl;

   611 

   612 	if ( L == 0 )

   613 		printf("%s <empty>", S);

   614 	else

   615 		printf("%s\n\t", S);

   616 	extra = L % 24;

   617 	if ( extra ) {

   618 		ctrb = 0;

   619 		for ( i=0; i<24-extra; i++ ) {

   620 			printf("  ");

   621 			if ( ++ctrb == 4) {

   622 				printf(" ");

   623 				ctrb = 0;

   624 			}

   625 		}

   626 

   627 		for ( i=0; i<extra; i++ ) {

   628 			printf("%02X", A[i]);

   629 			if ( ++ctrb == 4) {

   630 				printf(" ");

   631 				ctrb = 0;

   632 			}

   633 		}

   634 		printf("\n\t");

   635 	}

   636 

   637 	ctrb = ctrl = 0;

   638 	for ( i=extra; i<L; i++ ) {

   639 		printf("%02X", A[i]);

   640 		if ( ++ctrb == 4) {

   641 			ctrl++;

   642 			if ( ctrl == 6 ) {

   643 				printf("\n\t");

   644 				ctrl = 0;

   645 			}

   646 			else

   647 				printf(" ");

   648 			ctrb = 0;

   649 		}

   650 	}

   651 	printf("\n\n");

   652 }

   653 

   654 

   655 /**********************************************************************/

   656 /*  Performs byte reverse for PC based implementation (little endian) */

   657 /**********************************************************************/

   658 void byteReverse(ULONG *buffer, int byteCount)

   659 {

   660 	ULONG value;

   661 	int count;

   662 

   663 	byteCount /= sizeof( ULONG );

   664 	for( count = 0; count < byteCount; count++ ) {

   665 		value = ( buffer[ count ] << 16 ) | ( buffer[ count ] >> 16 );

   666 		buffer[ count ] = ( ( value & 0xFF00FF00L ) >> 8 ) | ( ( value & 0x00FF00FFL ) << 8 );

   667 	}

   668 }

   669 

   670 void

   671 ahtopb (char *ascii_hex, BYTE *p_binary, int bin_len)

   672 {

   673 	BYTE    nibble;

   674 	int     i; 

   675 	

   676 	for ( i=0; i<bin_len; i++ ) {

   677         nibble = ascii_hex[i * 2];

   678 	    if ( nibble > 'F' )

   679 	        nibble -= 0x20;   

   680 	    if ( nibble > '9' )

   681 	        nibble -= 7;      

   682 	    nibble -= '0';   

   683 	    p_binary[i] = (BYTE)(nibble << 4);

   684 		

   685 	    nibble = ascii_hex[i * 2 + 1];

   686 	    if ( nibble > 'F' )

   687 			nibble -= 0x20;

   688         if ( nibble > '9' )

   689             nibble -= 7;   

   690         nibble -= '0';

   691 		p_binary[i] = (BYTE)(p_binary[i] + nibble);

   692 	}

   693 }