fixed: constraints in asm operands for gcc

added: UInt::SetFromTable for 64bit code (now the support for 64bit
       platforms seems to be completed)
added: asin - arc sin, acos - arc cos


git-svn-id: svn://ttmath.org/publicrep/ttmath/trunk@16 e52654a7-88a9-db11-a3e9-0013d4bc506e
This commit is contained in:
Tomasz Sowa 2007-02-24 18:59:05 +00:00
parent d4b8bf2202
commit d04632ea74
7 changed files with 345 additions and 33 deletions

View File

@ -1,3 +1,7 @@
Version 0.7.0 (2007.02.24):
* finished support for 64bit platforms
* added ASin (arcsin), ACos (arccos) functions
Version 0.6.4 (2007.01.29):
* fixed the problem with a sign in the mathematical parser /-(1) was 1/
* added UInt::AddInt and UInt::SubInt

View File

@ -524,6 +524,218 @@ namespace ttmath
}
/*
*
* inverse trigonometric functions
*
*
*/
/*!
arcus sin
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;
}
/*!
arcus sin
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;
}
/*!
arc sin (x)
x is from <-1,1>
*/
template<class ValueType>
ValueType ASin(ValueType x, ErrorCode * err = 0)
{
ValueType one;
one.SetOne();
bool change_sign = false;
if( x.GreaterWithoutSignThan(one) )
{
if( err )
*err = err_improper_argument;
return one;
}
if( x.IsSign() )
{
change_sign = true;
x.Abs();
}
one.exponent.SubOne(); // =0.5
ValueType result;
// 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;
}
/*!
arc cos (x)
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

View File

@ -187,7 +187,7 @@ public:
*/
void SetPi()
{
static const uint temp_table[] = {
static const unsigned int temp_table[] = {
0xc90fdaa2, 0x2168c234, 0xc4c6628b, 0x80dc1cd1, 0x29024e08, 0x8a67cc74, 0x020bbea6, 0x3b139b22,
0x514a0879, 0x8e3404dd, 0xef9519b3, 0xcd3a431b, 0x302b0a6d, 0xf25f1437, 0x4fe1356d, 0x6d51c245,
0xe485b576, 0x625e7ec6, 0xf44c42e9, 0xa637ed6b, 0x0bff5cb6, 0xf406b7ed, 0xee386bfb, 0x5a899fa5,
@ -199,9 +199,11 @@ public:
0xecfb8504, 0x58dbef0a, 0x8aea7157, 0x5d060c7d, 0xb3970f85, 0xa6e1e4c7, 0xabf5ae8c, 0xdb0933d7,
0x1e8c94e0, 0x4a25619d, 0xcee3d226, 0x1ad2ee6b, 0xf0139f9d, 0x88e637cb
};
// 78 unsigned words
// 78 unsigned 32bits words
// this is static table which represents the value Pi (mantissa of its)
// (first is the highest word)
// we must define this table as 'unsigned int' because
// on 32bits and 64bits platforms this table is 32bits
mantissa.SetFromTable(temp_table, sizeof(temp_table) / sizeof(uint));
exponent = -sint(man)*sint(TTMATH_BITS_PER_UINT) + 2;
@ -235,7 +237,7 @@ public:
*/
void SetE()
{
static const uint temp_table[] = {
static const unsigned int temp_table[] = {
0xadf85458, 0xa2bb4a9a, 0xafdc5620, 0x273d3cf1, 0xd8b9c583, 0xce2d3695, 0xa9e13641, 0x146433fb,
0xcc939dce, 0x249b3ef9, 0x7d2fe363, 0x630c75d8, 0xf681b202, 0xaec4617a, 0xd3df1ed5, 0xd5fd6561,
0x2433f51f, 0x5f066ed0, 0x85636555, 0x3ded1af3, 0xb557135e, 0x7f57c935, 0x984f0c70, 0xe0e68b77,
@ -260,7 +262,7 @@ public:
*/
void SetLn2()
{
static const uint temp_table[] = {
static const unsigned int temp_table[] = {
0xb17217f7, 0xd1cf79ab, 0xc9e3b398, 0x03f2f6af, 0x40f34326, 0x7298b62d, 0x8a0d175b, 0x8baafa2b,
0xe7b87620, 0x6debac98, 0x559552fb, 0x4afa1b10, 0xed2eae35, 0xc1382144, 0x27573b29, 0x1169b825,
0x3e96ca16, 0x224ae8c5, 0x1acbda11, 0x317c387e, 0xb9ea9bc3, 0xb136603b, 0x256fa0ec, 0x7657f74b,
@ -738,15 +740,18 @@ public:
(we also can change 'Div' instead modifying this 'Pow' -- it'll be the same effect,
this error is only when we're using our mathematic parser)
gcc 3.4.6 (FreeBSD) with -O3 doesn't work perfect as well -- there must be
without the reference
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
#ifdef __GNUC__
uint Pow(const Big<exp, man> & pow)
#else
uint Pow(const Big<exp, man> pow)
#endif
//#ifdef __GNUC__
// uint Pow(const Big<exp, man> & pow)
//#else
uint Pow(Big<exp, man> pow)
//#endif
{
TTMATH_REFERENCE_ASSERT( pow )
// TTMATH_REFERENCE_ASSERT( pow )
if( IsZero() )
{

View File

@ -704,6 +704,32 @@ void Min(int sindex, int amount_of_args, ValueType & result)
}
void ASin(int sindex, int amount_of_args, ValueType & result)
{
if( amount_of_args != 1 )
Error( err_improper_amount_of_arguments );
ErrorCode err;
result = ttmath::ASin(stack[sindex].value, &err);
if(err != err_ok)
Error( err );
}
void ACos(int sindex, int amount_of_args, ValueType & result)
{
if( amount_of_args != 1 )
Error( err_improper_amount_of_arguments );
ErrorCode err;
result = ttmath::ACos(stack[sindex].value, &err);
if(err != err_ok)
Error( err );
}
/*!
this method returns the value from a user-defined function
@ -826,6 +852,8 @@ void CreateFunctionsTable()
InsertFunctionToTable(std::string("exp"), &Parser<ValueType>::Exp);
InsertFunctionToTable(std::string("max"), &Parser<ValueType>::Max);
InsertFunctionToTable(std::string("min"), &Parser<ValueType>::Min);
InsertFunctionToTable(std::string("asin"), &Parser<ValueType>::ASin);
InsertFunctionToTable(std::string("acos"), &Parser<ValueType>::ACos);
}

View File

@ -60,8 +60,8 @@
the version of the library
*/
#define TTMATH_MAJOR_VER 0
#define TTMATH_MINOR_VER 6
#define TTMATH_REVISION_VER 4
#define TTMATH_MINOR_VER 7
#define TTMATH_REVISION_VER 0
/*!

View File

@ -122,6 +122,7 @@ public:
SetZero();
}
#if !defined _M_X64 && !defined __x86_64__
/*!
this method copies the value stored in an another table
@ -135,7 +136,7 @@ public:
(this rounding isn't a perfect rounding -- look at the description below)
and if temp_table_len is smaller than value_size we'll clear the rest words
int table
in the table
*/
void SetFromTable(const uint * temp_table, uint temp_table_len)
{
@ -161,7 +162,7 @@ public:
can set a carry and then there'll be a small problem
for optimization
*/
if( table[0] != 0xffffffff )
if( table[0] != TTMATH_UINT_MAX_VALUE )
++table[0];
}
}
@ -178,7 +179,7 @@ public:
*
*/
#if !defined _M_X64 && !defined __x86_64__
/*!
@ -375,7 +376,7 @@ public:
"leal (%%ebx,%%edx,4), %%ebx \n"
"movl %4, %%edx \n"
"movl %%esi, %%edx \n"
"clc \n"
"1: \n"
@ -404,7 +405,7 @@ public:
"pop %%ebx \n"
: "=a" (c)
: "c" (b), "d" (index), "b" (p1), "q" (value)
: "c" (b), "d" (index), "b" (p1), "S" (value)
: "cc", "memory" );
#endif
@ -529,7 +530,7 @@ public:
"movl $0, %%edx \n"
"movl (%%ebx), %%eax \n"
"addl %4, %%eax \n"
"addl %%esi, %%eax \n"
"movl %%eax, (%%ebx) \n"
"inc %%ebx \n"
@ -538,7 +539,7 @@ public:
"inc %%ebx \n"
"movl (%%ebx), %%eax \n"
"adcl %5, %%eax \n"
"adcl %%edi, %%eax \n"
"movl %%eax, (%%ebx) \n"
"jnc 2f \n"
@ -570,7 +571,7 @@ public:
"pop %%ebx \n"
: "=a" (c)
: "c" (b), "d" (index), "b" (p1), "m" (x1), "m" (x2)
: "c" (b), "d" (index), "b" (p1), "S" (x1), "D" (x2)
: "cc", "memory" );
#endif
@ -773,7 +774,7 @@ public:
"leal (%%ebx,%%edx,4), %%ebx \n"
"movl %4, %%edx \n"
"movl %%esi, %%edx \n"
"clc \n"
"1: \n"
@ -802,7 +803,7 @@ public:
"pop %%ebx \n"
: "=a" (c)
: "c" (b), "d" (index), "b" (p1), "q" (value)
: "c" (b), "d" (index), "b" (p1), "S" (value)
: "cc", "memory" );
#endif
@ -1173,8 +1174,8 @@ public:
"movl $-1, %0 \n"
"1: \n"
: "=q" (result)
: "q" (x)
: "=R" (result)
: "R" (x)
: "cc" );
#endif
@ -1336,7 +1337,7 @@ public:
#ifdef __GNUC__
asm __volatile__(
__asm__ __volatile__(
"mull %%edx \n"
@ -3082,6 +3083,7 @@ public:
// these methods for 64bit processors are defined in 'ttmathuint64.h'
void SetFromTable(const unsigned int * temp_table, uint temp_table_len);
uint Add(const UInt<value_size> & ss2, uint c=0);
uint AddInt(uint value, uint index = 0);
uint AddTwoInts(uint x2, uint x1, uint index);

View File

@ -54,6 +54,67 @@ namespace ttmath
#if defined _M_X64 || defined __x86_64__
/*!
this method copies the value stored in an another table
(warning: first values in temp_table are the highest words -- it's different
from our table)
we copy as many words as it is possible
if temp_table_len is bigger than value_size we'll try to round
the lowest word from table depending on the last not used bit in temp_table
(this rounding isn't a perfect rounding -- look at the description below)
and if temp_table_len is smaller than value_size we'll clear the rest words
in the table
warning: we're using 'temp_table' as a pointer at 32bit words
*/
template<uint value_size>
void UInt<value_size>::SetFromTable(const unsigned int * temp_table, uint temp_table_len)
{
uint temp_table_index = 0;
sint i; // 'i' with a sign
for(i=value_size-1 ; i>=0 && temp_table_index<temp_table_len; --i, ++temp_table_index)
{
table[i] = uint(temp_table[ temp_table_index ]) << 32;
++temp_table_index;
if( temp_table_index<temp_table_len )
table[i] |= temp_table[ temp_table_index ];
}
// rounding mantissa
if( temp_table_index < temp_table_len )
{
if( (temp_table[temp_table_index] & TTMATH_UINT_HIGHEST_BIT) != 0 )
{
/*
very simply rounding
if the bit from not used last word from temp_table is set to one
we're rouding the lowest word in the table
in fact there should be a normal addition but
we don't use Add() or AddTwoInts() because these methods
can set a carry and then there'll be a small problem
for optimization
*/
if( table[0] != TTMATH_UINT_MAX_VALUE )
++table[0];
}
}
// cleaning the rest of the mantissa
for( ; i>=0 ; --i)
table[i] = 0;
}
/*!
this method adding ss2 to the this and adding carry if it's defined
(this = this + ss2 + c)
@ -187,7 +248,7 @@ namespace ttmath
"leaq (%%rbx,%%rdx,8), %%rbx \n"
"movq %4, %%rdx \n"
"movq %%rsi, %%rdx \n"
"clc \n"
"1: \n"
@ -220,7 +281,7 @@ namespace ttmath
"pop %%rbx \n"
: "=a" (c)
: "c" (b), "d" (index), "b" (p1), "q" (value)
: "c" (b), "d" (index), "b" (p1), "S" (value)
: "cc", "memory" );
#endif
@ -292,7 +353,7 @@ namespace ttmath
"movq $0, %%rdx \n"
"movq (%%rbx), %%rax \n"
"addq %4, %%rax \n"
"addq %%rsi, %%rax \n"
"movq %%rax, (%%rbx) \n"
"inc %%rbx \n"
@ -305,7 +366,7 @@ namespace ttmath
"inc %%rbx \n"
"movq (%%rbx), %%rax \n"
"adcq %5, %%rax \n"
"adcq %%rdi, %%rax \n"
"movq %%rax, (%%rbx) \n"
"jnc 2f \n"
@ -341,7 +402,7 @@ namespace ttmath
"pop %%rbx \n"
: "=a" (c)
: "c" (b), "d" (index), "b" (p1), "q" (x1), "q" (x2)
: "c" (b), "d" (index), "b" (p1), "S" (x1), "D" (x2)
: "cc", "memory" );
#endif
@ -481,7 +542,7 @@ namespace ttmath
"leaq (%%rbx,%%rdx,8), %%rbx \n"
"movq %4, %%rdx \n"
"movq %%rsi, %%rdx \n"
"clc \n"
"1: \n"
@ -514,7 +575,7 @@ namespace ttmath
"pop %%rbx \n"
: "=a" (c)
: "c" (b), "d" (index), "b" (p1), "q" (value)
: "c" (b), "d" (index), "b" (p1), "S" (value)
: "cc", "memory" );
#endif
@ -775,7 +836,7 @@ namespace ttmath
#ifdef __GNUC__
asm __volatile__(
__asm__ __volatile__(
"mulq %%rdx \n"