* added: global Gamma() function
* added: gamma() function to the parser * added: Big::IsInteger() method returns true if the value is integer * 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 git-svn-id: svn://ttmath.org/publicrep/ttmath/trunk@178 e52654a7-88a9-db11-a3e9-0013d4bc506e
This commit is contained in:
parent
d3a64b79ca
commit
53547cfab5
25
CHANGELOG
25
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)
|
||||
|
|
808
ttmath/ttmath.h
808
ttmath/ttmath.h
|
@ -73,7 +73,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
|
||||
|
@ -90,7 +90,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
|
||||
|
@ -222,7 +222,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)
|
||||
|
@ -263,7 +263,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)
|
||||
|
@ -304,7 +304,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)
|
||||
|
@ -1937,10 +1937,7 @@ namespace ttmath
|
|||
template<class ValueType>
|
||||
bool RootCheckIndexFrac(ValueType & x, const ValueType & index, ErrorCode * err)
|
||||
{
|
||||
ValueType indexfrac(index);
|
||||
indexfrac.RemainFraction();
|
||||
|
||||
if( !indexfrac.IsZero() )
|
||||
if( !index.IsInteger() )
|
||||
{
|
||||
// index must be integer
|
||||
if( err )
|
||||
|
@ -2072,154 +2069,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
|
||||
|
@ -2279,6 +2128,651 @@ 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
|
||||
}
|
||||
|
||||
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
|
||||
|
||||
|
||||
|
|
|
@ -1368,10 +1368,7 @@ public:
|
|||
|
||||
if( pow.exponent>-int(man*TTMATH_BITS_PER_UINT) && pow.exponent<=0 )
|
||||
{
|
||||
Big<exp, man> pow_frac( pow );
|
||||
pow_frac.RemainFraction();
|
||||
|
||||
if( pow_frac.IsZero() )
|
||||
if( pow.IsInteger() )
|
||||
return PowInt( pow );
|
||||
}
|
||||
|
||||
|
@ -4098,6 +4095,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
|
||||
|
|
|
@ -47,6 +47,7 @@
|
|||
#include "ttmathtypes.h"
|
||||
|
||||
#include <string>
|
||||
#include <vector>
|
||||
#include <list>
|
||||
#include <map>
|
||||
|
||||
|
@ -431,7 +432,7 @@ public:
|
|||
*/
|
||||
History()
|
||||
{
|
||||
buffer_max_size = 10;
|
||||
buffer_max_size = 15;
|
||||
}
|
||||
|
||||
|
||||
|
@ -488,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
|
||||
|
|
|
@ -137,7 +137,6 @@ namespace ttmath
|
|||
template<class ValueType>
|
||||
class Parser
|
||||
{
|
||||
|
||||
private:
|
||||
|
||||
/*!
|
||||
|
@ -427,10 +426,9 @@ VariablesTable variables_table;
|
|||
|
||||
|
||||
/*!
|
||||
you can't calculate the factorial if the argument is greater than 'factorial_max'
|
||||
default value is zero which means there are not any limitations
|
||||
some coefficients used when calculating the gamma (or factorial) function
|
||||
*/
|
||||
ValueType factorial_max;
|
||||
CGamma<ValueType> cgamma;
|
||||
|
||||
|
||||
/*!
|
||||
|
@ -676,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
|
||||
|
@ -686,11 +698,8 @@ void Factorial(int sindex, int amount_of_args, ValueType & result)
|
|||
Error( err_improper_amount_of_arguments );
|
||||
|
||||
ErrorCode err;
|
||||
|
||||
if( !factorial_max.IsZero() && stack[sindex].value > factorial_max )
|
||||
Error( err_too_big_factorial );
|
||||
|
||||
result = ttmath::Factorial(stack[sindex].value, &err, pstop_calculating);
|
||||
result = ttmath::Factorial(stack[sindex].value, cgamma, &err, pstop_calculating);
|
||||
|
||||
if(err != err_ok)
|
||||
Error( err );
|
||||
|
@ -1471,6 +1480,7 @@ void InsertVariableToTable(const tt_char * variable_name, pfunction_var pf)
|
|||
*/
|
||||
void CreateFunctionsTable()
|
||||
{
|
||||
InsertFunctionToTable(TTMATH_TEXT("gamma"), &Parser<ValueType>::Gamma);
|
||||
InsertFunctionToTable(TTMATH_TEXT("factorial"), &Parser<ValueType>::Factorial);
|
||||
InsertFunctionToTable(TTMATH_TEXT("abs"), &Parser<ValueType>::Abs);
|
||||
InsertFunctionToTable(TTMATH_TEXT("sin"), &Parser<ValueType>::Sin);
|
||||
|
@ -2419,7 +2429,6 @@ Parser(): default_stack_size(100)
|
|||
base = 10;
|
||||
deg_rad_grad = 1;
|
||||
error = err_ok;
|
||||
factorial_max.SetZero();
|
||||
|
||||
CreateFunctionsTable();
|
||||
CreateVariablesTable();
|
||||
|
@ -2439,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.
|
||||
|
@ -2521,17 +2529,6 @@ 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
|
||||
|
@ -2559,11 +2556,11 @@ return error;
|
|||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
} // namespace
|
||||
|
||||
|
||||
|
|
|
@ -54,7 +54,7 @@
|
|||
|
||||
#include <stdexcept>
|
||||
#include <sstream>
|
||||
|
||||
#include <vector>
|
||||
|
||||
/*!
|
||||
the version of the library
|
||||
|
@ -279,6 +279,19 @@ 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
|
||||
{
|
||||
|
||||
|
@ -312,7 +325,6 @@ namespace ttmath
|
|||
err_object_exists,
|
||||
err_unknown_object,
|
||||
err_still_calculating,
|
||||
err_too_big_factorial,
|
||||
err_in_short_form_used_function
|
||||
};
|
||||
|
||||
|
@ -490,7 +502,7 @@ namespace ttmath
|
|||
PrintLog(msg, std::cout);
|
||||
#endif
|
||||
|
||||
#define TTMATH_LOG(quote) TTMATH_LOG_HELPER(quote)
|
||||
#define TTMATH_LOG(msg) TTMATH_LOG_HELPER(msg)
|
||||
|
||||
#else
|
||||
|
||||
|
|
Loading…
Reference in New Issue