531 lines
10 KiB
C
531 lines
10 KiB
C
|
/*
|
||
|
* This file is part of TTMath Mathematical Library
|
||
|
* and is distributed under the (new) BSD licence.
|
||
|
* Author: Tomasz Sowa <t.sowa@slimaczek.pl>
|
||
|
*/
|
||
|
|
||
|
/*
|
||
|
* Copyright (c) 2006-2007, 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.
|
||
|
*/
|
||
|
|
||
|
#include "ttmathbig.h"
|
||
|
|
||
|
#include <string>
|
||
|
|
||
|
|
||
|
|
||
|
namespace ttmath
|
||
|
{
|
||
|
|
||
|
/*!
|
||
|
the factorial from given 'x'
|
||
|
e.g.
|
||
|
Factorial(4) = 4! = 1*2*3*4
|
||
|
*/
|
||
|
template<int exp,int man>
|
||
|
Big<exp, man> Factorial(const Big<exp, man> & x, ErrorCode * err = 0, const volatile StopCalculating * stop = 0)
|
||
|
{
|
||
|
Big<exp, man> result;
|
||
|
|
||
|
result.SetOne();
|
||
|
|
||
|
if( x.IsSign() )
|
||
|
{
|
||
|
if( err )
|
||
|
*err = err_improper_argument;
|
||
|
|
||
|
return result;
|
||
|
}
|
||
|
|
||
|
if( !x.exponent.IsSign() && !x.exponent.IsZero() )
|
||
|
{
|
||
|
// when x>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 our mantissa)
|
||
|
if( err )
|
||
|
*err = err_overflow;
|
||
|
|
||
|
return result;
|
||
|
}
|
||
|
|
||
|
Big<exp, man> multipler;
|
||
|
Big<exp, man> one;
|
||
|
uint carry = 0;
|
||
|
|
||
|
one = result; // =1
|
||
|
multipler = result; // =1
|
||
|
|
||
|
while( !carry && multipler < x )
|
||
|
{
|
||
|
if( stop && stop->WasStopSignal() )
|
||
|
{
|
||
|
if( err )
|
||
|
*err = err_interrupt;
|
||
|
|
||
|
return result;
|
||
|
}
|
||
|
|
||
|
carry += multipler.Add(one);
|
||
|
carry += result.Mul(multipler);
|
||
|
}
|
||
|
|
||
|
|
||
|
if( err )
|
||
|
{
|
||
|
if( carry )
|
||
|
*err = err_overflow;
|
||
|
else
|
||
|
*err = err_ok;
|
||
|
}
|
||
|
|
||
|
return result;
|
||
|
}
|
||
|
|
||
|
|
||
|
/*!
|
||
|
absolute value of x
|
||
|
e.g. -2 = 2
|
||
|
2 = 2
|
||
|
*/
|
||
|
template<int exp,int man>
|
||
|
Big<exp, man> Abs(const Big<exp, man> & x)
|
||
|
{
|
||
|
Big<exp, man> result( x );
|
||
|
result.Abs();
|
||
|
|
||
|
return result;
|
||
|
}
|
||
|
|
||
|
|
||
|
/*!
|
||
|
this method skips the fraction from x
|
||
|
e.g 2.2 = 2
|
||
|
2.7 = 2
|
||
|
-2.2 = 2
|
||
|
-2.7 = 2
|
||
|
*/
|
||
|
template<int exp,int man>
|
||
|
Big<exp, man> SkipFraction(const Big<exp, man> & x)
|
||
|
{
|
||
|
Big<exp, man> result( x );
|
||
|
result.SkipFraction();
|
||
|
|
||
|
return result;
|
||
|
}
|
||
|
|
||
|
|
||
|
/*!
|
||
|
this method rounds to the nearest integer value
|
||
|
e.g 2.2 = 2
|
||
|
2.7 = 3
|
||
|
-2.2 = -2
|
||
|
-2.7 = -3
|
||
|
*/
|
||
|
template<int exp,int man>
|
||
|
Big<exp, man> Round(const Big<exp, man> & x)
|
||
|
{
|
||
|
Big<exp, man> result( x );
|
||
|
result.Round();
|
||
|
|
||
|
return result;
|
||
|
}
|
||
|
|
||
|
|
||
|
/*!
|
||
|
this method calculates the natural logarithm (logarithm with the base 'e')
|
||
|
*/
|
||
|
template<int exp,int man>
|
||
|
Big<exp, man> Ln(const Big<exp, man> & x, ErrorCode * err = 0)
|
||
|
{
|
||
|
Big<exp, man> 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 method calculates the logarithm
|
||
|
*/
|
||
|
template<int exp,int man>
|
||
|
Big<exp, man> Log(const Big<exp, man> & base, const Big<exp, man> & x, ErrorCode * err = 0)
|
||
|
{
|
||
|
Big<exp, man> result;
|
||
|
|
||
|
uint state = result.Log(base,x);
|
||
|
|
||
|
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 method calculates the expression e^x
|
||
|
*/
|
||
|
template<int exp,int man>
|
||
|
Big<exp, man> Exp(const Big<exp, man> & x, ErrorCode * err = 0)
|
||
|
{
|
||
|
Big<exp, man> result;
|
||
|
|
||
|
uint state = result.Exp(x);
|
||
|
|
||
|
if( err )
|
||
|
if( state!=0 )
|
||
|
*err = err_overflow;
|
||
|
else
|
||
|
*err = err_ok;
|
||
|
|
||
|
return result;
|
||
|
}
|
||
|
|
||
|
|
||
|
/*!
|
||
|
*
|
||
|
* trigonometric functions
|
||
|
*
|
||
|
*/
|
||
|
|
||
|
|
||
|
/*!
|
||
|
an auxiliary function for calculating the Sin
|
||
|
(you don't have to call this function)
|
||
|
*/
|
||
|
template<int exp,int man>
|
||
|
void PrepareSin(Big<exp,man> & x, bool & change_sign)
|
||
|
{
|
||
|
Big<exp,man> 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 > temp )
|
||
|
{
|
||
|
x.Div( temp );
|
||
|
x.RemainFraction();
|
||
|
x.Mul( temp );
|
||
|
}
|
||
|
|
||
|
// 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;
|
||
|
}
|
||
|
}
|
||
|
|
||
|
|
||
|
/*!
|
||
|
an auxiliary function for calculating the Sin
|
||
|
(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<int exp,int man>
|
||
|
Big<exp,man> Sin0pi05(const Big<exp,man> & x)
|
||
|
{
|
||
|
Big<exp,man> result;
|
||
|
Big<exp,man> numerator, denominator;
|
||
|
Big<exp,man> d_numerator, d_denominator;
|
||
|
Big<exp,man> 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
|
||
|
Big<exp,man> pi05;
|
||
|
pi05.Set05Pi();
|
||
|
|
||
|
temp = x;
|
||
|
temp.Sub( pi05 );
|
||
|
d_numerator = temp;
|
||
|
d_numerator.Mul( temp );
|
||
|
|
||
|
d_denominator = one;
|
||
|
}
|
||
|
|
||
|
int c = 0;
|
||
|
bool addition = false;
|
||
|
|
||
|
old_result = result;
|
||
|
for(int i=1 ; i<5000 ; ++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 breaks 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;
|
||
|
}
|
||
|
|
||
|
|
||
|
/*!
|
||
|
this function calulates the Sin
|
||
|
*/
|
||
|
template<int exp,int man>
|
||
|
Big<exp,man> Sin(Big<exp,man> x)
|
||
|
{
|
||
|
Big<exp,man> one;
|
||
|
bool change_sign;
|
||
|
|
||
|
PrepareSin( x, change_sign );
|
||
|
Big<exp,man> result = Sin0pi05( x );
|
||
|
|
||
|
one.SetOne();
|
||
|
|
||
|
if( result > one )
|
||
|
result = one;
|
||
|
else
|
||
|
if( result.IsSign() )
|
||
|
result.SetZero();
|
||
|
|
||
|
if( change_sign )
|
||
|
result.ChangeSign();
|
||
|
|
||
|
return result;
|
||
|
}
|
||
|
|
||
|
|
||
|
/*!
|
||
|
this function calulates the Cos
|
||
|
we're using the formula cos(x) = sin(x + PI/2)
|
||
|
*/
|
||
|
template<int exp,int man>
|
||
|
Big<exp,man> Cos(Big<exp,man> x)
|
||
|
{
|
||
|
Big<exp,man> pi05;
|
||
|
pi05.Set05Pi();
|
||
|
|
||
|
x.Add( pi05 );
|
||
|
|
||
|
return Sin(x);
|
||
|
}
|
||
|
|
||
|
|
||
|
/*!
|
||
|
this function calulates the Tan
|
||
|
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 precise
|
||
|
because Tan receives its values from -infinity to +infinity
|
||
|
and when we calculate it from any series then we can make
|
||
|
a small mistake than calculating 'sin/cos'
|
||
|
*/
|
||
|
template<int exp,int man>
|
||
|
Big<exp,man> Tan(const Big<exp,man> & x, ErrorCode * err = 0)
|
||
|
{
|
||
|
Big<exp,man> result = Cos(x);
|
||
|
|
||
|
if( result.IsZero() )
|
||
|
{
|
||
|
if( err )
|
||
|
*err = err_improper_argument;
|
||
|
|
||
|
return result;
|
||
|
}
|
||
|
|
||
|
if( err )
|
||
|
*err = err_ok;
|
||
|
|
||
|
return Sin(x) / result;
|
||
|
}
|
||
|
|
||
|
|
||
|
/*!
|
||
|
this function calulates the CTan
|
||
|
we're using the formula tan(x) = cos(x) / sin(x)
|
||
|
|
||
|
(why do we make it in this way?
|
||
|
look at the info in Tan() function)
|
||
|
*/
|
||
|
template<int exp,int man>
|
||
|
Big<exp,man> CTan(const Big<exp,man> & x, ErrorCode * err = 0)
|
||
|
{
|
||
|
Big<exp,man> result = Sin(x);
|
||
|
|
||
|
if( result.IsZero() )
|
||
|
{
|
||
|
if( err )
|
||
|
*err = err_improper_argument;
|
||
|
|
||
|
return result;
|
||
|
}
|
||
|
|
||
|
if( err )
|
||
|
*err = err_ok;
|
||
|
|
||
|
return Cos(x) / result;
|
||
|
}
|
||
|
|
||
|
|
||
|
} // namespace
|
||
|
|
||
|
|
||
|
|
||
|
|
||
|
|
||
|
#endif
|