//
//  BigInt.m
//
//  Created by Josh Santomieri on 5/17/09.
//
// This is a port of C# code originally written by Chew Keong TAN
// See http://www.codeproject.com/csharp/biginteger.asp for more details.
// 
// Objective-C Big Integer Library

#import "BigInt.h"

#define _BI_DEBUG_ 0

static int primesBelow2000[] = {
2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97,
101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199,
211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293,
307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397,
401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499,
503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599,
601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691,
701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797,
809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887,
907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997,
1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097,
1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181, 1187, 1193,
1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291, 1297,
1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, 1381, 1399,
1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499,
1511, 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, 1597,
1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657, 1663, 1667, 1669, 1693, 1697, 1699,
1709, 1721, 1723, 1733, 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789,
1801, 1811, 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889,
1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997, 1999 };

@implementation BigInt

@synthesize dataLength;

-(id)init {
	if(self = [super init]) {
		dataLength = 1;
	}
	return self;
}

-(id)initWithInt:(int)value {
	return [[[BigInt alloc] initWithLong:(long)value] autorelease];
}

-(id)initWithLong:(long)value {
	if(self= [super init]) {
		
		bzero(data, sizeof(data));
		
		long long tempVal = (long long)value;
		long long tmp2 = (long long)value;
		
		// copy bytes from long to BigInteger without any assumption of
		// the length of the long datatype
		
		dataLength = 0;
		while (tmp2 != 0 && dataLength < MAX_LENGTH) {
			data[dataLength] = (uint)(tmp2 & 0xFFFFFFFF);
			tmp2 >>= 32;
			dataLength++;
		}
		
		if (tempVal > 0)         // overflow check for +ve value
		{
			if (tmp2 != 0 || (data[MAX_LENGTH - 1] & 0x80000000) != 0)
				@throw [NSException exceptionWithName:@"ArithmeticException" reason:@"Positive overflow in constructor." userInfo:nil];
		}
		else if (tempVal < 0)    // underflow check for -ve value
		{
			if (tmp2 != -1 || (data[MAX_LENGTH - 1] & 0x80000000) == 0)
				@throw [NSException exceptionWithName:@"ArithmeticException" reason:@"Negative underflow in constructor." userInfo:nil];
			
		}
		
		if (dataLength == 0)
			dataLength = 1;
	}
	return self;
}

-(id)initWithULong:(ulong)value {
	if(self = [super init]) {
		
		bzero(data, sizeof(data));
		
		ulong tempVal = (ulong)value;
		ulong tmp2 = (ulong)value;
		
#if _BI_DEBUG_
		NSLog(@"initWithULong:value = %qi", value);
#endif
		
		dataLength = 0;
		while (tmp2 != 0 && dataLength < MAX_LENGTH) {
			data[dataLength] = (uint)(tmp2 & 0xFFFFFFFF);
			tmp2 >>= 32;
			dataLength++;
		}
		
		if (tempVal > 0)         // overflow check for +ve value
		{
			if (tmp2 != 0 || (data[MAX_LENGTH - 1] & 0x80000000) != 0)
				@throw [NSException exceptionWithName:@"ArithmeticException" reason:@"Positive overflow in constructor." userInfo:nil];
		}
		else if (tempVal < 0)    // underflow check for -ve value
		{
			if (tmp2 != -1 || (data[MAX_LENGTH - 1] & 0x80000000) == 0)
				@throw [NSException exceptionWithName:@"ArithmeticException" reason:@"Negative underflow in constructor." userInfo:nil];
		}
		
		if (dataLength == 0)
			dataLength = 1;
	}
	return self;
}

-(id)initWithBigInt:(BigInt *)value {
	if(self = [super init]) {
		bzero(data, sizeof(data));
		
		dataLength = value.dataLength;		
		for (int i = 0; i < dataLength; i++)
			data[i] = [value getDataAtIndex:i];
	}
	return self;
}

-(id)initWithUIntArray:(uint *)value withSize:(int)size{
	if(self = [super init]) {
		bzero(data, sizeof(data));
		
		dataLength = size;
		
		if (dataLength > MAX_LENGTH)
			@throw [NSException exceptionWithName:@"ArithmaticException" reason:@"Byte overflow in constructor" userInfo:nil];
		
		for (int i = dataLength - 1, j = 0; i >= 0; i--, j++) {
			data[j] = (uint)(value[i]);
		}
		
		while (dataLength > 1 && data[dataLength - 1] == 0)
			dataLength--;
		
	}
	return self;
}

-(id)initWithString:(NSString *)value andRadix:(int)radix {
	if(self = [super init]) {
		bzero(data, sizeof(data));
		
		BigInt *multiplier = [BigInt createFromLong:1];
		BigInt *result = [BigInt create];
		
		value = [value uppercaseString];
		int limit = 0;
		
		if ([value characterAtIndex:0] == '-')
			limit = 1;
		
		for (int i = [value length] - 1; i >= limit; i--) {
			int posVal = (int)[value characterAtIndex:i];
			
			if (posVal >= '0' && posVal <= '9')
				posVal -= '0';
			else if (posVal >= 'A' && posVal <= 'Z')
				posVal = (posVal - 'A') + 10;
			else
				posVal = 9999999;       // arbitrary large
			
			if (posVal >= radix)
				@throw [NSException exceptionWithName:@"ArithmeticException" reason:@"Invalid string in constructor." userInfo:nil];
			else {
				if ([value characterAtIndex:0] == '-')
					posVal = -posVal;
				
				BigInt *mult = [multiplier multiply:[BigInt createFromInt:posVal]];
				result = [result add:mult];
				
				if ((i - 1) >= limit)
					multiplier = [multiplier multiply:[BigInt createFromInt:radix]];
			}
		}
		
		if ([value characterAtIndex:0] == '-')     // negative values
		{
			if (([result getDataAtIndex:(MAX_LENGTH - 1)] & 0x80000000) == 0)
				@throw [NSException exceptionWithName:@"ArithmeticException" reason:@"Negative underflow in constructor." userInfo:nil];
		}
		else    // positive values
		{
			if (([result getDataAtIndex:(MAX_LENGTH - 1)] & 0x80000000) != 0)
				@throw [NSException exceptionWithName:@"ArithmeticException" reason:@"Positive overflow in constructor." userInfo:nil];
		}
		
		for (int i = 0; i < result.dataLength; i++)
			[self setData:[result getDataAtIndex:i] atIndex:i];
		
		dataLength = result.dataLength;	
	}
	return self;
}

-(void)dealloc {
	
	[super dealloc];
}

+(BigInt *)create {
	return [BigInt createFromLong:0];
}

+(BigInt *)createFromInt:(int)value {
	return [[[BigInt alloc] initWithLong:(long)value] autorelease];
}

+(BigInt *)createFromLong:(long)value {
	return [[[BigInt alloc] initWithLong:value] autorelease];
}

+(BigInt *)createFromULong:(ulong)value {
	return [[[BigInt alloc] initWithULong:value] autorelease];
}

+(BigInt *)createFromBigInt:(BigInt *)value {
	return [[[BigInt alloc] initWithBigInt:value] autorelease];
}

+(BigInt *)createFromString:(NSString *)value andRadix:(int)radix {
	return [[[BigInt alloc] initWithString:value andRadix:radix] autorelease];	
}

-(uint *)getData {
	return data;
}

-(uint)getDataAtIndex:(int)index {
	return data[index];
}

-(void)setData:(uint)value atIndex:(int)index {
	data[index] = value;
}

-(BigInt *)add:(BigInt *)bi2 {
#if _BI_DEBUG_
	NSLog(@"add");
#endif
	
	BigInt *result = [BigInt create];
	
	result.dataLength = (self.dataLength > bi2.dataLength) ? self.dataLength : bi2.dataLength;
	
	ulong carry = 0;
	for (int i = 0; i < result.dataLength; i++) {
		ulong sum = (ulong)[self getDataAtIndex:i] + (ulong)[bi2 getDataAtIndex:i] + carry;
		carry = sum >> 32;
		[result setData:(uint)(sum & 0xFFFFFFFF) atIndex:i];
	}
	
	if (carry != 0 && result.dataLength < MAX_LENGTH) {
		[result setData:(uint)carry atIndex:result.dataLength];
		result.dataLength++;
	}
	
	while (result.dataLength > 1 && [result getDataAtIndex:(result.dataLength - 1)] == 0)
		result.dataLength--;
	
	// overflow check
	int lastPos = MAX_LENGTH - 1;
	if (([self getDataAtIndex:lastPos] & 0x80000000) == ([bi2 getDataAtIndex:lastPos] & 0x80000000) &&
		([result getDataAtIndex:lastPos] & 0x80000000) != ([self getDataAtIndex:lastPos] & 0x80000000)) {
		@throw [NSException exceptionWithName:@"ArithmeticException" reason:@"Overflow" userInfo:nil];
	}
	
	return result;
}

