- update to current root trunc's version

- update to root trunc's UNICODE support

git-svn-id: svn://ttmath.org/publicrep/ttmath/branches/chk@182 e52654a7-88a9-db11-a3e9-0013d4bc506e

CHANGELOG

@@ -1,3 +1,28 @@
+Version 0.9.0 prerelease (2009.07.16):
+    * added:   support for wide characters (wchar_t)
+               wide characters are used when macro TTMATH_USE_WCHAR is defined
+               this macro is defined automatically when there is macro UNICODE or _UNICODE defined
+               some types have been changed
+                 char               -> tt_char
+                 std::string        -> tt_string
+                 std::ostringstream -> tt_ostringstream
+                 std::ostream       -> tt_ostream
+                 std::istream       -> tt_istream
+               normally tt_char is equal char but when you are using wide characters then tt_char will be wchar_t (and so on)
+               (all typedef's are in ttmathtypes.h)
+    * added:   Big::IsInteger()
+               returns true if the value is integer (without fraction)
+               (NaN flag is not checked)
+    * added:   global Gamma() function
+    * added:   gamma() function to the parser
+    * added:   CGamma<ValueType> class
+               is used with Gamma() and Factorial() in multithreaded environment
+    * changed: Factorial() is using the Gamma() function now
+    * removed: Parser<>::SetFactorialMax() method
+               the factorial() is such a fast now that we don't need the method longer
+    * removed: ErrorCode::err_too_big_factorial
+
+
 Version 0.8.5 (2009.06.16):
     * fixed:   Big::Mod(x) didn't correctly return a carry
                and the result was sometimes very big (even greater than x)

ttmath/ttmath.h

@@ -1,7 +1,7 @@
 /*
  * This file is a part of TTMath Bignum Library
  * and is distributed under the (new) BSD licence.
- * Author: Tomasz Sowa <t.sowa@slimaczek.pl>
+ * Author: Tomasz Sowa <t.sowa@ttmath.org>
  */
 
 /* 
@@ -50,12 +50,9 @@
 #pragma warning( disable: 4127 )
 #endif
 
-
 #include "ttmathbig.h"
 #include "ttmathobjects.h"
 
-#include <string>
-
 
 namespace ttmath
 {
@@ -67,6 +64,7 @@ namespace ttmath
 	 */
 
 
+
 	/*
  	 *
 	 *  functions for rounding
@@ -76,7 +74,7 @@ namespace ttmath
 
 
 	/*!
-		this method skips the fraction from x
+		this function skips the fraction from x
 		e.g  2.2  = 2
 		     2.7  = 2
 			 -2.2 = 2
@@ -93,7 +91,7 @@ namespace ttmath
 
 
 	/*!
-		this method rounds to the nearest integer value
+		this function rounds to the nearest integer value
 		e.g  2.2  = 2
 		     2.7  = 3
 			 -2.2 = -2
@@ -225,7 +223,7 @@ namespace ttmath
 
 	
 	/*!
-		this method calculates the natural logarithm (logarithm with the base 'e')
+		this function calculates the natural logarithm (logarithm with the base 'e')
 	*/
 	template<class ValueType>
 	ValueType Ln(const ValueType & x, ErrorCode * err = 0)
@@ -266,7 +264,7 @@ namespace ttmath
 
 
 	/*!
-		this method calculates the logarithm
+		this function calculates the logarithm
 	*/
 	template<class ValueType>
 	ValueType Log(const ValueType & x, const ValueType & base, ErrorCode * err = 0)
@@ -307,7 +305,7 @@ namespace ttmath
 
 
 	/*!
-		this method calculates the expression e^x
+		this function calculates the expression e^x
 	*/
 	template<class ValueType>
 	ValueType Exp(const ValueType & x, ErrorCode * err = 0)
@@ -1940,10 +1938,7 @@ namespace ttmath
 		template<class ValueType>
 		bool RootCheckIndexFrac(ValueType & x, const ValueType & index, ErrorCode * err)
 		{
-			ValueType indexfrac(index);
+			if( !index.IsInteger() )
-			indexfrac.RemainFraction();
-
-			if( !indexfrac.IsZero() )
 			{
 				// index must be integer
 				if( err )
@@ -2075,154 +2070,6 @@ namespace ttmath
 
 
 
-
-	namespace auxiliaryfunctions
-	{
-	
-	template<class ValueType>
-	uint FactorialInt(	const ValueType & x, ErrorCode * err,
-						const volatile StopCalculating * stop,
-						ValueType & result)
-	{
-		uint maxvalue = TTMATH_UINT_MAX_VALUE;
-
-		if( x < TTMATH_UINT_MAX_VALUE )
-			x.ToUInt(maxvalue);
-
-		uint multipler = 1;
-		uint carry     = 0;
-
-		while( !carry && multipler<maxvalue  )
-		{
-			if( stop && (multipler & 127)==0 ) // it means 'stop && (multipler % 128)==0'
-			{
-				// after each 128 iterations we make a test
-				if( stop->WasStopSignal() )
-				{
-					if( err )
-						*err = err_interrupt;
-
-					return 2;
-				}
-			}
-			
-			++multipler;
-			carry += result.MulUInt(multipler);
-		}
-
-		if( err )
-			*err = carry ? err_overflow : err_ok;
-
-	return carry ? 1 : 0;
-	}
-
-
-	template<class ValueType>
-	int FactorialMore(	const ValueType & x, ErrorCode * err,
-						const volatile StopCalculating * stop,
-						ValueType & result)
-	{
-		ValueType multipler(TTMATH_UINT_MAX_VALUE);
-		ValueType one;
-
-		one.SetOne();
-		uint carry = 0;
-		uint iter = 1; // only for testing the stop object
-
-		while( !carry && multipler < x )
-		{
-			if( stop && (iter & 31)==0 ) // it means 'stop && (iter % 32)==0'
-			{
-				// after each 32 iterations we make a test
-				if( stop->WasStopSignal() )
-				{
-					if( err )
-						*err = err_interrupt;
-
-					return 2;
-				}
-			}
-			
-			carry += multipler.Add(one);
-			carry += result.Mul(multipler);
-			++iter;
-		}
-
-		if( err )
-			*err = carry ? err_overflow : err_ok;
-
-	return carry ? 1 : 0;
-	}
-
-
-	} // namespace
-	
-
-
-	/*!
-		the factorial from given 'x'
-		e.g.
-		Factorial(4) = 4! = 1*2*3*4
-	*/
-	template<class ValueType>
-	ValueType Factorial(const ValueType & x, ErrorCode * err = 0, const volatile StopCalculating * stop = 0)
-	{
-	using namespace auxiliaryfunctions;
-
-	static History<ValueType> history;
-	ValueType result;
-
-		if( x.IsNan() || x.IsSign() )
-		{
-			if( err )
-				*err = err_improper_argument;
-
-		return result; // NaN set by default
-		}
-
-		result.SetOne();
-
-		if( !x.exponent.IsSign() && !x.exponent.IsZero() )
-		{
-			// when x.exponent>0 there's no sense to calculate the formula
-			// (we can't add one into the x bacause
-			// we don't have enough bits in the mantissa)
-			if( err )
-				*err = err_overflow;
-
-			result.SetNan();
-
-		return result;
-		}
-
-		ErrorCode err_tmp;
-
-		if( history.Get(x, result, err_tmp) )
-		{
-			if( err )
-				*err = err_tmp;
-		
-			return result;
-		}
-
-		uint status = FactorialInt(x, err, stop, result);
-		if( status == 0 )
-			status = FactorialMore(x, err, stop, result);
-
-		if( status == 2 )
-		{
-			// the calculation has been interrupted
-			result.SetNan();
-			return result;
-		}
-
-		err_tmp = status==1 ? err_overflow : err_ok;
-		history.Add(x, result, err_tmp);
-
-	return result;
-	}
-
-
 	/*!
 		absolute value of x
 		e.g.  -2 = 2 
@@ -2283,6 +2130,653 @@ namespace ttmath
 
 
 
+	namespace auxiliaryfunctions
+	{
+	
+	/*!
+		this function is used to store factorials in a given container
+		'more' means how many values should be added at the end
+
+			e.g.
+			std::vector<ValueType> fact;
+			SetFactorialSequence(fact, 3);
+			// now the container has three values: 1  1  2
+
+			SetFactorialSequence(fact, 2);
+			// now the container has five values:  1  1  2  6  24
+	*/
+	template<class ValueType>
+	void SetFactorialSequence(std::vector<ValueType> & fact, uint more = 20)
+	{
+		if( more == 0 )
+			more = 1;
+
+		uint start = static_cast<uint>(fact.size());
+		fact.resize(fact.size() + more);
+
+		if( start == 0 )
+		{
+			fact[0] = 1;
+			++start;
+		}
+
+		for(uint i=start ; i<fact.size() ; ++i)
+		{
+			fact[i] = fact[i-1];
+			fact[i].MulInt(i);
+		}
+	}
+
+
+	/*!
+		an auxiliary function used to calculate Bernoulli numbers
+
+		this function returns a sum:
+		sum(m) = sum_{k=0}^{m-1} {2^k * (m k) * B(k)}    k in [0, m-1]   (m k) means binomial coefficient = (m! / (k! * (m-k)!))
+
+		you should have sufficient factorials in cgamma.fact
+		(cgamma.fact should have at least m items)
+
+		n_ should be equal 2
+	*/
+	template<class ValueType>
+	ValueType SetBernoulliNumbersSum(CGamma<ValueType> & cgamma, const ValueType & n_, uint m,
+									  const volatile StopCalculating * stop = 0)
+			{
+	ValueType k_, temp, temp2, temp3, sum;
+
+		sum.SetZero();
+		
+		for(uint k=0 ; k<m ; ++k)			// k<m means k<=m-1
+		{
+			if( stop && (k & 15)==0 )		// means: k % 16 == 0
+				if( stop->WasStopSignal() )
+					return ValueType();		// NaN
+
+			if( k>1 && (k & 1) == 1 )		// for that k the Bernoulli number is zero
+				continue;
+
+			k_ = k;
+
+			temp = n_;				// n_ is equal 2
+			temp.Pow(k_);
+			// temp = 2^k
+
+			temp2 = cgamma.fact[m];
+			temp3 = cgamma.fact[k];
+			temp3.Mul(cgamma.fact[m-k]);
+			temp2.Div(temp3);
+			// temp2 = (m k) = m! / ( k! * (m-k)! )
+
+			temp.Mul(temp2);
+			temp.Mul(cgamma.bern[k]);
+
+			sum.Add(temp);
+			// sum += 2^k * (m k) * B(k)
+
+			if( sum.IsNan() )
+				break;
+		}
+
+	return sum;
+	}
+
+
+	/*!
+		an auxiliary function used to calculate Bernoulli numbers
+		start is >= 2
+
+		we use the recurrence formula: 
+		   B(m) = 1 / (2*(1 - 2^m)) * sum(m)
+		   where sum(m) is calculated by SetBernoulliNumbersSum()
+	*/
+	template<class ValueType>
+	bool SetBernoulliNumbersMore(CGamma<ValueType> & cgamma, uint start, const volatile StopCalculating * stop = 0)
+				{
+	ValueType denominator, temp, temp2, temp3, m_, sum, sum2, n_, k_;
+
+		const uint n = 2;
+		n_ = n;
+
+		// start is >= 2
+		for(uint m=start ; m<cgamma.bern.size() ; ++m)
+		{
+			if( (m & 1) == 1 )
+			{
+				cgamma.bern[m].SetZero();
+			}
+			else
+			{
+				m_ = m;
+
+				temp = n_;				// n_ = 2
+				temp.Pow(m_);
+				// temp = 2^m
+
+				denominator.SetOne();
+				denominator.Sub(temp);
+				if( denominator.exponent.AddOne() ) // it means: denominator.MulInt(2)
+					denominator.SetNan();
+
+				// denominator = 2 * (1 - 2^m)
+
+				cgamma.bern[m] = SetBernoulliNumbersSum(cgamma, n_, m, stop);
+
+				if( stop && stop->WasStopSignal() )
+				{
+					cgamma.bern.resize(m);		// valid numbers are in [0, m-1]
+					return false;
+			}
+			
+				cgamma.bern[m].Div(denominator);
+		}
+		}
+
+	return true;
+	}
+
+
+	/*!
+		this function is used to calculate Bernoulli numbers,
+		returns false if there was a stop signal,
+		'more' means how many values should be added at the end
+
+			e.g.
+			typedef Big<1,2> MyBig;
+			CGamma<MyBig> cgamma;
+			SetBernoulliNumbers(cgamma, 3);
+			// now we have three first Bernoulli numbers:  1  -0.5  0.16667
+			
+			SetBernoulliNumbers(cgamma, 4);
+			// now we have 7 Bernoulli numbers:  1  -0.5  0.16667   0   -0.0333   0   0.0238
+	*/
+	template<class ValueType>
+	bool SetBernoulliNumbers(CGamma<ValueType> & cgamma, uint more = 20, const volatile StopCalculating * stop = 0)
+	{
+		if( more == 0 )
+			more = 1;
+
+		uint start = static_cast<uint>(cgamma.bern.size());
+		cgamma.bern.resize(cgamma.bern.size() + more);
+
+		if( start == 0 )
+		{
+			cgamma.bern[0].SetOne();
+			++start;
+	}
+
+		if( cgamma.bern.size() == 1 )
+			return true;
+
+		if( start == 1 )
+		{
+			cgamma.bern[1].Set05();
+			cgamma.bern[1].ChangeSign();
+			++start;
+		}
+
+		// we should have sufficient factorials in cgamma.fact
+		if( cgamma.fact.size() < cgamma.bern.size() )
+			SetFactorialSequence(cgamma.fact, static_cast<uint>(cgamma.bern.size() - cgamma.fact.size()));
+
+
+	return SetBernoulliNumbersMore(cgamma, start, stop);
+	}
+
+
+	/*!
+		an auxiliary function used to calculate the Gamma() function
+
+		we calculate a sum:
+		   sum(n) = sum_{m=2} { B(m) / ( (m^2 - m) * n^(m-1) )  } = 1/(12*n) - 1/(360*n^3) + 1/(1260*n^5) + ...
+	       B(m) means a mth Bernoulli number
+		   the sum starts from m=2, we calculate as long as the value will not change after adding a next part
+	*/
+	template<class ValueType>
+	ValueType GammaFactorialHighSum(const ValueType & n, CGamma<ValueType> & cgamma, ErrorCode & err,
+									const volatile StopCalculating * stop)
+	{
+	ValueType temp, temp2, denominator, sum, oldsum;
+
+		sum.SetZero();
+
+		for(uint m=2 ; m<TTMATH_ARITHMETIC_MAX_LOOP ; m+=2)
+		{
+			if( stop && (m & 3)==0 ) // (m & 3)==0 means: (m % 4)==0
+				if( stop->WasStopSignal() )
+				{
+					err = err_interrupt;
+					return ValueType(); // NaN
+				}
+
+			temp = (m-1);
+			denominator = n;
+			denominator.Pow(temp);
+			// denominator = n ^ (m-1)
+
+			temp = m;
+			temp2 = temp;
+			temp.Mul(temp2);
+			temp.Sub(temp2);
+			// temp = m^2 - m
+
+			denominator.Mul(temp);
+			// denominator = (m^2 - m) * n ^ (m-1)
+
+			if( m >= cgamma.bern.size() )
+			{
+				if( !SetBernoulliNumbers(cgamma, m - cgamma.bern.size() + 1 + 3, stop) ) // 3 more than needed
+				{
+					// there was the stop signal
+					err = err_interrupt;
+					return ValueType(); // NaN
+			}
+			}
+			
+			temp = cgamma.bern[m];
+			temp.Div(denominator);
+
+			oldsum = sum;
+			sum.Add(temp);
+
+			if( sum.IsNan() || oldsum==sum )
+				break;
+		}
+
+	return sum;
+	}
+
+
+	/*!
+		an auxiliary function used to calculate the Gamma() function
+
+		we calculate a helper function GammaFactorialHigh() by using Stirling's series:
+		   n! = (n/e)^n * sqrt(2*pi*n) * exp( sum(n) )
+		   where n is a real number (not only an integer) and is sufficient large (greater than TTMATH_GAMMA_BOUNDARY)
+		   and sum(n) is calculated by GammaFactorialHighSum()
+	*/
+	template<class ValueType>
+	ValueType GammaFactorialHigh(const ValueType & n, CGamma<ValueType> & cgamma, ErrorCode & err,
+								 const volatile StopCalculating * stop)
+	{
+	ValueType temp, temp2, temp3, denominator, sum;
+
+		temp.Set2Pi();
+		temp.Mul(n);
+		temp2 = Sqrt(temp);
+		// temp2 = sqrt(2*pi*n)
+
+		temp = n;
+		temp3.SetE();
+		temp.Div(temp3);
+		temp.Pow(n);
+		// temp = (n/e)^n
+
+		sum = GammaFactorialHighSum(n, cgamma, err, stop);
+		temp3.Exp(sum);
+		// temp3 = exp(sum)
+
+		temp.Mul(temp2);
+		temp.Mul(temp3);
+
+	return temp;
+	}
+
+
+	/*!
+		an auxiliary function used to calculate the Gamma() function
+	
+		Gamma(x) = GammaFactorialHigh(x-1)
+	*/
+	template<class ValueType>
+	ValueType GammaPlusHigh(ValueType n, CGamma<ValueType> & cgamma, ErrorCode & err, const volatile StopCalculating * stop)
+	{
+	ValueType one;
+
+		one.SetOne();
+		n.Sub(one);
+
+	return GammaFactorialHigh(n, cgamma, err, stop);
+	}
+
+
+	/*!
+		an auxiliary function used to calculate the Gamma() function
+	
+		we use this function when n is integer and a small value (from 0 to TTMATH_GAMMA_BOUNDARY]
+		we use the formula:
+		   gamma(n) = (n-1)! = 1 * 2 * 3 * ... * (n-1) 
+	*/
+	template<class ValueType>
+	ValueType GammaPlusLowIntegerInt(uint n, CGamma<ValueType> & cgamma)
+	{
+		TTMATH_ASSERT( n > 0 )
+
+		if( n - 1 < static_cast<uint>(cgamma.fact.size()) )
+			return cgamma.fact[n - 1];
+
+		ValueType res;
+		uint start = 2;
+
+		if( cgamma.fact.size() < 2 )
+		{
+			res.SetOne();
+		}
+		else
+		{
+			start = static_cast<uint>(cgamma.fact.size());
+			res   = cgamma.fact[start-1];
+		}
+
+		for(uint i=start ; i<n ; ++i)
+			res.MulInt(i);
+
+	return res;
+	}
+	
+
+	/*!
+		an auxiliary function used to calculate the Gamma() function
+
+		we use this function when n is integer and a small value (from 0 to TTMATH_GAMMA_BOUNDARY]
+	*/
+	template<class ValueType>
+	ValueType GammaPlusLowInteger(const ValueType & n, CGamma<ValueType> & cgamma)
+	{
+	sint n_;
+
+		n.ToInt(n_);
+
+	return GammaPlusLowIntegerInt(n_, cgamma);
+	}
+
+
+	/*!
+		an auxiliary function used to calculate the Gamma() function
+
+		we use this function when n is a small value (from 0 to TTMATH_GAMMA_BOUNDARY]
+		we use a recurrence formula:
+		   gamma(z+1) = z * gamma(z)
+		   then: gamma(z) = gamma(z+1) / z
+
+		e.g.
+		   gamma(3.89) = gamma(2001.89) / ( 3.89 * 4.89 * 5.89 * ... * 1999.89 * 2000.89 )
+	*/
+	template<class ValueType>
+	ValueType GammaPlusLow(ValueType n, CGamma<ValueType> & cgamma, ErrorCode & err, const volatile StopCalculating * stop)
+	{
+	ValueType one, denominator, temp, boundary;
+
+		if( n.IsInteger() )
+			return GammaPlusLowInteger(n, cgamma);
+
+		one.SetOne();
+		denominator = n;
+		boundary    = TTMATH_GAMMA_BOUNDARY;
+
+		while( n < boundary )
+		{
+			n.Add(one);
+			denominator.Mul(n);
+		}
+
+		n.Add(one);
+
+		// now n is sufficient big
+		temp = GammaPlusHigh(n, cgamma, err, stop);
+		temp.Div(denominator);
+
+	return temp;
+		}
+
+
+	/*!
+		an auxiliary function used to calculate the Gamma() function
+	*/
+	template<class ValueType>
+	ValueType GammaPlus(const ValueType & n, CGamma<ValueType> & cgamma, ErrorCode & err, const volatile StopCalculating * stop)
+		{
+		if( n > TTMATH_GAMMA_BOUNDARY )
+			return GammaPlusHigh(n, cgamma, err, stop);
+
+	return GammaPlusLow(n, cgamma, err, stop);
+	}
+
+
+	/*!
+		an auxiliary function used to calculate the Gamma() function
+
+		this function is used when n is negative
+		we use the reflection formula:
+		   gamma(1-z) * gamma(z) = pi / sin(pi*z)
+		   then: gamma(z) = pi / (sin(pi*z) * gamma(1-z))
+
+	*/
+	template<class ValueType>
+	ValueType GammaMinus(const ValueType & n, CGamma<ValueType> & cgamma, ErrorCode & err, const volatile StopCalculating * stop)
+	{
+	ValueType pi, denominator, temp, temp2;
+
+		if( n.IsInteger() )
+		{
+			// gamma function is not defined when n is negative and integer
+			err = err_improper_argument;
+			return temp; // NaN
+		}
+
+		pi.SetPi();
+
+		temp = pi;
+		temp.Mul(n);
+		temp2 = Sin(temp);
+		// temp2 = sin(pi * n)
+
+		temp.SetOne();
+		temp.Sub(n);
+		temp = GammaPlus(temp, cgamma, err, stop);
+		// temp = gamma(1 - n)
+
+		temp.Mul(temp2);
+		pi.Div(temp);
+
+	return pi;
+	}
+
+	} // namespace auxiliaryfunctions
+
+
+
+	/*!
+		this function calculates the Gamma function
+
+		it's multithread safe, you should create a CGamma<> object and use it whenever you call the Gamma()
+		e.g.
+			typedef Big<1,2> MyBig;
+			MyBig x=234, y=345.53;
+			CGamma<MyBig> cgamma;
+			std::cout << Gamma(x, cgamma) << std::endl;
+			std::cout << Gamma(y, cgamma) << std::endl;
+		in the CGamma<> object the function stores some coefficients (factorials, Bernoulli numbers),
+		and they will be reused in next calls to the function
+
+		each thread should have its own CGamma<> object, and you can use these objects with Factorial() function too
+	*/
+	template<class ValueType>
+	ValueType Gamma(const ValueType & n, CGamma<ValueType> & cgamma, ErrorCode * err = 0,
+					const volatile StopCalculating * stop = 0)
+	{
+	using namespace auxiliaryfunctions;
+
+	ValueType result;
+	ErrorCode err_tmp;
+
+		
+		if( n.IsNan() )
+		{
+			if( err )
+				*err = err_improper_argument;
+
+		return result; // NaN is set by default
+		}
+
+		TTMATH_USE_THREADSAFE_OBJ(cgamma.history);
+
+		if( cgamma.history.Get(n, result, err_tmp) )
+		{
+			if( err )
+				*err = err_tmp;
+
+			return result;
+		}
+
+		err_tmp = err_ok;
+
+		if( n.IsSign() )
+		{
+			result = GammaMinus(n, cgamma, err_tmp, stop);
+		}
+		else
+		if( n.IsZero() )
+		{
+			err_tmp = err_improper_argument;
+			result.SetNan();
+		}
+		else
+		{
+			result = GammaPlus(n, cgamma, err_tmp, stop);
+		}
+
+		if( result.IsNan() && err_tmp==err_ok )
+			err_tmp = err_overflow;
+
+		if( err )
+			*err = err_tmp;
+
+		if( stop && !stop->WasStopSignal() )
+			cgamma.history.Add(n, result, err_tmp);
+
+	return result;
+	}
+
+
+	/*!
+		this function calculates the Gamma function
+
+		note: this function should be used only in a single-thread environment
+	*/
+	template<class ValueType>
+	ValueType Gamma(const ValueType & n, ErrorCode * err = 0)
+	{
+	// warning: this static object is not thread safe
+	static CGamma<ValueType> cgamma;
+
+	return Gamma(n, cgamma, err);
+	}
+
+
+
+	namespace auxiliaryfunctions
+	{
+
+	/*!
+		an auxiliary function for calculating the factorial function
+
+		we use the formula:
+		   x! = gamma(x+1)
+	*/
+	template<class ValueType>
+	ValueType Factorial2(ValueType x, CGamma<ValueType> * cgamma = 0, ErrorCode * err = 0,
+						 const volatile StopCalculating * stop = 0)
+	{
+	ValueType result, one;
+
+		if( x.IsNan() || x.IsSign() || !x.IsInteger() )
+		{
+			if( err )
+				*err = err_improper_argument;
+
+		return result; // NaN set by default
+		}
+
+		one.SetOne();
+		x.Add(one);
+
+		if( cgamma )
+			return Gamma(x, *cgamma, err, stop);
+
+	return Gamma(x, err);
+	}
+	
+	} // namespace auxiliaryfunctions
+
+
+
+	/*!
+		the factorial from given 'x'
+		e.g.
+		Factorial(4) = 4! = 1*2*3*4
+
+		it's multithread safe, you should create a CGamma<> object and use it whenever you call the Factorial()
+		e.g.
+			typedef Big<1,2> MyBig;
+			MyBig x=234, y=345.53;
+			CGamma<MyBig> cgamma;
+			std::cout << Factorial(x, cgamma) << std::endl;
+			std::cout << Factorial(y, cgamma) << std::endl;
+		in the CGamma<> object the function stores some coefficients (factorials, Bernoulli numbers),
+		and they will be reused in next calls to the function
+
+		each thread should have its own CGamma<> object, and you can use these objects with Gamma() function too
+	*/
+	template<class ValueType>
+	ValueType Factorial(const ValueType & x, CGamma<ValueType> & cgamma, ErrorCode * err = 0,
+						const volatile StopCalculating * stop = 0)
+	{
+		return auxiliaryfunctions::Factorial2(x, &cgamma, err, stop);
+	}
+
+
+	/*!
+		the factorial from given 'x'
+		e.g.
+		Factorial(4) = 4! = 1*2*3*4
+
+		note: this function should be used only in a single-thread environment
+	*/
+	template<class ValueType>
+	ValueType Factorial(const ValueType & x, ErrorCode * err = 0)
+	{
+		return auxiliaryfunctions::Factorial2(x, 0, err, 0);
+	}
+
+
+	/*!
+		this method prepares some coefficients: factorials and Bernoulli numbers
+		stored in 'fact' and 'bern' objects
+
+		we're defining the method here because we're using Gamma() function which
+		is not available in ttmathobjects.h
+
+		read the doc info in ttmathobjects.h file where CGamma<> struct is declared
+	*/
+	template<class ValueType>
+	void CGamma<ValueType>::InitAll()
+	{
+		ValueType x = TTMATH_GAMMA_BOUNDARY + 1;
+		
+		// history.Remove(x) removes only one object
+		// we must be sure that there are not others objects with the key 'x'
+		while( history.Remove(x) )
+		{
+		}
+
+		// the simplest way to initialize is to call the Gamma function with (TTMATH_GAMMA_BOUNDARY + 1)
+		// when x is larger then less coefficients we need
+		Gamma(x, *this);
+	}
+
+
+
 } // namespace
 
 

ttmath/ttmathbig.h

@@ -1,7 +1,7 @@
 /*
  * This file is a part of TTMath Bignum Library
  * and is distributed under the (new) BSD licence.
- * Author: Tomasz Sowa <t.sowa@slimaczek.pl>
+ * Author: Tomasz Sowa <t.sowa@ttmath.org>
  */
 
 /* 
@@ -38,6 +38,8 @@
 #ifndef headerfilettmathbig
 #define headerfilettmathbig
 
+#include "ttmathconfig.h"
+
 /*!
 	\file ttmathbig.h
     \brief A Class for representing floating point numbers
@@ -47,6 +49,10 @@
 
 #include <iostream>
 
+#if defined(_MSC_VER)
+	#pragma warning(disable:4127) // conditional expression is constant
+#endif
+
 namespace ttmath
 {
 
@@ -81,7 +87,7 @@ public:
 
 Int<exp>	exponent;
 UInt<man>	mantissa;
-unsigned char info;
+tt_char		info;
 
 
 /*!
@@ -267,7 +273,7 @@ private:
 		// 3101 digits were taken from this website
 		//  (later the digits were compared with:
 		//   http://www.eveandersson.com/pi/digits/1000000 and http://www.geom.uiuc.edu/~huberty/math5337/groupe/digits.html )
-		// and they were set into Big<1,400> type (using operator=(const char*) on a 32bit platform)
+		// and they were set into Big<1,400> type (using operator=(const tt_char*) on a 32bit platform)
 		// and then the first 256 words were taken into this table
 		// (TTMATH_BUILTIN_VARIABLES_SIZE on 32bit platform should have the value 256,
 		// and on 64bit platform value 128 (256/2=128))
@@ -420,7 +426,7 @@ public:
 		// (TTMATH_BUILTIN_VARIABLES_SIZE on 32bit platform should have the value 256,
 		// and on 64bit platform value 128 (256/2=128))
 
-		mantissa.SetFromTable(temp_table, sizeof(temp_table) / sizeof(int));
+		mantissa.SetFromTable(temp_table, sizeof(temp_table) / sizeof(unsigned int));
 		exponent = -sint(man)*sint(TTMATH_BITS_PER_UINT);
 		info = 0;
 	}
@@ -1028,7 +1034,7 @@ public:
 
 	UInt<man*2> man1;
 	UInt<man*2> man2;
-	uint i,c = 0;
+	uint i,c;
 		
 		if( IsNan() || ss2.IsNan() || ss2.IsZero() )
 			return CheckCarry(1);
@@ -1049,9 +1055,7 @@ public:
 
 		i = man1.CompensationToLeft();
 
-		if( i )
+		c  = exponent.Sub(i);
-			c += exponent.Sub(i);
-
 		c += exponent.Sub(ss2.exponent);
 		
 		for(i=0 ; i<man ; ++i)
@@ -1073,7 +1077,7 @@ public:
 
 		e.g.
 		 12.6 mod  3 =  0.6   because 12.6 = 3*4 + 0.6
-		-12.6 mod  3 = -0.6   bacause -12.6 = 3*(-4) + (-0.6)
+		-12.6 mod  3 = -0.6
 		 12.6 mod -3 =  0.6
 		-12.6 mod -3 = -0.6
 
@@ -1107,7 +1111,6 @@ public:
 
 
 
-
 	/*!
 		power this = this ^ pow
 		(pow without a sign)
@@ -1368,10 +1371,7 @@ public:
 
 		if( pow.exponent>-int(man*TTMATH_BITS_PER_UINT) && pow.exponent<=0 )
 		{
-			Big<exp, man> pow_frac( pow );
+			if( pow.IsInteger() )
-			pow_frac.RemainFraction();
-
-			if( pow_frac.IsZero() )
 				return PowInt( pow );
 		}
 
@@ -1769,7 +1769,7 @@ public:
 		// MS Visual Express 2005 reports a warning (in the lines with 'uint man_diff = ...'):
 		// warning C4307: '*' : integral constant overflow
 		// but we're using 'if( man > another_man )' and 'if( man < another_man )' and there'll be no such situation here
-		#ifdef _MSC_VER
+		#ifndef __GNUC__
 		#pragma warning( disable: 4307 )
 		#endif
 
@@ -1785,7 +1785,7 @@ public:
 			c += exponent.AddInt(man_diff, 0);
 		}
 		
-		#ifdef _MSC_VER
+		#ifndef __GNUC__
 		#pragma warning( default: 4307 )
 		#endif
 
@@ -2010,7 +2010,7 @@ public:
 		// error but I leave it at the moment as is
 		TTMATH_ASSERT( sizeof(double) == 8 )
 
-		// I am not sure what will be on a platform which has 
+		// I am not sure what will be on a plaltform which has 
 		// a different endianness... but we use this library only
 		// on x86 and amd (intel) 64 bits (as there's a lot of assembler code)
 		union 
@@ -2175,7 +2175,7 @@ public:
 				// then V=(-1)**S * 2 ** (-1022) * (0.F)
 				// These are "unnormalized" values.
 
-				FromDouble_SetExpAndMan(bool(temp.u & 0x8000000000000000ul),
+				FromDouble_SetExpAndMan((temp.u & 0x8000000000000000ul) != 0,
 										e - 1022 - man*TTMATH_BITS_PER_UINT + 1, 0, m);
 				Standardizing();
 			}
@@ -2670,6 +2670,18 @@ public:
 		operator=(value);
 	}
 
+	class LogHistory
+		{
+		public:
+			Big<exp,man>	val[15];
+			
+							LogHistory()
+								{
+								for (int i = 0; i < 15; ++i)
+									val[i].SetZero();
+								}
+		TTMATH_IMPLEMENT_THREADSAFE_OBJ
+		};
 
 	/*!
 		a method for converting the value into a string with a base equal 'base'
@@ -2705,20 +2717,20 @@ public:
 
 		output:
 			return value:
-			0 - ok and 'result' will be an object of type std::string which holds the value
+			0 - ok and 'result' will be an object of type std::string (or std::wstring) which holds the value
 			1 - if there was a carry (shoudn't be in a normal situation - if is that means there
 			    is somewhere an error in the library)
 	*/
-	uint ToString(	std::string & result,
+	uint ToString(	tt_string & result,
 					uint base                  = 10,
 					bool always_scientific     = false,
 					sint when_scientific       = 15,
 					sint max_digit_after_comma = -1,
 					bool remove_trailing_zeroes = true,
-					char decimal_point = TTMATH_COMMA_CHARACTER_1 ) const
+					tt_char decimal_point = TTMATH_COMMA_CHARACTER_1 ) const
 	{										  
-		static char error_overflow_msg[] = "overflow";
+		static tt_char error_overflow_msg[] = TTMATH_TEXT("overflow");
-		static char error_nan_msg[]      = "NaN";
+		static tt_char error_nan_msg[]      = TTMATH_TEXT("NaN");
 		result.erase();
 		
 		if( IsNan() )
@@ -2735,7 +2747,7 @@ public:
 	
 		if( IsZero() )
 		{
-			result = "0";
+			result = '0';
 
 		return 0;
 		}
@@ -2860,7 +2872,7 @@ private:
 			but we need 'new'exp' as integer then we take:
 			new_exp = [log base (2^exponent)] + 1  <- where [x] means integer value from x
 	*/
-	uint ToString_CreateNewMantissaAndExponent(	std::string & new_man, uint base,
+	uint ToString_CreateNewMantissaAndExponent(	tt_string & new_man, uint base,
 												Int<exp+1> & new_exp) const
 	{
 	uint c = 0;
@@ -2964,11 +2976,12 @@ private:
 		// (LnSurrounding1() will return one immediately)
 		uint c = Ln(x);
 
-		// warning! this 'static' is not thread safe
-		static Big<exp,man> log_history[15] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
 		uint index = base - 2;
 
-		if( log_history[index].IsZero() )
+		static LogHistory	log_history;
+		TTMATH_USE_THREADSAFE_OBJ(log_history);
+		
+		if( log_history.val[index].IsZero() )
 		{	
 			// we don't have 'base' in 'log_history' then we calculate it now
 
@@ -2986,14 +2999,14 @@ private:
 
 			// the next time we'll get the 'Ln(base)' from the history,
 			// this 'log_history' can have (16-2+1) items max
-			log_history[index] = temp;
+			log_history.val[index] = temp;
 
 			c += Div(temp);
 		}
 		else
 		{
 			// we've calculated the 'Ln(base)' beforehand and we're getting it now
-			c += Div( log_history[index] );
+			c += Div( log_history.val[index] );
 		}
 
 	return (c==0)? 0 : 1;
@@ -3043,7 +3056,7 @@ private:
 		(we can make that speciality when the base is 4,8 or 16 as well 
 		but maybe in further time)
 	*/
-	uint ToString_CreateNewMantissaAndExponent_Base2(	std::string & new_man,
+	uint ToString_CreateNewMantissaAndExponent_Base2(	tt_string & new_man,
 														Int<exp+1> & new_exp     ) const
 	{
 		for( sint i=man-1 ; i>=0 ; --i )
@@ -3073,13 +3086,13 @@ private:
 		this method roundes the last character from the new mantissa
 		(it's used in systems where the base is different from 2)
 	*/
-	uint ToString_RoundMantissa(std::string & new_man, uint base, Int<exp+1> & new_exp, char decimal_point) const
+	uint ToString_RoundMantissa(tt_string & new_man, uint base, Int<exp+1> & new_exp, tt_char decimal_point) const
 	{
 		// we must have minimum two characters
 		if( new_man.length() < 2 )
 			return 0;
 
-		std::string::size_type i = new_man.length() - 1;
+		tt_string::size_type i = new_man.length() - 1;
 
 		// we're erasing the last character
 		uint digit = UInt<man>::CharToDigit( new_man[i] );
@@ -3100,7 +3113,7 @@ private:
 
 		this method addes one into the new mantissa
 	*/
-	void ToString_RoundMantissa_AddOneIntoMantissa(std::string & new_man, uint base, char decimal_point) const
+	void ToString_RoundMantissa_AddOneIntoMantissa(tt_string & new_man, uint base, tt_char decimal_point) const
 	{
 		if( new_man.empty() )
 			return;
@@ -3124,7 +3137,7 @@ private:
 			else
 				was_carry = false;
 
-			new_man[i] = static_cast<char>( UInt<man>::DigitToChar(digit) );
+			new_man[i] = static_cast<tt_char>( UInt<man>::DigitToChar(digit) );
 		}
 
 		if( i<0 && was_carry )
@@ -3138,13 +3151,13 @@ private:
 		this method sets the comma operator and/or puts the exponent
 		into the string
 	*/
-	uint ToString_SetCommaAndExponent(	std::string & new_man, uint base,
+	uint ToString_SetCommaAndExponent(	tt_string & new_man, uint base,
 										Int<exp+1> & new_exp,
 										bool always_scientific,
 										sint when_scientific,
 										sint max_digit_after_comma,
 										bool remove_trailing_zeroes,
-										char decimal_point) const
+										tt_char decimal_point) const
 	{
 	uint carry = 0;
 
@@ -3184,12 +3197,12 @@ private:
 		an auxiliary method for converting into the string
 	*/
 	void ToString_SetCommaAndExponent_Normal(
-											std::string & new_man,
+											tt_string & new_man,
 											uint base, 
 											Int<exp+1> & new_exp,
 											sint max_digit_after_comma,
 											bool remove_trailing_zeroes,
-											char decimal_point) const
+											tt_char decimal_point) const
 	{
 		if( !new_exp.IsSign() ) //if( new_exp >= 0 )
 			return ToString_SetCommaAndExponent_Normal_AddingZero(new_man, new_exp);
@@ -3201,7 +3214,7 @@ private:
 	/*!
 		an auxiliary method for converting into the string
 	*/
-	void ToString_SetCommaAndExponent_Normal_AddingZero(std::string & new_man,
+	void ToString_SetCommaAndExponent_Normal_AddingZero(tt_string & new_man,
 																Int<exp+1> & new_exp) const
 	{
 		// we're adding zero characters at the end
@@ -3221,12 +3234,12 @@ private:
 		an auxiliary method for converting into the string
 	*/
 	void ToString_SetCommaAndExponent_Normal_SetCommaInside(
-														std::string & new_man,
+														tt_string & new_man,
 														uint base,
 														Int<exp+1> & new_exp,
 														sint max_digit_after_comma,
 														bool remove_trailing_zeroes,
-														char decimal_point) const
+														tt_char decimal_point) const
 	{
 		// new_exp is < 0 
 
@@ -3245,7 +3258,7 @@ private:
 			// we're adding zero characters before the mantissa
 
 			uint how_many = e - new_man_len;
-			std::string man_temp(how_many+1, '0');
+			tt_string man_temp(how_many+1, '0');
 
 			man_temp.insert( man_temp.begin()+1, decimal_point);
 			new_man.insert(0, man_temp);
@@ -3258,12 +3271,12 @@ private:
 	/*!
 		an auxiliary method for converting into the string
 	*/
-	void ToString_SetCommaAndExponent_Scientific(	std::string & new_man,
+	void ToString_SetCommaAndExponent_Scientific(	tt_string & new_man,
 													uint base,
 													Int<exp+1> & scientific_exp,
 													sint max_digit_after_comma,
 													bool remove_trailing_zeroes,
-													char decimal_point) const
+													tt_char decimal_point) const
 	{
 		if( new_man.empty() )
 			return;
@@ -3277,16 +3290,16 @@ private:
 			new_man += 'e';
 
 			if( !scientific_exp.IsSign() )
-				new_man += "+";
+				new_man += '+';
 		}
 		else
 		{
 			// the 10 here is meant as the base 'base'
 			// (no matter which 'base' we're using there'll always be 10 here)
-			new_man += "*10^";
+			new_man += TTMATH_TEXT("*10^");
 		}
 
-		std::string temp_exp;
+		tt_string temp_exp;
 		scientific_exp.ToString( temp_exp, base );
 
 		new_man += temp_exp;
@@ -3296,11 +3309,11 @@ private:
 	/*!
 		an auxiliary method for converting into the string
 	*/
-	void ToString_CorrectDigitsAfterComma(	std::string & new_man,
+	void ToString_CorrectDigitsAfterComma(	tt_string & new_man,
 											uint base,
 											sint max_digit_after_comma,
 											bool remove_trailing_zeroes,
-											char decimal_point) const
+											tt_char decimal_point) const
 	{
 		if( max_digit_after_comma >= 0 )
 			ToString_CorrectDigitsAfterComma_Round(new_man, base, max_digit_after_comma, decimal_point);
@@ -3314,8 +3327,8 @@ private:
 		an auxiliary method for converting into the string
 	*/
 	void ToString_CorrectDigitsAfterComma_CutOffZeroCharacters(
-												std::string & new_man,
+												tt_string & new_man,
-												char decimal_point) const
+												tt_char decimal_point) const
 	{
 		// minimum two characters
 		if( new_man.length() < 2 )
@@ -3333,7 +3346,7 @@ private:
 		// we must have a comma 
 		// (the comma can be removed by ToString_CorrectDigitsAfterComma_Round
 		// which is called before)
-		if( new_man.find_last_of(decimal_point, i) == std::string::npos )
+		if( new_man.find_last_of(decimal_point, i) == tt_string::npos )
 			return;
 
 		// if directly before the first zero is the comma operator
@@ -3349,26 +3362,26 @@ private:
 		an auxiliary method for converting into the string
 	*/
 	void ToString_CorrectDigitsAfterComma_Round(
-											std::string & new_man,
+											tt_string & new_man,
 											uint base,
 											sint max_digit_after_comma,
-											char decimal_point) const
+											tt_char decimal_point) const
 	{
 		// first we're looking for the comma operator
-		std::string::size_type index = new_man.find(decimal_point, 0);
+		tt_string::size_type index = new_man.find(decimal_point, 0);
 
-		if( index == std::string::npos )
+		if( index == tt_string::npos )
 			// nothing was found (actually there can't be this situation)
 			return;
 
 		// we're calculating how many digits there are at the end (after the comma)
 		// 'after_comma' will be greater than zero because at the end
 		// we have at least one digit
-		std::string::size_type after_comma = new_man.length() - index - 1;
+		tt_string::size_type after_comma = new_man.length() - index - 1;
 
 		// if 'max_digit_after_comma' is greater than 'after_comma' (or equal)
 		// we don't have anything for cutting
-		if( std::string::size_type(max_digit_after_comma) >= after_comma )
+		if( tt_string::size_type(max_digit_after_comma) >= after_comma )
 			return;
 
 		uint last_digit = UInt<man>::CharToDigit( new_man[ index + max_digit_after_comma + 1 ], base );
@@ -3404,6 +3417,7 @@ public:
 		all digits after the comma we can ignore
 
 		'source' - pointer to the string for parsing
+		           'const char*' or 'const wchar_t*'
 
 		if 'after_source' is set that when this method finishes
 		it sets the pointer to the new first character after parsed value
@@ -3414,7 +3428,7 @@ public:
 		no value has been read (there are no digits)
 		on other words if 'value_read' is true -- there is at least one digit in the string
 	*/
-	uint FromString(const char * source, uint base = 10, const char ** after_source = 0, bool * value_read = 0)
+	uint FromString(const tt_char * source, uint base = 10, const tt_char ** after_source = 0, bool * value_read = 0)
 	{
 	bool is_sign;
 	bool value_read_temp = false;
@@ -3465,7 +3479,7 @@ private:
 
 		(this method is used from 'FromString_ReadPartScientific' too)
 	*/
-	void FromString_TestSign( const char * & source, bool & is_sign )
+	void FromString_TestSign( const tt_char * & source, bool & is_sign )
 	{
 		UInt<man>::SkipWhiteCharacters(source);
 
@@ -3487,7 +3501,7 @@ private:
 	/*!
 		we're testing whether there's a comma operator
 	*/
-	bool FromString_TestCommaOperator(const char * & source)
+	bool FromString_TestCommaOperator(const tt_char * & source)
 	{
 		if( (*source == TTMATH_COMMA_CHARACTER_1) || 
 			(*source == TTMATH_COMMA_CHARACTER_2 && TTMATH_COMMA_CHARACTER_2 != 0 ) )
@@ -3505,7 +3519,7 @@ private:
 		this method reads the first part of a string
 		(before the comma operator)
 	*/
-	uint FromString_ReadPartBeforeComma( const char * & source, uint base, bool & value_read )
+	uint FromString_ReadPartBeforeComma( const tt_char * & source, uint base, bool & value_read )
 	{
 		sint character;
 		Big<exp, man> temp;
@@ -3534,7 +3548,7 @@ private:
 		this method reads the second part of a string
 		(after the comma operator)
 	*/
-	uint FromString_ReadPartAfterComma( const char * & source, uint base, bool & value_read )
+	uint FromString_ReadPartAfterComma( const tt_char * & source, uint base, bool & value_read )
 	{
 	sint character;
 	uint c = 0, index = 1;
@@ -3592,12 +3606,12 @@ private:
 
 		it is called when the base is 10 and some digits were read before
 	*/
-	uint FromString_ReadScientificIfExists(const char * & source)
+	uint FromString_ReadScientificIfExists(const tt_char * & source)
 	{
 	uint c = 0;
 
 		bool scientific_read = false;
-		const char * before_scientific = source;
+		const tt_char * before_scientific = source;
 
 		if( FromString_TestScientific(source) )
 			c += FromString_ReadPartScientific( source, scientific_read );
@@ -3615,7 +3629,7 @@ private:
 
 		this character is only allowed when we're using the base equals 10
 	*/
-	bool FromString_TestScientific(const char * & source)
+	bool FromString_TestScientific(const tt_char * & source)
 	{
 		UInt<man>::SkipWhiteCharacters(source);
 
@@ -3634,7 +3648,7 @@ private:
 		this method reads the exponent (after 'e' character) when there's a scientific
 		format of value and only when we're using the base equals 10
 	*/
-	uint FromString_ReadPartScientific( const char * & source, bool & scientific_read )
+	uint FromString_ReadPartScientific( const tt_char * & source, bool & scientific_read )
 	{
 	uint c = 0;
 	Big<exp, man> new_exponent, temp;
@@ -3661,7 +3675,7 @@ private:
 		this method reads the value of the extra exponent when scientific format is used
 		(only when base == 10)
 	*/
-	uint FromString_ReadPartScientific_ReadExponent( const char * & source, Big<exp, man> & new_exponent, bool & scientific_read )
+	uint FromString_ReadPartScientific_ReadExponent( const tt_char * & source, Big<exp, man> & new_exponent, bool & scientific_read )
 	{
 	sint character;
 	Big<exp, man> base, temp;
@@ -3694,7 +3708,7 @@ public:
 	/*!
 		a method for converting a string into its value		
 	*/
-	uint FromString(const std::string & string, uint base = 10)
+	uint FromString(const tt_string & string, uint base = 10)
 	{
 		return FromString( string.c_str(), base );
 	}
@@ -3703,7 +3717,7 @@ public:
 	/*!
 		a constructor for converting a string into this class
 	*/
-	Big(const char * string)
+	Big(const tt_char * string)
 	{
 		FromString( string );
 	}
@@ -3712,7 +3726,7 @@ public:
 	/*!
 		a constructor for converting a string into this class
 	*/
-	Big(const std::string & string)
+	Big(const tt_string & string)
 	{
 		FromString( string.c_str() );
 	}
@@ -3721,7 +3735,7 @@ public:
 	/*!
 		an operator= for converting a string into its value
 	*/
-	Big<exp, man> & operator=(const char * string)
+	Big<exp, man> & operator=(const tt_char * string)
 	{
 		FromString( string );
 
@@ -3732,7 +3746,7 @@ public:
 	/*!
 		an operator= for converting a string into its value
 	*/
-	Big<exp, man> & operator=(const std::string & string)
+	Big<exp, man> & operator=(const tt_string & string)
 	{
 		FromString( string.c_str() );
 
@@ -3850,6 +3864,72 @@ public:
 	return false;
 	}
 
+	bool AboutEqual(const Big<exp,man> & ss2, int nBitsToIgnore = 4) const
+	{
+		// we should check the mantissas beforehand because sometimes we can have
+		// a mantissa set to zero but in the exponent something another value
+		// (maybe we've forgotten about calling CorrectZero() ?)
+		if( mantissa.IsZero())
+			{
+			if (ss2.mantissa.IsZero())
+				return true;
+			return(ss2.AboutEqual(*this,nBitsToIgnore));
+			}
+		
+		if (ss2.mantissa.IsZero())
+			{
+			return(this->exponent <= uint(2*(-sint(man*TTMATH_BITS_PER_UINT))+nBitsToIgnore));
+			}
+			
+		// exponents may not differ much!
+		ttmath::Int<exp>	expdiff(this->exponent - ss2.exponent);
+		
+		// they may differ one if for example mantissa1=0x80000000, mantissa2=0xffffffff
+		if (ttmath::Abs(expdiff) > 1)
+			return(false);		
+
+		// calculate the 'difference' mantissa		
+		ttmath::UInt<man>	man1(this->mantissa);
+		ttmath::UInt<man>	man2(ss2.mantissa);
+		ttmath::UInt<man>	mandiff;
+		
+		switch (expdiff.ToInt())
+			{
+			case +1:
+				man2.Rcr(1,0);
+				mandiff = man1;
+				mandiff.Sub(man2);
+				break;
+			case -1:
+				man1.Rcr(1,0);
+				mandiff = man2;
+				mandiff.Sub(man1);
+				break;
+			default:
+				if (man2 > man1)
+					{
+					mandiff = man2;
+					mandiff.Sub(man1);
+					}
+				  else
+					{
+					mandiff = man1;
+					mandiff.Sub(man2);
+					}
+				break;
+			}
+			
+		// faster to mask the bits!
+		ASSERT(nBitsToIgnore < TTMATH_BITS_PER_UINT);
+
+		for (int n = man-1; n > 0; --n)
+			{
+			if (mandiff.table[n] != 0)
+				return(false);
+			}
+		uint	nMask = ~((1 << nBitsToIgnore) - 1);
+		return((mandiff.table[0] & nMask) == 0);
+	}
 
 	bool operator<(const Big<exp,man> & ss2) const
 	{
@@ -3881,8 +3961,6 @@ public:
 	}
 
 
-
-
 	bool operator>(const Big<exp,man> & ss2) const
 	{
 		if( IsSign() && !ss2.IsSign() )
@@ -4097,6 +4175,47 @@ public:
 
 
 
+	/*!
+		this method returns true if the value is integer
+		(there is no a fraction)
+
+		(we don't check nan)
+	*/
+	bool IsInteger() const
+	{
+		if( IsZero() )
+			return true;
+
+		if( !exponent.IsSign() )
+			// exponent >=0 -- the value don't have any fractions
+			return true;
+
+		if( exponent <= -sint(man*TTMATH_BITS_PER_UINT) )
+			// the value is from (-1,1)
+			return false;
+
+		// exponent is in range (-man*TTMATH_BITS_PER_UINT, 0)
+		sint e = exponent.ToInt();
+		e = -e; // e means how many bits we must check
+
+		uint len  = e / TTMATH_BITS_PER_UINT;
+		uint rest = e % TTMATH_BITS_PER_UINT;
+		uint i    = 0;
+
+		for( ; i<len ; ++i )
+			if( mantissa.table[i] != 0 )
+				return false;
+
+		if( rest > 0 )
+		{
+			uint rest_mask = TTMATH_UINT_MAX_VALUE >> (TTMATH_BITS_PER_UINT - rest);
+			if( (mantissa.table[i] & rest_mask) != 0 )
+				return false;
+		}
+
+	return true;
+	}
+
 
 	/*!
 		this method rounds to the nearest integer value
@@ -4145,9 +4264,14 @@ public:
 	*
 	*/
 
-	friend std::ostream & operator<<(std::ostream & s, const Big<exp,man> & l)
+	/*!
+		output for standard streams
+
+		tt_ostream is either std::ostream or std::wostream
+	*/
+	friend tt_ostream & operator<<(tt_ostream & s, const Big<exp,man> & l)
 	{
-	std::string ss;
+	tt_string ss;
 
 		l.ToString(ss);
 		s << ss;
@@ -4156,12 +4280,17 @@ public:
 	}
 
 
-	friend std::istream & operator>>(std::istream & s, Big<exp,man> & l)
+	/*!
-	{
+		input from standard streams
-	std::string ss;
 
-	// 'char' for operator>>
+		tt_istream is either std::istream or std::wistream
-	unsigned char z;
+	*/
+	friend tt_istream & operator>>(tt_istream & s, Big<exp,man> & l)
+	{
+	tt_string ss;
+	
+	// 'tt_char' for operator>>
+	tt_char z;
 	bool was_comma = false;
 
 		// operator>> omits white characters if they're set for ommiting
@@ -4174,7 +4303,7 @@ public:
 		}
 
 		// we're reading only digits (base=10) and only one comma operator
-		for( ; s.good() ; z=s.get() )
+		for( ; s.good() ; z=static_cast<tt_char>(s.get()) )
 		{
 			if( z == TTMATH_COMMA_CHARACTER_1 ||
 			  ( z == TTMATH_COMMA_CHARACTER_2 && TTMATH_COMMA_CHARACTER_2 != 0 ) )
@@ -4204,6 +4333,9 @@ public:
 
 };
 
+#if defined(_MSC_VER)
+	#pragma warning(default:4127) // conditional expression is constant
+#endif
 
 } // namespace
 

ttmath/ttmathconfig.h

@@ -0,0 +1,110 @@
+/*
+ * This file is a part of TTMath Bignum Library
+ * and is distributed under the PNG licence.
+ * Author: Christian Kaiser <chk@online.de>
+ */
+
+/* 
+	Copyright (c) 2009 Christian Kaiser
+
+	This software is provided 'as-is', without any express or implied
+	warranty. In no event will the authors be held liable for any damages
+	arising from the use of this software.
+
+	Permission is granted to anyone to use this software for any purpose,
+	including commercial applications, and to alter it and redistribute it
+	freely, subject to the following restrictions:
+
+		1. The origin of this software must not be misrepresented; you must not
+		claim that you wrote the original software. If you use this software
+		in a product, an acknowledgment in the product documentation would be
+		appreciated but is not required.
+
+		2. Altered source versions must be plainly marked as such, and must not be
+		misrepresented as being the original software.
+
+		3. This notice may not be removed or altered from any source
+		distribution.
+ */
+
+#ifndef headerfilettmathconfig
+#define headerfilettmathconfig
+#pragma once
+
+#include <sstream>
+
+namespace ttmath
+{
+#if defined(WIN32)
+	#include <windows.h>
+	
+	#if defined(_MT)
+		class 	clsCrit
+			{
+			private:
+				mutable CRITICAL_SECTION 		_Crit;
+
+												clsCrit(const clsCrit&) // inhibit copy (easy mistake to do; use clsCritObj instead!!!)
+													{
+													}
+				clsCrit&						operator=(const clsCrit& rhs); // inhibit assignment
+			public:
+												clsCrit(void)
+													{
+													InitializeCriticalSection(&_Crit);
+													}
+				virtual							~clsCrit(void)
+													{
+													DeleteCriticalSection(&_Crit);
+													}
+
+				void							Enter(void) const
+													{
+													EnterCriticalSection(&_Crit);
+													}
+				void							Leave(void) const
+													{
+													LeaveCriticalSection(&_Crit);
+													}
+			};
+
+		class 	clsCritObj
+			{
+			private:
+				const clsCrit&					_Crit;
+				
+				clsCritObj&						operator=(const clsCritObj& rhs); // not applicable
+			public:
+												clsCritObj(const clsCrit& Sync)
+													: _Crit(Sync)
+													{
+													_Crit.Enter();
+													}
+												~clsCritObj(void)
+													{
+													_Crit.Leave();
+													}
+			};
+		#define TTMATH_IMPLEMENT_THREADSAFE_OBJ					\
+			private:											\
+				clsCrit CritSect;								\
+			public:												\
+				operator clsCrit&()								\
+				{												\
+					return(CritSect);							\
+				}
+		#define TTMATH_USE_THREADSAFE_OBJ(c)	clsCritObj	lock(c)
+	#endif
+#else // defined(WIN32) 
+	// not Windows world: no threading synchronization for now
+#endif
+
+#if !defined(TTMATH_IMPLEMENT_THREADSAFE_OBJ)
+	// if we don't know about serialization, make it a no-op
+	#define TTMATH_IMPLEMENT_THREADSAFE_OBJ		/* */
+	#define TTMATH_USE_THREADSAFE_OBJ(c)		/* */
+#endif
+
+} // namespace
+
+#endif // headerfilettmathconfig

ttmath/ttmathint.h

@@ -1,7 +1,7 @@
 /*
  * This file is a part of TTMath Bignum Library
  * and is distributed under the (new) BSD licence.
- * Author: Tomasz Sowa <t.sowa@slimaczek.pl>
+ * Author: Tomasz Sowa <t.sowa@ttmath.org>
  */
 
 /* 
@@ -47,6 +47,10 @@
 
 #include "ttmathuint.h"
 
+#if defined(_MSC_VER)
+	#pragma warning(disable:4127) // conditional expression is constant
+#endif
+
 namespace ttmath
 {
 
@@ -641,8 +645,14 @@ public:
 
 		// there can be a carry here when the size of this value is equal one word
 		// and the 'value' has the highest bit set
+		#if defined(_MSC_VER)
+			#pragma warning(disable:4127) // conditional expression is constant
+		#endif
 		if( value_size==1 && (value & TTMATH_UINT_HIGHEST_BIT)!=0 )
 			return 1;
+		#if defined(_MSC_VER)
+			#pragma warning(default:4127) // conditional expression is constant
+		#endif
 
 	return 0;
 	}
@@ -821,7 +831,7 @@ public:
 	/*!
 		a constructor for converting string to this class (with the base=10)
 	*/
-	Int(const char * s)
+	Int(const tt_char * s)
 	{
 		FromString(s);
 	}
@@ -830,7 +840,7 @@ public:
 	/*!
 		a constructor for converting a string to this class (with the base=10)
 	*/
-	Int(const std::string & s)
+	Int(const tt_string & s)
 	{
 		FromString( s.c_str() );
 	}
@@ -869,7 +879,7 @@ public:
 	/*!	
 		this method converts the value to a string with a base equal 'b'
 	*/
-	void ToString(std::string & result, uint b = 10) const
+	void ToString(tt_string & result, uint b = 10) const
 	{
 		if( IsSign() )
 		{
@@ -890,12 +900,14 @@ public:
 
 	/*!
 		this method converts a string into its value
+		string is given either as 'const char *' or 'const wchar_t *'
+
 		it returns carry=1 if the value will be too big or an incorrect base 'b' is given
 
 		string is ended with a non-digit value, for example:
 			"-12" will be translated to -12
 			as well as:
-			"- 12foo" will be translated to 12 too
+			"- 12foo" will be translated to -12 too
 
 		existing first white characters will be ommited
 		(between '-' and a first digit can be white characters too)
@@ -904,7 +916,7 @@ public:
 
 		value_read (if exists) tells whether something has actually been read (at least one digit)
 	*/
-	uint FromString(const char * s, uint b = 10, const char ** after_source = 0, bool * value_read = 0)
+	uint FromString(const tt_char * s, uint b = 10, const tt_char ** after_source = 0, bool * value_read = 0)
 	{
 	bool is_sign = false;
 	
@@ -961,16 +973,16 @@ public:
 		this method converts a string into its value
 		it returns carry=1 if the value will be too big or an incorrect base 'b' is given
 	*/
-	uint FromString(const std::string & s, uint b = 10)
+	uint FromString(const tt_string & s, uint b = 10)
 	{
-		return FromString( s.c_str() );
+		return FromString( s.c_str(), b );
 	}
 
 
 	/*!
 		this operator converts a string into its value (with base = 10)
 	*/
-	Int<value_size> & operator=(const char * s)
+	Int<value_size> & operator=(const tt_char * s)
 	{
 		FromString(s);
 
@@ -981,7 +993,7 @@ public:
 	/*!
 		this operator converts a string into its value (with base = 10)
 	*/
-	Int<value_size> & operator=(const std::string & s)
+	Int<value_size> & operator=(const tt_string & s)
 	{
 		FromString( s.c_str() );
 
@@ -1268,9 +1280,14 @@ public:
 	*
 	*/
 
-	friend std::ostream & operator<<(std::ostream & s, const Int<value_size> & l)
+	/*!
+		output for standard streams
+
+		tt_ostream is either std::ostream or std::wostream
+	*/
+	friend tt_ostream & operator<<(tt_ostream & s, const Int<value_size> & l)
 	{
-	std::string ss;
+	tt_string ss;
 
 		l.ToString(ss);
 		s << ss;
@@ -1279,13 +1296,17 @@ public:
 	}
 
 
+	/*!
+		input from standard streams
 
-	friend std::istream & operator>>(std::istream & s, Int<value_size> & l)
+		tt_istream is either std::istream or std::wistream
+	*/
+	friend tt_istream & operator>>(tt_istream & s, Int<value_size> & l)
 	{
-	std::string ss;
+	tt_string ss;
 	
-	// char for operator>>
+	// tt_char for operator>>
-	unsigned char z;
+	tt_char z;
 	
 		// operator>> omits white characters if they're set for ommiting
 		s >> z;
@@ -1300,7 +1321,7 @@ public:
 		while( s.good() && UInt<value_size>::CharToDigit(z, 10)>=0 )
 		{
 			ss += z;
-			z = s.get();
+			z = static_cast<tt_char>(s.get());
 		}
 
 		// we're leaving the last readed character
@@ -1316,4 +1337,9 @@ public:
 
 } // namespace
 
+#if defined(_MSC_VER)
+	#pragma warning(default:4127) // conditional expression is constant
+#endif
+
+
 #endif

ttmath/ttmathobjects.h

@@ -1,7 +1,7 @@
 /*
  * This file is a part of TTMath Mathematical Library
  * and is distributed under the (new) BSD licence.
- * Author: Tomasz Sowa <t.sowa@slimaczek.pl>
+ * Author: Tomasz Sowa <t.sowa@ttmath.org>
  */
 
 /* 
@@ -47,6 +47,7 @@
 #include "ttmathtypes.h"
 
 #include <string>
+#include <vector>
 #include <list>
 #include <map>
 
@@ -73,18 +74,19 @@ public:
 	struct Item
 	{
 		// name of a variable of a function
-		std::string value;
+		// (either std::string or std::wstring)
+		tt_string value;
 
 		// number of parameters required by the function
 		// (if there's a variable this 'param' is ignored)
 		int param;
 
 		Item() {}
-		Item(const std::string & v, int p) : value(v), param(p) {}
+		Item(const tt_string & v, int p) : value(v), param(p) {}
 	};
 
 	// 'Table' is the type of our table
-	typedef std::map<std::string, Item> Table;
+	typedef std::map<tt_string, Item> Table;
 	typedef	Table::iterator Iterator;
 	typedef	Table::const_iterator CIterator;
 
@@ -112,7 +114,7 @@ public:
 	/*!
 		this method returns true if the name can be as a name of an object
 	*/
-	static bool IsNameCorrect(const std::string & name)
+	static bool IsNameCorrect(const tt_string & name)
 	{
 		if( name.empty() )
 			return false;
@@ -120,7 +122,7 @@ public:
 		if( !CorrectCharacter(name[0], false) )
 			return false;
 
-		std::string::const_iterator i=name.begin();
+		tt_string::const_iterator i=name.begin();
 
 		for(++i ; i!=name.end() ; ++i)
 			if( !CorrectCharacter(*i, true) )
@@ -133,7 +135,7 @@ public:
 	/*!
 		this method returns true if such an object is defined (name exists)
 	*/
-	bool IsDefined(const std::string & name)
+	bool IsDefined(const tt_string & name)
 	{
 		Iterator i = table.find(name);
 
@@ -148,7 +150,7 @@ public:
 	/*!
 		this method adds one object (variable of function) into the table
 	*/
-	ErrorCode Add(const std::string & name, const std::string & value, int param = 0)
+	ErrorCode Add(const tt_string & name, const tt_string & value, int param = 0)
 	{
 		if( !IsNameCorrect(name) )
 			return err_incorrect_name;
@@ -205,7 +207,7 @@ public:
 	/*!
 		this method changes the value and the number of parameters for a specific object
 	*/
-	ErrorCode EditValue(const std::string & name, const std::string & value, int param = 0)
+	ErrorCode EditValue(const tt_string & name, const tt_string & value, int param = 0)
 	{
 		if( !IsNameCorrect(name) )
 			return err_incorrect_name;
@@ -225,7 +227,7 @@ public:
 	/*!
 		this method changes the name of a specific object
 	*/
-	ErrorCode EditName(const std::string & old_name, const std::string & new_name)
+	ErrorCode EditName(const tt_string & old_name, const tt_string & new_name)
 	{
 		if( !IsNameCorrect(old_name) || !IsNameCorrect(new_name) )
 			return err_incorrect_name;
@@ -256,7 +258,7 @@ public:
 	/*!
 		this method deletes an object
 	*/
-	ErrorCode Delete(const std::string & name)
+	ErrorCode Delete(const tt_string & name)
 	{
 		if( !IsNameCorrect(name) )
 			return err_incorrect_name;
@@ -275,7 +277,7 @@ public:
 	/*!
 		this method gets the value of a specific object
 	*/
-	ErrorCode GetValue(const std::string & name, std::string & value) const
+	ErrorCode GetValue(const tt_string & name, tt_string & value) const
 	{
 		if( !IsNameCorrect(name) )
 			return err_incorrect_name;
@@ -298,7 +300,7 @@ public:
 		this method gets the value of a specific object
 		(this version is used for not copying the whole string)
 	*/
-	ErrorCode GetValue(const std::string & name, const char ** value) const
+	ErrorCode GetValue(const tt_string & name, const tt_char ** value) const
 	{
 		if( !IsNameCorrect(name) )
 			return err_incorrect_name;
@@ -321,7 +323,7 @@ public:
 		this method gets the value and the number of parameters
 		of a specific object
 	*/
-	ErrorCode GetValueAndParam(const std::string & name, std::string & value, int * param) const
+	ErrorCode GetValueAndParam(const tt_string & name, tt_string & value, int * param) const
 	{
 		if( !IsNameCorrect(name) )
 			return err_incorrect_name;
@@ -347,7 +349,7 @@ public:
 		of a specific object
 		(this version is used for not copying the whole string)
 	*/
-	ErrorCode GetValueAndParam(const std::string & name, const char ** value, int * param) const
+	ErrorCode GetValueAndParam(const tt_string & name, const tt_char ** value, int * param) const
 	{
 		if( !IsNameCorrect(name) )
 			return err_incorrect_name;
@@ -430,7 +432,7 @@ public:
 	*/
 	History()
 	{
-		buffer_max_size = 10;
+		buffer_max_size = 15;
 	}
 
 
@@ -487,10 +489,118 @@ public:
 	return false;
 	}
 
+
+	/*!
+		this methods deletes an item
+
+		we assume that there is only one item with the 'key'
+		(this methods removes the first one)
+	*/
+	bool Remove(const ValueType & key)
+	{
+		typename buffer_type::iterator i = buffer.begin();
+
+		for( ; i != buffer.end() ; ++i )
+		{
+			if( i->key == key )
+			{
+				buffer.erase(i);
+				return true;
+			}
+		}
+
+	return false;
+	}
+
+
 }; // end of class History
 
 
 
+/*!
+	this is an auxiliary class used when calculating Gamma() or Factorial()
+
+	in multithreaded environment you can provide an object of this class to
+	the Gamma() or Factorial() function, e.g;
+		typedef Big<1, 3> MyBig;
+		MyBig x = 123456;
+		CGamma<MyBig> cgamma;
+		std::cout << Gamma(x, cgamma);
+	each thread should have its own CGamma<> object
+
+	in a single-thread environment a CGamma<> object is a static variable
+	in a second version of Gamma() and you don't have to explicitly use it, e.g.
+		typedef Big<1, 3> MyBig;
+		MyBig x = 123456;
+		std::cout << Gamma(x);
+*/
+template<class ValueType>
+struct CGamma
+{
+	/*!
+		this table holds factorials
+			1
+			1
+			2
+			6
+			24
+			120
+			720
+			.......
+	*/
+	std::vector<ValueType> fact;
+
+
+	/*!
+		this table holds Bernoulli numbers
+			1
+			-0.5
+			0.166666666666666666666666667
+			0
+			-0.0333333333333333333333333333
+			0
+			0.0238095238095238095238095238
+			0
+			-0.0333333333333333333333333333
+			0
+			0.075757575757575757575757576
+			.....
+	*/
+	std::vector<ValueType> bern;
+
+
+	/*!
+		here we store some calculated values
+		(this is for speeding up, if the next argument of Gamma() or Factorial()
+		is in the 'history' then the result we are not calculating but simply
+		return from the 'history' object)
+	*/
+	History<ValueType> history;
+
+
+	/*!
+		this method prepares some coefficients: factorials and Bernoulli numbers
+		stored in 'fact' and 'bern' objects
+		
+		how many values should be depends on the size of the mantissa - if
+		the mantissa is larger then we must calculate more values
+		    for a mantissa which consists of 256 bits (8 words on a 32bit platform)
+			we have to calculate about 30 values (the size of fact and bern will be 30),
+			and for a 2048 bits mantissa we have to calculate 306 coefficients
+
+		you don't have to call this method, these coefficients will be automatically calculated
+		when they are needed
+
+		you must note that calculating of the coefficients is a little time-consuming operation,
+		(especially when the mantissa is large) and first called to Gamma() or Factorial()
+		can take more time than next calls, and in the end this is the point when InitAll()
+		comes in handy: you can call this method somewhere at the beginning of your program
+	*/
+	void InitAll();
+	// definition is in ttmath.h
+};
+
+
 
 
 } // namespace

ttmath/ttmathparser.h

@@ -1,7 +1,7 @@
 /*
  * This file is a part of TTMath Bignum Library
  * and is distributed under the (new) BSD licence.
- * Author: Tomasz Sowa <t.sowa@slimaczek.pl>
+ * Author: Tomasz Sowa <t.sowa@ttmath.org>
  */
 
 /* 
@@ -45,7 +45,7 @@
     \brief A mathematical parser
 */
 
-#include <fstream>
+#include <cstdio>
 #include <vector>
 #include <map>
 #include <set>
@@ -137,7 +137,6 @@ namespace ttmath
 template<class ValueType>
 class Parser
 {
-
 private:
 
 /*!
@@ -256,7 +255,7 @@ public:
 		bool function;
 
 		// if function is true
-		std::string function_name;
+		tt_string function_name;
 
 		/*
 			the sign of value
@@ -311,10 +310,11 @@ ErrorCode error;
 
 /*!
 	pointer to the currently reading char
+	it's either char* or wchar_t*
 
 	when an error has occured it may be used to count the index of the wrong character
 */
-const char * pstring;
+const tt_char * pstring;
 
 
 /*!
@@ -351,7 +351,7 @@ const Objects * puser_variables;
 const Objects * puser_functions;
 
 
-typedef std::map<std::string, ValueType> FunctionLocalVariables;
+typedef std::map<tt_string, ValueType> FunctionLocalVariables;
 
 /*!
 	a pointer to the local variables of a function
@@ -362,13 +362,13 @@ const FunctionLocalVariables * pfunction_local_variables;
 /*!
 	a temporary set using during parsing user defined variables
 */
-std::set<std::string> visited_variables;
+std::set<tt_string> visited_variables;
 
 
 /*!
 	a temporary set using during parsing user defined functions
 */
-std::set<std::string> visited_functions;
+std::set<tt_string> visited_functions;
 
 
 
@@ -396,10 +396,10 @@ typedef void (ValueType::*pfunction_var)();
 	table of mathematic functions
 
 	this map consists of:
-		std::string - function's name
+		tt_string - function's name
 		pfunction - pointer to specific function
 */
-typedef std::map<std::string, pfunction> FunctionsTable;
+typedef std::map<tt_string, pfunction> FunctionsTable;
 FunctionsTable functions_table;
 
 
@@ -407,10 +407,10 @@ FunctionsTable functions_table;
 	table of mathematic operators
 
 	this map consists of:
-		std::string - operators's name
+		tt_string - operators's name
 		MatOperator::Type - type of the operator
 */
-typedef std::map<std::string, typename MatOperator::Type> OperatorsTable;
+typedef std::map<tt_string, typename MatOperator::Type> OperatorsTable;
 OperatorsTable operators_table;
 
 
@@ -418,18 +418,17 @@ OperatorsTable operators_table;
 	table of mathematic variables
 
 	this map consists of:
-		std::string   - variable's name
+		tt_string     - variable's name
 		pfunction_var - pointer to specific function which returns value of variable
 */
-typedef std::map<std::string, pfunction_var> VariablesTable;
+typedef std::map<tt_string, pfunction_var> VariablesTable;
 VariablesTable variables_table;
 
 
 /*!
-	you can't calculate the factorial if the argument is greater than 'factorial_max'
+	some coefficients used when calculating the gamma (or factorial) function
-	default value is zero which means there are not any limitations
 */
-ValueType factorial_max;
+CGamma<ValueType> cgamma;
 
 
 /*!
@@ -456,7 +455,7 @@ void SkipWhiteCharacters()
 /*!
 	an auxiliary method for RecurrenceParsingVariablesOrFunction(...)
 */
-void RecurrenceParsingVariablesOrFunction_CheckStopCondition(bool variable, const std::string & name)
+void RecurrenceParsingVariablesOrFunction_CheckStopCondition(bool variable, const tt_string & name)
 {
 	if( variable )
 	{
@@ -474,7 +473,7 @@ void RecurrenceParsingVariablesOrFunction_CheckStopCondition(bool variable, cons
 /*!
 	an auxiliary method for RecurrenceParsingVariablesOrFunction(...)
 */
-void RecurrenceParsingVariablesOrFunction_AddName(bool variable, const std::string & name)
+void RecurrenceParsingVariablesOrFunction_AddName(bool variable, const tt_string & name)
 {
 	if( variable )
 		visited_variables.insert( name );
@@ -486,7 +485,7 @@ void RecurrenceParsingVariablesOrFunction_AddName(bool variable, const std::stri
 /*!
 	an auxiliary method for RecurrenceParsingVariablesOrFunction(...)
 */
-void RecurrenceParsingVariablesOrFunction_DeleteName(bool variable, const std::string & name)
+void RecurrenceParsingVariablesOrFunction_DeleteName(bool variable, const tt_string & name)
 {
 	if( variable )
 		visited_variables.erase( name );
@@ -505,7 +504,8 @@ void RecurrenceParsingVariablesOrFunction_DeleteName(bool variable, const std::s
 	(there can be a recurrence here therefore we're using 'visited_variables'
 	and 'visited_functions' sets to make a stop condition)
 */
-ValueType RecurrenceParsingVariablesOrFunction(bool variable, const std::string & name, const char * new_string, FunctionLocalVariables * local_variables = 0)
+ValueType RecurrenceParsingVariablesOrFunction(bool variable, const tt_string & name, const tt_char * new_string,
+											   FunctionLocalVariables * local_variables = 0)
 {
 	RecurrenceParsingVariablesOrFunction_CheckStopCondition(variable, name);
 	RecurrenceParsingVariablesOrFunction_AddName(variable, name);
@@ -548,12 +548,12 @@ public:
 /*!
 	this method returns the user-defined value of a variable
 */
-bool GetValueOfUserDefinedVariable(const std::string & variable_name,ValueType & result)
+bool GetValueOfUserDefinedVariable(const tt_string & variable_name,ValueType & result)
 {
 	if( !puser_variables )
 		return false;
 
-	const char * string_value;
+	const tt_char * string_value;
 
 	if( puser_variables->GetValue(variable_name, &string_value) != err_ok )
 		return false;
@@ -567,7 +567,7 @@ return true;
 /*!
 	this method returns the value of a local variable of a function
 */
-bool GetValueOfFunctionLocalVariable(const std::string & variable_name, ValueType & result)
+bool GetValueOfFunctionLocalVariable(const tt_string & variable_name, ValueType & result)
 {
 	if( !pfunction_local_variables )
 		return false;
@@ -589,7 +589,7 @@ return true;
 	we make an object of type ValueType then call a method which 
 	sets the correct value in it and finally we'll return the object
 */
-ValueType GetValueOfVariable(const std::string & variable_name)
+ValueType GetValueOfVariable(const tt_string & variable_name)
 {
 ValueType result;
 
@@ -600,7 +600,7 @@ ValueType result;
 		return result;
 
 
-	typename std::map<std::string, pfunction_var>::iterator i =
+	typename std::map<tt_string, pfunction_var>::iterator i =
 													variables_table.find(variable_name);
 
 	if( i == variables_table.end() )
@@ -674,6 +674,20 @@ return result;
 }
 
 
+void Gamma(int sindex, int amount_of_args, ValueType & result)
+{
+	if( amount_of_args != 1 )
+		Error( err_improper_amount_of_arguments );
+
+	ErrorCode err;
+	
+	result = ttmath::Gamma(stack[sindex].value, cgamma, &err, pstop_calculating);
+
+	if(err != err_ok)
+		Error( err );
+}
+
+
 /*!
 	factorial
 	result = 1 * 2 * 3 * 4 * .... * x
@@ -685,10 +699,7 @@ void Factorial(int sindex, int amount_of_args, ValueType & result)
 
 	ErrorCode err;
 
-	if( !factorial_max.IsZero() && stack[sindex].value > factorial_max )
+	result = ttmath::Factorial(stack[sindex].value, cgamma, &err, pstop_calculating);
-		Error( err_too_big_factorial );
-
-	result = ttmath::Factorial(stack[sindex].value, &err, pstop_calculating);
 
 	if(err != err_ok)
 		Error( err );
@@ -1344,17 +1355,34 @@ void Avg(int sindex, int amount_of_args, ValueType & result)
 }	
 
 
+
+/*!
+	we use such a method because 'wvsprintf' is not everywhere defined
+*/
+void Sprintf(tt_char * buffer, int par)
+{
+char buf[30]; // char, not tt_char
+int i;
+
+	sprintf(buf, "%d", par);
+	for(i=0 ; buf[i] != 0 ; ++i)
+		buffer[i] = buf[i];
+
+	buffer[i] = 0;
+}
+
+
 /*!
 	this method returns the value from a user-defined function
 
 	(look at the description in 'CallFunction(...)')
 */
-bool GetValueOfUserDefinedFunction(const std::string & function_name, int amount_of_args, int sindex)
+bool GetValueOfUserDefinedFunction(const tt_string & function_name, int amount_of_args, int sindex)
 {
 	if( !puser_functions )
 		return false;
 
-	const char * string_value;
+	const tt_char * string_value;
 	int param;
 
 	if( puser_functions->GetValueAndParam(function_name, &string_value, &param) != err_ok )
@@ -1368,15 +1396,17 @@ bool GetValueOfUserDefinedFunction(const std::string & function_name, int amount
 
 	if( amount_of_args > 0 )
 	{
-		char buffer[20];
+		tt_char buffer[30];
 
 		// x = x1
-		sprintf(buffer,"x");
+		buffer[0] = 'x';
+		buffer[1] = 0;
 		local_variables.insert( std::make_pair(buffer, stack[sindex].value) );
 
 		for(int i=0 ; i<amount_of_args ; ++i)
 		{
-			sprintf(buffer,"x%d",i+1);
+			buffer[0] = 'x';
+			Sprintf(buffer+1, i+1);
 			local_variables.insert( std::make_pair(buffer, stack[sindex + i*2].value) );
 		}
 	}
@@ -1400,7 +1430,7 @@ return true;
 	result will be stored in 'stack[sindex-1].value'
 	(we don't have to set the correct type of this element, it'll be set later)
 */
-void CallFunction(const std::string & function_name, int amount_of_args, int sindex)
+void CallFunction(const tt_string & function_name, int amount_of_args, int sindex)
 {
 	if( GetValueOfUserDefinedFunction(function_name, amount_of_args, sindex) )
 		return;
@@ -1426,9 +1456,9 @@ void CallFunction(const std::string & function_name, int amount_of_args, int sin
 	function_name - name of the function
 	pf - pointer to the function (to the wrapper)
 */
-void InsertFunctionToTable(const char * function_name, pfunction pf)
+void InsertFunctionToTable(const tt_char * function_name, pfunction pf)
 {
-	functions_table.insert( std::make_pair(std::string(function_name), pf));
+	functions_table.insert( std::make_pair(tt_string(function_name), pf));
 }
 
 
@@ -1439,9 +1469,9 @@ void InsertFunctionToTable(const char * function_name, pfunction pf)
 	variable_name - name of the function
 	pf - pointer to the function
 */
-void InsertVariableToTable(const char * variable_name, pfunction_var pf)
+void InsertVariableToTable(const tt_char * variable_name, pfunction_var pf)
 {
-	variables_table.insert( std::make_pair(std::string(variable_name), pf));
+	variables_table.insert( std::make_pair(tt_string(variable_name), pf));
 }
 
 
@@ -1450,67 +1480,65 @@ void InsertVariableToTable(const char * variable_name, pfunction_var pf)
 */
 void CreateFunctionsTable()
 {
-	/*
+	InsertFunctionToTable(TTMATH_TEXT("gamma"),		&Parser<ValueType>::Gamma);
-		names of functions should consist of small letters
+	InsertFunctionToTable(TTMATH_TEXT("factorial"),	&Parser<ValueType>::Factorial);
-	*/
+	InsertFunctionToTable(TTMATH_TEXT("abs"),   	&Parser<ValueType>::Abs);
-	InsertFunctionToTable("factorial",	&Parser<ValueType>::Factorial);
+	InsertFunctionToTable(TTMATH_TEXT("sin"),   	&Parser<ValueType>::Sin);
-	InsertFunctionToTable("abs",   		&Parser<ValueType>::Abs);
+	InsertFunctionToTable(TTMATH_TEXT("cos"),   	&Parser<ValueType>::Cos);
-	InsertFunctionToTable("sin",   		&Parser<ValueType>::Sin);
+	InsertFunctionToTable(TTMATH_TEXT("tan"),   	&Parser<ValueType>::Tan);
-	InsertFunctionToTable("cos",   		&Parser<ValueType>::Cos);
+	InsertFunctionToTable(TTMATH_TEXT("tg"),		&Parser<ValueType>::Tan);
-	InsertFunctionToTable("tan",   		&Parser<ValueType>::Tan);
+	InsertFunctionToTable(TTMATH_TEXT("cot"),  		&Parser<ValueType>::Cot);
-	InsertFunctionToTable("tg",		  	&Parser<ValueType>::Tan);
+	InsertFunctionToTable(TTMATH_TEXT("ctg"),  		&Parser<ValueType>::Cot);
-	InsertFunctionToTable("cot",  		&Parser<ValueType>::Cot);
+	InsertFunctionToTable(TTMATH_TEXT("int"),	   	&Parser<ValueType>::Int);
-	InsertFunctionToTable("ctg",  		&Parser<ValueType>::Cot);
+	InsertFunctionToTable(TTMATH_TEXT("round"),	 	&Parser<ValueType>::Round);
-	InsertFunctionToTable("int",	   	&Parser<ValueType>::Int);
+	InsertFunctionToTable(TTMATH_TEXT("ln"),	    &Parser<ValueType>::Ln);
-	InsertFunctionToTable("round",	 	&Parser<ValueType>::Round);
+	InsertFunctionToTable(TTMATH_TEXT("log"),	   	&Parser<ValueType>::Log);
-	InsertFunctionToTable("ln",	    	&Parser<ValueType>::Ln);
+	InsertFunctionToTable(TTMATH_TEXT("exp"),	   	&Parser<ValueType>::Exp);
-	InsertFunctionToTable("log",	   	&Parser<ValueType>::Log);
+	InsertFunctionToTable(TTMATH_TEXT("max"),	   	&Parser<ValueType>::Max);
-	InsertFunctionToTable("exp",	   	&Parser<ValueType>::Exp);
+	InsertFunctionToTable(TTMATH_TEXT("min"),	   	&Parser<ValueType>::Min);
-	InsertFunctionToTable("max",	   	&Parser<ValueType>::Max);
+	InsertFunctionToTable(TTMATH_TEXT("asin"),   	&Parser<ValueType>::ASin);
-	InsertFunctionToTable("min",	   	&Parser<ValueType>::Min);
+	InsertFunctionToTable(TTMATH_TEXT("acos"),   	&Parser<ValueType>::ACos);
-	InsertFunctionToTable("asin",   	&Parser<ValueType>::ASin);
+	InsertFunctionToTable(TTMATH_TEXT("atan"),   	&Parser<ValueType>::ATan);
-	InsertFunctionToTable("acos",   	&Parser<ValueType>::ACos);
+	InsertFunctionToTable(TTMATH_TEXT("atg"),	   	&Parser<ValueType>::ATan);
-	InsertFunctionToTable("atan",   	&Parser<ValueType>::ATan);
+	InsertFunctionToTable(TTMATH_TEXT("acot"),   	&Parser<ValueType>::ACot);
-	InsertFunctionToTable("atg",	   	&Parser<ValueType>::ATan);
+	InsertFunctionToTable(TTMATH_TEXT("actg"),   	&Parser<ValueType>::ACot);
-	InsertFunctionToTable("acot",   	&Parser<ValueType>::ACot);
+	InsertFunctionToTable(TTMATH_TEXT("sgn"),   	&Parser<ValueType>::Sgn);
-	InsertFunctionToTable("actg",   	&Parser<ValueType>::ACot);
+	InsertFunctionToTable(TTMATH_TEXT("mod"),   	&Parser<ValueType>::Mod);
-	InsertFunctionToTable("sgn",   		&Parser<ValueType>::Sgn);
+	InsertFunctionToTable(TTMATH_TEXT("if"),   		&Parser<ValueType>::If);
-	InsertFunctionToTable("mod",   		&Parser<ValueType>::Mod);
+	InsertFunctionToTable(TTMATH_TEXT("or"),   		&Parser<ValueType>::Or);
-	InsertFunctionToTable("if",   		&Parser<ValueType>::If);
+	InsertFunctionToTable(TTMATH_TEXT("and"),  		&Parser<ValueType>::And);
-	InsertFunctionToTable("or",   		&Parser<ValueType>::Or);
+	InsertFunctionToTable(TTMATH_TEXT("not"),  		&Parser<ValueType>::Not);
-	InsertFunctionToTable("and",  		&Parser<ValueType>::And);
+	InsertFunctionToTable(TTMATH_TEXT("degtorad"),	&Parser<ValueType>::DegToRad);
-	InsertFunctionToTable("not",  		&Parser<ValueType>::Not);
+	InsertFunctionToTable(TTMATH_TEXT("radtodeg"),	&Parser<ValueType>::RadToDeg);
-	InsertFunctionToTable("degtorad",	&Parser<ValueType>::DegToRad);
+	InsertFunctionToTable(TTMATH_TEXT("degtodeg"),	&Parser<ValueType>::DegToDeg);
-	InsertFunctionToTable("radtodeg",	&Parser<ValueType>::RadToDeg);
+	InsertFunctionToTable(TTMATH_TEXT("gradtorad"),	&Parser<ValueType>::GradToRad);
-	InsertFunctionToTable("degtodeg",	&Parser<ValueType>::DegToDeg);
+	InsertFunctionToTable(TTMATH_TEXT("radtograd"),	&Parser<ValueType>::RadToGrad);
-	InsertFunctionToTable("gradtorad",	&Parser<ValueType>::GradToRad);
+	InsertFunctionToTable(TTMATH_TEXT("degtograd"),	&Parser<ValueType>::DegToGrad);
-	InsertFunctionToTable("radtograd",	&Parser<ValueType>::RadToGrad);
+	InsertFunctionToTable(TTMATH_TEXT("gradtodeg"),	&Parser<ValueType>::GradToDeg);
-	InsertFunctionToTable("degtograd",	&Parser<ValueType>::DegToGrad);
+	InsertFunctionToTable(TTMATH_TEXT("ceil"),		&Parser<ValueType>::Ceil);
-	InsertFunctionToTable("gradtodeg",	&Parser<ValueType>::GradToDeg);
+	InsertFunctionToTable(TTMATH_TEXT("floor"),		&Parser<ValueType>::Floor);
-	InsertFunctionToTable("ceil",		&Parser<ValueType>::Ceil);
+	InsertFunctionToTable(TTMATH_TEXT("sqrt"),		&Parser<ValueType>::Sqrt);
-	InsertFunctionToTable("floor",		&Parser<ValueType>::Floor);
+	InsertFunctionToTable(TTMATH_TEXT("sinh"),		&Parser<ValueType>::Sinh);
-	InsertFunctionToTable("sqrt",		&Parser<ValueType>::Sqrt);
+	InsertFunctionToTable(TTMATH_TEXT("cosh"),		&Parser<ValueType>::Cosh);
-	InsertFunctionToTable("sinh",		&Parser<ValueType>::Sinh);
+	InsertFunctionToTable(TTMATH_TEXT("tanh"),		&Parser<ValueType>::Tanh);
-	InsertFunctionToTable("cosh",		&Parser<ValueType>::Cosh);
+	InsertFunctionToTable(TTMATH_TEXT("tgh"),		&Parser<ValueType>::Tanh);
-	InsertFunctionToTable("tanh",		&Parser<ValueType>::Tanh);
+	InsertFunctionToTable(TTMATH_TEXT("coth"),		&Parser<ValueType>::Coth);
-	InsertFunctionToTable("tgh",		&Parser<ValueType>::Tanh);
+	InsertFunctionToTable(TTMATH_TEXT("ctgh"),		&Parser<ValueType>::Coth);
-	InsertFunctionToTable("coth",		&Parser<ValueType>::Coth);
+	InsertFunctionToTable(TTMATH_TEXT("root"),		&Parser<ValueType>::Root);
-	InsertFunctionToTable("ctgh",		&Parser<ValueType>::Coth);
+	InsertFunctionToTable(TTMATH_TEXT("asinh"),		&Parser<ValueType>::ASinh);
-	InsertFunctionToTable("root",		&Parser<ValueType>::Root);
+	InsertFunctionToTable(TTMATH_TEXT("acosh"),		&Parser<ValueType>::ACosh);
-	InsertFunctionToTable("asinh",		&Parser<ValueType>::ASinh);
+	InsertFunctionToTable(TTMATH_TEXT("atanh"),		&Parser<ValueType>::ATanh);
-	InsertFunctionToTable("acosh",		&Parser<ValueType>::ACosh);
+	InsertFunctionToTable(TTMATH_TEXT("atgh"),		&Parser<ValueType>::ATanh);
-	InsertFunctionToTable("atanh",		&Parser<ValueType>::ATanh);
+	InsertFunctionToTable(TTMATH_TEXT("acoth"),		&Parser<ValueType>::ACoth);
-	InsertFunctionToTable("atgh",		&Parser<ValueType>::ATanh);
+	InsertFunctionToTable(TTMATH_TEXT("actgh"),		&Parser<ValueType>::ACoth);
-	InsertFunctionToTable("acoth",		&Parser<ValueType>::ACoth);
+	InsertFunctionToTable(TTMATH_TEXT("bitand"),	&Parser<ValueType>::BitAnd);
-	InsertFunctionToTable("actgh",		&Parser<ValueType>::ACoth);
+	InsertFunctionToTable(TTMATH_TEXT("bitor"),		&Parser<ValueType>::BitOr);
-	InsertFunctionToTable("bitand",		&Parser<ValueType>::BitAnd);
+	InsertFunctionToTable(TTMATH_TEXT("bitxor"),	&Parser<ValueType>::BitXor);
-	InsertFunctionToTable("bitor",		&Parser<ValueType>::BitOr);
+	InsertFunctionToTable(TTMATH_TEXT("band"),		&Parser<ValueType>::BitAnd);
-	InsertFunctionToTable("bitxor",		&Parser<ValueType>::BitXor);
+	InsertFunctionToTable(TTMATH_TEXT("bor"),		&Parser<ValueType>::BitOr);
-	InsertFunctionToTable("band",		&Parser<ValueType>::BitAnd);
+	InsertFunctionToTable(TTMATH_TEXT("bxor"),		&Parser<ValueType>::BitXor);
-	InsertFunctionToTable("bor",		&Parser<ValueType>::BitOr);
+	InsertFunctionToTable(TTMATH_TEXT("sum"),		&Parser<ValueType>::Sum);
-	InsertFunctionToTable("bxor",		&Parser<ValueType>::BitXor);
+	InsertFunctionToTable(TTMATH_TEXT("avg"),		&Parser<ValueType>::Avg);
-	InsertFunctionToTable("sum",		&Parser<ValueType>::Sum);
-	InsertFunctionToTable("avg",		&Parser<ValueType>::Avg);
 }
 
 
@@ -1519,11 +1547,8 @@ void CreateFunctionsTable()
 */
 void CreateVariablesTable()
 {
-	/*
+	InsertVariableToTable(TTMATH_TEXT("pi"), &ValueType::SetPi);
-		names of variables should consist of small letters
+	InsertVariableToTable(TTMATH_TEXT("e"),  &ValueType::SetE);
-	*/
-	InsertVariableToTable("pi", &ValueType::SetPi);
-	InsertVariableToTable("e",  &ValueType::SetE);
 }
 
 
@@ -1549,7 +1574,7 @@ return c;
 	what should be returned is tested just by a '(' character that means if there's
 	a '(' character after a name that function returns 'true'
 */
-bool ReadName(std::string & result)
+bool ReadName(tt_string & result)
 {
 int character;
 
@@ -1566,7 +1591,7 @@ int character;
 
 	do
 	{
-		result   += static_cast<char>( character );
+		result   += static_cast<tt_char>( character );
 		character = * ++pstring;
 	}
 	while(	(character>='a' && character<='z') ||
@@ -1621,7 +1646,7 @@ return false;
 */
 bool ReadVariableOrFunction(Item & result)
 {
-std::string name;
+tt_string name;
 bool is_it_name_of_function = ReadName(name);
 
 	if( is_it_name_of_function )
@@ -1650,7 +1675,7 @@ return is_it_name_of_function;
 */
 void ReadValue(Item & result, int reading_base)
 {
-const char * new_stack_pointer;
+const tt_char * new_stack_pointer;
 bool value_read;
 
 	int carry = result.value.FromString(pstring, reading_base, &new_stack_pointer, &value_read);
@@ -1823,9 +1848,9 @@ return 0;
 }
 
 
-void InsertOperatorToTable(const std::string & name, typename MatOperator::Type type)
+void InsertOperatorToTable(const tt_char * name, typename MatOperator::Type type)
 {
-	operators_table.insert( std::make_pair(name, type) );
+	operators_table.insert( std::make_pair(tt_string(name), type) );
 }
 
 
@@ -1834,19 +1859,19 @@ void InsertOperatorToTable(const std::string & name, typename MatOperator::Type
 */
 void CreateMathematicalOperatorsTable()
 {
-	InsertOperatorToTable(std::string("||"), MatOperator::lor);
+	InsertOperatorToTable(TTMATH_TEXT("||"), MatOperator::lor);
-	InsertOperatorToTable(std::string("&&"), MatOperator::land);
+	InsertOperatorToTable(TTMATH_TEXT("&&"), MatOperator::land);
-	InsertOperatorToTable(std::string("!="), MatOperator::neq);
+	InsertOperatorToTable(TTMATH_TEXT("!="), MatOperator::neq);
-	InsertOperatorToTable(std::string("=="), MatOperator::eq);
+	InsertOperatorToTable(TTMATH_TEXT("=="), MatOperator::eq);
-	InsertOperatorToTable(std::string(">="), MatOperator::get);
+	InsertOperatorToTable(TTMATH_TEXT(">="), MatOperator::get);
-	InsertOperatorToTable(std::string("<="), MatOperator::let);
+	InsertOperatorToTable(TTMATH_TEXT("<="), MatOperator::let);
-	InsertOperatorToTable(std::string(">"),  MatOperator::gt);
+	InsertOperatorToTable(TTMATH_TEXT(">"),  MatOperator::gt);
-	InsertOperatorToTable(std::string("<"),  MatOperator::lt);
+	InsertOperatorToTable(TTMATH_TEXT("<"),  MatOperator::lt);
-	InsertOperatorToTable(std::string("-"),  MatOperator::sub);
+	InsertOperatorToTable(TTMATH_TEXT("-"),  MatOperator::sub);
-	InsertOperatorToTable(std::string("+"),  MatOperator::add);
+	InsertOperatorToTable(TTMATH_TEXT("+"),  MatOperator::add);
-	InsertOperatorToTable(std::string("/"),  MatOperator::div);
+	InsertOperatorToTable(TTMATH_TEXT("/"),  MatOperator::div);
-	InsertOperatorToTable(std::string("*"),  MatOperator::mul);
+	InsertOperatorToTable(TTMATH_TEXT("*"),  MatOperator::mul);
-	InsertOperatorToTable(std::string("^"),  MatOperator::pow);
+	InsertOperatorToTable(TTMATH_TEXT("^"),  MatOperator::pow);
 }
 
 
@@ -1856,12 +1881,12 @@ void CreateMathematicalOperatorsTable()
 	e.g.
 	true when str1="test" and str2="te"
 */
-bool IsSubstring(const std::string & str1, const std::string & str2)
+bool IsSubstring(const tt_string & str1, const tt_string & str2)
 {
 	if( str2.length() > str1.length() )
 		return false;
 
-	for(std::string::size_type i=0 ; i<str2.length() ; ++i)
+	for(tt_string::size_type i=0 ; i<str2.length() ; ++i)
 		if( str1[i] != str2[i] )
 			return false;
 
@@ -1874,7 +1899,7 @@ return true;
 */
 void ReadMathematicalOperator(Item & result)
 {
-std::string oper;
+tt_string oper;
 typename OperatorsTable::iterator iter_old, iter_new;
 
 	iter_old = operators_table.end();
@@ -2404,7 +2429,6 @@ Parser(): default_stack_size(100)
 	base = 10;
 	deg_rad_grad = 1;
 	error = err_ok;
-	factorial_max.SetZero();
 
 	CreateFunctionsTable();
 	CreateVariablesTable();
@@ -2424,7 +2448,6 @@ Parser<ValueType> & operator=(const Parser<ValueType> & p)
 	base              = p.base;
 	deg_rad_grad      = p.deg_rad_grad;
 	error             = err_ok;
-	factorial_max = p.factorial_max;
 
 	/*
 		we don't have to call 'CreateFunctionsTable()' etc.
@@ -2506,22 +2529,11 @@ void SetFunctions(const Objects * pf)
 }
 
 
-/*!
-	you will not be allowed to calculate the factorial 
-	if its argument is greater than 'm'
-	there'll be: ErrorCode::err_too_big_factorial
-	default 'factorial_max' is zero which means you can calculate what you want to
-*/
-void SetFactorialMax(const ValueType & m)
-{
-	factorial_max = m;
-}
-
 
 /*!
 	the main method using for parsing string
 */
-ErrorCode Parse(const char * str)
+ErrorCode Parse(const tt_char * str)
 {
 	stack_index  = 0;
 	pstring      = str;
@@ -2544,11 +2556,11 @@ return error;
 }
 
 
-
-
-
 };
 
+
+
+
 } // namespace
 
 

ttmath/ttmathtypes.h

@@ -1,7 +1,7 @@
 /*
  * This file is a part of TTMath Bignum Library
  * and is distributed under the (new) BSD licence.
- * Author: Tomasz Sowa <t.sowa@slimaczek.pl>
+ * Author: Tomasz Sowa <t.sowa@ttmath.org>
  */
 
 /* 
@@ -54,7 +54,7 @@
 
 #include <stdexcept>
 #include <sstream>
-
+#include <vector>
 
 /*!
 	the version of the library
@@ -63,9 +63,9 @@
 	if zero that means this is the release version of the library
 */
 #define TTMATH_MAJOR_VER		0
-#define TTMATH_MINOR_VER		8
+#define TTMATH_MINOR_VER		9
-#define TTMATH_REVISION_VER		5
+#define TTMATH_REVISION_VER		0
-#define TTMATH_PRERELEASE_VER	0
+#define TTMATH_PRERELEASE_VER	1
 
 
 /*!
@@ -120,6 +120,7 @@ namespace ttmath
 	typedef unsigned int uint;
 	typedef signed   int sint;
 
+
 	/*!
 		this type is twice bigger than uint
 		(64bit on a 32bit platforms)
@@ -128,43 +129,39 @@ namespace ttmath
 		but it is defined in C99 and in upcoming C++0x /3.9.1 (2)/ and many compilers support it
 
 		this type is used in UInt::MulTwoWords and UInt::DivTwoWords when macro TTMATH_NOASM is defined
-		but only on a 32bit platform
 	*/
-	#ifdef TTMATH_NOASM
 	typedef unsigned long long int ulint;
-	#endif
-
-	/*!
-		how many bits there are in the uint type
-	*/
-	#define TTMATH_BITS_PER_UINT 32u
 
 	/*!
 		the mask for the highest bit in the unsigned 32bit word (2^31)
 	*/
-	#define TTMATH_UINT_HIGHEST_BIT 2147483648u
+	const uint TTMATH_UINT_HIGHEST_BIT = 0x80000000ul;
 
 	/*!
 		the max value of the unsigned 32bit word (2^32 - 1)
 		(all bits equal one)
 	*/
-	#define TTMATH_UINT_MAX_VALUE 4294967295u
+	const uint TTMATH_UINT_MAX_VALUE = 0xfffffffful;
 
 	/*!
 		the number of words (32bit words on 32bit platform)
 		which are kept in built-in variables for a Big<> type
 		(these variables are defined in ttmathbig.h)
 	*/
-	#define TTMATH_BUILTIN_VARIABLES_SIZE 256u
+	const uint TTMATH_BUILTIN_VARIABLES_SIZE  = 256u;
 
 #else
 
 	/*!
 		on 64bit platforms one word (uint, sint) will be equal 64bits
 	*/
-	typedef unsigned long uint;
+	#if defined(_MSC_VER)
-	typedef signed   long sint;
+		typedef unsigned __int64 uint;
-
+		typedef signed __int64 sint;
+	#else
+		typedef unsigned long long uint;
+		typedef signed long long sint;
+	#endif
 	/*!
 		on 64bit platform we do not define ulint
 		sizeof(long long) is 8 (64bit) but we need 128bit
@@ -174,33 +171,61 @@ namespace ttmath
 	*/
 	//typedef unsigned long long int ulint;
 
-	/*!
-		how many bits there are in the uint type
-	*/
-	#define TTMATH_BITS_PER_UINT 64ul
-
 	/*!
 		the mask for the highest bit in the unsigned 64bit word (2^63)
 	*/
-	#define TTMATH_UINT_HIGHEST_BIT 9223372036854775808ul
+	const uint TTMATH_UINT_HIGHEST_BIT = 0x8000000000000000ul;
 
 	/*!
 		the max value of the unsigned 64bit word (2^64 - 1)
 		(all bits equal one)
 	*/
-	#define TTMATH_UINT_MAX_VALUE 18446744073709551615ul
+	const uint TTMATH_UINT_MAX_VALUE = 0xfffffffffffffffful;
 
 	/*!
 		the number of words (64bit words on 64bit platforms)
 		which are kept in built-in variables for a Big<> type
 		(these variables are defined in ttmathbig.h)
 	*/
-	#define TTMATH_BUILTIN_VARIABLES_SIZE 128ul
+	const uint TTMATH_BUILTIN_VARIABLES_SIZE = 128ul;
 
 #endif
+
+	const uint TTMATH_BITS_PER_UINT = (sizeof(uint)*8);
+
 }
 
 
+#if defined(UNICODE) || defined(_UNICODE)
+#define TTMATH_USE_WCHAR
+#endif
+
+
+#ifdef TTMATH_USE_WCHAR
+
+	typedef	wchar_t					tt_char;
+	typedef	std::wstring			tt_string;
+	typedef std::wostringstream		tt_ostringstream;
+	typedef std::wostream			tt_ostream;
+	typedef std::wistream			tt_istream;
+	#define TTMATH_TEXT_HELPER(txt)	L##txt
+
+#else
+
+	typedef	char					tt_char;
+	typedef	std::string				tt_string;
+	typedef std::ostringstream		tt_ostringstream;
+	typedef std::ostream			tt_ostream;
+	typedef std::istream			tt_istream;
+	#define TTMATH_TEXT_HELPER(txt)	txt
+
+#endif
+
+#define TTMATH_TEXT(txt) 	TTMATH_TEXT_HELPER(txt)
+
+
+
+
 
 /*!
 	characters which represent the comma operator
@@ -249,6 +274,18 @@ namespace ttmath
 #endif
 
 
+/*!
+	this is a special value used when calculating the Gamma(x) function
+	if x is greater than this value then the Gamma(x) will be calculated using
+	some kind of series
+
+	don't use smaller values than about 100
+*/
+#define TTMATH_GAMMA_BOUNDARY 2000
+
+
+
+
 namespace ttmath
 {
 
@@ -282,7 +319,6 @@ namespace ttmath
 		err_object_exists,
 		err_unknown_object,
 		err_still_calculating,
-		err_too_big_factorial,
 		err_in_short_form_used_function
 	};
 
@@ -313,20 +349,20 @@ namespace ttmath
 	*/
 	class ExceptionInfo
 	{
-	const char * file;
+	const tt_char * file;
 	int line;
 
 	public:
 		ExceptionInfo() : file(0), line(0) {}
-		ExceptionInfo(const char * f, int l) : file(f), line(l) {}
+		ExceptionInfo(const tt_char * f, int l) : file(f), line(l) {}
 
-		std::string Where() const
+		tt_string Where() const
 		{
 			if( !file )
-				return "unknown";
+				return TTMATH_TEXT("unknown");
 
-			std::ostringstream result;
+			tt_ostringstream result;
-			result << file << ":" << line;
+			result << file << TTMATH_TEXT(":") << line;
 
 		return result.str();
 		}
@@ -340,7 +376,7 @@ namespace ttmath
 		can throw an exception of this type
 
 		If you compile with gcc you can get a small benefit 
-		from using method Where() (it returns std::string with
+		from using method Where() (it returns std::string (or std::wstring) with
 		the name and the line of a file where the macro TTMATH_REFERENCE_ASSERT
 		was used)
 
@@ -368,12 +404,12 @@ namespace ttmath
 		{
 		}
 
-		ReferenceError(const char * f, int l) :
+		ReferenceError(const tt_char * f, int l) :
 							std::logic_error ("reference error"), ExceptionInfo(f,l)
 		{
 		}
 
-		std::string Where() const
+		tt_string Where() const
 		{
 			return ExceptionInfo::Where();
 		}
@@ -388,7 +424,7 @@ namespace ttmath
 		of this type
 
 		if you compile with gcc you can get a small benefit 
-		from using method Where() (it returns std::string with
+		from using method Where() (it returns std::string (or std::wstring) with
 		the name and the line of a file where the macro TTMATH_ASSERT
 		was used)
 	*/
@@ -400,12 +436,12 @@ namespace ttmath
 		{
 		}
 
-		RuntimeError(const char * f, int l) :
+		RuntimeError(const tt_char * f, int l) :
 						std::runtime_error ("internal error"), ExceptionInfo(f,l)
 		{
 		}
 
-		std::string Where() const
+		tt_string Where() const
 		{
 			return ExceptionInfo::Where();
 		}
@@ -420,11 +456,22 @@ namespace ttmath
 
 		#if defined(__FILE__) && defined(__LINE__)
 
+			#ifdef TTMATH_USE_WCHAR
+				#define TTMATH_FILE_HELPER2(arg)  L##arg
+				#define TTMATH_FILE_HELPER(x)     TTMATH_FILE_HELPER2(x)
+				#define TTMATH_FILE               TTMATH_FILE_HELPER(__FILE__)
+			#else
+				#define TTMATH_FILE               __FILE__
+			#endif
+
 			#define TTMATH_REFERENCE_ASSERT(expression) \
-				if( &(expression) == this ) throw ttmath::ReferenceError(__FILE__, __LINE__);
+				if( &(expression) == this ) throw ttmath::ReferenceError(TTMATH_FILE, __LINE__);
 
 			#define TTMATH_ASSERT(expression) \
-				if( !(expression) ) throw ttmath::RuntimeError(__FILE__, __LINE__);
+				if( !(expression) ) throw ttmath::RuntimeError(TTMATH_FILE, __LINE__);
+
+			#define TTMATH_VERIFY(expression) \
+				if( !(expression) ) throw ttmath::RuntimeError(TTMATH_TEXT(__FILE__), __LINE__);
 
 		#else
 
@@ -433,19 +480,32 @@ namespace ttmath
 
 			#define TTMATH_ASSERT(expression) \
 				if( !(expression) ) throw RuntimeError();
+
+			#define TTMATH_VERIFY(expression) \
+				if( !(expression) ) throw RuntimeError();
 		#endif
 
 	#else
 		#define TTMATH_REFERENCE_ASSERT(expression)
 		#define TTMATH_ASSERT(expression)
+		#define TTMATH_VERIFY(expression)	(void)(expression);
 	#endif
 
-
+	#if !defined(LOG_PRINTF)
+		#define LOG_PRINTF printf
+	#endif
 
 	#ifdef TTMATH_DEBUG_LOG
 	
-		#define TTMATH_LOG(msg) \
+		#ifdef TTMATH_USE_WCHAR
+			#define TTMATH_LOG_HELPER(msg) \
+				PrintLog(L##msg, std::wcout);
+		#else
+			#define TTMATH_LOG_HELPER(msg) \
 				PrintLog(msg, std::cout);
+		#endif
+
+		#define TTMATH_LOG(msg) TTMATH_LOG_HELPER(msg)
 
 	#else
 
@@ -454,7 +514,9 @@ namespace ttmath
 	#endif
 
 
+
 } // namespace
 
 
 #endif
+

ttmath/ttmathuint.h

@@ -1,7 +1,7 @@
 /*
  * This file is a part of TTMath Bignum Library
  * and is distributed under the (new) BSD licence.
- * Author: Tomasz Sowa <t.sowa@slimaczek.pl>
+ * Author: Tomasz Sowa <t.sowa@ttmath.org>
  */
 
 /*
@@ -52,7 +52,9 @@
 
 #include "ttmathtypes.h"
 
-
+#if defined(_MSC_VER)
+	#pragma warning(disable:4127) // conditional expression is constant
+#endif
 
 /*!
     \brief a namespace for the TTMath library
@@ -86,7 +88,10 @@ public:
 
 		it prints the table in a nice form of several columns
 	*/
-	void PrintTable(std::ostream & output) const
+	#ifndef TTMATH_USE_WCHAR
+	// gcc has a problem with std::setfill when wchar_t is used
+
+	void PrintTable(tt_ostream & output) const
 	{
 		// how many columns there'll be
 		const int columns = 8;
@@ -94,7 +99,7 @@ public:
 		int c = 1;
 		for(int i=value_size-1 ; i>=0 ; --i)
 		{
-			output	<< "0x" << std::setfill('0');
+			output << TTMATH_TEXT("0x") << std::setfill('0');
 
 			#ifdef TTMATH_PLATFORM32
 				output << std::setw(8);
@@ -106,7 +111,7 @@ public:
 
 			if( i>0 )
 			{
-				output << ", ";		
+				output << TTMATH_TEXT(", ");		
 
 				if( ++c > columns )
 				{
@@ -118,14 +123,14 @@ public:
 
 		output << std::dec << std::endl;
 	}
+	#endif
 
-
+	void PrintLog(const tt_char * msg, tt_ostream & output) const
-	void PrintLog(const char * msg, std::ostream & output) const
 	{
 		output << msg << std::endl;
 
 		for(uint i=0 ; i<value_size ; ++i)
-			output << " table[" << i << "]: " << table[i] << std::endl;
+			output << TTMATH_TEXT(" table[") << i << TTMATH_TEXT("]: ") << table[i] << std::endl;
 	}
 
 
@@ -1966,7 +1971,7 @@ private:
 		for(uint i = j+1 ; i<value_size ; ++i)
 			q.table[i] = 0;
 
-		while( true )
+		for (;;)
 		{
 			u1 = table[j+n-1];
 			u0 = table[j+n-2];
@@ -2283,7 +2288,7 @@ public:
 	*/
 	bool IsTheLowestBitSet() const
 	{
-		return (*table & 1) != 0;
+		return (table[0] & 1) != 0;
 	}
 
 
@@ -2367,7 +2372,7 @@ public:
 
 
 	/*!
-		this method converts a digit into a char
+		this method converts a digit into a tt_char
 		digit should be from <0,F>
 		(we don't have to get a base)
 
@@ -2377,12 +2382,12 @@ public:
 			10 -> A
 			15 -> F
 	*/
-	static uint DigitToChar(uint digit)
+	static tt_char DigitToChar(uint digit)
 	{
 		if( digit < 10 )
-			return digit + '0';
+			return (tt_char)(digit + '0');
 
-	return digit - 10 + 'A';
+	return((tt_char)(digit - 10 + 'A'));
 	}
 
 
@@ -2503,7 +2508,7 @@ public:
 		this constant 10 has the int type (signed int), if we don't give such
 		operators and constructors the compiler will not compile the program,
 		because it has to make a conversion and doesn't know into which type
-		(the UInt class has operator=(const char*), operator=(uint) etc.)
+		(the UInt class has operator=(const tt_char*), operator=(uint) etc.)
 	*/
 	UInt<value_size> & operator=(sint i)
 	{
@@ -2618,18 +2623,18 @@ public:
 	/*!
 		a constructor for converting a string to this class (with the base=10)
 	*/
-	UInt(const char * s)
+	UInt(const tt_char * s)
 	{
 		FromString(s);
 
-		TTMATH_LOG("UInt::UInt(const char *)")
+		TTMATH_LOG("UInt::UInt(const tt_char *)")
 	}
 
 
 	/*!
 		a constructor for converting a string to this class (with the base=10)
 	*/
-	UInt(const std::string & s)
+	UInt(const tt_string & s)
 	{
 		FromString( s.c_str() );
 	}
@@ -2644,6 +2649,7 @@ public:
 	{
 	}
 
+
 	/*!
 		a copy constructor
 	*/
@@ -2694,10 +2700,10 @@ public:
 	/*!
 		this method converts the value to a string with a base equal 'b'
 	*/
-	void ToString(std::string & result, uint b = 10) const
+	void ToString(tt_string & result, uint b = 10) const
 	{
 	UInt<value_size> temp( *this );
-	char character;
+	tt_char character;
 	uint rem;
 
 		result.clear();
@@ -2708,7 +2714,7 @@ public:
 		do
 		{
 			temp.DivInt(b, &rem);
-			character = static_cast<char>( DigitToChar(rem) );
+			character = static_cast<tt_char>( DigitToChar(rem) );
 			result.insert(result.begin(), character);
 		}
 		while( !temp.IsZero() );
@@ -2722,7 +2728,7 @@ public:
 	/*
 		this method's ommiting any white characters from the string
 	*/
-	static void SkipWhiteCharacters(const char * & c)
+	static void SkipWhiteCharacters(const tt_char * & c)
 	{
 		while( (*c==' ' ) || (*c=='\t') || (*c==13 ) || (*c=='\n') )
 			++c;
@@ -2746,7 +2752,7 @@ public:
 
 		value_read (if exists) tells whether something has actually been read (at least one digit)
 	*/
-	uint FromString(const char * s, uint b = 10, const char ** after_source = 0, bool * value_read = 0)
+	uint FromString(const tt_char * s, uint b = 10, const tt_char ** after_source = 0, bool * value_read = 0)
 	{
 	UInt<value_size> base( b );
 	UInt<value_size> temp;
@@ -2796,7 +2802,7 @@ public:
 
 		(it returns carry=1 if the value will be too big or an incorrect base 'b' is given)
 	*/
-	uint FromString(const std::string & s, uint b = 10)
+	uint FromString(const tt_string & s, uint b = 10)
 	{
 		return FromString( s.c_str(), b );
 	}
@@ -2806,11 +2812,11 @@ public:
 	/*!
 		this operator converts a string into its value (with base = 10)
 	*/
-	UInt<value_size> & operator=(const char * s)
+	UInt<value_size> & operator=(const tt_char * s)
 	{
 		FromString(s);
 
-		TTMATH_LOG("UInt::operator=(const char *)")
+		TTMATH_LOG("UInt::operator=(const tt_char *)")
 
 	return *this;
 	}
@@ -2819,7 +2825,7 @@ public:
 	/*!
 		this operator converts a string into its value (with base = 10)
 	*/
-	UInt<value_size> & operator=(const std::string & s)
+	UInt<value_size> & operator=(const tt_string & s)
 	{
 		FromString( s.c_str() );
 
@@ -3187,9 +3193,15 @@ public:
 	*
 	*/
 
-	friend std::ostream & operator<<(std::ostream & s, const UInt<value_size> & l)
+
+	/*!
+		output for standard streams
+
+		tt_ostream is either std::ostream or std::wostream
+	*/
+	friend tt_ostream & operator<<(tt_ostream & s, const UInt<value_size> & l)
 	{
-	std::string ss;
+	tt_string ss;
 
 		l.ToString(ss);
 		s << ss;
@@ -3199,12 +3211,17 @@ public:
 
 
 
-	friend std::istream & operator>>(std::istream & s, UInt<value_size> & l)
+	/*!
-	{
+		input from standard streams
-	std::string ss;
 
-	// char for operator>>
+		tt_istream is either std::istream or std::wistream
-	unsigned char z;
+	*/
+	friend tt_istream & operator>>(tt_istream & s, UInt<value_size> & l)
+	{
+	tt_string ss;
+
+	// tt_char for operator>>
+	tt_char z;
 
 		// operator>> omits white characters if they're set for ommiting
 		s >> z;
@@ -3213,10 +3230,10 @@ public:
 		while( s.good() && CharToDigit(z, 10)>=0 )
 		{
 			ss += z;
-			z = s.get();
+			z = static_cast<tt_char>(s.get());
 		}
 
-		// we're leaving the last readed character
+		// we're leaving the last read character
 		// (it's not belonging to the value)
 		s.unget();
 
@@ -3236,7 +3253,6 @@ public:
 			ttmathuint_noasm.h
 	*/
 
-#ifdef TTMATH_NOASM
 	static uint AddTwoWords(uint a, uint b, uint carry, uint * result);
 	static uint SubTwoWords(uint a, uint b, uint carry, uint * result);
 
@@ -3261,8 +3277,6 @@ public:
 	static void MultiplySubtract(uint_ & u_, unsigned int & u3, unsigned int & q, uint_ v_);
 
 #endif // TTMATH_PLATFORM64
-#endif // TTMATH_NOASM
-
 
 private:
 	uint Rcl2_one(uint c);
@@ -3304,6 +3318,10 @@ public:
 
 } //namespace
 
+#if defined(_MSC_VER)
+	#pragma warning(default:4127) // conditional expression is constant
+#endif
+
 
 #include "ttmathuint_x86.h"
 #include "ttmathuint_x86_64.h"

ttmath/ttmathuint_noasm.h

@@ -1,7 +1,7 @@
 /*
  * This file is a part of TTMath Bignum Library
  * and is distributed under the (new) BSD licence.
- * Author: Tomasz Sowa <t.sowa@slimaczek.pl>
+ * Author: Tomasz Sowa <t.sowa@ttmath.org>
  */
 
 /* 
@@ -41,6 +41,8 @@
 
 #ifdef TTMATH_NOASM
 
+#pragma message("TTMATH_NOASM")
+
 /*!
 	\file ttmathuint_noasm.h
     \brief template class UInt<uint> with methods without any assembler code

ttmath/ttmathuint_x86.h

@@ -1,7 +1,7 @@
 /*
  * This file is a part of TTMath Bignum Library
  * and is distributed under the (new) BSD licence.
- * Author: Tomasz Sowa <t.sowa@slimaczek.pl>
+ * Author: Tomasz Sowa <t.sowa@ttmath.org>
  */
 
 /* 
@@ -36,14 +36,13 @@
  */
 
 
-
 #ifndef headerfilettmathuint_x86
 #define headerfilettmathuint_x86
 
-
 #ifndef TTMATH_NOASM
 #ifdef TTMATH_PLATFORM32
 
+#pragma message("TTMATH_ASM32")
 
 /*!
 	\file ttmathuint_x86.h
@@ -66,7 +65,6 @@ namespace ttmath
 	*
 	*/
 	
-
 	/*!
 		adding ss2 to the this and adding carry if it's defined
 		(this = this + ss2 + c)
@@ -85,46 +83,31 @@ namespace ttmath
 		// this algorithm doesn't require it
 
 		#ifndef __GNUC__
-			
 			//	this part might be compiled with for example visual c
 
 			__asm
 			{
-				push eax
+				xor eax,eax  // eax=0
-				push ebx
+				xor edx,edx  // edx=0
-				push ecx
-				push edx
-				push esi
-
 				mov ecx,[b]
-				
 				mov ebx,[p1]
 				mov esi,[p2]
 
-				xor edx,edx          // edx=0
+				sub eax,[c]  // CF=c
-				mov eax,[c]
-				neg eax              // CF=1 if rax!=0 , CF=0 if rax==0
 
+				ALIGN 16
 			ttmath_loop:				
-				mov eax,[esi+edx*4]
+				mov eax,[esi+edx*4+0]
-				adc [ebx+edx*4],eax
+				adc [ebx+edx*4+0],eax
 
-				inc edx
+				lea edx, [edx+1] // inc edx, but faster (no flags dependencies)
 				dec ecx
 			jnz ttmath_loop
 
-				adc ecx, ecx
+				setc	al
-				mov [c], ecx
+				movzx	eax, al
-
+				mov		[c], eax
-				pop esi
-				pop edx
-				pop ecx
-				pop ebx
-				pop eax
 			}
-
-
-
 		#endif		
 			
 
@@ -188,14 +171,8 @@ namespace ttmath
 		TTMATH_ASSERT( index < value_size )
 
 		#ifndef __GNUC__
-
 			__asm
 			{
-				push eax
-				push ebx
-				push ecx
-				push edx
-
 				mov ecx, [b]
 				sub ecx, [index]				
 
@@ -204,26 +181,21 @@ namespace ttmath
 
 				mov eax, [value]
 
+				ALIGN 16
 			ttmath_loop:
 				add [ebx+edx*4], eax
 			jnc ttmath_end
 
 				mov eax, 1
-				inc edx
+				lea edx, [edx+1] // inc edx, but faster (no flags dependencies)
 				dec ecx
 			jnz ttmath_loop
 
 			ttmath_end:
 				setc	al
-				movzx edx, al
+				movzx	eax, al
-				mov [c], edx
+				mov		[c], eax
-
-				pop edx
-				pop ecx
-				pop ebx
-				pop eax
 			}
-
 		#endif		
 			
 
@@ -303,52 +275,40 @@ namespace ttmath
 		#ifndef __GNUC__
 			__asm
 			{
-				push eax
-				push ebx
-				push ecx
-				push edx
-
 				mov ecx, [b]
-				sub ecx, [index]				
-
-				mov ebx, [p1]
 				mov edx, [index]
-
+				mov ebx, [p1]
 				mov eax, [x1]
+				sub ecx, edx	// max uints to add (value_size - index)
+
 				add [ebx+edx*4], eax
-				inc edx
+				lea ecx, [ecx-1]
-				dec ecx
 
 				mov eax, [x2]
 			
+				ALIGN 16
 			ttmath_loop:
-				adc [ebx+edx*4], eax
+				adc [ebx+edx*4+4], eax
 			jnc ttmath_end
 
 				mov eax, 0
-				inc edx
+				lea edx, [edx+1] // inc edx, but faster (no flags dependencies)
 				dec ecx
 			jnz ttmath_loop
 
 			ttmath_end:
 				setc	al
-				movzx edx, al
+				movzx	eax, al
-				mov [c], edx
+				mov		[c], eax
-				
-				pop edx
-				pop ecx
-				pop ebx
-				pop eax
-
 			}
 		#endif		
 
-
 		#ifdef __GNUC__
-		uint dummy, dummy2;
-
 			__asm__ __volatile__(
 			
+				"push %%ecx						\n"
+				"push %%edx						\n"
+
 				"subl %%edx, %%ecx 				\n"
 				
 				"addl %%esi, (%%ebx,%%edx,4) 	\n"
@@ -414,21 +374,19 @@ namespace ttmath
 			//	this part might be compiled with for example visual c
 			__asm
 			{
-				pushad
-
 				mov ecx, [ss2_size]
 				xor edx, edx               // edx = 0, cf = 0
 
 				mov esi, [ss1]
 				mov ebx, [ss2]
 				mov edi, [result]
-
+				ALIGN 16
 			ttmath_loop:
 				mov eax, [esi+edx*4]
 				adc eax, [ebx+edx*4]
 				mov [edi+edx*4], eax
 
-				inc edx
+				lea edx, [edx+1] // inc edx, but faster (no flags dependencies)
 				dec ecx
 			jnz ttmath_loop
 
@@ -447,7 +405,7 @@ namespace ttmath
 				adc eax, ebx 
 				mov [edi+edx*4], eax
 
-				inc edx
+				lea edx, [edx+1] // inc edx, but faster (no flags dependencies)
 				dec ecx
 			jnz ttmath_loop2
 
@@ -455,8 +413,6 @@ namespace ttmath
 
 			ttmath_end:
 				mov [c], ecx
-
-				popad
 			}
 
 		#endif		
@@ -506,8 +462,6 @@ namespace ttmath
 
 		#endif
 
-		TTMATH_LOG("UInt::AddVector")
-
 	return c;
 	}
 
@@ -531,40 +485,30 @@ namespace ttmath
 		// this algorithm doesn't require it
 
 		#ifndef __GNUC__
-
 			__asm
 			{
-				push eax
-				push ebx
-				push ecx
-				push edx
-				push esi
-
 				mov ecx,[b]
 				
 				mov ebx,[p1]
 				mov esi,[p2]
 
-				xor edx,edx          // edx=0
+				xor eax, eax
-				mov eax,[c]
+				mov edx, eax
-				neg eax              // CF=1 if rax!=0 , CF=0 if rax==0
 
+				sub eax, [c]
+
+				ALIGN 16
 			ttmath_loop:
 				mov eax, [esi+edx*4]
 				sbb [ebx+edx*4], eax
 
-				inc edx
+				lea edx, [edx+1] // inc edx, but faster (no flags dependencies)
 				dec ecx
 			jnz ttmath_loop
 
-				adc ecx, ecx
+				setc	al
-				mov [c], ecx
+				movzx	eax, al
-
+				mov		[c], eax
-				pop esi
-				pop edx
-				pop ecx
-				pop ebx
-				pop eax
 			}
 
 		#endif
@@ -631,14 +575,8 @@ namespace ttmath
 		TTMATH_ASSERT( index < value_size )
 
 		#ifndef __GNUC__
-
 			__asm
 			{
-				push eax
-				push ebx
-				push ecx
-				push edx
-
 				mov ecx, [b]
 				sub ecx, [index]				
 
@@ -647,26 +585,21 @@ namespace ttmath
 
 				mov eax, [value]
 
+				ALIGN 16
 			ttmath_loop:
 				sub [ebx+edx*4], eax
 			jnc ttmath_end
 
 				mov eax, 1
-				inc edx
+				lea edx, [edx+1] // inc edx, but faster (no flags dependencies)
 				dec ecx
 			jnz ttmath_loop
 
 			ttmath_end:
 				setc	al
-				movzx edx, al
+				movzx	eax, al
-				mov [c], edx
+				mov		[c], eax
-
-				pop edx
-				pop ecx
-				pop ebx
-				pop eax
 			}
-
 		#endif		
 			
 
@@ -742,7 +675,6 @@ namespace ttmath
 			*/
 			__asm
 			{
-				pushad
 
 				mov ecx, [ss2_size]
 				xor edx, edx               // edx = 0, cf = 0
@@ -756,7 +688,7 @@ namespace ttmath
 				sbb eax, [ebx+edx*4]
 				mov [edi+edx*4], eax
 
-				inc edx
+				lea edx, [edx+1]
 				dec ecx
 			jnz ttmath_loop
 
@@ -775,7 +707,7 @@ namespace ttmath
 				sbb eax, ebx 
 				mov [edi+edx*4], eax
 
-				inc edx
+				lea edx, [edx+1]
 				dec ecx
 			jnz ttmath_loop2
 
@@ -783,8 +715,6 @@ namespace ttmath
 
 			ttmath_end:
 				mov [c], ecx
-
-				popad
 			}
 
 		#endif		
@@ -834,8 +764,6 @@ namespace ttmath
 
 		#endif
 
-		TTMATH_LOG("UInt::SubVector")
-
 	return c;
 	}
 
@@ -862,29 +790,25 @@ namespace ttmath
 		#ifndef __GNUC__
 			__asm
 			{
-				push ebx
-				push ecx
-				push edx
-
 				mov ebx, [p1]
+
 				xor edx, edx
-				mov ecx, [c]
+				mov ecx, edx
-				neg ecx
+				sub ecx, [c]
+
 				mov ecx, [b]
 
+				ALIGN 16
 			ttmath_loop:
 				rcl dword ptr [ebx+edx*4], 1
 				
-				inc edx
+				lea edx, [edx+1] // inc edx, but faster (no flags dependencies)
 				dec ecx
 			jnz ttmath_loop
 
-				adc ecx, ecx
+				setc	al
-				mov [c], ecx
+				movzx	eax, al
-				
+				mov		[c], eax
-				pop edx
-				pop ecx
-				pop ebx
 			}
 		#endif
 
@@ -940,25 +864,22 @@ namespace ttmath
 		#ifndef __GNUC__
 			__asm
 			{
-				push ebx
+				xor ecx, ecx
-				push ecx
+				sub ecx, [c]
 
 				mov ebx, [p1]
-				mov ecx, [c]
-				neg ecx
 				mov ecx, [b]
 
+				ALIGN 16
 			ttmath_loop:
 				rcr dword ptr [ebx+ecx*4-4], 1
 				
 				dec ecx
 			jnz ttmath_loop
 
-				adc ecx, ecx
+				setc	al
-				mov [c], ecx
+				movzx	eax, al
-
+				mov		[c], eax
-				pop ecx
-				pop ebx
 			}
 		#endif
 
@@ -991,13 +912,6 @@ namespace ttmath
 
 
 
-#ifdef _MSC_VER
-#pragma warning (disable : 4731)
-//warning C4731: frame pointer register 'ebp' modified by inline assembly code
-#endif
-	
-
-
 	/*!
 		this method moves all bits into the left hand side
 		return value <- this <- c
@@ -1015,62 +929,47 @@ namespace ttmath
 	{
 	TTMATH_ASSERT( bits>0 && bits<TTMATH_BITS_PER_UINT )
 		
-	uint b = value_size;
+	register sint b = value_size;
-	uint * p1 = table;
+	register uint * p1 = table;
+	register uint mask;
 
 		#ifndef __GNUC__
 			__asm
 			{
-				push eax
-				push ebx
-				push ecx
-				push edx
-				push esi
-				push edi
-				push ebp
-
 				mov edi, [b]
 
 				mov ecx, 32
 				sub ecx, [bits]
 				mov edx, -1
 				shr edx, cl
+				mov [mask], edx
 
 				mov ecx, [bits]
 				mov ebx, [p1]
-				mov eax, [c]
-
-				mov ebp, edx         // ebp = mask (modified ebp - don't read/write to variables)
 
 				xor edx, edx   // edx = 0
-				mov esi, edx
+				mov esi, edx   // old value = 0 
-				or eax, eax
-				cmovnz esi, ebp      // if(c) esi=mask else esi=0
 
+				mov eax, [c]
+				or eax, eax
+				cmovnz esi, [mask] // if c then old value = mask
+
+				ALIGN 16
 			ttmath_loop:
 				rol dword ptr [ebx+edx*4], cl
 				
 				mov eax, [ebx+edx*4]
-				and eax, ebp
+				and eax, [mask] 
 				xor [ebx+edx*4], eax // clearing bits
 				or [ebx+edx*4], esi  // saving old value
 				mov esi, eax
 
-				inc edx
+				lea edx, [edx+1] // inc edx, but faster (no flags dependencies)
 				dec edi
 			jnz ttmath_loop
 
-				pop ebp              // restoring ebp
-
 				and eax, 1
-				mov [c], eax
+				mov	dword ptr [c], eax
-
-				pop edi
-				pop esi
-				pop edx
-				pop ecx
-				pop ebx
-				pop eax
 			}
 		#endif
 
@@ -1145,43 +1044,37 @@ namespace ttmath
 
 	uint b = value_size;
 	uint * p1 = table;
+	uint mask;
 
 		#ifndef __GNUC__
 			__asm
 			{
-				push eax
-				push ebx
-				push ecx
-				push edx
-				push esi
-				push edi
-				push ebp
-
 				mov edi, [b]
 
 				mov ecx, 32
 				sub ecx, [bits]
 				mov edx, -1
 				shl edx, cl
+				mov [mask], edx
 
 				mov ecx, [bits]
 				mov ebx, [p1]
-				mov eax, [c]
-
-				mov ebp, edx         // ebp = mask (modified ebp - don't read/write to variables)
 
 				xor edx, edx   // edx = 0
-				mov esi, edx
+				mov esi, edx   // old value = 0 
 				add edx, edi   
-				dec edx              // edx is pointing at the end of the table (on last word)
+				dec edx        // edx - is pointing at the last word
-				or eax, eax
-				cmovnz esi, ebp      // if(c) esi=mask else esi=0
 
+				mov eax, [c]
+				or eax, eax
+				cmovnz esi, [mask] // if c then old value = mask
+
+				ALIGN 16
 			ttmath_loop:
 				ror dword ptr [ebx+edx*4], cl
 				
 				mov eax, [ebx+edx*4]
-				and eax, ebp 
+				and eax, [mask] 
 				xor [ebx+edx*4], eax // clearing bits
 				or [ebx+edx*4], esi  // saving old value
 				mov esi, eax
@@ -1190,18 +1083,10 @@ namespace ttmath
 				dec edi
 			jnz ttmath_loop
 
-				pop ebp              // restoring ebp
+				rol eax, 1    // bit 31 will be bit 0
-
-				rol eax, 1           // 31bit will be first
 				and eax, 1  
-				mov [c], eax
 				
-				pop edi
+				mov	dword ptr [c], eax
-				pop esi
-				pop edx
-				pop ecx
-				pop ebx
-				pop eax
 			}
 		#endif
 
@@ -1258,10 +1143,6 @@ namespace ttmath
 	}
 
 
-#ifdef _MSC_VER
-#pragma warning (default : 4731)
-#endif
-
 
 	/*
 		this method returns the number of the highest set bit in one 32-bit word
@@ -1275,16 +1156,11 @@ namespace ttmath
 		#ifndef __GNUC__
 			__asm
 			{
-				push eax
-				push edx
-
 				mov edx,-1
 				bsr eax,[x]
 				cmovz eax,edx
-				mov [result], eax
 				
-				pop edx
+				mov [result], eax
-				pop eax
 			}
 		#endif
 
@@ -1332,9 +1208,6 @@ namespace ttmath
 		#ifndef __GNUC__
 			__asm
 			{
-			push ebx
-			push eax
-
 				mov eax, [v]
 				mov ebx, [bit]
 				bts eax, ebx
@@ -1343,9 +1216,6 @@ namespace ttmath
 				setc bl
 				movzx ebx, bl
 				mov [old_bit], ebx
-
-			pop eax
-			pop ebx
 			}
 		#endif
 
@@ -1364,7 +1234,6 @@ namespace ttmath
 		#endif
 
 		value = v;
-
 		return old_bit;
 	}
 
@@ -1396,17 +1265,11 @@ namespace ttmath
 
 			__asm
 			{
-			push eax
-			push edx
-
 				mov eax, [a]
 				mul dword ptr [b]
 
 				mov [result2_], edx
 				mov [result1_], eax
-
-			pop edx
-			pop eax
 			}
 
 		#endif
@@ -1469,18 +1332,12 @@ namespace ttmath
 		#ifndef __GNUC__
 			__asm
 			{
-				push eax
-				push edx
-
 				mov edx, [a]
 				mov eax, [b]
 				div dword ptr [c]
 
 				mov [r_], eax
 				mov [rest_], edx
-
-				pop edx
-				pop eax
 			}
 		#endif
 

ttmath/ttmathuint_x86_64.h

@@ -1,7 +1,7 @@
 /*
  * This file is a part of TTMath Bignum Library
  * and is distributed under the (new) BSD licence.
- * Author: Tomasz Sowa <t.sowa@slimaczek.pl>
+ * Author: Tomasz Sowa <t.sowa@ttmath.org>
  */
 
 /*
@@ -39,11 +39,10 @@
 #ifndef headerfilettmathuint_x86_64
 #define headerfilettmathuint_x86_64
 
-
 #ifndef TTMATH_NOASM
 #ifdef TTMATH_PLATFORM64
 
-
+#pragma message("TTMATH_ASM64")
 /*!
 	\file ttmathuint_x86_64.h
     \brief template class UInt<uint> with assembler code for 64bit x86_64 processors
@@ -51,10 +50,31 @@
 	this file is included at the end of ttmathuint.h
 */
 
+#define WIN32_LEAN_AND_MEAN
+#include <windows.h>
+
 
 namespace ttmath
 {
 
+	#if defined(_M_X64)
+		#include <intrin.h>
+
+		extern "C"
+			{
+			uint	__fastcall	adc_x64(uint* p1, const uint* p2, uint nSize, uint c);
+			uint	__fastcall	addindexed_x64(uint* p1, uint nSize, uint nPos, uint nValue);
+			uint	__fastcall	addindexed2_x64(uint* p1, uint nSize, uint nPos, uint nValue1, uint nValue2);
+			uint	__fastcall	sbb_x64(uint* p1, const uint* p2, uint nSize, uint c);
+			uint	__fastcall	subindexed_x64(uint* p1, uint nSize, uint nPos, uint nValue);
+			uint	__fastcall	rcl_x64(uint* p1, uint nSize, uint nLowestBit);
+			uint	__fastcall	rcr_x64(uint* p1, uint nSize, uint nLowestBit);
+			uint	__fastcall	div_x64(uint* pnValHi, uint* pnValLo, uint nDiv);
+			uint	__fastcall	rcl2_x64(uint* p1, uint nSize, uint nBits, uint c);
+			uint	__fastcall	rcr2_x64(uint* p1, uint nSize, uint nBits, uint c);
+			};
+	#endif
+
 	/*!
 	*
 	*	basic mathematic functions
@@ -83,12 +103,15 @@ namespace ttmath
 		// this algorithm doesn't require it
 
 		#ifndef __GNUC__
+			#if defined(_M_X64)
+				c = adc_x64(p1,p2,b,c);
+			#else
 				#error "another compiler than GCC is currently not supported in 64bit mode"
 			#endif
+		#endif
 
 		#ifdef __GNUC__
 			uint dummy, dummy2;
-
 			/*
 				this part should be compiled with gcc
 			*/
@@ -150,8 +173,12 @@ namespace ttmath
 		TTMATH_ASSERT( index < value_size )
 
 		#ifndef __GNUC__
+			#if defined(_M_X64)
+				c = addindexed_x64(p1,b,index,value);
+			#else
 				#error "another compiler than GCC is currently not supported in 64bit mode"
 			#endif
+		#endif
 
 		#ifdef __GNUC__
 			uint dummy, dummy2;
@@ -220,6 +247,30 @@ namespace ttmath
 	*/
 	template<uint value_size>
 	uint UInt<value_size>::AddTwoInts(uint x2, uint x1, uint index)
+	#if 0
+	{
+	uint i, c;
+
+		TTMATH_ASSERT( index < value_size )
+
+		printf("add %Id + %Id\n",x1,x2);
+		for(int i=index ; i<value_size ; ++i)
+			printf("%d: %Id\n",i,table[i]);
+
+		c = AddTwoWords(table[index],   x1, 0, &table[index]);
+		c = AddTwoWords(table[index+1], x2, c, &table[index+1]);
+
+		for(i=index+2 ; i<value_size && c ; ++i)
+			c = AddTwoWords(table[i], 0, c, &table[i]);
+		for(i=index ; i<value_size ; ++i)
+			printf("%d: %Id\n",i,table[i]);
+		printf(" -> %d\n",c);
+
+		TTMATH_LOG("UInt::AddTwoInts")
+
+	return c;
+	}
+	#else
 	{
 	uint b = value_size;
 	uint * p1 = table;
@@ -228,8 +279,19 @@ namespace ttmath
 		TTMATH_ASSERT( index < value_size - 1 )
 
 		#ifndef __GNUC__
+			#if defined(_M_X64)
+				//printf("add %Id + %Id\n",x1,x2);
+				//for(int i=index ; i<value_size ; ++i)
+				//	printf("%d: %Id\n",i,table[i]);
+				//if (table[0] == 1265784741359897913) DebugBreak();
+				c = addindexed2_x64(p1,b,index,x1,x2);
+				//for(int i=index ; i<value_size ; ++i)
+				//	printf("%d: %Id\n",i,table[i]);
+				//printf(" -> %d\n",c);
+			#else
 				#error "another compiler than GCC is currently not supported in 64bit mode"
 			#endif
+		#endif
 
 		#ifdef __GNUC__
 			uint dummy, dummy2;
@@ -265,89 +327,8 @@ namespace ttmath
 
 	return c;
 	}
-
-
-
-	/*!
-		this static method addes one vector to the other
-		'ss1' is larger in size or equal to 'ss2'
-
-		ss1 points to the first (larger) vector
-		ss2 points to the second vector
-		ss1_size - size of the ss1 (and size of the result too)
-		ss2_size - size of the ss2
-		result - is the result vector (which has size the same as ss1: ss1_size)
-
-		Example:  ss1_size is 5, ss2_size is 3
-		ss1:      ss2:   result (output):
-		  5        1         5+1
-		  4        3         4+3
-		  2        7         2+7
-		  6                  6
-		  9                  9
-	  of course the carry is propagated and will be returned from the last item
-	  (this method is used by the Karatsuba multiplication algorithm)
-	*/
-	template<uint value_size>
-	uint UInt<value_size>::AddVector(const uint * ss1, const uint * ss2, uint ss1_size, uint ss2_size, uint * result)
-	{
-		TTMATH_ASSERT( ss1_size >= ss2_size )
-
-		uint rest = ss1_size - ss2_size;
-		uint c;
-
-		#ifndef __GNUC__
-			#error "another compiler than GCC is currently not supported in 64bit mode"
 	#endif
 
-		#ifdef __GNUC__
-		uint dummy1, dummy2, dummy3;	
-			
-			//	this part should be compiled with gcc
-		
-			__asm__ __volatile__(
-				"mov %%rdx, %%r8					\n"
-				"xor %%rdx, %%rdx					\n"   // rdx = 0, cf = 0
-			"1:										\n"
-				"mov (%%rsi,%%rdx,8), %%rax			\n"
-				"adc (%%rbx,%%rdx,8), %%rax			\n"
-				"mov %%rax, (%%rdi,%%rdx,8)			\n"
-
-				"inc %%rdx							\n"
-				"dec %%rcx							\n"
-			"jnz 1b									\n"
-
-				"adc %%rcx, %%rcx					\n"   // rcx has the cf state
-
-				"or %%r8, %%r8						\n"
-				"jz 3f								\n"
-				
-				"xor %%rbx, %%rbx					\n"   // ebx = 0
-				"neg %%rcx							\n"   // setting cf from rcx
-				"mov %%r8, %%rcx					\n"   // rcx=rest and is != 0
-			"2:										\n"
-				"mov (%%rsi, %%rdx, 8), %%rax		\n"
-				"adc %%rbx, %%rax 					\n"
-				"mov %%rax, (%%rdi, %%rdx, 8)		\n"
-
-				"inc %%rdx							\n"
-				"dec %%rcx							\n"
-			"jnz 2b									\n"
-
-				"adc %%rcx, %%rcx					\n"
-			"3:										\n"
-
-				: "=a" (dummy1), "=b" (dummy2), "=c" (c),       "=d" (dummy3)
-				:                "1" (ss2),     "2" (ss2_size), "3" (rest),   "S" (ss1),  "D" (result)
-				: "%r8", "cc", "memory" );
-
-		#endif
-
-		TTMATH_LOG("UInt::AddVector")
-
-	return c;
-	}
-
 
 
 
@@ -369,13 +350,16 @@ namespace ttmath
 	uint * p1 = table;
 	const uint * p2 = ss2.table;
 
-
 		// we don't have to use TTMATH_REFERENCE_ASSERT here
 		// this algorithm doesn't require it
 
 		#ifndef __GNUC__
+			#if defined(_M_X64)
+				c = sbb_x64(p1,p2,b,c);
+			#else
 				#error "another compiler than GCC is currently not supported in 64bit mode"
 			#endif
+		#endif
 
 		#ifdef __GNUC__
 			uint dummy, dummy2;
@@ -399,6 +383,7 @@ namespace ttmath
 				: "0" (b), "1" (c), "b" (p1), "S" (p2)
 				: "cc", "memory" );
 
+
 		#endif
 
 		TTMATH_LOG("UInt::Sub")
@@ -432,15 +417,20 @@ namespace ttmath
 	uint b = value_size;
 	uint * p1 = table;
 	uint c;
-	uint dummy, dummy2;
 
 		TTMATH_ASSERT( index < value_size )
 
 		#ifndef __GNUC__
+			#if defined(_M_X64)
+				c = subindexed_x64(p1,b,index,value);
+			#else
 				#error "another compiler than GCC is currently not supported in 64bit mode"
 			#endif
+		#endif
 
 		#ifdef __GNUC__
+			uint dummy, dummy2;
+
 			__asm__ __volatile__(
 
 				"subq %%rdx, %%rcx 				\n"
@@ -464,100 +454,12 @@ namespace ttmath
 
 		#endif
 
-		TTMATH_LOG("UInt::SubInt")
+		TTMATH_LOG("UInt64::SubInt")
 
 	return c;
 	}
 
 
-
-
-	/*!
-		this static method subtractes one vector from the other
-		'ss1' is larger in size or equal to 'ss2'
-
-		ss1 points to the first (larger) vector
-		ss2 points to the second vector
-		ss1_size - size of the ss1 (and size of the result too)
-		ss2_size - size of the ss2
-		result - is the result vector (which has size the same as ss1: ss1_size)
-
-		Example:  ss1_size is 5, ss2_size is 3
-		ss1:      ss2:   result (output):
-		  5        1         5-1
-		  4        3         4-3
-		  2        7         2-7
-		  6                  6-1  (the borrow from previous item)
-		  9                  9
-		               return (carry): 0
-	  of course the carry (borrow) is propagated and will be returned from the last item
-	  (this method is used by the Karatsuba multiplication algorithm)
-	*/
-	template<uint value_size>
-	uint UInt<value_size>::SubVector(const uint * ss1, const uint * ss2, uint ss1_size, uint ss2_size, uint * result)
-	{
-		TTMATH_ASSERT( ss1_size >= ss2_size )
-
-		uint rest = ss1_size - ss2_size;
-		uint c;
-
-		#ifndef __GNUC__
-			#error "another compiler than GCC is currently not supported in 64bit mode"
-		#endif
-
-		#ifdef __GNUC__
-
-		/*
-			the asm code is nearly the same as in AddVector
-			only two instructions 'adc' are changed to 'sbb'
-		*/
-		uint dummy1, dummy2, dummy3;
-
-			__asm__ __volatile__(
-				"mov %%rdx, %%r8					\n"
-				"xor %%rdx, %%rdx					\n"   // rdx = 0, cf = 0
-			"1:										\n"
-				"mov (%%rsi,%%rdx,8), %%rax			\n"
-				"sbb (%%rbx,%%rdx,8), %%rax			\n"
-				"mov %%rax, (%%rdi,%%rdx,8)			\n"
-
-				"inc %%rdx							\n"
-				"dec %%rcx							\n"
-			"jnz 1b									\n"
-
-				"adc %%rcx, %%rcx					\n"   // rcx has the cf state
-
-				"or %%r8, %%r8						\n"
-				"jz 3f								\n"
-				
-				"xor %%rbx, %%rbx					\n"   // ebx = 0
-				"neg %%rcx							\n"   // setting cf from rcx
-				"mov %%r8, %%rcx					\n"   // rcx=rest and is != 0
-			"2:										\n"
-				"mov (%%rsi, %%rdx, 8), %%rax		\n"
-				"sbb %%rbx, %%rax 					\n"
-				"mov %%rax, (%%rdi, %%rdx, 8)		\n"
-
-				"inc %%rdx							\n"
-				"dec %%rcx							\n"
-			"jnz 2b									\n"
-
-				"adc %%rcx, %%rcx					\n"
-			"3:										\n"
-
-				: "=a" (dummy1), "=b" (dummy2), "=c" (c),       "=d" (dummy3)
-				:                "1" (ss2),     "2" (ss2_size), "3" (rest),   "S" (ss1),  "D" (result)
-				: "%r8", "cc", "memory" );
-
-		#endif
-
-		TTMATH_LOG("UInt::SubVector")
-
-	return c;
-	}
-
-
-
 	/*!
 		this method moves all bits into the left hand side
 		return value <- this <- c
@@ -578,10 +480,13 @@ namespace ttmath
 	sint b = value_size;
 	uint * p1 = table;
 
-
 		#ifndef __GNUC__
+			#if defined(_M_X64)
+				c = rcl_x64(p1,b,c);
+			#else
 				#error "another compiler than GCC is currently not supported in 64bit mode"
 			#endif
+		#endif
 
 		#ifdef __GNUC__
 		uint dummy, dummy2;
@@ -632,10 +537,13 @@ namespace ttmath
 	sint b = value_size;
 	uint * p1 = table;
 
-
 		#ifndef __GNUC__
+			#if defined(_M_X64)
+				c = rcr_x64(p1,b,c);
+			#else
 				#error "another compiler than GCC is currently not supported in 64bit mode"
 			#endif
+		#endif
 
 		#ifdef __GNUC__
 		uint dummy;
@@ -687,10 +595,13 @@ namespace ttmath
 	uint b = value_size;
 	uint * p1 = table;
 
-
 		#ifndef __GNUC__
+			#if defined(_M_X64)
+				c = rcl2_x64(p1,b,bits,c);
+			#else
 				#error "another compiler than GCC is currently not supported in 64bit mode"
 			#endif
+		#endif
 
 		#ifdef __GNUC__
 		uint dummy, dummy2, dummy3;
@@ -707,6 +618,7 @@ namespace ttmath
 
 			"xorq %%rdx, %%rdx				\n"
 			"movq %%rdx, %%rsi				\n"
+
 			"orq %%rax, %%rax				\n"
 			"cmovnz %%r8, %%rsi				\n"
 
@@ -758,14 +670,19 @@ namespace ttmath
 
 	sint b = value_size;
 	uint * p1 = table;
-	uint dummy, dummy2, dummy3;
 
 		#ifndef __GNUC__
+			#if defined(_M_X64)
+				c = rcr2_x64(p1,b,bits,c);
+			#else
 				#error "another compiler than GCC is currently not supported in 64bit mode"
 			#endif
+		#endif
 
 
 		#ifdef __GNUC__
+			uint dummy, dummy2, dummy3;
+
 			__asm__  __volatile__(
 
 			"movq %%rcx, %%rsi				\n"
@@ -780,6 +697,7 @@ namespace ttmath
 			"movq %%rdx, %%rsi				\n"
 			"addq %%rdi, %%rdx				\n"
 			"decq %%rdx						\n"
+
 			"orq %%rax, %%rax				\n"
 			"cmovnz %%R8, %%rsi				\n"
 
@@ -820,12 +738,20 @@ namespace ttmath
 	template<uint value_size>
 	sint UInt<value_size>::FindLeadingBitInWord(uint x)
 	{
-	sint result;
+	register sint result;
-
 
 		#ifndef __GNUC__
+			#if defined(_MSC_VER)
+				unsigned long	nIndex(0);
+
+				if (_BitScanReverse64(&nIndex,x) == 0)
+					result = -1;
+				  else
+					result = nIndex;
+			#else
 				#error "another compiler than GCC is currently not supported in 64bit mode"
 			#endif
+		#endif
 
 		#ifdef __GNUC__
 		uint dummy;
@@ -870,8 +796,16 @@ namespace ttmath
 
 
 		#ifndef __GNUC__
+			#if defined(_MSC_VER)
+				#if defined(TTMATH_PLATFORM64)
+					old_bit = _bittestandset64((__int64*)&value,bit) != 0;
+				  #else
+					old_bit = _bittestandset((long*)&value,bit) != 0;
+				#endif
+			#else
 				#error "another compiler than GCC is currently not supported in 64bit mode"
 			#endif
+		#endif
 
 		#ifdef __GNUC__
 
@@ -925,8 +859,12 @@ namespace ttmath
 	uint result2_;
 
 		#ifndef __GNUC__
+			#if defined(_MSC_VER)
+				result1_ = _umul128(a,b,&result2_);
+			#else
 				#error "another compiler than GCC is currently not supported in 64bit mode"
 			#endif
+		#endif
 
 		#ifdef __GNUC__
 
@@ -955,6 +893,7 @@ namespace ttmath
 	 *
 	*/
 
+	#ifndef __GNUC__
 
 	/*!
 		this method calculates 64bits word a:b / 32bits c (a higher, b lower word)
@@ -982,8 +921,14 @@ namespace ttmath
 		TTMATH_ASSERT( c != 0 )
 
 		#ifndef __GNUC__
+			#if defined(_MSC_VER)
+				div_x64(&a,&b,c);
+				r_ = a;
+				rest_ = b;
+			#else
 				#error "another compiler than GCC is currently not supported in 64bit mode"
 			#endif
+		#endif
 
 		#ifdef __GNUC__
 		
@@ -1002,6 +947,132 @@ namespace ttmath
 		*rest = rest_;
 	}
 
+	template<uint value_size>
+	uint UInt<value_size>::AddTwoWords(uint a, uint b, uint carry, uint * result)
+	{
+	uint temp;
+
+		if( carry == 0 )
+		{
+			temp = a + b;
+
+			if( temp < a )
+				carry = 1;
+		}
+		else
+		{
+			carry = 1;
+			temp  = a + b + carry;
+
+			if( temp > a ) // !(temp<=a)
+				carry = 0;
+		}
+
+		*result = temp;
+
+	return carry;
+	}
+
+	template<uint value_size>
+	uint UInt<value_size>::SubTwoWords(uint a, uint b, uint carry, uint * result)
+	{
+		if( carry == 0 )
+		{
+			*result = a - b;
+
+			if( a < b )
+				carry = 1;
+		}
+		else
+		{
+			carry   = 1;
+			*result = a - b - carry;
+
+			if( a > b ) // !(a <= b )
+				carry = 0;
+		}
+
+	return carry;
+	}
+
+	/*!
+		this static method addes one vector to the other
+		'ss1' is larger in size or equal to 'ss2'
+
+		ss1 points to the first (larger) vector
+		ss2 points to the second vector
+		ss1_size - size of the ss1 (and size of the result too)
+		ss2_size - size of the ss2
+		result - is the result vector (which has size the same as ss1: ss1_size)
+
+		Example:  ss1_size is 5, ss2_size is 3
+		ss1:      ss2:   result (output):
+		  5        1         5+1
+		  4        3         4+3
+		  2        7         2+7
+		  6                  6
+		  9                  9
+	  of course the carry is propagated and will be returned from the last item
+	  (this method is used by the Karatsuba multiplication algorithm)
+	*/
+	template<uint value_size>
+	uint UInt<value_size>::AddVector(const uint * ss1, const uint * ss2, uint ss1_size, uint ss2_size, uint * result)
+	{
+	uint i, c = 0;
+
+		TTMATH_ASSERT( ss1_size >= ss2_size )
+
+		for(i=0 ; i<ss2_size ; ++i)
+			c = AddTwoWords(ss1[i], ss2[i], c, &result[i]);
+
+		for( ; i<ss1_size ; ++i)
+			c = AddTwoWords(ss1[i], 0, c, &result[i]);
+
+		TTMATH_LOG("UInt::AddVector")
+
+	return c;
+	}
+
+	/*!
+		this static method subtractes one vector from the other
+		'ss1' is larger in size or equal to 'ss2'
+
+		ss1 points to the first (larger) vector
+		ss2 points to the second vector
+		ss1_size - size of the ss1 (and size of the result too)
+		ss2_size - size of the ss2
+		result - is the result vector (which has size the same as ss1: ss1_size)
+
+		Example:  ss1_size is 5, ss2_size is 3
+		ss1:      ss2:   result (output):
+		  5        1         5-1
+		  4        3         4-3
+		  2        7         2-7
+		  6                  6-1  (the borrow from previous item)
+		  9                  9
+		                 return (carry): 0
+	  of course the carry (borrow) is propagated and will be returned from the last item
+	  (this method is used by the Karatsuba multiplication algorithm)
+	*/
+	template<uint value_size>
+	uint UInt<value_size>::SubVector(const uint * ss1, const uint * ss2, uint ss1_size, uint ss2_size, uint * result)
+	{
+	uint i, c = 0;
+
+		TTMATH_ASSERT( ss1_size >= ss2_size )
+
+		for(i=0 ; i<ss2_size ; ++i)
+			c = SubTwoWords(ss1[i], ss2[i], c, &result[i]);
+
+		for( ; i<ss1_size ; ++i)
+			c = SubTwoWords(ss1[i], 0, c, &result[i]);
+
+		TTMATH_LOG("UInt::SubVector")
+
+	return c;
+	}
+
+	#endif // #ifndef __GNUC__
 
 } //namespace
 

ttmath/ttmathuint_x86_amd64_msvc.asm

@@ -0,0 +1,386 @@
+PUBLIC	adc_x64
+PUBLIC	addindexed_x64
+PUBLIC	addindexed2_x64
+
+PUBLIC	sbb_x64
+PUBLIC	subindexed_x64
+
+PUBLIC	rcl_x64
+PUBLIC	rcr_x64
+
+PUBLIC	rcl2_x64
+PUBLIC	rcr2_x64
+
+PUBLIC	div_x64
+
+;
+;	"rax, rcx, rdx, r8-r11 are volatile."
+;	"rbx, rbp, rdi, rsi, r12-r15 are nonvolatile."
+;
+
+.CODE
+
+        ALIGN       8
+
+;----------------------------------------
+
+adc_x64				PROC
+        ; rcx = p1
+        ; rdx = p2
+        ; r8 = nSize
+        ; r9 = nCarry
+
+        xor		rax, rax
+        xor		r11, r11
+        sub		rax, r9		; sets CARRY if r9 != 0
+
+		ALIGN 16
+ loop1:
+		mov		rax,qword ptr [rdx + r11 * 8]
+		adc		qword ptr [rcx + r11 * 8], rax
+		lea		r11, [r11+1]
+		dec		r8
+		jnz		loop1
+
+		setc	al
+		movzx	rax, al
+
+		ret
+
+adc_x64				ENDP
+
+;----------------------------------------
+
+        ALIGN       8
+
+;----------------------------------------
+
+addindexed_x64	PROC
+
+        ; rcx = p1
+        ; rdx = nSize
+        ; r8 = nPos
+        ; r9 = nValue
+
+		xor		rax, rax			; rax = result
+		sub		rdx, r8				; rdx = remaining count of uints
+
+		add		qword ptr [rcx + r8 * 8], r9
+		jc		next1
+
+		ret
+
+next1:
+		mov		r9, 1
+
+		ALIGN 16
+loop1:
+		dec		rdx
+		jz		done_with_cy
+		lea		r8, [r8+1]
+		add		qword ptr [rcx + r8 * 8], r9
+		jc		loop1
+
+		ret
+
+done_with_cy:
+		lea		rax, [rax+1]		; rax = 1
+
+		ret
+
+addindexed_x64	ENDP
+
+;----------------------------------------
+
+        ALIGN       8
+
+;----------------------------------------
+
+addindexed2_x64	PROC
+
+        ; rcx = p1 (pointer)
+        ; rdx = b  (value size)
+        ; r8 = nPos
+        ; r9 = nValue1
+        ; [esp+0x28] = nValue2
+
+		xor		rax, rax			; return value
+		mov		r11, rcx			; table
+		sub		rdx, r8				; rdx = remaining count of uints
+		mov		r10, [esp+028h]		; r10 = nValue2
+
+		add		qword ptr [r11 + r8 * 8], r9
+		lea		r8, [r8+1]
+		lea		rdx, [rdx-1]
+		adc		qword ptr [r11 + r8 * 8], r10
+		jc		next
+		ret
+
+		ALIGN 16
+loop1:
+		lea		r8, [r8+1]
+		add		qword ptr [r11 + r8 * 8], 1
+		jc		next
+		ret
+
+next:
+		dec		rdx					; does not modify CY too...
+		jnz		loop1
+		lea		rax, [rax+1]
+		ret
+
+addindexed2_x64	ENDP
+
+;----------------------------------------
+
+        ALIGN       8
+
+;----------------------------------------
+
+sbb_x64				PROC
+
+        ; rcx = p1
+        ; rdx = p2
+        ; r8 = nCount
+        ; r9 = nCarry
+
+        xor		rax, rax
+        xor		r11, r11
+        sub		rax, r9				; sets CARRY if r9 != 0
+
+		ALIGN 16
+ loop1:
+		mov		rax,qword ptr [rdx + r11 * 8]
+		sbb		qword ptr [rcx + r11 * 8], rax
+		lea		r11, [r11+1]
+		dec		r8
+		jnz		loop1
+
+		setc	al
+		movzx	rax, al
+
+		ret
+
+sbb_x64				ENDP
+
+;----------------------------------------
+
+        ALIGN       8
+
+;----------------------------------------
+
+subindexed_x64	PROC
+        ; rcx = p1
+        ; rdx = nSize
+        ; r8 = nPos
+        ; r9 = nValue
+
+		sub		rdx, r8				; rdx = remaining count of uints
+
+		ALIGN 16
+loop1:
+		sub		qword ptr [rcx + r8 * 8], r9
+		jnc		done
+
+		lea		r8, [r8+1]
+		mov		r9, 1
+		dec		rdx
+		jnz		loop1
+		jc		return_1	; most of the times, there will be NO carry (I hope)
+
+done:
+		xor		rax, rax
+		ret
+
+  return_1:
+		mov		rax, 1
+		ret
+
+subindexed_x64	ENDP
+
+;----------------------------------------
+
+        ALIGN       8
+
+;----------------------------------------
+
+rcl_x64	PROC
+        ; rcx = p1
+        ; rdx = b
+        ; r8 = nLowestBit
+
+		mov		r11, rcx			; table
+		xor		r10, r10
+		neg		r8					; CY set if r8 <> 0
+
+		ALIGN 16
+loop1:
+		rcl		qword ptr [r11 + r10 * 8], 1
+		lea		r10, [r10+1]
+		dec		rdx
+		jnz		loop1
+
+		setc	al
+		movzx	rax, al
+
+        ret
+
+rcl_x64	ENDP
+
+;----------------------------------------
+
+        ALIGN       8
+
+;----------------------------------------
+
+rcr_x64	PROC
+        ; rcx = p1
+        ; rdx = nSize
+        ; r8 = nLowestBit
+
+		xor		r10, r10
+		neg		r8					; CY set if r8 <> 0
+
+		ALIGN 16
+loop1:
+		rcr		qword ptr -8[rcx + rdx * 8], 1
+		dec		rdx
+		jnz		loop1
+
+		setc	al
+		movzx	rax, al
+
+        ret
+
+rcr_x64	ENDP
+
+;----------------------------------------
+
+        ALIGN       8
+
+;----------------------------------------
+
+div_x64	PROC
+
+        ; rcx = &Hi
+        ; rdx = &Lo
+        ; r8 = nDiv
+
+        mov		r11, rcx
+        mov		r10, rdx
+
+        mov		rdx, qword ptr [r11]
+        mov		rax, qword ptr [r10]
+        div		r8
+        mov		qword ptr [r10], rdx ; remainder
+        mov		qword ptr [r11], rax ; value
+
+        ret
+
+div_x64	ENDP
+
+;----------------------------------------
+
+        ALIGN       8
+
+;----------------------------------------
+
+rcl2_x64	PROC
+        ; rcx = p1
+        ; rdx = nSize
+        ; r8 = bits
+        ; r9 = c
+
+        push	rbx
+
+        mov		r10, rcx	; r10 = p1
+        xor		rax, rax
+
+        mov		rcx, 64
+        sub		rcx, r8
+
+        mov		r11, -1
+        shr		r11, cl		; r11 = mask
+
+		mov		rcx, r8		; rcx = count of bits
+
+		mov		rbx, rax	; rbx = old value = 0
+		or		r9, r9
+		cmovnz	rbx, r11	; if (c) then old value = mask
+
+        mov		r9, rax		; r9 = index (0..nSize-1)
+
+		ALIGN 16
+loop1:
+		rol		qword ptr [r10+r9*8], cl
+		mov		rax, qword ptr [r10+r9*8]
+		and		rax, r11
+		xor		qword ptr [r10+r9*8], rax
+		or		qword ptr [r10+r9*8], rbx
+		mov		rbx, rax
+
+		lea		r9, [r9+1]
+		dec		rdx
+
+		jnz		loop1
+
+		and		rax, 1
+		pop		rbx
+        ret
+
+rcl2_x64	ENDP
+
+;----------------------------------------
+
+        ALIGN       8
+
+;----------------------------------------
+
+rcr2_x64	PROC
+        ; rcx = p1
+        ; rdx = nSize
+        ; r8 = bits
+        ; r9 = c
+
+        push	rbx
+        mov		r10, rcx	; r10 = p1
+        xor		rax, rax
+
+        mov		rcx, 64
+        sub		rcx, r8
+
+        mov		r11, -1
+        shl		r11, cl		; r11 = mask
+
+		mov		rcx, r8		; rcx = count of bits
+
+		mov		rbx, rax	; rbx = old value = 0
+		or		r9, r9
+		cmovnz	rbx, r11	; if (c) then old value = mask
+
+        mov		r9, rdx		; r9 = index (0..nSize-1)
+		lea		r9, [r9-1]
+
+		ALIGN 16
+loop1:
+		ror		qword ptr [r10+r9*8], cl
+		mov		rax, qword ptr [r10+r9*8]
+		and		rax, r11
+		xor		qword ptr [r10+r9*8], rax
+		or		qword ptr [r10+r9*8], rbx
+		mov		rbx, rax
+
+		lea		r9, [r9-1]
+		dec		rdx
+
+		jnz		loop1
+
+		rol		rax, 1
+		and		rax, 1
+		pop		rbx
+
+        ret
+
+rcr2_x64	ENDP
+
+END