tomasz.sowa
/
ttmath
1
0
Fork 0
Code Issues Pull Requests Releases Activity
ttmath/ttmath/ttmath.h

2796 lines
55 KiB
C++
Raw Blame History

/*
* This file is a part of TTMath Bignum Library
* and is distributed under the (new) BSD licence.
* Author: Tomasz Sowa <t.sowa@ttmath.org>
*/
/*
* Copyright (c) 2006-2009, Tomasz Sowa
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* * Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
*
* * Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* * Neither the name Tomasz Sowa nor the names of contributors to this
* project may be used to endorse or promote products derived
* from this software without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
* THE POSSIBILITY OF SUCH DAMAGE.
*/
#ifndef headerfilettmathmathtt
#define headerfilettmathmathtt
/*!
\file ttmath.h
\brief Mathematics functions.
*/
#ifdef _MSC_VER
//warning C4127: conditional expression is constant
#pragma warning( disable: 4127 )
#endif
#include "ttmathbig.h"
#include "ttmathobjects.h"
namespace ttmath
{
/*
*
* functions defined here are used only with Big<> types
*
*
*/
/*
*
* functions for rounding
*
*
*/
/*!
this function skips the fraction from x
e.g 2.2 = 2
2.7 = 2
-2.2 = 2
-2.7 = 2
*/
template<class ValueType>
ValueType SkipFraction(const ValueType & x)
{
ValueType result( x );
result.SkipFraction();
return result;
}
/*!
this function rounds to the nearest integer value
e.g 2.2 = 2
2.7 = 3
-2.2 = -2
-2.7 = -3
*/
template<class ValueType>
ValueType Round(const ValueType & x, ErrorCode * err = 0)
{
if( x.IsNan() )
{
if( err )
*err = err_improper_argument;
return x; // NaN
}
ValueType result( x );
uint c = result.Round();
if( err )
*err = c ? err_overflow : err_ok;
return result;
}
/*!
this function returns a value representing the smallest integer
that is greater than or equal to x
Ceil(-3.7) = -3
Ceil(-3.1) = -3
Ceil(-3.0) = -3
Ceil(4.0) = 4
Ceil(4.2) = 5
Ceil(4.8) = 5
*/
template<class ValueType>
ValueType Ceil(const ValueType & x, ErrorCode * err = 0)
{
if( x.IsNan() )
{
if( err )
*err = err_improper_argument;
return x; // NaN
}
ValueType result(x);
uint c = 0;
result.SkipFraction();
if( result != x )
{
// x is with fraction
// if x is negative we don't have to do anything
if( !x.IsSign() )
{
ValueType one;
one.SetOne();
c += result.Add(one);
}
}
if( err )
*err = c ? err_overflow : err_ok;
return result;
}
/*!
this function returns a value representing the largest integer
that is less than or equal to x
Floor(-3.6) = -4
Floor(-3.1) = -4
Floor(-3) = -3
Floor(2) = 2
Floor(2.3) = 2
Floor(2.8) = 2
*/
template<class ValueType>
ValueType Floor(const ValueType & x, ErrorCode * err = 0)
{
if( x.IsNan() )
{
if( err )
*err = err_improper_argument;
return x; // NaN
}
ValueType result(x);
uint c = 0;
result.SkipFraction();
if( result != x )
{
// x is with fraction
// if x is positive we don't have to do anything
if( x.IsSign() )
{
ValueType one;
one.SetOne();
c += result.Sub(one);
}
}
if( err )
*err = c ? err_overflow : err_ok;
return result;
}
/*
*
* logarithms and the exponent
*
*
*/
/*!
this function calculates the natural logarithm (logarithm with the base 'e')
*/
template<class ValueType>
ValueType Ln(const ValueType & x, ErrorCode * err = 0)
{
if( x.IsNan() )
{
if( err )
*err = err_improper_argument;
return x; // NaN
}
ValueType result;
uint state = result.Ln(x);
if( err )
{
switch( state )
{
case 0:
*err = err_ok;
break;
case 1:
*err = err_overflow;
break;
case 2:
*err = err_improper_argument;
break;
default:
*err = err_internal_error;
break;
}
}
return result;
}
/*!
this function calculates the logarithm
*/
template<class ValueType>
ValueType Log(const ValueType & x, const ValueType & base, ErrorCode * err = 0)
{
if( x.IsNan() || base.IsNan() )
{
if( err )
*err = err_improper_argument;
return ValueType(); // default NaN
}
ValueType result;
uint state = result.Log(x, base);
if( err )
{
switch( state )
{
case 0:
*err = err_ok;
break;
case 1:
*err = err_overflow;
break;
case 2:
case 3:
*err = err_improper_argument;
break;
default:
*err = err_internal_error;
break;
}
}
return result;
}
/*!
this function calculates the expression e^x
*/
template<class ValueType>
ValueType Exp(const ValueType & x, ErrorCode * err = 0)
{
if( x.IsNan() )
{
if( err )
*err = err_improper_argument;
return x; // NaN
}
ValueType result;
uint c = result.Exp(x);
if( err )
*err = c ? err_overflow : err_ok;
return result;
}
/*!
*
* trigonometric functions
*
*/
/*
this namespace consists of auxiliary functions
(something like 'private' in a class)
*/
namespace auxiliaryfunctions
{
/*!
an auxiliary function for calculating the Sine
(you don't have to call this function)
*/
template<class ValueType>
uint PrepareSin(ValueType & x, bool & change_sign)
{
ValueType temp;
change_sign = false;
if( x.IsSign() )
{
// we're using the formula 'sin(-x) = -sin(x)'
change_sign = !change_sign;
x.ChangeSign();
}
// we're reducing the period 2*PI
// (for big values there'll always be zero)
temp.Set2Pi();
if( x.Mod(temp) )
return 1;
// we're setting 'x' as being in the range of <0, 0.5PI>
temp.SetPi();
if( x > temp )
{
// x is in (pi, 2*pi>
x.Sub( temp );
change_sign = !change_sign;
}
temp.Set05Pi();
if( x > temp )
{
// x is in (0.5pi, pi>
x.Sub( temp );
x = temp - x;
}
return 0;
}
/*!
an auxiliary function for calculating the Sine
(you don't have to call this function)
it returns Sin(x) where 'x' is from <0, PI/2>
we're calculating the Sin with using Taylor series in zero or PI/2
(depending on which point of these two points is nearer to the 'x')
Taylor series:
sin(x) = sin(a) + cos(a)*(x-a)/(1!)
- sin(a)*((x-a)^2)/(2!) - cos(a)*((x-a)^3)/(3!)
+ sin(a)*((x-a)^4)/(4!) + ...
when a=0 it'll be:
sin(x) = (x)/(1!) - (x^3)/(3!) + (x^5)/(5!) - (x^7)/(7!) + (x^9)/(9!) ...
and when a=PI/2:
sin(x) = 1 - ((x-PI/2)^2)/(2!) + ((x-PI/2)^4)/(4!) - ((x-PI/2)^6)/(6!) ...
*/
template<class ValueType>
ValueType Sin0pi05(const ValueType & x)
{
ValueType result;
ValueType numerator, denominator;
ValueType d_numerator, d_denominator;
ValueType one, temp, old_result;
// temp = pi/4
temp.Set05Pi();
temp.exponent.SubOne();
one.SetOne();
if( x < temp )
{
// we're using the Taylor series with a=0
result = x;
numerator = x;
denominator = one;
// d_numerator = x^2
d_numerator = x;
d_numerator.Mul(x);
d_denominator = 2;
}
else
{
// we're using the Taylor series with a=PI/2
result = one;
numerator = one;
denominator = one;
// d_numerator = (x-pi/2)^2
ValueType pi05;
pi05.Set05Pi();
temp = x;
temp.Sub( pi05 );
d_numerator = temp;
d_numerator.Mul( temp );
d_denominator = one;
}
uint c = 0;
bool addition = false;
old_result = result;
for(uint i=1 ; i<=TTMATH_ARITHMETIC_MAX_LOOP ; ++i)
{
// we're starting from a second part of the formula
c += numerator. Mul( d_numerator );
c += denominator. Mul( d_denominator );
c += d_denominator.Add( one );
c += denominator. Mul( d_denominator );
c += d_denominator.Add( one );
temp = numerator;
c += temp.Div(denominator);
if( c )
// Sin is from <-1,1> and cannot make an overflow
// but the carry can be from the Taylor series
// (then we only break our calculations)
break;
if( addition )
result.Add( temp );
else
result.Sub( temp );
addition = !addition;
// we're testing whether the result has changed after adding
// the next part of the Taylor formula, if not we end the loop
// (it means 'x' is zero or 'x' is PI/2 or this part of the formula
// is too small)
if( result == old_result )
break;
old_result = result;
}
return result;
}
} // namespace auxiliaryfunctions
/*!
this function calculates the Sine
*/
template<class ValueType>
ValueType Sin(ValueType x, ErrorCode * err = 0)
{
using namespace auxiliaryfunctions;
ValueType one, result;
bool change_sign;
if( x.IsNan() )
{
if( err )
*err = err_improper_argument;
return result; // NaN is set by default
}
if( err )
*err = err_ok;
if( PrepareSin( x, change_sign ) )
{
// x is too big, we cannnot reduce the 2*PI period
// prior to version 0.8.5 the result was zero
// result has NaN flag set by default
if( err )
*err = err_overflow; // maybe another error code? err_improper_argument?
return result; // NaN is set by default
}
result = Sin0pi05( x );
one.SetOne();
// after calculations there can be small distortions in the result
if( result > one )
result = one;
else
if( result.IsSign() )
// we've calculated the sin from <0, pi/2> and the result
// should be positive
result.SetZero();
if( change_sign )
result.ChangeSign();
return result;
}
/*!
this function calulates the Cosine
we're using the formula cos(x) = sin(x + PI/2)
*/
template<class ValueType>
ValueType Cos(ValueType x, ErrorCode * err = 0)
{
if( x.IsNan() )
{
if( err )
*err = err_improper_argument;
return x; // NaN
}
ValueType pi05;
pi05.Set05Pi();
uint c = x.Add( pi05 );
if( c )
{
if( err )
*err = err_overflow;
return ValueType(); // result is undefined (NaN is set by default)
}
return Sin(x, err);
}
/*!
this function calulates the Tangent
we're using the formula tan(x) = sin(x) / cos(x)
it takes more time than calculating the Tan directly
from for example Taylor series but should be a bit preciser
because Tan receives its values from -infinity to +infinity
and when we calculate it from any series then we can make
a greater mistake than calculating 'sin/cos'
*/
template<class ValueType>
ValueType Tan(const ValueType & x, ErrorCode * err = 0)
{
ValueType result = Cos(x, err);
if( err && *err != err_ok )
return result;
if( result.IsZero() )
{
if( err )
*err = err_improper_argument;
result.SetNan();
return result;
}
return Sin(x, err) / result;
}
/*!
this function calulates the Tangent
look at the description of Tan(...)
(the abbreviation of Tangent can be 'tg' as well)
*/
template<class ValueType>
ValueType Tg(const ValueType & x, ErrorCode * err = 0)
{
return Tan(x, err);
}
/*!
this function calulates the Cotangent
we're using the formula tan(x) = cos(x) / sin(x)
(why do we make it in this way?
look at information in Tan() function)
*/
template<class ValueType>
ValueType Cot(const ValueType & x, ErrorCode * err = 0)
{
ValueType result = Sin(x, err);
if( err && *err != err_ok )
return result;
if( result.IsZero() )
{
if( err )
*err = err_improper_argument;
result.SetNan();
return result;
}
return Cos(x, err) / result;
}
/*!
this function calulates the Cotangent
look at the description of Cot(...)
(the abbreviation of Cotangent can be 'ctg' as well)
*/
template<class ValueType>
ValueType Ctg(const ValueType & x, ErrorCode * err = 0)
{
return Cot(x, err);
}
/*
*
* inverse trigonometric functions
*
*
*/
namespace auxiliaryfunctions
{
/*!
an auxiliary function for calculating the Arc Sine
we're calculating asin from the following formula:
asin(x) = x + (1*x^3)/(2*3) + (1*3*x^5)/(2*4*5) + (1*3*5*x^7)/(2*4*6*7) + ...
where abs(x) <= 1
we're using this formula when x is from <0, 1/2>
*/
template<class ValueType>
ValueType ASin_0(const ValueType & x)
{
ValueType nominator, denominator, nominator_add, nominator_x, denominator_add, denominator_x;
ValueType two, result(x), x2(x);
ValueType nominator_temp, denominator_temp, old_result = result;
uint c = 0;
x2.Mul(x);
two = 2;
nominator.SetOne();
denominator = two;
nominator_add = nominator;
denominator_add = denominator;
nominator_x = x;
denominator_x = 3;
for(uint i=1 ; i<=TTMATH_ARITHMETIC_MAX_LOOP ; ++i)
{
c += nominator_x.Mul(x2);
nominator_temp = nominator_x;
c += nominator_temp.Mul(nominator);
denominator_temp = denominator;
c += denominator_temp.Mul(denominator_x);
c += nominator_temp.Div(denominator_temp);
// if there is a carry somewhere we only break the calculating
// the result should be ok -- it's from <-pi/2, pi/2>
if( c )
break;
result.Add(nominator_temp);
if( result == old_result )
// there's no sense to calculate more
break;
old_result = result;
c += nominator_add.Add(two);
c += denominator_add.Add(two);
c += nominator.Mul(nominator_add);
c += denominator.Mul(denominator_add);
c += denominator_x.Add(two);
}
return result;
}
/*!
an auxiliary function for calculating the Arc Sine
we're calculating asin from the following formula:
asin(x) = pi/2 - sqrt(2)*sqrt(1-x) * asin_temp
asin_temp = 1 + (1*(1-x))/((2*3)*(2)) + (1*3*(1-x)^2)/((2*4*5)*(4)) + (1*3*5*(1-x)^3)/((2*4*6*7)*(8)) + ...
where abs(x) <= 1
we're using this formula when x is from (1/2, 1>
*/
template<class ValueType>
ValueType ASin_1(const ValueType & x)
{
ValueType nominator, denominator, nominator_add, nominator_x, nominator_x_add, denominator_add, denominator_x;
ValueType denominator2;
ValueType one, two, result;
ValueType nominator_temp, denominator_temp, old_result;
uint c = 0;
two = 2;
one.SetOne();
nominator = one;
result = one;
old_result = result;
denominator = two;
nominator_add = nominator;
denominator_add = denominator;
nominator_x = one;
nominator_x.Sub(x);
nominator_x_add = nominator_x;
denominator_x = 3;
denominator2 = two;
for(uint i=1 ; i<=TTMATH_ARITHMETIC_MAX_LOOP ; ++i)
{
nominator_temp = nominator_x;
c += nominator_temp.Mul(nominator);
denominator_temp = denominator;
c += denominator_temp.Mul(denominator_x);
c += denominator_temp.Mul(denominator2);
c += nominator_temp.Div(denominator_temp);
// if there is a carry somewhere we only break the calculating
// the result should be ok -- it's from <-pi/2, pi/2>
if( c )
break;
result.Add(nominator_temp);
if( result == old_result )
// there's no sense to calculate more
break;
old_result = result;
c += nominator_x.Mul(nominator_x_add);
c += nominator_add.Add(two);
c += denominator_add.Add(two);
c += nominator.Mul(nominator_add);
c += denominator.Mul(denominator_add);
c += denominator_x.Add(two);
c += denominator2.Mul(two);
}
nominator_x_add.exponent.AddOne(); // *2
one.exponent.SubOne(); // =0.5
nominator_x_add.Pow(one); // =sqrt(nominator_x_add)
result.Mul(nominator_x_add);
one.Set05Pi();
one.Sub(result);
return one;
}
} // namespace auxiliaryfunctions
/*!
this function calculates the Arc Sine
x is from <-1,1>
*/
template<class ValueType>
ValueType ASin(ValueType x, ErrorCode * err = 0)
{
using namespace auxiliaryfunctions;
ValueType result, one;
one.SetOne();
bool change_sign = false;
if( x.IsNan() )
{
if( err )
*err = err_improper_argument;
return result; // NaN is set by default
}
if( x.GreaterWithoutSignThan(one) )
{
if( err )
*err = err_improper_argument;
return result; // NaN is set by default
}
if( x.IsSign() )
{
change_sign = true;
x.Abs();
}
one.exponent.SubOne(); // =0.5
// asin(-x) = -asin(x)
if( x.GreaterWithoutSignThan(one) )
result = ASin_1(x);
else
result = ASin_0(x);
if( change_sign )
result.ChangeSign();
if( err )
*err = err_ok;
return result;
}
/*!
this function calculates the Arc Cosine
we're using the formula:
acos(x) = pi/2 - asin(x)
*/
template<class ValueType>
ValueType ACos(const ValueType & x, ErrorCode * err = 0)
{
ValueType temp;
temp.Set05Pi();
temp.Sub(ASin(x,err));
return temp;
}
namespace auxiliaryfunctions
{
/*!
an auxiliary function for calculating the Arc Tangent
arc tan (x) where x is in <0; 0.5)
(x can be in (-0.5 ; 0.5) too)
we're using the Taylor series expanded in zero:
atan(x) = x - (x^3)/3 + (x^5)/5 - (x^7)/7 + ...
*/
template<class ValueType>
ValueType ATan0(const ValueType & x)
{
ValueType nominator, denominator, nominator_add, denominator_add, temp;
ValueType result, old_result;
bool adding = false;
uint c = 0;
result = x;
old_result = result;
nominator = x;
nominator_add = x;
nominator_add.Mul(x);
denominator.SetOne();
denominator_add = 2;
for(uint i=1 ; i<=TTMATH_ARITHMETIC_MAX_LOOP ; ++i)
{
c += nominator.Mul(nominator_add);
c += denominator.Add(denominator_add);
temp = nominator;
c += temp.Div(denominator);
if( c )
// the result should be ok
break;
if( adding )
result.Add(temp);
else
result.Sub(temp);
if( result == old_result )
// there's no sense to calculate more
break;
old_result = result;
adding = !adding;
}
return result;
}
/*!
an auxiliary function for calculating the Arc Tangent
where x is in <0 ; 1>
*/
template<class ValueType>
ValueType ATan01(const ValueType & x)
{
ValueType half;
half.Set05();
/*
it would be better if we chose about sqrt(2)-1=0.41... instead of 0.5 here
because as you can see below:
when x = sqrt(2)-1
abs(x) = abs( (x-1)/(1+x) )
so when we're calculating values around x
then they will be better converged to each other
for example if we have x=0.4999 then during calculating ATan0(0.4999)
we have to make about 141 iterations but when we have x=0.5
then during calculating ATan0( (x-1)/(1+x) ) we have to make
only about 89 iterations (both for Big<3,9>)
in the future this 0.5 can be changed
*/
if( x.SmallerWithoutSignThan(half) )
return ATan0(x);
/*
x>=0.5 and x<=1
(x can be even smaller than 0.5)
y = atac(x)
x = tan(y)
tan(y-b) = (tan(y)-tab(b)) / (1+tan(y)*tan(b))
y-b = atan( (tan(y)-tab(b)) / (1+tan(y)*tan(b)) )
y = b + atan( (x-tab(b)) / (1+x*tan(b)) )
let b = pi/4
tan(b) = tan(pi/4) = 1
y = pi/4 + atan( (x-1)/(1+x) )
so
atac(x) = pi/4 + atan( (x-1)/(1+x) )
when x->1 (x converges to 1) the (x-1)/(1+x) -> 0
and we can use ATan0() function here
*/
ValueType n(x),d(x),one,result;
one.SetOne();
n.Sub(one);
d.Add(one);
n.Div(d);
result = ATan0(n);
n.Set05Pi();
n.exponent.SubOne(); // =pi/4
result.Add(n);
return result;
}
/*!
an auxiliary function for calculating the Arc Tangent
where x > 1
we're using the formula:
atan(x) = pi/2 - atan(1/x) for x>0
*/
template<class ValueType>
ValueType ATanGreaterThanPlusOne(const ValueType & x)
{
ValueType temp, atan;
temp.SetOne();
if( temp.Div(x) )
{
// if there was a carry here that means x is very big
// and atan(1/x) fast converged to 0
atan.SetZero();
}
else
atan = ATan01(temp);
temp.Set05Pi();
temp.Sub(atan);
return temp;
}
} // namespace auxiliaryfunctions
/*!
this function calculates the Arc Tangent
*/
template<class ValueType>
ValueType ATan(ValueType x)
{
using namespace auxiliaryfunctions;
ValueType one, result;
one.SetOne();
bool change_sign = false;
if( x.IsNan() )
return result; // NaN is set by default
// if x is negative we're using the formula:
// atan(-x) = -atan(x)
if( x.IsSign() )
{
change_sign = true;
x.Abs();
}
if( x.GreaterWithoutSignThan(one) )
result = ATanGreaterThanPlusOne(x);
else
result = ATan01(x);
if( change_sign )
result.ChangeSign();
return result;
}
/*!
this function calculates the Arc Tangent
look at the description of ATan(...)
(the abbreviation of Arc Tangent can be 'atg' as well)
*/
template<class ValueType>
ValueType ATg(const ValueType & x)
{
return ATan(x);
}
/*!
this function calculates the Arc Cotangent
we're using the formula:
actan(x) = pi/2 - atan(x)
*/
template<class ValueType>
ValueType ACot(const ValueType & x)
{
ValueType result;
result.Set05Pi();
result.Sub(ATan(x));
return result;
}
/*!
this function calculates the Arc Cotangent
look at the description of ACot(...)
(the abbreviation of Arc Cotangent can be 'actg' as well)
*/
template<class ValueType>
ValueType ACtg(const ValueType & x)
{
return ACot(x);
}
/*
*
* hyperbolic functions
*
*
*/
/*!
this function calculates the Hyperbolic Sine
we're using the formula sinh(x)= ( e^x - e^(-x) ) / 2
*/
template<class ValueType>
ValueType Sinh(const ValueType & x, ErrorCode * err = 0)
{
if( x.IsNan() )
{
if( err )
*err = err_improper_argument;
return x; // NaN
}
ValueType ex, emx;
uint c = 0;
c += ex.Exp(x);
c += emx.Exp(-x);
c += ex.Sub(emx);
c += ex.exponent.SubOne();
if( err )
*err = c ? err_overflow : err_ok;
return ex;
}
/*!
this function calculates the Hyperbolic Cosine
we're using the formula cosh(x)= ( e^x + e^(-x) ) / 2
*/
template<class ValueType>
ValueType Cosh(const ValueType & x, ErrorCode * err = 0)
{
if( x.IsNan() )
{
if( err )
*err = err_improper_argument;
return x; // NaN
}
ValueType ex, emx;
uint c = 0;
c += ex.Exp(x);
c += emx.Exp(-x);
c += ex.Add(emx);
c += ex.exponent.SubOne();
if( err )
*err = c ? err_overflow : err_ok;
return ex;
}
/*!
this function calculates the Hyperbolic Tangent
we're using the formula tanh(x)= ( e^x - e^(-x) ) / ( e^x + e^(-x) )
*/
template<class ValueType>
ValueType Tanh(const ValueType & x, ErrorCode * err = 0)
{
if( x.IsNan() )
{
if( err )
*err = err_improper_argument;
return x; // NaN
}
ValueType ex, emx, nominator, denominator;
uint c = 0;
c += ex.Exp(x);
c += emx.Exp(-x);
nominator = ex;
c += nominator.Sub(emx);
denominator = ex;
c += denominator.Add(emx);
c += nominator.Div(denominator);
if( err )
*err = c ? err_overflow : err_ok;
return nominator;
}
/*!
this function calculates the Hyperbolic Tangent
look at the description of Tanh(...)
(the abbreviation of Hyperbolic Tangent can be 'tgh' as well)
*/
template<class ValueType>
ValueType Tgh(const ValueType & x, ErrorCode * err = 0)
{
return Tanh(x, err);
}
/*!
this function calculates the Hyperbolic Cotangent
we're using the formula coth(x)= ( e^x + e^(-x) ) / ( e^x - e^(-x) )
*/
template<class ValueType>
ValueType Coth(const ValueType & x, ErrorCode * err = 0)
{
if( x.IsNan() )
{
if( err )
*err = err_improper_argument;
return x; // NaN
}
if( x.IsZero() )
{
if( err )
*err = err_improper_argument;
return ValueType(); // NaN is set by default
}
ValueType ex, emx, nominator, denominator;
uint c = 0;
c += ex.Exp(x);
c += emx.Exp(-x);
nominator = ex;
c += nominator.Add(emx);
denominator = ex;
c += denominator.Sub(emx);
c += nominator.Div(denominator);
if( err )
*err = c ? err_overflow : err_ok;
return nominator;
}
/*!
this function calculates the Hyperbolic Cotangent
look at the description of Coth(...)
(the abbreviation of Hyperbolic Cotangent can be 'ctgh' as well)
*/
template<class ValueType>
ValueType Ctgh(const ValueType & x, ErrorCode * err = 0)
{
return Coth(x, err);
}
/*
*
* inverse hyperbolic functions
*
*
*/
/*!
inverse hyperbolic sine
asinh(x) = ln( x + sqrt(x^2 + 1) )
*/
template<class ValueType>
ValueType ASinh(const ValueType & x, ErrorCode * err = 0)
{
if( x.IsNan() )
{
if( err )
*err = err_improper_argument;
return x; // NaN
}
ValueType xx(x), one, result;
uint c = 0;
one.SetOne();
c += xx.Mul(x);
c += xx.Add(one);
one.exponent.SubOne(); // one=0.5
// xx is >= 1
c += xx.PowFrac(one); // xx=sqrt(xx)
c += xx.Add(x);
c += result.Ln(xx); // xx > 0
// here can only be a carry
if( err )
*err = c ? err_overflow : err_ok;
return result;
}
/*!
inverse hyperbolic cosine
acosh(x) = ln( x + sqrt(x^2 - 1) ) x in <1, infinity)
*/
template<class ValueType>
ValueType ACosh(const ValueType & x, ErrorCode * err = 0)
{
if( x.IsNan() )
{
if( err )
*err = err_improper_argument;
return x; // NaN
}
ValueType xx(x), one, result;
uint c = 0;
one.SetOne();
if( x < one )
{
if( err )
*err = err_improper_argument;
return result; // NaN is set by default
}
c += xx.Mul(x);
c += xx.Sub(one);
// xx is >= 0
// we can't call a PowFrac when the 'x' is zero
// if x is 0 the sqrt(0) is 0
if( !xx.IsZero() )
{
one.exponent.SubOne(); // one=0.5
c += xx.PowFrac(one); // xx=sqrt(xx)
}
c += xx.Add(x);
c += result.Ln(xx); // xx >= 1
// here can only be a carry
if( err )
*err = c ? err_overflow : err_ok;
return result;
}
/*!
inverse hyperbolic tangent
atanh(x) = 0.5 * ln( (1+x) / (1-x) ) x in (-1, 1)
*/
template<class ValueType>
ValueType ATanh(const ValueType & x, ErrorCode * err = 0)
{
if( x.IsNan() )
{
if( err )
*err = err_improper_argument;
return x; // NaN
}
ValueType nominator(x), denominator, one, result;
uint c = 0;
one.SetOne();
if( !x.SmallerWithoutSignThan(one) )
{
if( err )
*err = err_improper_argument;
return result; // NaN is set by default
}
c += nominator.Add(one);
denominator = one;
c += denominator.Sub(x);
c += nominator.Div(denominator);
c += result.Ln(nominator);
c += result.exponent.SubOne();
// here can only be a carry
if( err )
*err = c ? err_overflow : err_ok;
return result;
}
/*!
inverse hyperbolic tantent
*/
template<class ValueType>
ValueType ATgh(const ValueType & x, ErrorCode * err = 0)
{
return ATanh(x, err);
}
/*!
inverse hyperbolic cotangent
acoth(x) = 0.5 * ln( (x+1) / (x-1) ) x in (-infinity, -1) or (1, infinity)
*/
template<class ValueType>
ValueType ACoth(const ValueType & x, ErrorCode * err = 0)
{
if( x.IsNan() )
{
if( err )
*err = err_improper_argument;
return x; // NaN
}
ValueType nominator(x), denominator(x), one, result;
uint c = 0;
one.SetOne();
if( !x.GreaterWithoutSignThan(one) )
{
if( err )
*err = err_improper_argument;
return result; // NaN is set by default
}
c += nominator.Add(one);
c += denominator.Sub(one);
c += nominator.Div(denominator);
c += result.Ln(nominator);
c += result.exponent.SubOne();
// here can only be a carry
if( err )
*err = c ? err_overflow : err_ok;
return result;
}
/*!
inverse hyperbolic cotantent
*/
template<class ValueType>
ValueType ACtgh(const ValueType & x, ErrorCode * err = 0)
{
return ACoth(x, err);
}
/*
*
* functions for converting between degrees, radians and gradians
*
*
*/
/*!
this function converts degrees to radians
it returns: x * pi / 180
*/
template<class ValueType>
ValueType DegToRad(const ValueType & x, ErrorCode * err = 0)
{
ValueType result, temp;
uint c = 0;
if( x.IsNan() )
{
if( err )
*err = err_improper_argument;
return result; // NaN is set by default
}
result = x;
// it is better to make division first and then multiplication
// the result is more accurate especially when x is: 90,180,270 or 360
temp = 180;
c += result.Div(temp);
temp.SetPi();
c += result.Mul(temp);
if( err )
*err = c ? err_overflow : err_ok;
return result;
}
/*!
this function converts radians to degrees
it returns: x * 180 / pi
*/
template<class ValueType>
ValueType RadToDeg(const ValueType & x, ErrorCode * err = 0)
{
ValueType result, delimiter;
uint c = 0;
if( x.IsNan() )
{
if( err )
*err = err_improper_argument;
return result; // NaN is set by default
}
result = 180;
c += result.Mul(x);
delimiter.SetPi();
c += result.Div(delimiter);
if( err )
*err = c ? err_overflow : err_ok;
return result;
}
/*!
this function converts degrees in the long format into one value
long format: (degrees, minutes, seconds)
minutes and seconds must be greater than or equal zero
result:
if d>=0 : result= d + ((s/60)+m)/60
if d<0 : result= d - ((s/60)+m)/60
((s/60)+m)/60 = (s+60*m)/3600 (second version is faster because
there's only one division)
for example:
DegToDeg(10, 30, 0) = 10.5
DegToDeg(10, 24, 35.6)=10.4098(8)
*/
template<class ValueType>
ValueType DegToDeg( const ValueType & d, const ValueType & m, const ValueType & s,
ErrorCode * err = 0)
{
ValueType delimiter, multipler;
uint c = 0;
if( d.IsNan() || m.IsNan() || s.IsNan() || m.IsSign() || s.IsSign() )
{
if( err )
*err = err_improper_argument;
return delimiter; // NaN is set by default
}
multipler = 60;
delimiter = 3600;
c += multipler.Mul(m);
c += multipler.Add(s);
c += multipler.Div(delimiter);
if( d.IsSign() )
multipler.ChangeSign();
c += multipler.Add(d);
if( err )
*err = c ? err_overflow : err_ok;
return multipler;
}
/*!
this function converts degrees in the long format to radians
*/
template<class ValueType>
ValueType DegToRad( const ValueType & d, const ValueType & m, const ValueType & s,
ErrorCode * err = 0)
{
ValueType temp_deg = DegToDeg(d,m,s,err);
if( err && *err!=err_ok )
return temp_deg;
return DegToRad(temp_deg, err);
}
/*!
this function converts gradians to radians
it returns: x * pi / 200
*/
template<class ValueType>
ValueType GradToRad(const ValueType & x, ErrorCode * err = 0)
{
ValueType result, temp;
uint c = 0;
if( x.IsNan() )
{
if( err )
*err = err_improper_argument;
return result; // NaN is set by default
}
result = x;
// it is better to make division first and then multiplication
// the result is more accurate especially when x is: 100,200,300 or 400
temp = 200;
c += result.Div(temp);
temp.SetPi();
c += result.Mul(temp);
if( err )
*err = c ? err_overflow : err_ok;
return result;
}
/*!
this function converts radians to gradians
it returns: x * 200 / pi
*/
template<class ValueType>
ValueType RadToGrad(const ValueType & x, ErrorCode * err = 0)
{
ValueType result, delimiter;
uint c = 0;
if( x.IsNan() )
{
if( err )
*err = err_improper_argument;
return result; // NaN is set by default
}
result = 200;
c += result.Mul(x);
delimiter.SetPi();
c += result.Div(delimiter);
if( err )
*err = c ? err_overflow : err_ok;
return result;
}
/*!
this function converts degrees to gradians
it returns: x * 200 / 180
*/
template<class ValueType>
ValueType DegToGrad(const ValueType & x, ErrorCode * err = 0)
{
ValueType result, temp;
uint c = 0;
if( x.IsNan() )
{
if( err )
*err = err_improper_argument;
return result; // NaN is set by default
}
result = x;
temp = 200;
c += result.Mul(temp);
temp = 180;
c += result.Div(temp);
if( err )
*err = c ? err_overflow : err_ok;
return result;
}
/*!
this function converts degrees in the long format to gradians
*/
template<class ValueType>
ValueType DegToGrad( const ValueType & d, const ValueType & m, const ValueType & s,
ErrorCode * err = 0)
{
ValueType temp_deg = DegToDeg(d,m,s,err);
if( err && *err!=err_ok )
return temp_deg;
return DegToGrad(temp_deg, err);
}
/*!
this function converts degrees to gradians
it returns: x * 180 / 200
*/
template<class ValueType>
ValueType GradToDeg(const ValueType & x, ErrorCode * err = 0)
{
ValueType result, temp;
uint c = 0;
if( x.IsNan() )
{
if( err )
*err = err_improper_argument;
return result; // NaN is set by default
}
result = x;
temp = 180;
c += result.Mul(temp);
temp = 200;
c += result.Div(temp);
if( err )
*err = c ? err_overflow : err_ok;
return result;
}
/*
*
* another functions
*
*
*/
/*!
this function calculates the square root
Sqrt(9) = 3
*/
template<class ValueType>
ValueType Sqrt(ValueType x, ErrorCode * err = 0)
{
if( x.IsNan() || x.IsSign() )
{
if( err )
*err = err_improper_argument;
return ValueType(); // NaN is set by default
}
if( x.IsZero() )
{
// Sqrt(0) = 0
if( err )
*err = err_ok;
return x;
}
ValueType pow;
pow.Set05();
// PowFrac can return only a carry because x is greater than zero
uint c = x.PowFrac(pow);
if( err )
*err = c ? err_overflow : err_ok;
return x;
}
namespace auxiliaryfunctions
{
template<class ValueType>
bool RootCheckIndexSign(ValueType & x, const ValueType & index, ErrorCode * err)
{
if( index.IsSign() )
{
// index cannot be negative
if( err )
*err = err_improper_argument;
x.SetNan();
return true;
}
return false;
}
template<class ValueType>
bool RootCheckIndexZero(ValueType & x, const ValueType & index, ErrorCode * err)
{
if( index.IsZero() )
{
if( x.IsZero() )
{
// there isn't root(0;0) - we assume it's not defined
if( err )
*err = err_improper_argument;
x.SetNan();
return true;
}
// root(x;0) is 1 (if x!=0)
x.SetOne();
if( err )
*err = err_ok;
return true;
}
return false;
}
template<class ValueType>
bool RootCheckIndexOne(const ValueType & index, ErrorCode * err)
{
ValueType one;
one.SetOne();
if( index == one )
{
//root(x;1) is x
// we do it because if we used the PowFrac function
// we would lose the precision
if( err )
*err = err_ok;
return true;
}
return false;
}
template<class ValueType>
bool RootCheckIndexFrac(ValueType & x, const ValueType & index, ErrorCode * err)
{
if( !index.IsInteger() )
{
// index must be integer
if( err )
*err = err_improper_argument;
x.SetNan();
return true;
}
return false;
}
template<class ValueType>
bool RootCheckXZero(ValueType & x, ErrorCode * err)
{
if( x.IsZero() )
{
// root(0;index) is zero (if index!=0)
// RootCheckIndexZero() must be called beforehand
x.SetZero();
if( err )
*err = err_ok;
return true;
}
return false;
}
template<class ValueType>
bool RootCheckIndex(ValueType & x, const ValueType & index, ErrorCode * err, bool * change_sign)
{
*change_sign = false;
if( index.Mod2() )
{
// index is odd (1,3,5...)
if( x.IsSign() )
{
*change_sign = true;
x.Abs();
}
}
else
{
// index is even
// x cannot be negative
if( x.IsSign() )
{
if( err )
*err = err_improper_argument;
x.SetNan();
return true;
}
}
return false;
}
}
/*!
indexth Root of x
index must be integer and not negative <0;1;2;3....)
if index==0 the result is one
if x==0 the result is zero and we assume root(0;0) is not defined
if index is even (2;4;6...) the result is x^(1/index) and x>0
if index is odd (1;2;3;...) the result is either
-(abs(x)^(1/index)) if x<0 or
x^(1/index)) if x>0
(for index==1 the result is equal x)
*/
template<class ValueType>
ValueType Root(ValueType x, const ValueType & index, ErrorCode * err = 0)
{
using namespace auxiliaryfunctions;
if( x.IsNan() || index.IsNan() )
{
if( err )
*err = err_improper_argument;
return ValueType(); // NaN is set by default
}
if( RootCheckIndexSign(x, index, err) ) return x;
if( RootCheckIndexZero(x, index, err) ) return x;
if( RootCheckIndexOne ( index, err) ) return x;
if( RootCheckIndexFrac(x, index, err) ) return x;
if( RootCheckXZero (x, err) ) return x;
// index integer and index!=0
// x!=0
uint c = 0;
bool change_sign;
if( RootCheckIndex(x, index, err, &change_sign ) ) return x;
ValueType newindex;
newindex.SetOne();
c += newindex.Div(index);
c += x.PowFrac(newindex); // here can only be a carry
if( change_sign )
{
// the value of x should be different from zero
// (x is actually tested by RootCheckXZero)
TTMATH_ASSERT( x.IsZero() == false )
x.SetSign();
}
if( err )
*err = c ? err_overflow : err_ok;
return x;
}
/*!
absolute value of x
e.g. -2 = 2
2 = 2
*/
template<class ValueType>
ValueType Abs(const ValueType & x)
{
ValueType result( x );
result.Abs();
return result;
}
/*!
it returns the sign of the value
e.g. -2 = -1
0 = 0
10 = 1
*/
template<class ValueType>
ValueType Sgn(ValueType x)
{
x.Sgn();
return x;
}
/*!
the remainder from a division
e.g.
mod( 12.6 ; 3) = 0.6 because 12.6 = 3*4 + 0.6
mod(-12.6 ; 3) = -0.6 bacause -12.6 = 3*(-4) + (-0.6)
mod( 12.6 ; -3) = 0.6
mod(-12.6 ; -3) = -0.6
*/
template<class ValueType>
ValueType Mod(ValueType a, const ValueType & b, ErrorCode * err = 0)
{
if( a.IsNan() || b.IsNan() )
{
if( err )
*err = err_improper_argument;
return ValueType(); // NaN is set by default
}
uint c = a.Mod(b);
if( err )
*err = c ? err_overflow : err_ok;
return a;
}
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
/*!
this is for convenience for the user
he can only use '#include <ttmath/ttmath.h>' even if he uses the parser
*/
#include "ttmathparser.h"
#ifdef _MSC_VER
#pragma warning( default: 4127 )
//warning C4127: conditional expression is constant
#endif
#endif
Reference in New Issue View Git Blame Copy Permalink