-(BigInt *)subtract:(BigInt *)bi2 {
#if _BI_DEBUG_
	NSLog(@"subtract");
#endif
	
	BigInt *result = [[BigInt alloc] init];
	
	result.dataLength = (self.dataLength > bi2.dataLength) ? self.dataLength : bi2.dataLength;
	
	long carryIn = 0;
	for (int i = 0; i < result.dataLength; i++) {
		long long diff = (((long long)[self getDataAtIndex:i] - (long long)[bi2 getDataAtIndex:i]) - carryIn);
		[result setData:(uint)(diff & 0xFFFFFFFF) atIndex:i];
		
		if (diff < 0)
			carryIn = 1;
		else
			carryIn = 0;
	}
	
	// roll over to negative
	if (carryIn != 0) {
		for (int i = result.dataLength; i < MAX_LENGTH; i++) {
			[result setData:0xFFFFFFFF atIndex:i];
		}
		result.dataLength = MAX_LENGTH;
	}
	
	// fixed in v1.03 to give correct datalength for a - (-b)
	while (result.dataLength > 1 && [result getDataAtIndex:(result.dataLength - 1)] == 0)
		result.dataLength--;
	
	// overflow check
	
	int lastPos = MAX_LENGTH - 1;
	if (([self getDataAtIndex:lastPos] & 0x80000000) != ([bi2 getDataAtIndex:lastPos] & 0x80000000) &&
		([result getDataAtIndex:lastPos] & 0x80000000) != ([self getDataAtIndex:lastPos] & 0x80000000)) {
		@throw [NSException exceptionWithName:@"ArithmeticException" reason:@"Overflow" userInfo:nil];
	}
	
	return [result autorelease];
}

-(BigInt *)multiply:(BigInt *)bi2 {
#if _BI_DEBUG_	
	NSLog(@"multiply");
#endif
	
	NSAutoreleasePool *pool = [[NSAutoreleasePool alloc] init];
	
	BigInt *result = [BigInt create];
	int lastPos = MAX_LENGTH - 1;
	bool bi1Neg = false, bi2Neg = false;
	BigInt *bi1 = [BigInt createFromBigInt:self];
	
	// take the absolute value of the inputs
	@try {
		if (([bi1 getDataAtIndex:lastPos] & 0x80000000) != 0)     // bi1 negative
		{
			bi1Neg = true; bi1 = [bi1 negate];
		}
		if (([bi2 getDataAtIndex:lastPos] & 0x80000000) != 0)     // bi2 negative
		{
			bi2Neg = true; bi2 = [bi2 negate];
		}
	}
	@catch (NSException *ex) { }
	
	// multiply the absolute values
	@try {
		for (int i = 0; i < bi1.dataLength; i++) {
			if ([bi1 getDataAtIndex:i] == 0) continue;
			
			ulong mcarry = 0;
			for (int j = 0, k = i; j < bi2.dataLength; j++, k++) {
				// k = i + j
				
				ulong bi1_val = (ulong)[bi1 getDataAtIndex:i];
				ulong bi2_val = (ulong)[bi2 getDataAtIndex:j];
				ulong res_val = (ulong)[result getDataAtIndex:k];
				
#if _BI_DEBUG_
				NSLog(@"bi1_val = %qi", bi1_val);
#endif
				
				ulong val = (bi1_val * bi2_val) + res_val + mcarry;
				//ulong val = ((ulong)([bi1 getDataAtIndex:i] * [bi2 getDataAtIndex:j]) + (ulong)[result getDataAtIndex:k] + mcarry);
				
				[result setData: (uint)(val & 0xFFFFFFFF) atIndex:k];
				mcarry = val >> 32;
#if _BI_DEBUG_
				NSLog(@"mcarry=%qi", mcarry);
#endif
			}
			
			if (mcarry != 0) {
				[result setData:(uint)mcarry atIndex:(i + bi2.dataLength)];
				
#if _BI_DEBUG_
				NSLog(@"mcarry != 0; mcarry=%qi", mcarry);
#endif
			}
		}
	}
	@catch (NSException *ex) {
		@throw [NSException exceptionWithName:@"ArithmeticException" reason:@"Multiplication overflow." userInfo:nil];
	}
	
	result.dataLength = self.dataLength + bi2.dataLength;
	if (result.dataLength > MAX_LENGTH) {
		result.dataLength = MAX_LENGTH;
	}
	
	while (result.dataLength > 1 && [result getDataAtIndex:(result.dataLength - 1)] == 0) {
		result.dataLength--;
	}
	
	// overflow check (result is -ve)
	if (([result getDataAtIndex:lastPos] & 0x80000000) != 0) {
		if (bi1Neg != bi2Neg && [result getDataAtIndex:lastPos] == 0x80000000)    // different sign
		{
			// handle the special case where multiplication produces
			// a max negative number in 2's complement.
			
			if (result.dataLength == 1) {
				[result retain];
				[pool release];
				return [result autorelease];
			} else {
				bool isMaxNeg = true;
				for (int i = 0; i < result.dataLength - 1 && isMaxNeg; i++) {
					if ([result getDataAtIndex:i] != 0)
						isMaxNeg = false;
				}
				
				if (isMaxNeg) {
					[result retain];
					[pool release];
					return [result autorelease];
				}
			}
		}
		
		@throw [NSException exceptionWithName:@"ArithmeticException" reason:@"overflow." userInfo:nil];
	}
	
	// if input has different signs, then result is -ve
	if (bi1Neg != bi2Neg)
		result = [result negate];
	
	[result retain];
	[pool release];
	
	return [result autorelease];	
}

-(BigInt *)negate {
#if _BI_DEBUG_
	NSLog(@"negagte");
#endif
	
	if (self.dataLength == 1 && [self getDataAtIndex:0] == 0)
		return [BigInt create];
	
	NSAutoreleasePool *pool = [[NSAutoreleasePool alloc] init];
	
	BigInt *result = [BigInt createFromBigInt:self];
	
	// 1's complement
	for (int i = 0; i < MAX_LENGTH; i++)
		[result setData: (uint)(~([self getDataAtIndex:i])) atIndex: i];
	
	// add one to result of 1's complement
	ulong val, carry = 1;
	int index = 0;
	
	while (carry != 0 && index < MAX_LENGTH) {
		val = (long long)([result getDataAtIndex:index]);
		val++;
		
		[result setData:(uint)(val & 0xFFFFFFFF) atIndex: index];
		carry = val >> 32;
		
		index++;
	}
	
	if (([self getDataAtIndex:(MAX_LENGTH - 1)] & 0x80000000) == ([result getDataAtIndex:(MAX_LENGTH - 1)] & 0x80000000))
		@throw [NSException exceptionWithName:@"ArithmeticException" reason:@"Overflow in negation." userInfo:nil];
	
	result.dataLength = MAX_LENGTH;
	
	while (result.dataLength > 1 && [result getDataAtIndex:(result.dataLength - 1)] == 0)
		result.dataLength--;
	
	[result retain];
	
	[pool release];
	
	return [result autorelease];
}

-(BigInt *)divide:(BigInt *)bi2 {
#if _BI_DEBUG_
	NSLog(@"divide");
#endif
	
	NSAutoreleasePool *pool = [[NSAutoreleasePool alloc] init];
	
	BigInt *quotient = [BigInt createFromLong:0];
	BigInt *remainder = [BigInt createFromLong:0];
	BigInt *bi1 = [BigInt createFromBigInt:self];
	
	int lastPos = MAX_LENGTH - 1;
	bool divisorNeg = false, dividendNeg = false;
	
	if (([bi1 getDataAtIndex:lastPos] & 0x80000000) != 0)     // bi1 negative
	{
		bi1 = [bi1 negate];
		dividendNeg = true;
	}
	if (([bi2 getDataAtIndex:lastPos] & 0x80000000) != 0)     // bi2 negative
	{
		bi2 = [bi2 negate];
		divisorNeg = true;
	}
	
	if ([bi1 lessThan: bi2]) {
		//return quotient;
	} else {
		if (bi2.dataLength == 1)
			[BigInt singleByteDivide:bi1 bi2:bi2 outQuotient:quotient outRemainder:remainder];
		else
			[BigInt multiByteDivide:bi1 bi2:bi2 outQuotient:quotient outRemainder:remainder];
		
		if (dividendNeg != divisorNeg)
			quotient = [quotient negate];
	}	
	
	[quotient retain];
	
	[pool release];
	
	return [quotient autorelease];
}

-(BigInt *)not {
#if _BI_DEBUG_
	NSLog(@"not");
#endif	
	
	BigInt *result = [[BigInt alloc] initWithBigInt:self];
	
	for (int i = 0; i < MAX_LENGTH; i++)
		[result setData:(uint)(~([self getDataAtIndex:i])) atIndex:i];
	
	result.dataLength = MAX_LENGTH;
	
	while (result.dataLength > 1 && [result getDataAtIndex:(result.dataLength - 1)] == 0)
		result.dataLength--;
	
	return [result autorelease];	
}


