crypto/weakcrypto/source/bigint/primes.cpp
changeset 72 de46a57f75fb
equal deleted inserted replaced
65:970c0057d9bc 72:de46a57f75fb
       
     1 /*
       
     2 * Copyright (c) 2003-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 the License "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 *
       
    16 */
       
    17 
       
    18 
       
    19 #include <bigint.h>
       
    20 #include <e32std.h>
       
    21 #include <securityerr.h>
       
    22 #include "words.h"
       
    23 #include "primes.h"
       
    24 #include "algorithms.h"
       
    25 #include "mont.h"
       
    26 #include "stackinteger.h"
       
    27 
       
    28 static inline void EliminateComposites(TUint* aS, TUint aPrime, TUint aJ, 
       
    29 	TUint aMaxIndex)
       
    30 	{
       
    31 	for(; aJ<aMaxIndex; aJ+=aPrime)
       
    32 		ArraySetBit(aS, aJ);
       
    33 	}
       
    34 
       
    35 static inline TInt FindLeastSignificantZero(TUint aX)
       
    36 	{
       
    37 	aX = ~aX;
       
    38 	int i = 0;
       
    39 	if( aX << 16 == 0 ) aX>>=16, i+=16;
       
    40 	if( aX << 24 == 0 ) aX>>=8, i+=8;
       
    41 	if( aX << 28 == 0 ) aX>>=4, i+=4;
       
    42 	if( aX << 30 == 0 ) aX>>=2, i+=2;
       
    43 	if( aX << 31 == 0 ) ++i;
       
    44 	return i;
       
    45 	}
       
    46 
       
    47 static inline TInt FindFirstPrimeCandidate(TUint* aS, TUint aBitLength)
       
    48 	{
       
    49 	assert(aBitLength % WORD_BITS == 0);
       
    50 	TUint i=0;
       
    51 	//The empty statement at the end of this is stop warnings in all compilers
       
    52 	for(; aS[i] == KMaxTUint && i<BitsToWords(aBitLength); i++) {;}
       
    53 
       
    54 	if(i == BitsToWords(aBitLength))
       
    55 		return -1;
       
    56 	else
       
    57 		{
       
    58 		assert( FindLeastSignificantZero((TUint)(aS[i])) >= 0 );
       
    59 		assert( FindLeastSignificantZero((TUint)(aS[i])) <= 31 );
       
    60 		return i*WORD_BITS + FindLeastSignificantZero((TUint32)(aS[i]));
       
    61 		}
       
    62 	}
       
    63 
       
    64 static inline TUint FindSmallestIndex(TUint aPrime, TUint aRemainder)
       
    65 	{
       
    66 	TUint& j = aRemainder;
       
    67 	if(j)
       
    68 		{
       
    69 		j = aPrime - aRemainder;
       
    70 		if( j & 0x1L )
       
    71 			{
       
    72 			//if j is odd then this + j is even so we actually want 
       
    73 			//the next number for which (this + j % p == 0) st this + j is odd
       
    74 			//that is: this + j + p == 0 mod p
       
    75 			j += aPrime;
       
    76 			}
       
    77 		//Turn j into an index for a bit array representing odd numbers only
       
    78 		j>>=1;
       
    79 		}
       
    80 	return j;
       
    81 	}
       
    82 
       
    83 static TBool IsSmallPrime(TUint aK) 
       
    84 	{
       
    85 	//This is just a binary search of our small prime table.
       
    86 	TUint l = 0;
       
    87 	TUint u = KPrimeTableSize;
       
    88 	while( u > l )
       
    89 		{
       
    90 		TUint m = (l+u)>>1;
       
    91 		TUint p = KPrimeTable[m];
       
    92 		if(aK < p)
       
    93 			u = m;
       
    94 		else if (aK > p)
       
    95 			l = m + 1;
       
    96 		else
       
    97 			return ETrue;
       
    98 		}
       
    99 	return EFalse;
       
   100 	}
       
   101 
       
   102 static inline TUint RabinMillerRounds(TUint aBits) 
       
   103 	{
       
   104 	//See HAC Table 4.4
       
   105 	if(aBits > 1300)
       
   106 		return 2;
       
   107 	if (aBits > 850)
       
   108 		return 3;
       
   109 	if (aBits > 650)
       
   110 		return 4;
       
   111 	if (aBits > 550)
       
   112 		return 5;
       
   113 	if (aBits > 450)
       
   114 		return 6;
       
   115 	if (aBits > 400)
       
   116 		return 7;
       
   117 	if (aBits > 350)
       
   118 		return 8;
       
   119 	if (aBits > 300)
       
   120 		return 9;
       
   121 	if (aBits > 250)
       
   122 		return 12;
       
   123 	if (aBits > 200)
       
   124 		return 15;
       
   125 	if (aBits > 150)
       
   126 		return 18;
       
   127 	if (aBits > 100)
       
   128 		return 27;
       
   129 	//All of the above are optimisations on the worst case.  The worst case
       
   130 	//chance of odd composite integers being declared prime by Rabin-Miller is
       
   131 	//(1/4)^t where t is the number of rounds.  Thus, t = 40 means that the
       
   132 	//chance of declaring a composite integer prime is less than 2^(-80).  See
       
   133 	//HAC Fact 4.25 and most of chapter 4 for more details.
       
   134 	return 40;
       
   135 	}
       
   136 
       
   137 static TBool HasSmallDivisorL(const TInteger& aPossiblePrime)
       
   138 	{
       
   139 	assert(aPossiblePrime.IsOdd());
       
   140 	//Start checking at the first odd prime, whether it is even should have
       
   141 	//already been checked
       
   142 	for( TUint i=1; i < KPrimeTableSize; i++ )
       
   143 		{
       
   144 		if( aPossiblePrime.ModuloL(KPrimeTable[i]) == 0 )
       
   145 			{
       
   146 			return ETrue;
       
   147 			}
       
   148 		}
       
   149 	return EFalse;
       
   150 	}
       
   151 
       
   152 static TBool RabinMillerIterationL(const CMontgomeryStructure& aMont, 
       
   153 	const TInteger& aProbablePrime, const TInteger& aBase)
       
   154 	{
       
   155 	//see HAC 4.24
       
   156 	const TInteger& n = aProbablePrime;
       
   157 	assert(n > KLastSmallPrimeSquared);
       
   158 	assert(n.IsOdd());
       
   159 	assert(aBase > TInteger::One());
       
   160 
       
   161 	RInteger nminus1 = n.MinusL(TInteger::One());
       
   162 	CleanupStack::PushL(nminus1);
       
   163 	assert(aBase < nminus1);
       
   164 
       
   165 	// 1) find (s | 2^s*r == n-1) where r is odd
       
   166 	// we want the largest power of 2 that divides n-1
       
   167 	TUint s=0;
       
   168 	for(;;s++)
       
   169 		{
       
   170 		if(nminus1.Bit(s))
       
   171 			{
       
   172 			break;
       
   173 			}
       
   174 		}
       
   175 	// (r = (n-1) / 2^s) which is equiv to (n-1 >>= s)
       
   176 	RInteger r = RInteger::NewL(nminus1);
       
   177 	CleanupStack::PushL(r);
       
   178 	r >>= s;
       
   179 
       
   180 	//At no point do we own y, aMont owns it
       
   181 	const TInteger* y = &(aMont.ExponentiateL(aBase, r));
       
   182 
       
   183 	TBool probablePrime = EFalse;
       
   184 	
       
   185 	TUint j=1;
       
   186 	if( *y == TInteger::One() || *y == nminus1 )
       
   187 		{
       
   188 		probablePrime = ETrue;
       
   189 		}
       
   190 	else
       
   191 		{
       
   192 		for(j=1; j<s; j++)
       
   193 			{
       
   194 			y = &(aMont.SquareL(*y));
       
   195 			if(*y == nminus1)
       
   196 				{
       
   197 				probablePrime = ETrue;
       
   198 				break;
       
   199 				}
       
   200 			if(*y == TInteger::One())
       
   201 				{
       
   202 				assert(probablePrime == EFalse);
       
   203 				break;
       
   204 				}
       
   205 			}
       
   206 		}
       
   207 	CleanupStack::PopAndDestroy(&r);
       
   208 	CleanupStack::PopAndDestroy(&nminus1);//y,r,nminus1
       
   209 	return probablePrime;
       
   210 	}
       
   211 
       
   212 static TBool RabinMillerTestL(const CMontgomeryStructure& aMont, 
       
   213 	const TInteger& aProbablePrime, TUint aRounds) 
       
   214 	{
       
   215 	const TInteger& n = aProbablePrime;
       
   216 	assert(n > KLastSmallPrimeSquared);
       
   217 	
       
   218 	RInteger nminus2 = n.MinusL(TInteger::Two());
       
   219 	CleanupStack::PushL(nminus2);
       
   220 
       
   221 	for(TUint i=0; i<aRounds; i++)
       
   222 		{
       
   223 		RInteger base = RInteger::NewRandomL(TInteger::Two(), nminus2);
       
   224 		CleanupStack::PushL(base);
       
   225 		if(!RabinMillerIterationL(aMont, n, base))
       
   226 			{
       
   227 			CleanupStack::PopAndDestroy(2, &nminus2);//base, nminus2
       
   228 			return EFalse;
       
   229 			}
       
   230 		CleanupStack::PopAndDestroy(&base);
       
   231 		}
       
   232 	CleanupStack::PopAndDestroy(&nminus2);
       
   233 	return ETrue;
       
   234 	}
       
   235 
       
   236 static TBool IsStrongProbablePrimeL(const TInteger& aPrime) 
       
   237 	{
       
   238 	CMontgomeryStructure* mont = CMontgomeryStructure::NewLC(aPrime);
       
   239 	//This should be using short circuit evaluation
       
   240 	TBool probablePrime = RabinMillerIterationL(*mont, aPrime, TInteger::Two())
       
   241 		&& RabinMillerTestL(*mont, aPrime,RabinMillerRounds(aPrime.BitCount()));
       
   242 	CleanupStack::PopAndDestroy(mont);
       
   243 	return probablePrime;
       
   244 	}
       
   245 
       
   246 /* In the _vast_ majority of cases this simply checks that your chosen random
       
   247  * number is >= KLastSmallPrimeSquared and return EFalse and lets the normal
       
   248  * prime generation routines handle the situation.  In the case where it is
       
   249  * smaller, it generates a provable prime and returns ETrue.  The algorithm for
       
   250  * finding a provable prime < KLastPrimeSquared is not the most efficient in the
       
   251  * world, but two points come to mind
       
   252  * 1) The two if statements hardly _ever_ evaluate to ETrue in real life.
       
   253  * 2) Even when it is, the distribution of primes < KLastPrimeSquared is pretty
       
   254  * dense, so you aren't going to have check many.
       
   255  * This function is essentially here for two reasons:
       
   256  * 1) Ensures that it is possible to generate primes < KLastPrimeSquared (the
       
   257  * test code does this)
       
   258  * 2) Ensures that if you request a prime of a large bit size that there is an
       
   259  * even probability distribution across all integers < 2^aBits
       
   260  */
       
   261 TBool TInteger::SmallPrimeRandomizeL(void)
       
   262 	{
       
   263 	TBool foundPrime = EFalse;
       
   264 	//If the random number we've chosen is less than KLastSmallPrime,
       
   265 	//testing for primality is easy.
       
   266 	if(*this <= KLastSmallPrime)
       
   267 		{
       
   268 		//If Zero or One, or two, next prime number is two
       
   269 		if(IsZero() || *this == One() || *this == Two())
       
   270 			{
       
   271 			CopyL(TInteger::Two());
       
   272 			foundPrime = ETrue;
       
   273 			}
       
   274 		else
       
   275 			{
       
   276 			//Make sure any number we bother testing is at least odd
       
   277 			SetBit(0);
       
   278 			//Binary search the small primes.
       
   279 			while(!IsSmallPrime(ConvertToUnsignedLong()))
       
   280 				{
       
   281 				//If not prime, add two and try the next odd number.
       
   282 
       
   283 				//will never carry as the minimum size of an RInteger is 2
       
   284 				//words.  Much bigger than KLastSmallPrime on 32bit
       
   285 				//architectures.
       
   286 				IncrementNoCarry(Ptr(), Size(), 2);
       
   287 				}
       
   288 			assert(IsSmallPrime(ConvertToUnsignedLong()));
       
   289 			foundPrime = ETrue;
       
   290 			}
       
   291 		}
       
   292 	else if(*this <= KLastSmallPrimeSquared)
       
   293 		{
       
   294 		//Make sure any number we bother testing is at least odd
       
   295 		SetBit(0);
       
   296 
       
   297 		while(HasSmallDivisorL(*this) && *this <= KLastSmallPrimeSquared)
       
   298 			{
       
   299 			//If not prime, add two and try the next odd number.
       
   300 
       
   301 			//will never carry as the minimum size of an RInteger is 2
       
   302 			//words.  Much bigger than KLastSmallPrime on 32bit
       
   303 			//architectures.
       
   304 			IncrementNoCarry(Ptr(), Size(), 2);
       
   305 			}
       
   306 		//If we exited while loop because it had no small divisor then it is
       
   307 		//prime.  Otherwise, we've exceeded the limit of what we can provably
       
   308 		//generate.  Therefore the normal prime gen routines will be run on it
       
   309 		//now.
       
   310 		if(*this < KLastSmallPrimeSquared)
       
   311 			{
       
   312 			foundPrime = ETrue;
       
   313 			}
       
   314 		else
       
   315 			{
       
   316 			assert(foundPrime == EFalse);
       
   317 			}
       
   318 		}
       
   319 	//This doesn't mean there is no such prime, simply means that the number
       
   320 	//wasn't less than KSmallPrimeSquared and needs to be handled by the normal
       
   321 	//prime generation routines.
       
   322 	return foundPrime;
       
   323 	}
       
   324 
       
   325 void TInteger::PrimeRandomizeL(TUint aBits, TRandomAttribute aAttr)
       
   326 	{
       
   327 	assert(aBits > 1); 
       
   328 	
       
   329 	//"this" is "empty" currently.  Consists of Size() words of 0's.  This is just
       
   330 	//checking that sign flag is positive as we don't set it later.
       
   331 	assert(NotNegative());
       
   332 
       
   333 	//Flag for the whole function saying if we've found a prime
       
   334 	TBool foundProbablePrime = EFalse;
       
   335 
       
   336 	//Find 2^aBits + 1 -- any prime we find must be less than this.
       
   337 	RInteger max = RInteger::NewEmptyL(BitsToWords(aBits)+1);
       
   338 	CleanupStack::PushL(max);
       
   339 	max.SetBit(aBits);
       
   340 	assert(max.BitCount()-1 == aBits);
       
   341 
       
   342 	// aBits 	| approx number of odd numbers you must try to have a 50% 
       
   343 	//			chance of finding a prime
       
   344 	//---------------------------------------------------------
       
   345 	// 512		| 122		
       
   346 	// 1024		| 245
       
   347 	// 2048		| 1023
       
   348 	//Therefore if we are generating larger than 1024 bit numbers we'll use a
       
   349 	//bigger bit array to have a better chance of avoiding re-generating it.
       
   350 	TUint sLength = aBits > 1024 ? 1024 : 512;
       
   351 	RInteger S = RInteger::NewEmptyL(BitsToWords(sLength));
       
   352 	CleanupStack::PushL(S);
       
   353 
       
   354 	while(!foundProbablePrime)
       
   355 		{
       
   356 		//Randomly choose aBits
       
   357 		RandomizeL(aBits, aAttr);
       
   358 
       
   359 		//If the random number chosen is less than KSmallPrimeSquared, we have a
       
   360 		//special set of routines.
       
   361 		if(SmallPrimeRandomizeL())
       
   362 			{
       
   363 			foundProbablePrime = ETrue;
       
   364 			}
       
   365 		else
       
   366 			{
       
   367 			//if it was <= KLastSmallPrimeSquared then it would have been
       
   368 			//handled by SmallPrimeRandomizeL()
       
   369 			assert(*this > KLastSmallPrimeSquared);
       
   370 
       
   371 			//Make sure any number we bother testing is at least odd
       
   372 			SetBit(0);
       
   373 
       
   374 			//Ensure that this + 2*sLength < max
       
   375 			RInteger temp = max.MinusL(*this);
       
   376 			CleanupStack::PushL(temp);
       
   377 			++temp;
       
   378 			temp >>=1;
       
   379 			if(temp < sLength)
       
   380 				{
       
   381 				//if this + 2*sLength >= max then we use a smaller sLength to
       
   382 				//ensure we don't find a number that is outside of our bounds
       
   383 				//(and bigger than our allocated memory for this)
       
   384 
       
   385 				//temp must be less than KMaxTUint as sLength is a TUint 
       
   386 				sLength = temp.ConvertToUnsignedLong();	
       
   387 				}
       
   388 			CleanupStack::PopAndDestroy(&temp);
       
   389 
       
   390 			//Start at 1 as no point in checking against 2 (all odd numbers)
       
   391 			for(TUint i=1; i<KPrimeTableSize; i++)
       
   392 				{
       
   393 				//no need to call ModuloL as we know KPrimeTable[i] is not 0
       
   394 				TUint remainder = Modulo(*this, KPrimeTable[i]);
       
   395 				TUint index = FindSmallestIndex(KPrimeTable[i], remainder);
       
   396 				EliminateComposites(S.Ptr(), KPrimeTable[i], index, sLength);
       
   397 				}
       
   398 			TInt j = FindFirstPrimeCandidate(S.Ptr(), sLength);
       
   399 			TInt prev = 0;
       
   400 			for(; j>=0; j=FindFirstPrimeCandidate(S.Ptr(), sLength))
       
   401 				{
       
   402 				ArraySetBit(S.Ptr(), j);
       
   403 
       
   404 				//should never carry as we earlier made sure that 2*j + this < max
       
   405 				//where max is 1 bit more than we asked for.
       
   406 				IncrementNoCarry(Ptr(), Size(), 2*(j-prev));
       
   407 
       
   408 				assert(*this < max);
       
   409 				assert(!HasSmallDivisorL(*this));
       
   410 
       
   411 				prev = j;
       
   412 
       
   413 				if( IsStrongProbablePrimeL(*this) )
       
   414 					{
       
   415 					foundProbablePrime = ETrue;
       
   416 					break;
       
   417 					}
       
   418 				}
       
   419 			//This clears the memory
       
   420 			S.CopyL(0, EFalse);
       
   421 			}
       
   422 		}
       
   423 	CleanupStack::PopAndDestroy(2, &max);
       
   424 	}
       
   425 
       
   426 TBool TInteger::IsPrimeL(void) const
       
   427 	{
       
   428 	if( NotPositive() )
       
   429 		{
       
   430 		return EFalse;
       
   431 		}
       
   432 	else if( IsEven() )
       
   433 		{
       
   434 		return *this == Two();
       
   435 		}
       
   436 	else if( *this <= KLastSmallPrime )
       
   437 		{
       
   438 		assert(KLastSmallPrime < KMaxTUint);
       
   439 		return IsSmallPrime(this->ConvertToUnsignedLong());
       
   440 		}
       
   441 	else if( *this <= KLastSmallPrimeSquared )
       
   442 		{
       
   443 		return !HasSmallDivisorL(*this);
       
   444 		}
       
   445 	else 
       
   446 		{
       
   447 		return !HasSmallDivisorL(*this) && IsStrongProbablePrimeL(*this);
       
   448 		}
       
   449 	}