From e14e65002b62afe05e73100b009070b78c70243c Mon Sep 17 00:00:00 2001 From: Tomasz Sowa Date: Tue, 27 Feb 2007 20:18:33 +0000 Subject: [PATCH] fixed: removed 'const' from some methods (operators: += -= *= /=) in the Big class fixed: bad sizes in tables in some 'Set...' methods in the Big class fixed: Big::FromInt(Int value) - the sign must be set at the end because SetSign checks whether there is zero and depends on it sets the sign or not (this was the stupid error which causes sometimes the errors 'overflow during printing') fixed: Big::SetMin - the sign must be set at the and changed: Big::Pow can use the reference now (the problem was actually with the Big::FromInt) added: a namespace 'auxiliaryfunctions' (in ttmath.h) added: ATan - arc tan, ACTan - arc ctan git-svn-id: svn://ttmath.org/publicrep/ttmath/trunk@17 e52654a7-88a9-db11-a3e9-0013d4bc506e --- CHANGELOG | 6 ++ TODO | 1 + ttmath/ttmath.h | 224 +++++++++++++++++++++++++++++++++++++++++- ttmath/ttmathbig.h | 148 ++++++++++++---------------- ttmath/ttmathint.h | 7 +- ttmath/ttmathparser.h | 20 ++++ ttmath/ttmathtypes.h | 2 +- ttmath/ttmathuint.h | 6 +- ttmath/ttmathuint64.h | 2 +- 9 files changed, 322 insertions(+), 94 deletions(-) diff --git a/CHANGELOG b/CHANGELOG index f7a9856..0aacdd1 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -1,3 +1,9 @@ +Version 0.7.1 (2007.02.27): + * fixed the error 'overflow during printing' which was caused by Big::FromInt(Int value) + (the sign has to be set at the end) + * fixed many small errors + * added ATan (arctan), ACTan (arc ctan) functions + Version 0.7.0 (2007.02.24): * finished support for 64bit platforms * added ASin (arcsin), ACos (arccos) functions diff --git a/TODO b/TODO index 4683f6c..f7e18ca 100644 --- a/TODO +++ b/TODO @@ -4,3 +4,4 @@ TODO TTMath Library * to add the method Mod (a remainder from a division) for the Big type * to add operators (or functions) and, or, xor * to add functions for generating random values +* to add constructors (UInt, Int) for 64bit platforms (constractors which take int and unsigned) diff --git a/ttmath/ttmath.h b/ttmath/ttmath.h index 989b9f4..38db0dd 100644 --- a/ttmath/ttmath.h +++ b/ttmath/ttmath.h @@ -261,7 +261,9 @@ namespace ttmath * trigonometric functions * */ - + + namespace auxiliaryfunctions + { /*! an auxiliary function for calculating the Sin @@ -420,6 +422,7 @@ namespace ttmath return result; } + } // namespace auxiliaryfunctions /*! this function calulates the Sin @@ -427,6 +430,8 @@ namespace ttmath template ValueType Sin(ValueType x) { + using namespace auxiliaryfunctions; + ValueType one; bool change_sign; @@ -531,6 +536,9 @@ namespace ttmath * */ + namespace auxiliaryfunctions + { + /*! arcus sin @@ -673,6 +681,9 @@ namespace ttmath } + } // namespace auxiliaryfunctions + + /*! arc sin (x) x is from <-1,1> @@ -680,6 +691,8 @@ namespace ttmath template ValueType ASin(ValueType x, ErrorCode * err = 0) { + using namespace auxiliaryfunctions; + ValueType one; one.SetOne(); bool change_sign = false; @@ -736,10 +749,213 @@ namespace ttmath } + + namespace auxiliaryfunctions + { + + + /*! + 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 + 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; + } + + + /*! + arc tan (x) where x is in <0 ; 1> + */ + template + ValueType ATan01(const ValueType & x) + { + ValueType half; + half.SetDotOne(); + + /* + 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; + } + + + /*! + arc tan (x) where x > 1 + + we're using the formula: + atan(x) = pi/2 - atan(1/x) for x>0 + */ + template + 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 + + + /*! + arc tan + */ + template + ValueType ATan(ValueType x) + { + using namespace auxiliaryfunctions; + + ValueType one, result; + one.SetOne(); + bool change_sign = false; + + // 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; + } + + + /*! + arc ctan + + we're using the formula: + actan(x) = pi/2 - atan(x) + */ + template + ValueType ACTan(const ValueType & x) + { + ValueType result; + + result.Set05Pi(); + result.Sub(ATan(x)); + + return result; + } + + } // namespace - - - #endif diff --git a/ttmath/ttmathbig.h b/ttmath/ttmathbig.h index 875b042..6eadebd 100644 --- a/ttmath/ttmathbig.h +++ b/ttmath/ttmathbig.h @@ -45,6 +45,8 @@ #include "ttmathint.h" + + namespace ttmath { @@ -59,8 +61,8 @@ class Big /* value = mantissa * 2^exponent - exponent - integer value with a sign - mantissa - integer value without a sing + exponent - an integer value with a sign + mantissa - an integer value without a sing mantissa must be pushed into the left side that is the highest bit from mantissa must be one (of course if there's another value than zero) -- this job @@ -205,7 +207,7 @@ public: // 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)); + mantissa.SetFromTable(temp_table, sizeof(temp_table) / sizeof(int)); exponent = -sint(man)*sint(TTMATH_BITS_PER_UINT) + 2; info = 0; } @@ -250,7 +252,7 @@ public: 0xb4130c93, 0xbc437944, 0xf4fd4452, 0xe2d74dd3, 0x645b2194, 0x41468794 }; - mantissa.SetFromTable(temp_table, sizeof(temp_table) / sizeof(uint)); + mantissa.SetFromTable(temp_table, sizeof(temp_table) / sizeof(int)); exponent = -sint(man)*sint(TTMATH_BITS_PER_UINT) + 2; info = 0; } @@ -275,7 +277,7 @@ public: 0x44a02554, 0x731cdc8e, 0xa17293d1, 0x228a4ef8, 0x6e1adf84, 0x08689fa8 }; - mantissa.SetFromTable(temp_table, sizeof(temp_table) / sizeof(uint)); + mantissa.SetFromTable(temp_table, sizeof(temp_table) / sizeof(int)); exponent = -sint(man)*sint(TTMATH_BITS_PER_UINT); info = 0; } @@ -301,10 +303,10 @@ public: void SetMin() { info = 0; - SetSign(); mantissa.SetMaxValue(); exponent.SetMaxValue(); + SetSign(); // we don't have to use 'Standardizing()' because the last bit from // the mantissa is set @@ -334,7 +336,7 @@ public: /*! it clears the sign - (there'll be a absolute value) + (there'll be an absolute value) e.g. -1 -> 1 @@ -397,6 +399,7 @@ public: { Int exp_offset( exponent ); Int mantissa_size_in_bits( man * TTMATH_BITS_PER_UINT ); + uint c = 0; exp_offset.Sub( ss2.exponent ); @@ -410,7 +413,8 @@ public: ss2 = *this; *this = temp; } - + + if( exp_offset > mantissa_size_in_bits ) { // the second value is too short for taking into consideration in the sum @@ -729,29 +733,9 @@ public: 1 - carry 2 - incorrect argument ('this' or 'pow') */ - /* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - there should be 'sint Pow(const Big & pow)' - but vc2005express doesn't want to compile it perfect, that means - when using 'Maximize Speed /O2' the result of compilation doesn't work property - for example 10^(1/2) is a big value - i don't know where is the problem, with this source code or in the compilator - (it is when we're using 'Big<3,10>' or bigger values in parsing) - /gcc 3.4.2 works perfect (with -O3 optimalization flag)/ - - (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(Big pow) -//#endif + uint Pow(const Big & pow) { -// TTMATH_REFERENCE_ASSERT( pow ) + TTMATH_REFERENCE_ASSERT( pow ) if( IsZero() ) { @@ -767,6 +751,7 @@ public: Big pow_frac( pow ); pow_frac.RemainFraction(); + if( pow_frac.IsZero() ) return PowBInt( pow ); @@ -781,7 +766,7 @@ public: uint c = temp.Ln(*this); c += temp.Mul(pow); c += Exp(temp); - + return (c==0)? 0 : 1; } @@ -1038,6 +1023,7 @@ public: // m will be the value of the mantissa in range <1,2) Big m(x); m.exponent = -sint(man*TTMATH_BITS_PER_UINT - 1); + LnSurrounding1(m); Big exponent_temp; @@ -1048,7 +1034,6 @@ public: Big ln2; ln2.SetLn2(); - c += exponent_temp.Mul(ln2); c += Add(exponent_temp); @@ -1282,6 +1267,7 @@ public: FromInt( sint(value) ); } + #if defined _M_X64 || defined __x86_64__ /*! a constructor for converting 'int' to this class @@ -1301,17 +1287,18 @@ public: void FromInt(Int value) { info = 0; + bool is_sign = false; if( value.IsSign() ) { value.ChangeSign(); - SetSign(); + is_sign = true; } uint minimum_size = (int_size < man)? int_size : man; - sint compensation = (sint)value.CompensationToLeft(); + sint compensation = (sint)value.CompensationToLeft(); exponent = (sint(int_size)-sint(man)) * sint(TTMATH_BITS_PER_UINT) - compensation; - + // copying the highest words uint i; for(i=1 ; i<=minimum_size ; ++i) @@ -1321,6 +1308,9 @@ public: for( ; i<=man ; ++i) mantissa.table[man-i] = 0; + if( is_sign ) + SetSign(); + } @@ -1367,6 +1357,7 @@ public: /*! the default assignment operator */ + Big & operator=(const Big & value) { info = value.info; @@ -1375,26 +1366,27 @@ public: return *this; } - + /*! a constructor for copying from another object of this class */ + Big(const Big & value) { operator=(value); } - + /*! a method for converting the value into the string with a base equal 'base' input: - base - the base on which the value will be showed + base - the base on which the value will be shown - if 'always_scientific' is true the result will be showed in 'scientific' mode + if 'always_scientific' is true the result will be shown in 'scientific' mode - if 'always_scientific' is false the result will be showed + if 'always_scientific' is false the result will be shown either as 'scientific' or 'normal' mode, it depends whether the abs(exponent) is greater than 'when_scientific' or not, if it's greater the value will be printed as 'scientific' @@ -1407,8 +1399,8 @@ public: (zero characters at the end -- after the comma operator) if 'max_digit_after_comma' is equal or greater than zero - that only 'max_digit_after_comma' after the comma operator will be showed - (if 'max_digit_after_comma' is equal zero there'll be showed only + that only 'max_digit_after_comma' after the comma operator will be shown + (if 'max_digit_after_comma' is equal zero there'll be shown only integer value without the comma) for example when the value is: 12.345678 and max_digit_after_comma is 4 @@ -1431,7 +1423,7 @@ public: { static char error_overflow_msg[] = "overflow"; result.erase(); - + if(base<2 || base>16) { result = error_overflow_msg; @@ -1612,19 +1604,17 @@ private: // the sign don't interest us here temp.mantissa = mantissa; temp.exponent = exponent; - c += temp.Div( base_ ); // moving all bits of the mantissa into the right // (how many times to move depend on the exponent) c += temp.ToString_MoveMantissaIntoRight(); - // on account of we take 'new_exp' as small as it's + // because we took 'new_exp' as small as it was // possible ([log base (2^exponent)] + 1) that after the division // (temp.Div( base_ )) the value of exponent should be equal zero or // minimum smaller than zero then we've got the mantissa which has // maximum valid bits - temp.mantissa.ToString(new_man, base); // because we had used a bigger type for calculating I think we @@ -1651,6 +1641,7 @@ private: uint ToString_Log(const Big & x, uint base) { TTMATH_REFERENCE_ASSERT( x ) + TTMATH_ASSERT( base>=2 && base<=16 ) Big temp; temp.SetOne(); @@ -1669,41 +1660,29 @@ private: // (LnSurrounding1() will return one immediately) uint c = Ln(x); + static Big log_history[15] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; + uint index = base - 2; + + if( log_history[index].IsZero() ) + { + // we don't have 'base' in 'log_history' then we calculate it now + Big base_(base); + c += temp.Ln(base_); + c += Div(temp); + + // the next time we'll get the 'Ln(base)' from the history, + // this 'log_history' can have (16-2+1) items max + log_history[index] = temp; + } + else + { + // we've calculated the 'Ln(base)' beforehand and we're getting it now + c += Div( log_history[index] ); + } + + return (c==0)? 0 : 1; + } - /* - * - * this is only temporarily (for testing) - * - */ - -// static Big log_history[15] = { 0l,0l,0l,0l,0l,0l,0l,0l,0l,0l, -// 0l,0l,0l,0l,0l }; - static Big log_history[15] = { 0,0,0,0,0,0,0,0,0,0, - 0,0,0,0,0 }; - - Big * log_value = log_history + base - 2; - - if( log_value->IsZero() ) - { - // we don't have 'base' in 'log_history' then we calculate it now - Big base_(base); - c += temp.Ln(base_); - c += Div(temp); - - // the next time we'll get the 'Ln(base)' from the history, - // this 'log_history' can have (16-2+1) items max - - *log_value = temp; - } - else - { - // we've calculated the 'Ln(base)' beforehand and we're using it now - - c += Div( *log_value ); - } - - return (c==0)? 0 : 1; - } /*! @@ -1722,6 +1701,7 @@ private: if( !exponent.IsSign() ) return 1; + if( exponent <= -sint(man*TTMATH_BITS_PER_UINT) ) // if 'exponent' is <= than '-sint(man*TTMATH_BITS_PER_UINT)' // it means that we must cut the whole mantissa @@ -2613,7 +2593,7 @@ public: return temp; } - Big & operator-=(const Big & ss2) const + Big & operator-=(const Big & ss2) { Sub(ss2); @@ -2631,7 +2611,7 @@ public: } - Big & operator+=(const Big & ss2) const + Big & operator+=(const Big & ss2) { Add(ss2); @@ -2649,7 +2629,7 @@ public: } - Big & operator*=(const Big & ss2) const + Big & operator*=(const Big & ss2) { Mul(ss2); @@ -2667,7 +2647,7 @@ public: } - Big & operator/=(const Big & ss2) const + Big & operator/=(const Big & ss2) { Div(ss2); diff --git a/ttmath/ttmathint.h b/ttmath/ttmathint.h index c250d27..6cb2800 100644 --- a/ttmath/ttmathint.h +++ b/ttmath/ttmathint.h @@ -493,8 +493,11 @@ public: UInt::table[i] = p.table[i]; - if( i < value_size ) +// if( i < value_size ) + + if( value_size > argument_size ) { + // 'this' is longer than 'p' uint fill = (p.table[argument_size-1] & TTMATH_UINT_HIGHEST_BIT)? TTMATH_UINT_MAX_VALUE : 0; for( ; i type to this class + this operator converts an UInt type to this class it doesn't return a carry */ diff --git a/ttmath/ttmathparser.h b/ttmath/ttmathparser.h index 49c9b04..7cf02b8 100644 --- a/ttmath/ttmathparser.h +++ b/ttmath/ttmathparser.h @@ -730,6 +730,24 @@ void ACos(int sindex, int amount_of_args, ValueType & result) } +void ATan(int sindex, int amount_of_args, ValueType & result) +{ + if( amount_of_args != 1 ) + Error( err_improper_amount_of_arguments ); + + result = ttmath::ATan(stack[sindex].value); +} + + +void ACTan(int sindex, int amount_of_args, ValueType & result) +{ + if( amount_of_args != 1 ) + Error( err_improper_amount_of_arguments ); + + result = ttmath::ACTan(stack[sindex].value); +} + + /*! this method returns the value from a user-defined function @@ -854,6 +872,8 @@ void CreateFunctionsTable() InsertFunctionToTable(std::string("min"), &Parser::Min); InsertFunctionToTable(std::string("asin"), &Parser::ASin); InsertFunctionToTable(std::string("acos"), &Parser::ACos); + InsertFunctionToTable(std::string("atan"), &Parser::ATan); + InsertFunctionToTable(std::string("actan"), &Parser::ACTan); } diff --git a/ttmath/ttmathtypes.h b/ttmath/ttmathtypes.h index f7d11e0..aea5774 100644 --- a/ttmath/ttmathtypes.h +++ b/ttmath/ttmathtypes.h @@ -61,7 +61,7 @@ */ #define TTMATH_MAJOR_VER 0 #define TTMATH_MINOR_VER 7 -#define TTMATH_REVISION_VER 0 +#define TTMATH_REVISION_VER 1 /*! diff --git a/ttmath/ttmathuint.h b/ttmath/ttmathuint.h index a99a0e1..4b43954 100644 --- a/ttmath/ttmathuint.h +++ b/ttmath/ttmathuint.h @@ -2519,7 +2519,7 @@ public: /*! - this operator convert an UInt type to this class + this operator converts an UInt type to this class it doesn't return a carry */ @@ -2589,7 +2589,6 @@ public: { } - /*! a copy constructor */ @@ -2599,6 +2598,7 @@ public: } + /*! a template for producting constructors for copying from another types */ @@ -2610,6 +2610,8 @@ public: } + + /*! a destructor */ diff --git a/ttmath/ttmathuint64.h b/ttmath/ttmathuint64.h index 9687f6a..123ddac 100644 --- a/ttmath/ttmathuint64.h +++ b/ttmath/ttmathuint64.h @@ -109,7 +109,7 @@ namespace ttmath } // cleaning the rest of the mantissa - for( ; i>=0 ; --i) + for( ; i >= 0 ; --i) table[i] = 0; }