-(BigInt *)shiftLeft:(int)shiftVal {
#if _BI_DEBUG_
	NSLog(@"shiftLeft");
#endif
	
	BigInt *result = [[BigInt alloc] initWithBigInt:self];
	
	/*
	uint tmp[MAX_LENGTH];
	bzero(tmp, sizeof(tmp));
	
	for(int i = 0; i < result.dataLength; i++) {
		tmp[i] = [result getDataAtIndex:i];
	}
	 */
	
	result.dataLength = [BigInt shiftLeft:[result getData] withSizeOf:result.dataLength bits:shiftVal];
	
	/*
	for(int i = 0; i < result.dataLength; i++) {
		[result setData: tmp[i] atIndex: i];
	}
	 */
	
	return [result autorelease];
}


-(BigInt *)shiftRight:(int)shiftVal {
#if _BI_DEBUG_
	NSLog(@"shiftRight");
#endif
	
	BigInt *result = [[BigInt alloc] initWithBigInt:self];
	result.dataLength = [BigInt shiftRight:[result getData] withSizeOf:result.dataLength bits:shiftVal];
	return [result autorelease];
}

-(BOOL)equals:(BigInt *)bi {
#if _BI_DEBUG_
	NSLog(@"equals");
#endif
	
	if (self.dataLength != bi.dataLength)
		return false;
	
	for (int i = 0; i < self.dataLength; i++) {
		if ([self getDataAtIndex:i] != [bi getDataAtIndex:i])
			return false;
	}
	return true;	
}

-(BOOL)greaterThan:(BigInt *)bi2 {
#if _BI_DEBUG_
	NSLog(@"greaterThan");
#endif
	
	int pos = MAX_LENGTH - 1;
	
	// bi1 is negative, bi2 is positive
	if (([self getDataAtIndex:pos] & 0x80000000) != 0 && ([bi2 getDataAtIndex:pos] & 0x80000000) == 0)
		return false;
	
	// bi1 is positive, bi2 is negative
	else if (([self getDataAtIndex:pos] & 0x80000000) == 0 && ([bi2 getDataAtIndex:pos] & 0x80000000) != 0)
		return true;
	
	// same sign
	int len = (self.dataLength > bi2.dataLength) ? self.dataLength : bi2.dataLength;
	for (pos = len - 1; pos >= 0 && [self getDataAtIndex:pos] == [bi2 getDataAtIndex:pos]; pos--) ;
	
	if (pos >= 0) {
		if ([self getDataAtIndex:pos] > [bi2 getDataAtIndex:pos])
			return true;
		return false;
	}
	return false;
}

-(BOOL)lessThan:(BigInt *)bi2 {
#if _BI_DEBUG_
	NSLog(@"lessThan");
#endif
	
	int pos = MAX_LENGTH - 1;
	
	// bi1 is negative, bi2 is positive
	if (([self getDataAtIndex:pos] & 0x80000000) != 0 && ([bi2 getDataAtIndex:pos] & 0x80000000) == 0)
		return true;
	
	// bi1 is positive, bi2 is negative
	else if (([self getDataAtIndex:pos] & 0x80000000) == 0 && ([bi2 getDataAtIndex:pos] & 0x80000000) != 0)
		return false;
	
	// same sign
	int len = (self.dataLength > bi2.dataLength) ? self.dataLength : bi2.dataLength;
	for (pos = len - 1; pos >= 0 && [self getDataAtIndex:pos] == [bi2 getDataAtIndex:pos]; pos--) ;
	
	if (pos >= 0) {
		if ([self getDataAtIndex:pos] < [bi2 getDataAtIndex:pos])
			return true;
		return false;
	}
	return false;
}

-(BOOL)lessThanOrEqualTo:(BigInt *)bi2 {
	return ([self lessThan:bi2] || [self equals:bi2]);
}

-(BOOL)greaterThanOrEqualTo:(BigInt *)bi2 {
	return ([self greaterThan:bi2] || [self equals:bi2]);
}

//***********************************************************************
// Returns the absolute value
//***********************************************************************
-(BigInt *)abs {
#if _BI_DEBUG_
	NSLog(@"abs");
#endif
	
	if (([self getDataAtIndex:(MAX_LENGTH - 1)] & 0x80000000) != 0)
		return ([self negate]);
	else
		return [BigInt createFromBigInt:self];
}

-(BigInt *)sqrt {
#if _BI_DEBUG_
	NSLog(@"sqrt");
#endif
	
	NSAutoreleasePool *pool = [[NSAutoreleasePool alloc] init];
	
	uint numBits = (uint)[self bitCount];
	
	if ((numBits & 0x1) != 0)        // odd number of bits
		numBits = (numBits >> 1) + 1;
	else
		numBits = (numBits >> 1);
	
	uint bytePos = numBits >> 5;
	int bitPos = (int)(numBits & 0x1F);
	
	uint mask;
	
	BigInt *result = [BigInt create];
	if (bitPos == 0)
		mask = 0x80000000;
	else {
		mask = (uint)1 << bitPos;
		bytePos++;
	}
	result.dataLength = (int)bytePos;
	
	for (int i = (int)bytePos - 1; i >= 0; i--) {
		while (mask != 0) {
			// guess
			[result setData:([result getDataAtIndex:i] ^ mask) atIndex:i];
			
			// undo the guess if its square is larger than this
			if ([[result multiply: result] greaterThan: self])
				[result setData:([result getDataAtIndex:i] ^ mask) atIndex:i];
			
			mask >>= 1;
		}
		mask = 0x80000000;
	}
	
	[result retain];
	[pool release];
	return result;	
}

-(BigInt *)gcd:(BigInt *)bi {
#if _BI_DEBUG_
	NSLog(@"gcd");
#endif
	
	NSAutoreleasePool *pool = [[NSAutoreleasePool alloc] init];
	
	BigInt *x;
	BigInt *y;
	
	if (([self getDataAtIndex:(MAX_LENGTH - 1)] & 0x80000000) != 0)     // negative
		x = [self negate];
	else
		x = [BigInt createFromBigInt:self];
	
	if (([bi getDataAtIndex:(MAX_LENGTH - 1)] & 0x80000000) != 0)     // negative
		y = [bi negate];
	else
		y = [BigInt createFromBigInt:bi];
	
	BigInt *g = y;
	
	while (x.dataLength > 1 || (x.dataLength == 1 && [x getDataAtIndex:0] != 0)) {
		g = x;
		x = [y mod: x];
		y = g;
	}
	
	[g retain];
	[pool release];
	
	return g;	
}

//***********************************************************************
// Returns a string representing the BigInt in sign-and-magnitude
// format in the specified radix.
//
// Example
// -------
// If the value of BigInt is -255 in base 10, then
// [self toStringWithRadix:16] returns "-FF"
//
//***********************************************************************
-(NSString *)toStringWithRadix:(int)radix {
#if _BI_DEBUG_
	NSLog(@"toStringWithRadix");
#endif
	
	if (radix < 2 || radix > 36)
		@throw [NSException exceptionWithName:@"ArgumentException" reason:@"Radix must be >= 2 and <= 36" userInfo:nil];
	
	NSAutoreleasePool *pool = [[NSAutoreleasePool alloc] init];
	
	NSString *charSet = @"ABCDEFGHIJKLMNOPQRSTUVWXYZ";
	NSMutableString *result = [[NSMutableString alloc] init];
	
	BigInt *a = self;
	
	bool negative = false;
	if (([a getDataAtIndex:(MAX_LENGTH - 1)] & 0x80000000) != 0) {
		negative = true;
		@try {
			a = [a negate];
		}
		@catch (NSException *ex) { }
	}
	
	BigInt *quotient = [BigInt create];
	BigInt *remainder = [BigInt create];
	BigInt *biRadix = [BigInt createFromLong:radix];
	
	if (a.dataLength == 1 && [a getDataAtIndex:0] == 0)
		[result appendString:@"0"];
	else {
		while (a.dataLength > 1 || (a.dataLength == 1 && [a getDataAtIndex:0] != 0)) {
			[BigInt singleByteDivide:a bi2:biRadix outQuotient:quotient outRemainder:remainder];
			
			if ([remainder getDataAtIndex:0] < 10) {
				[result insertString:[NSString stringWithFormat:@"%d", [remainder getDataAtIndex:0]] atIndex:0];
			} else {
				//NSLog(@"Value = %d", ([remainder getDataAtIndex:0] - 10));
				[result insertString:[NSString stringWithFormat:@"%c", [charSet characterAtIndex:([remainder getDataAtIndex:0] - 10)]] atIndex:0];
			}
			
			a = quotient;
		}
		if (negative)
			[result insertString:@"-" atIndex:0];
	}
	
	NSString *sOut = [result copy];
	[result release];
	
	[sOut retain];
	[pool release];
	
	return [sOut autorelease];	
}

