diff --git a/CHANGELOG b/CHANGELOG index 4d2d516..f7a9856 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -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 diff --git a/ttmath/ttmath.h b/ttmath/ttmath.h index c39fc73..989b9f4 100644 --- a/ttmath/ttmath.h +++ b/ttmath/ttmath.h @@ -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 + 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 + 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 + 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 + ValueType ACos(const ValueType & x, ErrorCode * err = 0) + { + ValueType temp; + + temp.Set05Pi(); + temp.Sub(ASin(x,err)); + + return temp; + } + + } // namespace diff --git a/ttmath/ttmathbig.h b/ttmath/ttmathbig.h index b78da54..875b042 100644 --- a/ttmath/ttmathbig.h +++ b/ttmath/ttmathbig.h @@ -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 & pow) -#else - uint Pow(const Big pow) -#endif +//#ifdef __GNUC__ +// uint Pow(const Big & pow) +//#else + uint Pow(Big pow) +//#endif { - TTMATH_REFERENCE_ASSERT( pow ) +// TTMATH_REFERENCE_ASSERT( pow ) if( IsZero() ) { diff --git a/ttmath/ttmathparser.h b/ttmath/ttmathparser.h index c7f93af..49c9b04 100644 --- a/ttmath/ttmathparser.h +++ b/ttmath/ttmathparser.h @@ -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::Exp); InsertFunctionToTable(std::string("max"), &Parser::Max); InsertFunctionToTable(std::string("min"), &Parser::Min); + InsertFunctionToTable(std::string("asin"), &Parser::ASin); + InsertFunctionToTable(std::string("acos"), &Parser::ACos); } diff --git a/ttmath/ttmathtypes.h b/ttmath/ttmathtypes.h index 07b8fb3..f7d11e0 100644 --- a/ttmath/ttmathtypes.h +++ b/ttmath/ttmathtypes.h @@ -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 /*! diff --git a/ttmath/ttmathuint.h b/ttmath/ttmathuint.h index ab75fe0..a99a0e1 100644 --- a/ttmath/ttmathuint.h +++ b/ttmath/ttmathuint.h @@ -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 & ss2, uint c=0); uint AddInt(uint value, uint index = 0); uint AddTwoInts(uint x2, uint x1, uint index); diff --git a/ttmath/ttmathuint64.h b/ttmath/ttmathuint64.h index 664ed9f..9687f6a 100644 --- a/ttmath/ttmathuint64.h +++ b/ttmath/ttmathuint64.h @@ -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 + void UInt::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=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"