+(int)shiftLeft:(uint *)buffer withSizeOf:(int)bufferSize bits:(int)shiftVal {
#if _BI_DEBUG_
	NSLog(@"shiftLeft");

	
	for(int i = 0; i < bufferSize; i++)
		NSLog(@"buffer[%d] = %ld", i, buffer[i]);
#endif
	
	int shiftAmount = 32;
	int bufLen = bufferSize;
	
	while (bufLen > 1 && buffer[bufLen - 1] == 0)
		bufLen--;
	
	for (int count = shiftVal; count > 0; ) {
		if (count < shiftAmount)
			shiftAmount = count;
		
		//Console.WriteLine("shiftAmount = {0}", shiftAmount);
		
		ulong carry = 0;
		for (int i = 0; i < bufLen; i++) {
			ulong val = ((ulong)buffer[i]) << shiftAmount;
			val |= carry;
			
			buffer[i] = (uint)(val & 0xFFFFFFFF);
			carry = val >> 32;
		}
		
		if (carry != 0) {
			if (bufLen + 1 <= bufferSize) {
				buffer[bufLen] = (uint)carry;
				bufLen++;
			}
		}
		count -= shiftAmount;
	}
	return bufLen;	
}

+(int)shiftRight:(uint *)buffer withSizeOf:(int)bufferSize bits:(int)shiftVal {
#if _BI_DEBUG_
	NSLog(@"shiftRight");
#endif
	
	int shiftAmount = 32;
	int invShift = 0;
	int bufLen = bufferSize;
	
	while (bufLen > 1 && buffer[bufLen - 1] == 0)
		bufLen--;
	
	for (int count = shiftVal; count > 0; ) {
		if (count < shiftAmount) {
			shiftAmount = count;
			invShift = 32 - shiftAmount;
		}
		
		ulong carry = 0;
		for (int i = bufLen - 1; i >= 0; i--) {
			ulong val = ((ulong)buffer[i]) >> shiftAmount;
			val |= carry;
			
			carry = ((ulong)buffer[i]) << invShift;
			buffer[i] = (uint)(val);
		}
		
		count -= shiftAmount;
	}
	
	while (bufLen > 1 && buffer[bufLen - 1] == 0)
		bufLen--;
	
	return bufLen;
}

+(void)singleByteDivide:(BigInt *)bi1 bi2:(BigInt *)bi2 outQuotient:(BigInt *)outQuotient outRemainder:(BigInt *)outRemainder {
#if _BI_DEBUG_
	NSLog(@"singleByteDivide");
#endif
	
	uint result[MAX_LENGTH];
	int resultPos = 0;
	
	// copy dividend to reminder
	for (int i = 0; i < MAX_LENGTH; i++)
		[outRemainder setData:[bi1 getDataAtIndex:i] atIndex: i];
	(outRemainder).dataLength = bi1.dataLength;
	
	while ((outRemainder).dataLength > 1 && [outRemainder getDataAtIndex:((outRemainder).dataLength - 1)] == 0)
		(outRemainder).dataLength--;
	
	ulong divisor = (ulong)[bi2 getDataAtIndex:0];
	int pos = (outRemainder).dataLength - 1;
	ulong dividend = (ulong)[outRemainder getDataAtIndex:pos];
	
	//NSLog(@"divisor = %d dividend = %d", divisor, dividend);
	
	if (dividend >= divisor) {
		ulong quotient = dividend / divisor;
		result[resultPos++] = (uint)quotient;
		
		[outRemainder setData:(uint)(dividend % divisor) atIndex:pos];
	}
	pos--;
	
	while (pos >= 0) {
		dividend = ((ulong)[outRemainder getDataAtIndex:pos + 1] << 32) + (ulong)[outRemainder getDataAtIndex:pos];
		ulong quotient = dividend / divisor;
		result[resultPos++] = (uint)quotient;
		
		[outRemainder setData:0 atIndex:pos + 1];
		[outRemainder setData:(uint)(dividend % divisor) atIndex:pos--];
	}
	
	(*outQuotient).dataLength = resultPos;
	int j = 0;
	for (int i =(*outQuotient).dataLength - 1; i >= 0; i--, j++)
		[outQuotient setData:result[i] atIndex:j];
	
	for (; j < MAX_LENGTH; j++)
		[outQuotient setData:0 atIndex:j];
	
	while ((outQuotient).dataLength > 1 && [outQuotient getDataAtIndex:(outQuotient).dataLength - 1] == 0)
		(outQuotient).dataLength--;
	
	if ((outQuotient).dataLength == 0)
		(outQuotient).dataLength = 1;
	
	while ((outRemainder).dataLength > 1 && [outRemainder getDataAtIndex:(outRemainder).dataLength - 1] == 0)
		(outRemainder).dataLength--;
}

+(void)multiByteDivide:(BigInt *)bi1 bi2:(BigInt *)bi2 outQuotient:(BigInt *)outQuotient outRemainder:(BigInt *)outRemainder {
#if _BI_DEBUG_	
	NSLog(@"multiByteDivide");
#endif
	
	NSAutoreleasePool *pool = [[NSAutoreleasePool alloc] init];
	
	uint result[MAX_LENGTH];
	bzero(result, sizeof(result));
	
	int remainderLen = bi1.dataLength + 1;
	uint remainder[remainderLen];
	bzero(remainder, sizeof(remainder));
	
	uint mask = 0x80000000;
	uint val = [bi2 getDataAtIndex:(bi2.dataLength - 1)];
	int shift = 0, resultPos = 0;
	
	while (mask != 0 && (val & mask) == 0) {
		shift++; mask >>= 1;
	}
	
	for (int i = 0; i < bi1.dataLength; i++)
		remainder[i] = [bi1 getDataAtIndex:i];
	
#if _BI_DEBUG_
	for(int i = 0; i < remainderLen; i++)
		NSLog(@"before: remainder[%d] = %u", i, remainder[i]);
#endif
	
	[BigInt shiftLeft:(uint *)&remainder withSizeOf:remainderLen bits:shift];
	
#if _BI_DEBUG_
	for(int i = 0; i < remainderLen; i++)
		NSLog(@"after: remainder[%d] = %u", i, remainder[i]);
#endif
	
	bi2 = [bi2 shiftLeft: shift];
	
	int j = remainderLen - bi2.dataLength;
	int pos = remainderLen - 1;
	
	ulong firstDivisorByte = [bi2 getDataAtIndex:(bi2.dataLength - 1)];
	ulong secondDivisorByte = [bi2 getDataAtIndex:(bi2.dataLength - 2)];
	
	int divisorLen = bi2.dataLength + 1;
	uint dividendPart[divisorLen];	
	bzero(dividendPart, sizeof(dividendPart));
	
	while (j > 0) {
		ulong dividend = ((ulong)remainder[pos] << 32) + (ulong)remainder[pos - 1];
		
		ulong q_hat = dividend / firstDivisorByte;
		ulong r_hat = dividend % firstDivisorByte;
		
		bool done = false;
		while (!done) {
			done = true;
			
			if (q_hat == 0x10000000 ||
				(q_hat * secondDivisorByte) > ((r_hat << 32) + (ulong)remainder[pos - 2])) {
				q_hat--;
				r_hat += firstDivisorByte;
				
				if (r_hat < 0x10000000)
					done = false;
			}
		}
		
		for (int h = 0; h < divisorLen; h++) {
#if _BI_DEBUG_
			NSLog(@"dividendPart[%d] = %u", h, remainder[pos - h]);
#endif
			
			dividendPart[h] = remainder[pos - h];
		}
		
#if _BI_DEBUG_
		NSLog(@"q_hat = %qi", q_hat);
#endif
		
		BigInt *kk = [[BigInt alloc] initWithUIntArray:dividendPart withSize:divisorLen];
		BigInt *ss = [bi2 multiply:[BigInt createFromULong:(ulong)q_hat]];
		
		while ([ss greaterThan: kk]) {
			q_hat--;
			
			/*
			BigInt *tmp = [ss divide: bi2];
			if([tmp getDataAtIndex:0] == 0) {
				tmp = [tmp add:[BigInt createFromInt:1]];
			}
			
			ss = [ss subtract:[bi2 multiply:tmp]];
			*/
			ss = [ss subtract:bi2];
		}
		BigInt *yy = [kk subtract:ss];
		
		for (int h = 0; h < divisorLen; h++) {
			remainder[pos - h] = [yy getDataAtIndex:(bi2.dataLength - h)];
		}
		
#if _BI_DEBUG_
		for(int i = 0; i < remainderLen; i++) {
			NSLog(@"remainder[%d] = %u", i, remainder[i]);
		}
#endif
		
		
		result[resultPos++] = (uint)q_hat;
		
		pos--;
		j--;
		
		[kk release];
	}
	
	outQuotient.dataLength = resultPos;
	int y = 0;
	
	for (int x = outQuotient.dataLength - 1; x >= 0; x--, y++) {
		[outQuotient setData:result[x] atIndex:y];
	}
	
	for (; y < MAX_LENGTH; y++) {
		[outQuotient setData:0 atIndex:y];
	}
	
	while (outQuotient.dataLength > 1 && [outQuotient getDataAtIndex:(outQuotient.dataLength - 1)] == 0) {
		outQuotient.dataLength--;
	}
	
	if (outQuotient.dataLength == 0) {
		outQuotient.dataLength = 1;
	}
	
	outRemainder.dataLength = [BigInt shiftRight:remainder withSizeOf:remainderLen bits:shift];
	
	for (y = 0; y < outRemainder.dataLength; y++) {
		[outRemainder setData:remainder[y] atIndex:y];
	}
	
	for (; y < MAX_LENGTH; y++) {
		[outRemainder setData:0 atIndex:y];
	}
	
	[pool release];
}

//***********************************************************************
// Modulo
//***********************************************************************
-(BigInt *)mod:(BigInt *)bi2 {
	
#if _BI_DEBUG_
	NSLog(@"mod");
#endif
	
	NSAutoreleasePool *pool = [[NSAutoreleasePool alloc] init];
	
	BigInt *quotient = [BigInt create];
	BigInt *remainder =	[BigInt createFromBigInt:self];
	BigInt *bi1 = [BigInt createFromBigInt:self];
	
	int lastPos = MAX_LENGTH - 1;
	bool dividendNeg = false;
	
	if (([bi1 getDataAtIndex:lastPos] & 0x80000000) != 0)     // bi1 negative
	{
		bi1 = [bi1 negate];
		dividendNeg = true;
	}
	if (([bi2 getDataAtIndex:lastPos] & 0x80000000) != 0)     // bi2 negative
		bi2 = [bi2 negate];
	
	if ([bi1 lessThan: bi2]) {
		//return remainder;
	} else {
		if (bi2.dataLength == 1)
			[BigInt singleByteDivide:bi1 bi2:bi2 outQuotient:quotient outRemainder:remainder];
		else
			[BigInt multiByteDivide:bi1 bi2:bi2 outQuotient:quotient outRemainder:remainder];
		if (dividendNeg)
			remainder = [remainder negate];
	}	
	
	[remainder retain];
	[pool release];
	
	return remainder;
}

//***********************************************************************
// bitwise AND
//***********************************************************************
-(BigInt *)and:(BigInt *)bi2 {
#if _BI_DEBUG_
	NSLog(@"and");
#endif
	
	BigInt *result = [BigInt createFromLong:0];
	
	int len = (self.dataLength > bi2.dataLength) ? self.dataLength : bi2.dataLength;
	
	for (int i = 0; i < len; i++) {
		uint sum = (uint)([self getDataAtIndex:i] & [bi2 getDataAtIndex:i]);
		[result setData:sum atIndex:i];
	}
	
	result.dataLength = MAX_LENGTH;
	
	while (result.dataLength > 1 && [result getDataAtIndex:(result.dataLength - 1)] == 0)
		result.dataLength--;
	
	return result;	
}

//***********************************************************************
// bitwise OR
//***********************************************************************
-(BigInt *)or:(BigInt *)bi2 {
#if _BI_DEBUG_
	NSLog(@"or");
#endif
	
	BigInt *result = [BigInt createFromLong:0];
	
	int len = (self.dataLength > bi2.dataLength) ? self.dataLength : bi2.dataLength;
	
	for (int i = 0; i < len; i++) {
		uint sum = (uint)([self getDataAtIndex:i] | [bi2 getDataAtIndex:i]);
		[result setData:sum atIndex:i];
	}
	
	result.dataLength = MAX_LENGTH;
	
	while (result.dataLength > 1 && [result getDataAtIndex:(result.dataLength - 1)] == 0)
		result.dataLength--;
	
	return result;	
}

//***********************************************************************
// bitwise XOR
//***********************************************************************
-(BigInt *)xOr:(BigInt *)bi2 {
#if _BI_DEBUG_
	NSLog(@"xOr");
#endif
	
	BigInt *result = [BigInt createFromLong:0];
	
	int len = (self.dataLength > bi2.dataLength) ? self.dataLength : bi2.dataLength;
	
	for (int i = 0; i < len; i++) {
		uint sum = (uint)([self getDataAtIndex:i] ^ [bi2 getDataAtIndex:i]);
		[result setData:sum atIndex:i];
	}
	
	result.dataLength = MAX_LENGTH;
	
	while (result.dataLength > 1 && [result getDataAtIndex:(result.dataLength - 1)] == 0)
		result.dataLength--;
	
	return result;	
}

//***********************************************************************
// Returns the position of the most significant bit in the BigInt.
//
// Eg.  The result is 0, if the value of BigInt is 0...0000 0000
//      The result is 1, if the value of BigInt is 0...0000 0001
//      The result is 2, if the value of BigInt is 0...0000 0010
//      The result is 2, if the value of BigInt is 0...0000 0011
//
//***********************************************************************
-(int)bitCount {
#if _BI_DEBUG_
	NSLog(@"bitCount");
#endif
	
	while (self.dataLength > 1 && [self getDataAtIndex:(self.dataLength - 1)] == 0)
		self.dataLength--;
	
	uint value = [self getDataAtIndex:(self.dataLength - 1)];
	uint mask = 0x80000000;
	int bits = 32;
	
	while (bits > 0 && (value & mask) == 0) {
		bits--;
		mask >>= 1;
	}
	bits += ((self.dataLength - 1) << 5);
	
	return bits;
}

//***********************************************************************
// Modulo Exponentiation
//***********************************************************************
-(BigInt *)modPow:(BigInt *)exp withMod:(BigInt *)mod {
#if _BI_DEBUG_
	NSLog(@"modPow");
#endif
	
	if (([exp getDataAtIndex:(MAX_LENGTH - 1)] & 0x80000000) != 0)
		@throw [NSException exceptionWithName:@"ArithmeticException" reason:@"Positive exponents only." userInfo:nil];
	
	NSAutoreleasePool *pool = [[NSAutoreleasePool alloc] init];
	
	BigInt *resultNum = [BigInt createFromLong:1];
	BigInt *tempNum;
	bool thisNegative = false;
	
	if (([self getDataAtIndex:(MAX_LENGTH - 1)] & 0x80000000) != 0)   // negative this
	{
		tempNum = [[self negate] mod: mod];
		thisNegative = true;
	}
	else
		tempNum = [self mod: mod];  // ensures (tempNum * tempNum) < b^(2k)
	
	if (([mod getDataAtIndex:(MAX_LENGTH - 1)] & 0x80000000) != 0)   // negative n
		mod = [mod negate];
	
	// calculate constant = b^(2k) / m
	BigInt *constant = [BigInt create];
	
	int i = mod.dataLength << 1;
	[constant setData:0x00000001 atIndex:i];
	constant.dataLength = i + 1;
	
	constant = [constant divide: mod];
	int totalBits = [exp bitCount];
	int count = 0;
	
	// perform squaring and multiply exponentiation
	for (int pos = 0; pos < exp.dataLength; pos++) {
		uint mask = 0x01;
		
		//NSLog(@"Pos = %d of %d", pos, exp.dataLength);
		
		for (int index = 0; index < 32; index++) {
			
			//NSLog(@"Pos = %d; Index = %d", pos, index);
			
			if (([exp getDataAtIndex:pos] & mask) != 0)
				resultNum = [BigInt barrettReduction:[resultNum multiply: tempNum] andN:mod andConstant:constant];
			
			mask <<= 1;
			
			tempNum = [BigInt barrettReduction:[tempNum multiply:tempNum] andN:mod andConstant:constant];
			
			if (tempNum.dataLength == 1 && [tempNum getDataAtIndex:0] == 1) {
				if (thisNegative && ([exp getDataAtIndex:0] & 0x1) != 0)  {   //odd exp
					resultNum = [resultNum negate];
				}
				[resultNum retain];
				[pool release];
				return [resultNum autorelease];
			}
			count++;
			if (count == totalBits)
				break;
		}
	}
	
	if (thisNegative && ([exp getDataAtIndex:0] & 0x1) != 0)    //odd exp
		resultNum = [resultNum negate];
	
	[resultNum retain];
	[pool release];
	
	return [resultNum autorelease];
}


//***********************************************************************
// Fast calculation of modular reduction using Barrett's reduction.
// Requires x < b^(2k), where b is the base.  In this case, base is
// 2^32 (uint).
//
// Reference [4]
//***********************************************************************
+(BigInt *)barrettReduction:(BigInt *)x andN:(BigInt *)n andConstant:(BigInt *)constant {
#if _BI_DEBUG_
	NSLog(@"barrettReduction: x:%@ n:%@ c:%@", [x toStringWithRadix:10], [n toStringWithRadix:10], [constant toStringWithRadix:10]);
#endif
	
	NSAutoreleasePool *pool = [[NSAutoreleasePool alloc] init];
	
	int k = n.dataLength;
	int kPlusOne = k + 1;
	int kMinusOne = k - 1;
	
	BigInt *q1 = [BigInt create];	
	
	// q1 = x / b^(k-1)
	for (int i = kMinusOne, j = 0; i < x.dataLength; i++, j++) {
		[q1 setData:[x getDataAtIndex:i] atIndex:j];
	}
	
	q1.dataLength = x.dataLength - kMinusOne;
	if (q1.dataLength <= 0) {
		q1.dataLength = 1;
	}
	
	BigInt *q2 = [q1 multiply: constant];
	BigInt *q3 = [BigInt create];
	
	// q3 = q2 / b^(k+1)
	for (int i = kPlusOne, j = 0; i < q2.dataLength; i++, j++) {
		[q3 setData:[q2 getDataAtIndex:i] atIndex:j];
	}
	
	q3.dataLength = q2.dataLength - kPlusOne;
	if (q3.dataLength <= 0) {
		q3.dataLength = 1;
	}
	
	// r1 = x mod b^(k+1)
	// i.e. keep the lowest (k+1) words
	BigInt *r1 = [BigInt create];
	int lengthToCopy = (x.dataLength > kPlusOne) ? kPlusOne : x.dataLength;
	for (int i = 0; i < lengthToCopy; i++)
		[r1 setData:[x getDataAtIndex:i] atIndex:i];
	
	r1.dataLength = lengthToCopy;
	
	// r2 = (q3 * n) mod b^(k+1)
	// partial multiplication of q3 and n
	
	BigInt *r2 = [BigInt create];
	for (int i = 0; i < q3.dataLength; i++) {
		if ([q3 getDataAtIndex:i] == 0) continue;
		
		ulong mcarry = 0;
		int t = i;
		for (int j = 0; j < n.dataLength && t < kPlusOne; j++, t++) {
			// t = i + j
			ulong val = ((ulong)[q3 getDataAtIndex:i] * (ulong)[n getDataAtIndex:j]) + (ulong)[r2 getDataAtIndex:t] + mcarry;
			
			[r2 setData:(uint)(val & 0xFFFFFFFF) atIndex:t];
			mcarry = (val >> 32);
		}
		
		if (t < kPlusOne) {
			[r2 setData:(uint)mcarry atIndex:t];
		}
	}
	
	r2.dataLength = kPlusOne;
	while (r2.dataLength > 1 && [r2 getDataAtIndex:(r2.dataLength - 1)] == 0) {
		r2.dataLength--;
	}
	
	r1 = [r1 subtract:r2];
	
	if (([r1 getDataAtIndex:(MAX_LENGTH - 1)] & 0x80000000) != 0)        // negative
	{
		BigInt *val = [BigInt create];
		[val setData:0x00000001 atIndex:kPlusOne];
		val.dataLength = kPlusOne + 1;
		r1 = [r1 add: val];
	}
	
#if _BI_DEBUG_	
	long cnt = 0;
	
	NSLog(@"r1=%@; n=%@", [r1 toStringWithRadix:10], [n toStringWithRadix:10]);
#endif
	
	while ([r1 greaterThanOrEqualTo:n]) {
#if _BI_DEBUG_
		NSLog(@"loop count %ld", ++cnt);
#endif
		r1 = [r1 subtract:n];
		
		/*BigInt *tmp = [r1 divide:n];
		if([tmp getDataAtIndex:0] == 0)
			tmp = [tmp add:[BigInt createFromInt:1]];
		r1 = [r1 subtract:[n multiply:tmp]];*/
	}
	
	[r1 retain];
	[pool release];
	return [r1 autorelease];
}

//***********************************************************************
// Populates "this" with the specified amount of random bits
//***********************************************************************
-(void)getRandomBits:(int)bits {
#if _BI_DEBUG_
	NSLog(@"getRandomBits");
#endif
	
	int dwords = bits >> 5;
	int remBits = bits & 0x1F;
	
	if (remBits != 0)
		dwords++;
	
	if (dwords > MAX_LENGTH)
		@throw [NSException exceptionWithName:@"ArithmeticException" reason:@"Number of required bits > maxLength." userInfo:nil];
	
	for (int i = 0; i < dwords; i++) {
		double d1 = ((double)(1 + arc4random())) / (double)10000000000;
		data[i] = (uint)((long long)(d1 * 0x100000000) & 0xFFFFFFFF);
	}
	
	for (int i = dwords; i < MAX_LENGTH; i++)
		data[i] = 0;
	
	if (remBits != 0) {
		uint mask = (uint)(0x01 << (remBits - 1));
		data[dwords - 1] |= mask;
		
		mask = (uint)(0xFFFFFFFF >> (32 - remBits));
		data[dwords - 1] &= mask;
	}
	else
		data[dwords - 1] |= 0x80000000;
	
	self.dataLength = dwords;
	
	if (dataLength == 0)
		dataLength = 1;
}

//***********************************************************************
// Returns the lowest 4 bytes of the BigInt as an int.
//***********************************************************************
-(int)intValue {
	return (int)data[0];
}

//***********************************************************************
// Determines whether a number is probably prime, using the Rabin-Miller's
// test.  Before applying the test, the number is tested for divisibility
// by primes < 2000
//
// Returns true if number is probably prime.
//***********************************************************************
-(BOOL)isProbablePrimeWithConfidence:(int)confidence {
#if _BI_DEBUG_
	NSLog(@"isProbablePrimeWithConfidence");
#endif
	
	NSAutoreleasePool *pool = [[NSAutoreleasePool alloc] init];
	
	BigInt *thisVal;
	if (([self getDataAtIndex:(MAX_LENGTH - 1)] & 0x80000000) != 0)        // negative
		thisVal = [self negate];
	else
		thisVal = self;
	
	// test for divisibility by primes < 2000
	for (int p = 0; p < (sizeof(primesBelow2000) / sizeof(int)); p++) {
		BigInt *divisor = [BigInt createFromLong:primesBelow2000[p]];
		
#if _BI_DEBUG_
		NSLog(@"Testing prime[%d]", p);
#endif
		
		if ([divisor greaterThanOrEqualTo: thisVal])
			break;
		
		BigInt *resultNum = [thisVal mod: divisor];
		if ([resultNum intValue] == 0) {
			[pool release];
			return false;
		}
	}
	
	BOOL result = FALSE;
	
	if ([thisVal rabinMillerTestWithConfidence:confidence]) {
		result = true;
	} else {
		result = false;
	}
	
	[pool release];
	
	return result;
}

//***********************************************************************
// Determines whether this BigInteger is probably prime using a
// combination of base 2 strong pseudoprime test and Lucas strong
// pseudoprime test.
//
// The sequence of the primality test is as follows,
//
// 1) Trial divisions are carried out using prime numbers below 2000.
//    if any of the primes divides this BigInteger, then it is not prime.
//
// 2) Perform base 2 strong pseudoprime test.  If this BigInteger is a
//    base 2 strong pseudoprime, proceed on to the next step.
//
// 3) Perform strong Lucas pseudoprime test.
//
// Returns True if this BigInteger is both a base 2 strong pseudoprime
// and a strong Lucas pseudoprime.
//
// For a detailed discussion of this primality test, see [6].
//
//***********************************************************************
-(BOOL)isProbablePrime {
#if _BI_DEBUG_
	NSLog(@"isProbablePrime");
#endif
	
	BigInt *thisVal;
	if (([self getDataAtIndex:(MAX_LENGTH - 1)] & 0x80000000) != 0)        // negative
		thisVal = [self negate];
	else
		thisVal = self;
	
	if (thisVal.dataLength == 1) {
		// test small numbers
		if ([thisVal getDataAtIndex:0] == 0 || [thisVal getDataAtIndex:0] == 1)
			return false;
		else if ([thisVal getDataAtIndex:0] == 2 || [thisVal getDataAtIndex:0] == 3)
			return true;
	}
	
	if (([thisVal getDataAtIndex:0] & 0x1) == 0)     // even numbers
		return false;
	
	// test for divisibility by primes < 2000
	for (int p = 0; p < (sizeof(primesBelow2000)/sizeof(int)); p++) {
		BigInt *divisor = [BigInt createFromLong:primesBelow2000[p]];
		
		if ([divisor greaterThanOrEqualTo: thisVal])
			break;
		
		BigInt *resultNum = [thisVal  mod: divisor];
		if ([resultNum intValue] == 0) {
			//Console.WriteLine("Not prime!  Divisible by {0}\n",
			//                  primesBelow2000[p]);
			
			return false;
		}
	}
	
	// Perform BASE 2 Rabin-Miller Test
	
	// calculate values of s and t
	BigInt *p_sub1 = [thisVal subtract:[BigInt createFromLong:1]];
	int s = 0;
	
	for (int index = 0; index < p_sub1.dataLength; index++) {
		uint mask = 0x01;
		
		for (int i = 0; i < 32; i++) {
			if (([p_sub1 getDataAtIndex:index] & mask) != 0) {
				index = p_sub1.dataLength;      // to break the outer loop
				break;
			}
			mask <<= 1;
			s++;
		}
	}
	
	BigInt *t = [p_sub1 shiftRight: s];
	
	//int bits = [thisVal bitCount];
	BigInt *a = [BigInt createFromLong: 2];
	
	// b = a^t mod p
	BigInt *b = [a modPow:t withMod:thisVal];
	bool result = false;
	
	if (b.dataLength == 1 && [b getDataAtIndex:0] == 1)         // a^t mod p = 1
		result = true;
	
	for (int j = 0; result == false && j < s; j++) {
		if ([b equals: p_sub1])         // a^((2^j)*t) mod p = p-1 for some 0 <= j <= s-1
		{
			result = true;
			break;
		}
		
		b = [[b multiply: b] mod: thisVal];
	}
	
	// if number is strong pseudoprime to base 2, then do a strong lucas test
	if (result)
		result = [BigInt lucasStrongTest:thisVal];
	
	return result;	
}

//***********************************************************************
// Generates a positive BigInt that is probably prime.
//***********************************************************************
+(BigInt *)generatePseudoPrimeWithBits:(int)bits andConfidence:(int)confidence {
	
#if _BI_DEBUG_
	NSLog(@"generatePseudoPrimeWithBits");
#endif
	
	NSAutoreleasePool *pool = [[NSAutoreleasePool alloc] init];
	BigInt *result = [BigInt create];
	bool done = false;
	
	while (!done) {
		NSAutoreleasePool *whilePool = [[NSAutoreleasePool alloc] init];
		
		[result getRandomBits: bits];
		[result setData:([result getDataAtIndex:0] | 0x01) atIndex:0];		// make it odd
		
		// prime test
		done = [result isProbablePrimeWithConfidence:confidence];
		
#if _BI_DEBUG_
		NSLog(@"Is Prime: %s (%s)", [[result toStringWithRadix:16] UTF8String], [done ? @"Yes" : @"No" UTF8String]);
#endif
		
		[whilePool release];
	}
	
	[result retain];
	[pool release];
	return [result autorelease];	
}


//***********************************************************************
// Performs the calculation of the kth term in the Lucas Sequence.
// For details of the algorithm, see reference [9].
//
// k must be odd.  i.e LSB == 1
//***********************************************************************
+(NSMutableArray *)lucasSequence:(BigInt *)P andQ:(BigInt *)Q andk:(BigInt *)k andn:(BigInt *)n andConstant:(BigInt *)constant ands:(int)s {
	
#if _BI_DEBUG_	
	NSLog(@"lucasSequence: P:%@ Q:%@ k:%@ n:%@ constant:%@ s:%d",
	[P toStringWithRadix:10],
	[Q toStringWithRadix:10],
	[k toStringWithRadix:10],
	[n toStringWithRadix:10],
	[constant toStringWithRadix:10],
	s);
#endif
	
	NSAutoreleasePool *pool = [[NSAutoreleasePool alloc] init];
	
	BigInt *result[3];
	
	result[0] = [BigInt create];
	result[1] = [BigInt create];
	result[2] = [BigInt create];
	
	if (([k getDataAtIndex:0] & 0x00000001) == 0)
		@throw [NSException exceptionWithName:@"ArgumentException" reason:@"Argument k must be odd." userInfo:nil];
	
	int numbits = [k bitCount];
	uint mask = (uint)0x1 << ((numbits & 0x1F) - 1);
	
	// v = v0, v1 = v1, u1 = u1, Q_k = Q^0
	
	BigInt *v = [[BigInt createFromLong:2] mod: n];
	BigInt *Q_k = [[BigInt createFromLong:1] mod: n];
	BigInt *v1 = [P mod: n]; 
	BigInt *u1 = [BigInt createFromBigInt:Q_k];
	
	bool flag = true;
	
	for (int i = k.dataLength - 1; i >= 0; i--)     // iterate on the binary expansion of k
	{
		//Console.WriteLine("round");
		while (mask != 0) {
			if (i == 0 && mask == 0x00000001)        // last bit
				break;
			
			if (([k getDataAtIndex:i] & mask) != 0)             // bit is set
			{
				// index doubling with addition
				
				u1 = [[u1 multiply: v1] mod: n];
				
				v = [[[v multiply: v1] subtract: [P multiply: Q_k]] mod: n];

				v1 = [BigInt barrettReduction:[v1 multiply:v1] andN:n andConstant:constant];
				//v1 = [n barrettReduction: [v1 multiply: v1] , n, constant);
				
				v1 = [[v1 subtract:[[Q_k multiply:Q] shiftLeft:1]] mod:n];
				//v1 = (v1 - ((Q_k * Q) << 1)) % n;
				
				if (flag)
					flag = false;
				else
					Q_k = [BigInt barrettReduction:[Q_k multiply: Q_k] andN:n andConstant:constant];
				
				Q_k = [[Q_k multiply: Q] mod: n];
			}
			else {
				// index doubling
				u1 = [[[u1 multiply: v] subtract: Q_k] mod: n];
				
				v1 = [[[v multiply: v1] subtract: [P multiply: Q_k]] mod: n];
				v = [BigInt barrettReduction:[v multiply: v] andN: n andConstant: constant];
				v = [[v subtract: [Q_k shiftLeft: 1]] mod: n];
				
				if (flag) {
					Q_k = [Q mod: n];
					flag = false;
				}
				else
					Q_k = [BigInt barrettReduction:[Q_k multiply: Q_k] andN: n andConstant: constant];
			}
			
			mask >>= 1;
		}
		mask = 0x80000000;
	}
	
	// at this point u1 = u(n+1) and v = v(n)
	// since the last bit always 1, we need to transform u1 to u(2n+1) and v to v(2n+1)
	
	u1 = [[[u1 multiply: v] subtract: Q_k] mod: n];
	v = [[[v multiply: v1] subtract: [P multiply: Q_k]] mod: n];
	
	if (flag)
		flag = false;
	else
		Q_k = [BigInt barrettReduction:[Q_k multiply: Q_k] andN: n andConstant: constant];
	
	Q_k = [[Q_k multiply: Q] mod:n];
	
	for (int i = 0; i < s; i++) {
		// index doubling
		u1 = [[u1 multiply: v] mod: n];
		v = [[[v multiply: v] subtract:[Q_k shiftLeft: 1]] mod: n];
		
		if (flag) {
			Q_k = [Q mod: n];
			flag = false;
		}
		else
			Q_k = [BigInt barrettReduction:[Q_k multiply: Q_k] andN: n andConstant: constant];
	}
	
	result[0] = u1;
	result[1] = v;
	result[2] = Q_k;
	
	NSMutableArray *aRet = [[NSMutableArray alloc] init];
	[aRet addObject:result[0]];
	[aRet addObject:result[1]];
	[aRet addObject:result[2]];
	[aRet retain];
	
	[pool release];
	
	return [aRet autorelease];
}


//***********************************************************************
// Implementation of the Lucas Strong Pseudo Prime test.
//
// Let n be an odd number with gcd(n,D) = 1, and n - J(D, n) = 2^s * d
// with d odd and s >= 0.
//
// If Ud mod n = 0 or V2^r*d mod n = 0 for some 0 <= r < s, then n
// is a strong Lucas pseudoprime with parameters (P, Q).  We select
// P and Q based on Selfridge.
//
// Returns True if number is a strong Lucus pseudo prime.
// Otherwise, returns False indicating that number is composite.
//***********************************************************************
+(BOOL)lucasStrongTest:(BigInt *) thisVal {
	
#if _BI_DEBUG_
	NSLog(@"lucasStrongTest");
#endif
	
	// Do the test (selects D based on Selfridge)
	// Let D be the first element of the sequence
	// 5, -7, 9, -11, 13, ... for which J(D,n) = -1
	// Let P = 1, Q = (1-D) / 4
	
	NSAutoreleasePool *pool = [[NSAutoreleasePool alloc] init];
	
	if (([thisVal getDataAtIndex:(MAX_LENGTH - 1)] & 0x80000000) != 0)        // negative
		thisVal = [thisVal negate];
	
	if (thisVal.dataLength == 1) {
		// test small numbers
		if ([thisVal getDataAtIndex:0] == 0 || [thisVal getDataAtIndex:0] == 1)
			return false;
		else if ([thisVal getDataAtIndex:0] == 2 || [thisVal getDataAtIndex:0] == 3)
			return true;
	}
	
	if (([thisVal getDataAtIndex:0] & 0x1) == 0)     // even numbers
		return false;
	
	long D = 5, sign = -1, dCount = 0;
	bool done = false;
	
	while (!done) {
		int Jresult = [BigInt jacobi:[BigInt createFromLong:D] andB:thisVal];
		
		if (Jresult == -1)
			done = true;    // J(D, this) = 1
		else {
			if (Jresult == 0 && [[BigInt createFromLong:(long)abs(D)] lessThan: thisVal])       // divisor found
				return false;
			
			if (dCount == 20) {
				// check for square
				BigInt *root = [thisVal sqrt];
				if ([[root multiply: root] equals: thisVal])
					return false;
			}
			
			//Console.WriteLine(D);
			D = ((long)abs(D) + 2) * sign;
			sign = -sign;
		}
		dCount++;
	}
	
	long Q = (1 - D) >> 2;
	
	BigInt *p_add1 = [thisVal add:[BigInt createFromLong: 1]];
	int s = 0;
	
	for (int index = 0; index < p_add1.dataLength; index++) {
		uint mask = 0x01;
		
		for (int i = 0; i < 32; i++) {
			if (([p_add1 getDataAtIndex:index] & mask) != 0) {
				index = p_add1.dataLength;      // to break the outer loop
				break;
			}
			mask <<= 1;
			s++;
		}
	}
	
	BigInt *t = [p_add1 shiftRight: s];
	
	// calculate constant = b^(2k) / m
	// for Barrett Reduction
	BigInt *constant = [BigInt create];
	
	int nLen = thisVal.dataLength << 1;
	[constant setData:0x00000001 atIndex:nLen];
	constant.dataLength = nLen + 1;
	
	constant = [constant divide: thisVal];
	
	NSMutableArray *aLucus = [BigInt lucasSequence:[BigInt createFromLong:1] andQ:[BigInt createFromLong:Q] andk:t andn:thisVal andConstant:constant ands:0];
	BigInt *lucas[3];

	for(int j = 0; j < [aLucus count] && j < 3; j++) 
		lucas[j] = [aLucus objectAtIndex:j];
	
	bool isPrime = false;
	
	if ((lucas[0].dataLength == 1 && [lucas[0] getDataAtIndex:0] == 0) ||
		(lucas[1].dataLength == 1 && [lucas[1] getDataAtIndex:0] == 0)) {
		// u(t) = 0 or V(t) = 0
		isPrime = true;
	}
	
	for (int i = 1; i < s; i++) {
		if (!isPrime) {
			// doubling of index
			lucas[1] = [BigInt barrettReduction: [lucas[1] multiply: lucas[1]] andN: thisVal andConstant: constant];
			lucas[1] = [[lucas[1] subtract: [lucas[2] shiftLeft: 1]] mod: thisVal];
			
			//lucas[1] = ((lucas[1] * lucas[1]) - (lucas[2] << 1)) % thisVal;
			
			if ((lucas[1].dataLength == 1 && [lucas[1] getDataAtIndex:0] == 0))
				isPrime = true;
		}
	
		lucas[2] = [BigInt barrettReduction:[lucas[2] multiply: lucas[2]] andN: thisVal andConstant: constant];     //Q^k
	}
	
	
	if (isPrime)     // additional checks for composite numbers
	{
		// If n is prime and gcd(n, Q) == 1, then
		// Q^((n+1)/2) = Q * Q^((n-1)/2) is congruent to (Q * J(Q, n)) mod n
		
		BigInt *g = [thisVal gcd:[BigInt createFromLong:Q]];
		if (g.dataLength == 1 && [g getDataAtIndex:0] == 1)         // gcd(this, Q) == 1
		{
			if (([lucas[2] getDataAtIndex:(MAX_LENGTH - 1)] & 0x80000000) != 0)
				lucas[2] = [lucas[2] add: thisVal];
			
			BigInt *temp = [[[BigInt createFromLong: Q] multiply: [BigInt createFromLong: [BigInt jacobi:[BigInt createFromLong:Q] andB:thisVal]]] mod: thisVal];
			if (([temp getDataAtIndex:(MAX_LENGTH - 1)] & 0x80000000) != 0)
				temp = [temp add:thisVal];
			
			if (![lucas[2] equals: temp])
				isPrime = false;
		}
	}
	
	[pool release];
	
	return isPrime;
}

//***********************************************************************
// Computes the Jacobi Symbol for a and b.
// Algorithm adapted from [3] and [4] with some optimizations
//***********************************************************************
+(int)jacobi:(BigInt *)a andB:(BigInt *)b {
	
#if _BI_DEBUG_
	NSLog(@"jacobi");
#endif
	
	// Jacobi defined only for odd integers
	if (([b getDataAtIndex:0] & 0x1) == 0)
		@throw [NSException exceptionWithName:@"ArgumentException" reason:@"Jacobi defined only for odd integers." userInfo:nil];
	
	if ([a greaterThanOrEqualTo: b]) a = [a mod: b];
	if (a.dataLength == 1 && [a getDataAtIndex:0] == 0) return 0;  // a == 0
	if (a.dataLength == 1 && [a getDataAtIndex:0] == 1) return 1;  // a == 1
	
	if ([a lessThan:[BigInt createFromLong: 0]]) {
		if ((([[b subtract:[BigInt createFromLong:1]] getDataAtIndex:0]) & 0x2) == 0)       //if( (((b-1) >> 1).data[0] & 0x1) == 0)
			return [BigInt jacobi:[a negate] andB: b];
		else
			return -[BigInt jacobi:[a negate] andB:b];
	}
	
	int e = 0;
	for (int index = 0; index < a.dataLength; index++) {
		uint mask = 0x01;
		
		for (int i = 0; i < 32; i++) {
			if (([a getDataAtIndex:index] & mask) != 0) {
				index = a.dataLength;      // to break the outer loop
				break;
			}
			mask <<= 1;
			e++;
		}
	}
	
	BigInt *a1 = [a shiftRight: e];
	
	int s = 1;
	if ((e & 0x1) != 0 && (([b getDataAtIndex:0] & 0x7) == 3 || ([b getDataAtIndex:0] & 0x7) == 5))
		s = -1;
	
	if (([b getDataAtIndex:0] & 0x3) == 3 && ([a1 getDataAtIndex:0] & 0x3) == 3)
		s = -s;
	
	if (a1.dataLength == 1 && [a1 getDataAtIndex:0] == 1)
		return s;
	else
		return (s * [BigInt jacobi:[b mod: a1] andB: a1]);
}

//***********************************************************************
// Probabilistic prime test based on Rabin-Miller's
//
// for any p > 0 with p - 1 = 2^s * t
//
// p is probably prime (strong pseudoprime) if for any a < p,
// 1) a^t mod p = 1 or
// 2) a^((2^j)*t) mod p = p-1 for some 0 <= j <= s-1
//
// Otherwise, p is composite.
//
// Returns
// -------
// True if "this" is a strong pseudoprime to randomly chosen
// bases.  The number of chosen bases is given by the "confidence"
// parameter.
//
// False if "this" is definitely NOT prime.
//
//***********************************************************************
-(BOOL)rabinMillerTestWithConfidence:(int)confidence {
	
#if _BI_DEBUG_
	NSLog(@"rabinMillerTestWithConfidence");
#endif
	
	NSAutoreleasePool *pool = [[NSAutoreleasePool alloc] init];
	
	BigInt *thisVal;
	if (([self getDataAtIndex:(MAX_LENGTH - 1)] & 0x80000000) != 0)        // negative
		thisVal = [self negate];
	else
		thisVal = self;
	
	if (thisVal.dataLength == 1) {
		// test small numbers
		if ([thisVal getDataAtIndex:0] == 0 || [thisVal getDataAtIndex:0] == 1) {
			[pool release];
			return false;
		} else if ([thisVal getDataAtIndex:0] == 2 || [thisVal getDataAtIndex:0] == 3) {
			[pool release];
			return true;
		}
	}
	
	if (([thisVal getDataAtIndex:0] & 0x1) == 0) {    // even numbers
		[pool release];
		return false;
	}
	
	
	// calculate values of s and t
	BigInt *p_sub1 = [thisVal subtract: [BigInt createFromLong: 1]];
	int s = 0;
	
	for (int index = 0; index < p_sub1.dataLength; index++) {
		uint mask = 0x01;
		
		for (int i = 0; i < 32; i++) {
			if (([p_sub1 getDataAtIndex:index] & mask) != 0) {
				index = p_sub1.dataLength;      // to break the outer loop
				break;
			}
			mask <<= 1;
			s++;
		}
	}
	
	BigInt *t = [p_sub1 shiftRight: s];
	
	int bits = [thisVal bitCount];
	BigInt *a = [BigInt create];
	
	for (int round = 0; round < confidence; round++) {
		bool done = false;
		
		while (!done)		// generate a < n
		{
			int testBits = 0;
			
			// make sure "a" has at least 2 bits
			while (testBits < 2) {
				double d1 = (((double)(1 + arc4random())) / (double)10000000000);
				
				testBits = (int)(((ulong)(d1 * bits)) & 0xFFFF);
				//testBits = (int)(((1 + arc4random()) * bits) & MAX_LENGTH);
			}
			
			[a getRandomBits:testBits];
			//a.genRandomBits(testBits, rand);
			
			int byteLen = a.dataLength;
			
			// make sure "a" is not 0
			if (byteLen > 1 || (byteLen == 1 && [a getDataAtIndex:0] != 1))
				done = true;
			
#if _BI_DEBUG_			
			NSLog(@"Rabin Miller Test: TestBits: %d; ByteLen: %d", testBits, byteLen);
#endif
		}
		
		// check whether a factor exists (fix for version 1.03)
		BigInt *gcdTest = [a gcd:thisVal];
		if (gcdTest.dataLength == 1 && [gcdTest getDataAtIndex:0] != 1) {
			[pool release];
			return false;
		}
		
		BigInt *b = [a modPow:t withMod:thisVal];
			
		bool result = false;
		
		if (b.dataLength == 1 && [b getDataAtIndex:0] == 1)         // a^t mod p = 1
			result = true;
		
		for (int j = 0; result == false && j < s; j++) {
			if ([b equals: p_sub1])         // a^((2^j)*t) mod p = p-1 for some 0 <= j <= s-1
			{
				result = true;
				break;
			}
			
			b = [[b multiply: b] mod: thisVal];
		}
		
		if (result == false) {
			[pool release];
			return false;
		}
	}
	[pool release];
	return true;
}

@end
