Stephen Paulger
Library for big numbers from http://www.ttmath.org/
Dependents: PIDHeater82 Conceptcontroller_v_1_0 AlarmClockApp COG4050_adxl355_tilt ... more
Embed:
(wiki syntax)
Show/hide line numbers
ttmath.h
Go to the documentation of this file.
00001 /* 00002 * This file is a part of TTMath Bignum Library 00003 * and is distributed under the (new) BSD licence. 00004 * Author: Tomasz Sowa <t.sowa@ttmath.org> 00005 */ 00006 00007 /* 00008 * Copyright (c) 2006-2012, Tomasz Sowa 00009 * All rights reserved. 00010 * 00011 * Redistribution and use in source and binary forms, with or without 00012 * modification, are permitted provided that the following conditions are met: 00013 * 00014 * * Redistributions of source code must retain the above copyright notice, 00015 * this list of conditions and the following disclaimer. 00016 * 00017 * * Redistributions in binary form must reproduce the above copyright 00018 * notice, this list of conditions and the following disclaimer in the 00019 * documentation and/or other materials provided with the distribution. 00020 * 00021 * * Neither the name Tomasz Sowa nor the names of contributors to this 00022 * project may be used to endorse or promote products derived 00023 * from this software without specific prior written permission. 00024 * 00025 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 00026 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 00027 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 00028 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 00029 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 00030 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 00031 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 00032 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 00033 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 00034 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF 00035 * THE POSSIBILITY OF SUCH DAMAGE. 00036 */ 00037 00038 00039 00040 #ifndef headerfilettmathmathtt 00041 #define headerfilettmathmathtt 00042 00043 /*! 00044 \file ttmath.h 00045 \brief Mathematics functions. 00046 */ 00047 00048 #ifdef _MSC_VER 00049 //warning C4127: conditional expression is constant 00050 #pragma warning( disable: 4127 ) 00051 //warning C4702: unreachable code 00052 #pragma warning( disable: 4702 ) 00053 //warning C4800: forcing value to bool 'true' or 'false' (performance warning) 00054 #pragma warning( disable: 4800 ) 00055 #endif 00056 00057 00058 #include "ttmathbig.h" 00059 #include "ttmathobjects.h" 00060 00061 00062 namespace ttmath 00063 { 00064 /* 00065 * 00066 * functions defined here are used only with Big<> types 00067 * 00068 * 00069 */ 00070 00071 00072 /* 00073 * 00074 * functions for rounding 00075 * 00076 * 00077 */ 00078 00079 00080 /*! 00081 this function skips the fraction from x 00082 e.g 2.2 = 2 00083 2.7 = 2 00084 -2.2 = 2 00085 -2.7 = 2 00086 */ 00087 template<class ValueType> 00088 ValueType SkipFraction (const ValueType & x) 00089 { 00090 ValueType result( x ); 00091 result.SkipFraction(); 00092 00093 return result; 00094 } 00095 00096 00097 /*! 00098 this function rounds to the nearest integer value 00099 e.g 2.2 = 2 00100 2.7 = 3 00101 -2.2 = -2 00102 -2.7 = -3 00103 */ 00104 template<class ValueType> 00105 ValueType Round (const ValueType & x, ErrorCode * err = 0) 00106 { 00107 if( x.IsNan() ) 00108 { 00109 if( err ) 00110 *err = err_improper_argument; 00111 00112 return x; // NaN 00113 } 00114 00115 ValueType result( x ); 00116 uint c = result.Round(); 00117 00118 if( err ) 00119 *err = c ? err_overflow : err_ok; 00120 00121 return result; 00122 } 00123 00124 00125 00126 /*! 00127 this function returns a value representing the smallest integer 00128 that is greater than or equal to x 00129 00130 Ceil(-3.7) = -3 00131 Ceil(-3.1) = -3 00132 Ceil(-3.0) = -3 00133 Ceil(4.0) = 4 00134 Ceil(4.2) = 5 00135 Ceil(4.8) = 5 00136 */ 00137 template<class ValueType> 00138 ValueType Ceil (const ValueType & x, ErrorCode * err = 0) 00139 { 00140 if( x.IsNan() ) 00141 { 00142 if( err ) 00143 *err = err_improper_argument; 00144 00145 return x; // NaN 00146 } 00147 00148 ValueType result(x); 00149 uint c = 0; 00150 00151 result.SkipFraction(); 00152 00153 if( result != x ) 00154 { 00155 // x is with fraction 00156 // if x is negative we don't have to do anything 00157 if( !x.IsSign() ) 00158 { 00159 ValueType one; 00160 one.SetOne(); 00161 00162 c += result.Add(one); 00163 } 00164 } 00165 00166 if( err ) 00167 *err = c ? err_overflow : err_ok; 00168 00169 return result; 00170 } 00171 00172 00173 /*! 00174 this function returns a value representing the largest integer 00175 that is less than or equal to x 00176 00177 Floor(-3.6) = -4 00178 Floor(-3.1) = -4 00179 Floor(-3) = -3 00180 Floor(2) = 2 00181 Floor(2.3) = 2 00182 Floor(2.8) = 2 00183 */ 00184 template<class ValueType> 00185 ValueType Floor (const ValueType & x, ErrorCode * err = 0) 00186 { 00187 if( x.IsNan() ) 00188 { 00189 if( err ) 00190 *err = err_improper_argument; 00191 00192 return x; // NaN 00193 } 00194 00195 ValueType result(x); 00196 uint c = 0; 00197 00198 result.SkipFraction(); 00199 00200 if( result != x ) 00201 { 00202 // x is with fraction 00203 // if x is positive we don't have to do anything 00204 if( x.IsSign() ) 00205 { 00206 ValueType one; 00207 one.SetOne(); 00208 00209 c += result.Sub(one); 00210 } 00211 } 00212 00213 if( err ) 00214 *err = c ? err_overflow : err_ok; 00215 00216 return result; 00217 } 00218 00219 00220 00221 /* 00222 * 00223 * logarithms and the exponent 00224 * 00225 * 00226 */ 00227 00228 00229 /*! 00230 this function calculates the natural logarithm (logarithm with the base 'e') 00231 */ 00232 template<class ValueType> 00233 ValueType Ln (const ValueType & x, ErrorCode * err = 0) 00234 { 00235 if( x.IsNan() ) 00236 { 00237 if( err ) 00238 *err = err_improper_argument; 00239 00240 return x; // NaN 00241 } 00242 00243 ValueType result; 00244 uint state = result.Ln(x); 00245 00246 if( err ) 00247 { 00248 switch( state ) 00249 { 00250 case 0: 00251 *err = err_ok; 00252 break; 00253 case 1: 00254 *err = err_overflow; 00255 break; 00256 case 2: 00257 *err = err_improper_argument; 00258 break; 00259 default: 00260 *err = err_internal_error; 00261 break; 00262 } 00263 } 00264 00265 00266 return result; 00267 } 00268 00269 00270 /*! 00271 this function calculates the logarithm 00272 */ 00273 template<class ValueType> 00274 ValueType Log (const ValueType & x, const ValueType & base, ErrorCode * err = 0) 00275 { 00276 if( x.IsNan() ) 00277 { 00278 if( err ) *err = err_improper_argument; 00279 return x; 00280 } 00281 00282 if( base.IsNan() ) 00283 { 00284 if( err ) *err = err_improper_argument; 00285 return base; 00286 } 00287 00288 ValueType result; 00289 uint state = result.Log(x, base); 00290 00291 if( err ) 00292 { 00293 switch( state ) 00294 { 00295 case 0: 00296 *err = err_ok; 00297 break; 00298 case 1: 00299 *err = err_overflow; 00300 break; 00301 case 2: 00302 case 3: 00303 *err = err_improper_argument; 00304 break; 00305 default: 00306 *err = err_internal_error; 00307 break; 00308 } 00309 } 00310 00311 return result; 00312 } 00313 00314 00315 /*! 00316 this function calculates the expression e^x 00317 */ 00318 template<class ValueType> 00319 ValueType Exp (const ValueType & x, ErrorCode * err = 0) 00320 { 00321 if( x.IsNan() ) 00322 { 00323 if( err ) 00324 *err = err_improper_argument; 00325 00326 return x; // NaN 00327 } 00328 00329 ValueType result; 00330 uint c = result.Exp(x); 00331 00332 if( err ) 00333 *err = c ? err_overflow : err_ok; 00334 00335 return result; 00336 } 00337 00338 00339 /*! 00340 * 00341 * trigonometric functions 00342 * 00343 */ 00344 00345 00346 /* 00347 this namespace consists of auxiliary functions 00348 (something like 'private' in a class) 00349 */ 00350 namespace auxiliaryfunctions 00351 { 00352 00353 /*! 00354 an auxiliary function for calculating the Sine 00355 (you don't have to call this function) 00356 */ 00357 template<class ValueType> 00358 uint PrepareSin (ValueType & x, bool & change_sign) 00359 { 00360 ValueType temp; 00361 00362 change_sign = false; 00363 00364 if( x.IsSign() ) 00365 { 00366 // we're using the formula 'sin(-x) = -sin(x)' 00367 change_sign = !change_sign; 00368 x.ChangeSign(); 00369 } 00370 00371 // we're reducing the period 2*PI 00372 // (for big values there'll always be zero) 00373 temp.Set2Pi(); 00374 00375 if( x.Mod(temp) ) 00376 return 1; 00377 00378 00379 // we're setting 'x' as being in the range of <0, 0.5PI> 00380 00381 temp.SetPi(); 00382 00383 if( x > temp ) 00384 { 00385 // x is in (pi, 2*pi> 00386 x.Sub( temp ); 00387 change_sign = !change_sign; 00388 } 00389 00390 temp.Set05Pi(); 00391 00392 if( x > temp ) 00393 { 00394 // x is in (0.5pi, pi> 00395 x.Sub( temp ); 00396 x = temp - x; 00397 } 00398 00399 return 0; 00400 } 00401 00402 00403 /*! 00404 an auxiliary function for calculating the Sine 00405 (you don't have to call this function) 00406 00407 it returns Sin(x) where 'x' is from <0, PI/2> 00408 we're calculating the Sin with using Taylor series in zero or PI/2 00409 (depending on which point of these two points is nearer to the 'x') 00410 00411 Taylor series: 00412 sin(x) = sin(a) + cos(a)*(x-a)/(1!) 00413 - sin(a)*((x-a)^2)/(2!) - cos(a)*((x-a)^3)/(3!) 00414 + sin(a)*((x-a)^4)/(4!) + ... 00415 00416 when a=0 it'll be: 00417 sin(x) = (x)/(1!) - (x^3)/(3!) + (x^5)/(5!) - (x^7)/(7!) + (x^9)/(9!) ... 00418 00419 and when a=PI/2: 00420 sin(x) = 1 - ((x-PI/2)^2)/(2!) + ((x-PI/2)^4)/(4!) - ((x-PI/2)^6)/(6!) ... 00421 */ 00422 template<class ValueType> 00423 ValueType Sin0pi05 (const ValueType & x) 00424 { 00425 ValueType result; 00426 ValueType numerator, denominator; 00427 ValueType d_numerator, d_denominator; 00428 ValueType one, temp, old_result; 00429 00430 // temp = pi/4 00431 temp.Set05Pi(); 00432 temp.exponent.SubOne(); 00433 00434 one.SetOne(); 00435 00436 if( x < temp ) 00437 { 00438 // we're using the Taylor series with a=0 00439 result = x; 00440 numerator = x; 00441 denominator = one; 00442 00443 // d_numerator = x^2 00444 d_numerator = x; 00445 d_numerator.Mul(x); 00446 00447 d_denominator = 2; 00448 } 00449 else 00450 { 00451 // we're using the Taylor series with a=PI/2 00452 result = one; 00453 numerator = one; 00454 denominator = one; 00455 00456 // d_numerator = (x-pi/2)^2 00457 ValueType pi05; 00458 pi05.Set05Pi(); 00459 00460 temp = x; 00461 temp.Sub( pi05 ); 00462 d_numerator = temp; 00463 d_numerator.Mul( temp ); 00464 00465 d_denominator = one; 00466 } 00467 00468 uint c = 0; 00469 bool addition = false; 00470 00471 old_result = result; 00472 for(uint i=1 ; i<=TTMATH_ARITHMETIC_MAX_LOOP ; ++i) 00473 { 00474 // we're starting from a second part of the formula 00475 c += numerator. Mul( d_numerator ); 00476 c += denominator. Mul( d_denominator ); 00477 c += d_denominator.Add( one ); 00478 c += denominator. Mul( d_denominator ); 00479 c += d_denominator.Add( one ); 00480 temp = numerator; 00481 c += temp.Div(denominator); 00482 00483 if( c ) 00484 // Sin is from <-1,1> and cannot make an overflow 00485 // but the carry can be from the Taylor series 00486 // (then we only break our calculations) 00487 break; 00488 00489 if( addition ) 00490 result.Add( temp ); 00491 else 00492 result.Sub( temp ); 00493 00494 00495 addition = !addition; 00496 00497 // we're testing whether the result has changed after adding 00498 // the next part of the Taylor formula, if not we end the loop 00499 // (it means 'x' is zero or 'x' is PI/2 or this part of the formula 00500 // is too small) 00501 if( result == old_result ) 00502 break; 00503 00504 old_result = result; 00505 } 00506 00507 return result; 00508 } 00509 00510 } // namespace auxiliaryfunctions 00511 00512 00513 00514 /*! 00515 this function calculates the Sine 00516 */ 00517 template<class ValueType> 00518 ValueType Sin (ValueType x, ErrorCode * err = 0) 00519 { 00520 using namespace auxiliaryfunctions; 00521 00522 ValueType one, result; 00523 bool change_sign; 00524 00525 if( x.IsNan() ) 00526 { 00527 if( err ) 00528 *err = err_improper_argument; 00529 00530 return x; 00531 } 00532 00533 if( err ) 00534 *err = err_ok; 00535 00536 if( PrepareSin ( x, change_sign ) ) 00537 { 00538 // x is too big, we cannnot reduce the 2*PI period 00539 // prior to version 0.8.5 the result was zero 00540 00541 // result has NaN flag set by default 00542 00543 if( err ) 00544 *err = err_overflow; // maybe another error code? err_improper_argument? 00545 00546 return result; // NaN is set by default 00547 } 00548 00549 result = Sin0pi05 ( x ); 00550 00551 one.SetOne(); 00552 00553 // after calculations there can be small distortions in the result 00554 if( result > one ) 00555 result = one; 00556 else 00557 if( result.IsSign() ) 00558 // we've calculated the sin from <0, pi/2> and the result 00559 // should be positive 00560 result.SetZero(); 00561 00562 if( change_sign ) 00563 result.ChangeSign(); 00564 00565 return result; 00566 } 00567 00568 00569 /*! 00570 this function calulates the Cosine 00571 we're using the formula cos(x) = sin(x + PI/2) 00572 */ 00573 template<class ValueType> 00574 ValueType Cos (ValueType x, ErrorCode * err = 0) 00575 { 00576 if( x.IsNan() ) 00577 { 00578 if( err ) 00579 *err = err_improper_argument; 00580 00581 return x; // NaN 00582 } 00583 00584 ValueType pi05; 00585 pi05.Set05Pi(); 00586 00587 uint c = x.Add( pi05 ); 00588 00589 if( c ) 00590 { 00591 if( err ) 00592 *err = err_overflow; 00593 00594 return ValueType(); // result is undefined (NaN is set by default) 00595 } 00596 00597 return Sin (x, err); 00598 } 00599 00600 00601 /*! 00602 this function calulates the Tangent 00603 we're using the formula tan(x) = sin(x) / cos(x) 00604 00605 it takes more time than calculating the Tan directly 00606 from for example Taylor series but should be a bit preciser 00607 because Tan receives its values from -infinity to +infinity 00608 and when we calculate it from any series then we can make 00609 a greater mistake than calculating 'sin/cos' 00610 */ 00611 template<class ValueType> 00612 ValueType Tan (const ValueType & x, ErrorCode * err = 0) 00613 { 00614 ValueType result = Cos (x, err); 00615 00616 if( err && *err != err_ok ) 00617 return result; 00618 00619 if( result.IsZero() ) 00620 { 00621 if( err ) 00622 *err = err_improper_argument; 00623 00624 result.SetNan(); 00625 00626 return result; 00627 } 00628 00629 return Sin (x, err) / result; 00630 } 00631 00632 00633 /*! 00634 this function calulates the Tangent 00635 look at the description of Tan(...) 00636 00637 (the abbreviation of Tangent can be 'tg' as well) 00638 */ 00639 template<class ValueType> 00640 ValueType Tg (const ValueType & x, ErrorCode * err = 0) 00641 { 00642 return Tan (x, err); 00643 } 00644 00645 00646 /*! 00647 this function calulates the Cotangent 00648 we're using the formula tan(x) = cos(x) / sin(x) 00649 00650 (why do we make it in this way? 00651 look at information in Tan() function) 00652 */ 00653 template<class ValueType> 00654 ValueType Cot (const ValueType & x, ErrorCode * err = 0) 00655 { 00656 ValueType result = Sin (x, err); 00657 00658 if( err && *err != err_ok ) 00659 return result; 00660 00661 if( result.IsZero() ) 00662 { 00663 if( err ) 00664 *err = err_improper_argument; 00665 00666 result.SetNan(); 00667 00668 return result; 00669 } 00670 00671 return Cos (x, err) / result; 00672 } 00673 00674 00675 /*! 00676 this function calulates the Cotangent 00677 look at the description of Cot(...) 00678 00679 (the abbreviation of Cotangent can be 'ctg' as well) 00680 */ 00681 template<class ValueType> 00682 ValueType Ctg (const ValueType & x, ErrorCode * err = 0) 00683 { 00684 return Cot (x, err); 00685 } 00686 00687 00688 /* 00689 * 00690 * inverse trigonometric functions 00691 * 00692 * 00693 */ 00694 00695 namespace auxiliaryfunctions 00696 { 00697 00698 /*! 00699 an auxiliary function for calculating the Arc Sine 00700 00701 we're calculating asin from the following formula: 00702 asin(x) = x + (1*x^3)/(2*3) + (1*3*x^5)/(2*4*5) + (1*3*5*x^7)/(2*4*6*7) + ... 00703 where abs(x) <= 1 00704 00705 we're using this formula when x is from <0, 1/2> 00706 */ 00707 template<class ValueType> 00708 ValueType ASin_0 (const ValueType & x) 00709 { 00710 ValueType nominator, denominator, nominator_add, nominator_x, denominator_add, denominator_x; 00711 ValueType two, result(x), x2(x); 00712 ValueType nominator_temp, denominator_temp, old_result = result; 00713 uint c = 0; 00714 00715 x2.Mul(x); 00716 two = 2; 00717 00718 nominator.SetOne(); 00719 denominator = two; 00720 nominator_add = nominator; 00721 denominator_add = denominator; 00722 nominator_x = x; 00723 denominator_x = 3; 00724 00725 for(uint i=1 ; i<=TTMATH_ARITHMETIC_MAX_LOOP ; ++i) 00726 { 00727 c += nominator_x.Mul(x2); 00728 nominator_temp = nominator_x; 00729 c += nominator_temp.Mul(nominator); 00730 denominator_temp = denominator; 00731 c += denominator_temp.Mul(denominator_x); 00732 c += nominator_temp.Div(denominator_temp); 00733 00734 // if there is a carry somewhere we only break the calculating 00735 // the result should be ok -- it's from <-pi/2, pi/2> 00736 if( c ) 00737 break; 00738 00739 result.Add(nominator_temp); 00740 00741 if( result == old_result ) 00742 // there's no sense to calculate more 00743 break; 00744 00745 old_result = result; 00746 00747 00748 c += nominator_add.Add(two); 00749 c += denominator_add.Add(two); 00750 c += nominator.Mul(nominator_add); 00751 c += denominator.Mul(denominator_add); 00752 c += denominator_x.Add(two); 00753 } 00754 00755 return result; 00756 } 00757 00758 00759 00760 /*! 00761 an auxiliary function for calculating the Arc Sine 00762 00763 we're calculating asin from the following formula: 00764 asin(x) = pi/2 - sqrt(2)*sqrt(1-x) * asin_temp 00765 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)) + ... 00766 00767 where abs(x) <= 1 00768 00769 we're using this formula when x is from (1/2, 1> 00770 */ 00771 template<class ValueType> 00772 ValueType ASin_1 (const ValueType & x) 00773 { 00774 ValueType nominator, denominator, nominator_add, nominator_x, nominator_x_add, denominator_add, denominator_x; 00775 ValueType denominator2; 00776 ValueType one, two, result; 00777 ValueType nominator_temp, denominator_temp, old_result; 00778 uint c = 0; 00779 00780 two = 2; 00781 00782 one.SetOne(); 00783 nominator = one; 00784 result = one; 00785 old_result = result; 00786 denominator = two; 00787 nominator_add = nominator; 00788 denominator_add = denominator; 00789 nominator_x = one; 00790 nominator_x.Sub(x); 00791 nominator_x_add = nominator_x; 00792 denominator_x = 3; 00793 denominator2 = two; 00794 00795 00796 for(uint i=1 ; i<=TTMATH_ARITHMETIC_MAX_LOOP ; ++i) 00797 { 00798 nominator_temp = nominator_x; 00799 c += nominator_temp.Mul(nominator); 00800 denominator_temp = denominator; 00801 c += denominator_temp.Mul(denominator_x); 00802 c += denominator_temp.Mul(denominator2); 00803 c += nominator_temp.Div(denominator_temp); 00804 00805 // if there is a carry somewhere we only break the calculating 00806 // the result should be ok -- it's from <-pi/2, pi/2> 00807 if( c ) 00808 break; 00809 00810 result.Add(nominator_temp); 00811 00812 if( result == old_result ) 00813 // there's no sense to calculate more 00814 break; 00815 00816 old_result = result; 00817 00818 c += nominator_x.Mul(nominator_x_add); 00819 c += nominator_add.Add(two); 00820 c += denominator_add.Add(two); 00821 c += nominator.Mul(nominator_add); 00822 c += denominator.Mul(denominator_add); 00823 c += denominator_x.Add(two); 00824 c += denominator2.Mul(two); 00825 } 00826 00827 00828 nominator_x_add.exponent.AddOne(); // *2 00829 one.exponent.SubOne(); // =0.5 00830 nominator_x_add.Pow(one); // =sqrt(nominator_x_add) 00831 result.Mul(nominator_x_add); 00832 00833 one.Set05Pi(); 00834 one.Sub(result); 00835 00836 return one; 00837 } 00838 00839 00840 } // namespace auxiliaryfunctions 00841 00842 00843 /*! 00844 this function calculates the Arc Sine 00845 x is from <-1,1> 00846 */ 00847 template<class ValueType> 00848 ValueType ASin (ValueType x, ErrorCode * err = 0) 00849 { 00850 using namespace auxiliaryfunctions; 00851 00852 ValueType result, one; 00853 one.SetOne(); 00854 bool change_sign = false; 00855 00856 if( x.IsNan() ) 00857 { 00858 if( err ) 00859 *err = err_improper_argument; 00860 00861 return x; 00862 } 00863 00864 if( x.GreaterWithoutSignThan(one) ) 00865 { 00866 if( err ) 00867 *err = err_improper_argument; 00868 00869 return result; // NaN is set by default 00870 } 00871 00872 if( x.IsSign() ) 00873 { 00874 change_sign = true; 00875 x.Abs(); 00876 } 00877 00878 one.exponent.SubOne(); // =0.5 00879 00880 // asin(-x) = -asin(x) 00881 if( x.GreaterWithoutSignThan(one) ) 00882 result = ASin_1 (x); 00883 else 00884 result = ASin_0 (x); 00885 00886 if( change_sign ) 00887 result.ChangeSign(); 00888 00889 if( err ) 00890 *err = err_ok; 00891 00892 return result; 00893 } 00894 00895 00896 /*! 00897 this function calculates the Arc Cosine 00898 00899 we're using the formula: 00900 acos(x) = pi/2 - asin(x) 00901 */ 00902 template<class ValueType> 00903 ValueType ACos (const ValueType & x, ErrorCode * err = 0) 00904 { 00905 ValueType temp; 00906 00907 temp.Set05Pi(); 00908 temp.Sub(ASin (x, err)); 00909 00910 return temp; 00911 } 00912 00913 00914 00915 namespace auxiliaryfunctions 00916 { 00917 00918 /*! 00919 an auxiliary function for calculating the Arc Tangent 00920 00921 arc tan (x) where x is in <0; 0.5) 00922 (x can be in (-0.5 ; 0.5) too) 00923 00924 we're using the Taylor series expanded in zero: 00925 atan(x) = x - (x^3)/3 + (x^5)/5 - (x^7)/7 + ... 00926 */ 00927 template<class ValueType> 00928 ValueType ATan0 (const ValueType & x) 00929 { 00930 ValueType nominator, denominator, nominator_add, denominator_add, temp; 00931 ValueType result, old_result; 00932 bool adding = false; 00933 uint c = 0; 00934 00935 result = x; 00936 old_result = result; 00937 nominator = x; 00938 nominator_add = x; 00939 nominator_add.Mul(x); 00940 00941 denominator.SetOne(); 00942 denominator_add = 2; 00943 00944 for(uint i=1 ; i<=TTMATH_ARITHMETIC_MAX_LOOP ; ++i) 00945 { 00946 c += nominator.Mul(nominator_add); 00947 c += denominator.Add(denominator_add); 00948 00949 temp = nominator; 00950 c += temp.Div(denominator); 00951 00952 if( c ) 00953 // the result should be ok 00954 break; 00955 00956 if( adding ) 00957 result.Add(temp); 00958 else 00959 result.Sub(temp); 00960 00961 if( result == old_result ) 00962 // there's no sense to calculate more 00963 break; 00964 00965 old_result = result; 00966 adding = !adding; 00967 } 00968 00969 return result; 00970 } 00971 00972 00973 /*! 00974 an auxiliary function for calculating the Arc Tangent 00975 00976 where x is in <0 ; 1> 00977 */ 00978 template<class ValueType> 00979 ValueType ATan01 (const ValueType & x) 00980 { 00981 ValueType half; 00982 half.Set05(); 00983 00984 /* 00985 it would be better if we chose about sqrt(2)-1=0.41... instead of 0.5 here 00986 00987 because as you can see below: 00988 when x = sqrt(2)-1 00989 abs(x) = abs( (x-1)/(1+x) ) 00990 so when we're calculating values around x 00991 then they will be better converged to each other 00992 00993 for example if we have x=0.4999 then during calculating ATan0(0.4999) 00994 we have to make about 141 iterations but when we have x=0.5 00995 then during calculating ATan0( (x-1)/(1+x) ) we have to make 00996 only about 89 iterations (both for Big<3,9>) 00997 00998 in the future this 0.5 can be changed 00999 */ 01000 if( x.SmallerWithoutSignThan(half) ) 01001 return ATan0 (x); 01002 01003 01004 /* 01005 x>=0.5 and x<=1 01006 (x can be even smaller than 0.5) 01007 01008 y = atac(x) 01009 x = tan(y) 01010 01011 tan(y-b) = (tan(y)-tab(b)) / (1+tan(y)*tan(b)) 01012 y-b = atan( (tan(y)-tab(b)) / (1+tan(y)*tan(b)) ) 01013 y = b + atan( (x-tab(b)) / (1+x*tan(b)) ) 01014 01015 let b = pi/4 01016 tan(b) = tan(pi/4) = 1 01017 y = pi/4 + atan( (x-1)/(1+x) ) 01018 01019 so 01020 atac(x) = pi/4 + atan( (x-1)/(1+x) ) 01021 when x->1 (x converges to 1) the (x-1)/(1+x) -> 0 01022 and we can use ATan0() function here 01023 */ 01024 01025 ValueType n(x),d(x),one,result; 01026 01027 one.SetOne(); 01028 n.Sub(one); 01029 d.Add(one); 01030 n.Div(d); 01031 01032 result = ATan0 (n); 01033 01034 n.Set05Pi(); 01035 n.exponent.SubOne(); // =pi/4 01036 result.Add(n); 01037 01038 return result; 01039 } 01040 01041 01042 /*! 01043 an auxiliary function for calculating the Arc Tangent 01044 where x > 1 01045 01046 we're using the formula: 01047 atan(x) = pi/2 - atan(1/x) for x>0 01048 */ 01049 template<class ValueType> 01050 ValueType ATanGreaterThanPlusOne (const ValueType & x) 01051 { 01052 ValueType temp, atan; 01053 01054 temp.SetOne(); 01055 01056 if( temp.Div(x) ) 01057 { 01058 // if there was a carry here that means x is very big 01059 // and atan(1/x) fast converged to 0 01060 atan.SetZero(); 01061 } 01062 else 01063 atan = ATan01 (temp); 01064 01065 temp.Set05Pi(); 01066 temp.Sub(atan); 01067 01068 return temp; 01069 } 01070 01071 } // namespace auxiliaryfunctions 01072 01073 01074 /*! 01075 this function calculates the Arc Tangent 01076 */ 01077 template<class ValueType> 01078 ValueType ATan (ValueType x) 01079 { 01080 using namespace auxiliaryfunctions; 01081 01082 ValueType one, result; 01083 one.SetOne(); 01084 bool change_sign = false; 01085 01086 if( x.IsNan() ) 01087 return x; 01088 01089 // if x is negative we're using the formula: 01090 // atan(-x) = -atan(x) 01091 if( x.IsSign() ) 01092 { 01093 change_sign = true; 01094 x.Abs(); 01095 } 01096 01097 if( x.GreaterWithoutSignThan(one) ) 01098 result = ATanGreaterThanPlusOne (x); 01099 else 01100 result = ATan01 (x); 01101 01102 if( change_sign ) 01103 result.ChangeSign(); 01104 01105 return result; 01106 } 01107 01108 01109 /*! 01110 this function calculates the Arc Tangent 01111 look at the description of ATan(...) 01112 01113 (the abbreviation of Arc Tangent can be 'atg' as well) 01114 */ 01115 template<class ValueType> 01116 ValueType ATg (const ValueType & x) 01117 { 01118 return ATan (x); 01119 } 01120 01121 01122 /*! 01123 this function calculates the Arc Cotangent 01124 01125 we're using the formula: 01126 actan(x) = pi/2 - atan(x) 01127 */ 01128 template<class ValueType> 01129 ValueType ACot (const ValueType & x) 01130 { 01131 ValueType result; 01132 01133 result.Set05Pi(); 01134 result.Sub(ATan (x)); 01135 01136 return result; 01137 } 01138 01139 01140 /*! 01141 this function calculates the Arc Cotangent 01142 look at the description of ACot(...) 01143 01144 (the abbreviation of Arc Cotangent can be 'actg' as well) 01145 */ 01146 template<class ValueType> 01147 ValueType ACtg (const ValueType & x) 01148 { 01149 return ACot (x); 01150 } 01151 01152 01153 /* 01154 * 01155 * hyperbolic functions 01156 * 01157 * 01158 */ 01159 01160 01161 /*! 01162 this function calculates the Hyperbolic Sine 01163 01164 we're using the formula sinh(x)= ( e^x - e^(-x) ) / 2 01165 */ 01166 template<class ValueType> 01167 ValueType Sinh (const ValueType & x, ErrorCode * err = 0) 01168 { 01169 if( x.IsNan() ) 01170 { 01171 if( err ) 01172 *err = err_improper_argument; 01173 01174 return x; // NaN 01175 } 01176 01177 ValueType ex, emx; 01178 uint c = 0; 01179 01180 c += ex.Exp(x); 01181 c += emx.Exp(-x); 01182 01183 c += ex.Sub(emx); 01184 c += ex.exponent.SubOne(); 01185 01186 if( err ) 01187 *err = c ? err_overflow : err_ok; 01188 01189 return ex; 01190 } 01191 01192 01193 /*! 01194 this function calculates the Hyperbolic Cosine 01195 01196 we're using the formula cosh(x)= ( e^x + e^(-x) ) / 2 01197 */ 01198 template<class ValueType> 01199 ValueType Cosh (const ValueType & x, ErrorCode * err = 0) 01200 { 01201 if( x.IsNan() ) 01202 { 01203 if( err ) 01204 *err = err_improper_argument; 01205 01206 return x; // NaN 01207 } 01208 01209 ValueType ex, emx; 01210 uint c = 0; 01211 01212 c += ex.Exp(x); 01213 c += emx.Exp(-x); 01214 01215 c += ex.Add(emx); 01216 c += ex.exponent.SubOne(); 01217 01218 if( err ) 01219 *err = c ? err_overflow : err_ok; 01220 01221 return ex; 01222 } 01223 01224 01225 /*! 01226 this function calculates the Hyperbolic Tangent 01227 01228 we're using the formula tanh(x)= ( e^x - e^(-x) ) / ( e^x + e^(-x) ) 01229 */ 01230 template<class ValueType> 01231 ValueType Tanh (const ValueType & x, ErrorCode * err = 0) 01232 { 01233 if( x.IsNan() ) 01234 { 01235 if( err ) 01236 *err = err_improper_argument; 01237 01238 return x; // NaN 01239 } 01240 01241 ValueType ex, emx, nominator, denominator; 01242 uint c = 0; 01243 01244 c += ex.Exp(x); 01245 c += emx.Exp(-x); 01246 01247 nominator = ex; 01248 c += nominator.Sub(emx); 01249 denominator = ex; 01250 c += denominator.Add(emx); 01251 01252 c += nominator.Div(denominator); 01253 01254 if( err ) 01255 *err = c ? err_overflow : err_ok; 01256 01257 return nominator; 01258 } 01259 01260 01261 /*! 01262 this function calculates the Hyperbolic Tangent 01263 look at the description of Tanh(...) 01264 01265 (the abbreviation of Hyperbolic Tangent can be 'tgh' as well) 01266 */ 01267 template<class ValueType> 01268 ValueType Tgh (const ValueType & x, ErrorCode * err = 0) 01269 { 01270 return Tanh (x, err); 01271 } 01272 01273 /*! 01274 this function calculates the Hyperbolic Cotangent 01275 01276 we're using the formula coth(x)= ( e^x + e^(-x) ) / ( e^x - e^(-x) ) 01277 */ 01278 template<class ValueType> 01279 ValueType Coth (const ValueType & x, ErrorCode * err = 0) 01280 { 01281 if( x.IsNan() ) 01282 { 01283 if( err ) 01284 *err = err_improper_argument; 01285 01286 return x; // NaN 01287 } 01288 01289 if( x.IsZero() ) 01290 { 01291 if( err ) 01292 *err = err_improper_argument; 01293 01294 return ValueType(); // NaN is set by default 01295 } 01296 01297 ValueType ex, emx, nominator, denominator; 01298 uint c = 0; 01299 01300 c += ex.Exp(x); 01301 c += emx.Exp(-x); 01302 01303 nominator = ex; 01304 c += nominator.Add(emx); 01305 denominator = ex; 01306 c += denominator.Sub(emx); 01307 01308 c += nominator.Div(denominator); 01309 01310 if( err ) 01311 *err = c ? err_overflow : err_ok; 01312 01313 return nominator; 01314 } 01315 01316 01317 /*! 01318 this function calculates the Hyperbolic Cotangent 01319 look at the description of Coth(...) 01320 01321 (the abbreviation of Hyperbolic Cotangent can be 'ctgh' as well) 01322 */ 01323 template<class ValueType> 01324 ValueType Ctgh (const ValueType & x, ErrorCode * err = 0) 01325 { 01326 return Coth (x, err); 01327 } 01328 01329 01330 /* 01331 * 01332 * inverse hyperbolic functions 01333 * 01334 * 01335 */ 01336 01337 01338 /*! 01339 inverse hyperbolic sine 01340 01341 asinh(x) = ln( x + sqrt(x^2 + 1) ) 01342 */ 01343 template<class ValueType> 01344 ValueType ASinh (const ValueType & x, ErrorCode * err = 0) 01345 { 01346 if( x.IsNan() ) 01347 { 01348 if( err ) 01349 *err = err_improper_argument; 01350 01351 return x; // NaN 01352 } 01353 01354 ValueType xx(x), one, result; 01355 uint c = 0; 01356 one.SetOne(); 01357 01358 c += xx.Mul(x); 01359 c += xx.Add(one); 01360 one.exponent.SubOne(); // one=0.5 01361 // xx is >= 1 01362 c += xx.PowFrac(one); // xx=sqrt(xx) 01363 c += xx.Add(x); 01364 c += result.Ln(xx); // xx > 0 01365 01366 // here can only be a carry 01367 if( err ) 01368 *err = c ? err_overflow : err_ok; 01369 01370 return result; 01371 } 01372 01373 01374 /*! 01375 inverse hyperbolic cosine 01376 01377 acosh(x) = ln( x + sqrt(x^2 - 1) ) x in <1, infinity) 01378 */ 01379 template<class ValueType> 01380 ValueType ACosh (const ValueType & x, ErrorCode * err = 0) 01381 { 01382 if( x.IsNan() ) 01383 { 01384 if( err ) 01385 *err = err_improper_argument; 01386 01387 return x; // NaN 01388 } 01389 01390 ValueType xx(x), one, result; 01391 uint c = 0; 01392 one.SetOne(); 01393 01394 if( x < one ) 01395 { 01396 if( err ) 01397 *err = err_improper_argument; 01398 01399 return result; // NaN is set by default 01400 } 01401 01402 c += xx.Mul(x); 01403 c += xx.Sub(one); 01404 // xx is >= 0 01405 // we can't call a PowFrac when the 'x' is zero 01406 // if x is 0 the sqrt(0) is 0 01407 if( !xx.IsZero() ) 01408 { 01409 one.exponent.SubOne(); // one=0.5 01410 c += xx.PowFrac(one); // xx=sqrt(xx) 01411 } 01412 c += xx.Add(x); 01413 c += result.Ln(xx); // xx >= 1 01414 01415 // here can only be a carry 01416 if( err ) 01417 *err = c ? err_overflow : err_ok; 01418 01419 return result; 01420 } 01421 01422 01423 /*! 01424 inverse hyperbolic tangent 01425 01426 atanh(x) = 0.5 * ln( (1+x) / (1-x) ) x in (-1, 1) 01427 */ 01428 template<class ValueType> 01429 ValueType ATanh (const ValueType & x, ErrorCode * err = 0) 01430 { 01431 if( x.IsNan() ) 01432 { 01433 if( err ) 01434 *err = err_improper_argument; 01435 01436 return x; // NaN 01437 } 01438 01439 ValueType nominator(x), denominator, one, result; 01440 uint c = 0; 01441 one.SetOne(); 01442 01443 if( !x.SmallerWithoutSignThan(one) ) 01444 { 01445 if( err ) 01446 *err = err_improper_argument; 01447 01448 return result; // NaN is set by default 01449 } 01450 01451 c += nominator.Add(one); 01452 denominator = one; 01453 c += denominator.Sub(x); 01454 c += nominator.Div(denominator); 01455 c += result.Ln(nominator); 01456 c += result.exponent.SubOne(); 01457 01458 // here can only be a carry 01459 if( err ) 01460 *err = c ? err_overflow : err_ok; 01461 01462 return result; 01463 } 01464 01465 01466 /*! 01467 inverse hyperbolic tantent 01468 */ 01469 template<class ValueType> 01470 ValueType ATgh (const ValueType & x, ErrorCode * err = 0) 01471 { 01472 return ATanh (x, err); 01473 } 01474 01475 01476 /*! 01477 inverse hyperbolic cotangent 01478 01479 acoth(x) = 0.5 * ln( (x+1) / (x-1) ) x in (-infinity, -1) or (1, infinity) 01480 */ 01481 template<class ValueType> 01482 ValueType ACoth (const ValueType & x, ErrorCode * err = 0) 01483 { 01484 if( x.IsNan() ) 01485 { 01486 if( err ) 01487 *err = err_improper_argument; 01488 01489 return x; // NaN 01490 } 01491 01492 ValueType nominator(x), denominator(x), one, result; 01493 uint c = 0; 01494 one.SetOne(); 01495 01496 if( !x.GreaterWithoutSignThan(one) ) 01497 { 01498 if( err ) 01499 *err = err_improper_argument; 01500 01501 return result; // NaN is set by default 01502 } 01503 01504 c += nominator.Add(one); 01505 c += denominator.Sub(one); 01506 c += nominator.Div(denominator); 01507 c += result.Ln(nominator); 01508 c += result.exponent.SubOne(); 01509 01510 // here can only be a carry 01511 if( err ) 01512 *err = c ? err_overflow : err_ok; 01513 01514 return result; 01515 } 01516 01517 01518 /*! 01519 inverse hyperbolic cotantent 01520 */ 01521 template<class ValueType> 01522 ValueType ACtgh (const ValueType & x, ErrorCode * err = 0) 01523 { 01524 return ACoth (x, err); 01525 } 01526 01527 01528 01529 01530 01531 /* 01532 * 01533 * functions for converting between degrees, radians and gradians 01534 * 01535 * 01536 */ 01537 01538 01539 /*! 01540 this function converts degrees to radians 01541 01542 it returns: x * pi / 180 01543 */ 01544 template<class ValueType> 01545 ValueType DegToRad (const ValueType & x, ErrorCode * err = 0) 01546 { 01547 ValueType result, temp; 01548 uint c = 0; 01549 01550 if( x.IsNan() ) 01551 { 01552 if( err ) 01553 *err = err_improper_argument; 01554 01555 return x; 01556 } 01557 01558 result = x; 01559 01560 // it is better to make division first and then multiplication 01561 // the result is more accurate especially when x is: 90,180,270 or 360 01562 temp = 180; 01563 c += result.Div(temp); 01564 01565 temp.SetPi(); 01566 c += result.Mul(temp); 01567 01568 if( err ) 01569 *err = c ? err_overflow : err_ok; 01570 01571 return result; 01572 } 01573 01574 01575 /*! 01576 this function converts radians to degrees 01577 01578 it returns: x * 180 / pi 01579 */ 01580 template<class ValueType> 01581 ValueType RadToDeg (const ValueType & x, ErrorCode * err = 0) 01582 { 01583 ValueType result, delimiter; 01584 uint c = 0; 01585 01586 if( x.IsNan() ) 01587 { 01588 if( err ) 01589 *err = err_improper_argument; 01590 01591 return x; 01592 } 01593 01594 result = 180; 01595 c += result.Mul(x); 01596 01597 delimiter.SetPi(); 01598 c += result.Div(delimiter); 01599 01600 if( err ) 01601 *err = c ? err_overflow : err_ok; 01602 01603 return result; 01604 } 01605 01606 01607 /*! 01608 this function converts degrees in the long format into one value 01609 01610 long format: (degrees, minutes, seconds) 01611 minutes and seconds must be greater than or equal zero 01612 01613 result: 01614 if d>=0 : result= d + ((s/60)+m)/60 01615 if d<0 : result= d - ((s/60)+m)/60 01616 01617 ((s/60)+m)/60 = (s+60*m)/3600 (second version is faster because 01618 there's only one division) 01619 01620 for example: 01621 DegToDeg(10, 30, 0) = 10.5 01622 DegToDeg(10, 24, 35.6)=10.4098(8) 01623 */ 01624 template<class ValueType> 01625 ValueType DegToDeg ( const ValueType & d, const ValueType & m, const ValueType & s, 01626 ErrorCode * err = 0) 01627 { 01628 ValueType delimiter, multipler; 01629 uint c = 0; 01630 01631 if( d.IsNan() || m.IsNan() || s.IsNan() || m.IsSign() || s.IsSign() ) 01632 { 01633 if( err ) 01634 *err = err_improper_argument; 01635 01636 delimiter.SetZeroNan(); // not needed, only to get rid of GCC warning about an uninitialized variable 01637 01638 return delimiter; 01639 } 01640 01641 multipler = 60; 01642 delimiter = 3600; 01643 01644 c += multipler.Mul(m); 01645 c += multipler.Add(s); 01646 c += multipler.Div(delimiter); 01647 01648 if( d.IsSign() ) 01649 multipler.ChangeSign(); 01650 01651 c += multipler.Add(d); 01652 01653 if( err ) 01654 *err = c ? err_overflow : err_ok; 01655 01656 return multipler; 01657 } 01658 01659 01660 /*! 01661 this function converts degrees in the long format to radians 01662 */ 01663 template<class ValueType> 01664 ValueType DegToRad ( const ValueType & d, const ValueType & m, const ValueType & s, 01665 ErrorCode * err = 0) 01666 { 01667 ValueType temp_deg = DegToDeg (d,m,s,err); 01668 01669 if( err && *err!=err_ok ) 01670 return temp_deg; 01671 01672 return DegToRad (temp_deg, err); 01673 } 01674 01675 01676 /*! 01677 this function converts gradians to radians 01678 01679 it returns: x * pi / 200 01680 */ 01681 template<class ValueType> 01682 ValueType GradToRad (const ValueType & x, ErrorCode * err = 0) 01683 { 01684 ValueType result, temp; 01685 uint c = 0; 01686 01687 if( x.IsNan() ) 01688 { 01689 if( err ) 01690 *err = err_improper_argument; 01691 01692 return x; 01693 } 01694 01695 result = x; 01696 01697 // it is better to make division first and then multiplication 01698 // the result is more accurate especially when x is: 100,200,300 or 400 01699 temp = 200; 01700 c += result.Div(temp); 01701 01702 temp.SetPi(); 01703 c += result.Mul(temp); 01704 01705 if( err ) 01706 *err = c ? err_overflow : err_ok; 01707 01708 return result; 01709 } 01710 01711 01712 /*! 01713 this function converts radians to gradians 01714 01715 it returns: x * 200 / pi 01716 */ 01717 template<class ValueType> 01718 ValueType RadToGrad (const ValueType & x, ErrorCode * err = 0) 01719 { 01720 ValueType result, delimiter; 01721 uint c = 0; 01722 01723 if( x.IsNan() ) 01724 { 01725 if( err ) 01726 *err = err_improper_argument; 01727 01728 return x; 01729 } 01730 01731 result = 200; 01732 c += result.Mul(x); 01733 01734 delimiter.SetPi(); 01735 c += result.Div(delimiter); 01736 01737 if( err ) 01738 *err = c ? err_overflow : err_ok; 01739 01740 return result; 01741 } 01742 01743 01744 /*! 01745 this function converts degrees to gradians 01746 01747 it returns: x * 200 / 180 01748 */ 01749 template<class ValueType> 01750 ValueType DegToGrad (const ValueType & x, ErrorCode * err = 0) 01751 { 01752 ValueType result, temp; 01753 uint c = 0; 01754 01755 if( x.IsNan() ) 01756 { 01757 if( err ) 01758 *err = err_improper_argument; 01759 01760 return x; 01761 } 01762 01763 result = x; 01764 01765 temp = 200; 01766 c += result.Mul(temp); 01767 01768 temp = 180; 01769 c += result.Div(temp); 01770 01771 if( err ) 01772 *err = c ? err_overflow : err_ok; 01773 01774 return result; 01775 } 01776 01777 01778 /*! 01779 this function converts degrees in the long format to gradians 01780 */ 01781 template<class ValueType> 01782 ValueType DegToGrad ( const ValueType & d, const ValueType & m, const ValueType & s, 01783 ErrorCode * err = 0) 01784 { 01785 ValueType temp_deg = DegToDeg (d,m,s,err); 01786 01787 if( err && *err!=err_ok ) 01788 return temp_deg; 01789 01790 return DegToGrad (temp_deg, err); 01791 } 01792 01793 01794 /*! 01795 this function converts degrees to gradians 01796 01797 it returns: x * 180 / 200 01798 */ 01799 template<class ValueType> 01800 ValueType GradToDeg (const ValueType & x, ErrorCode * err = 0) 01801 { 01802 ValueType result, temp; 01803 uint c = 0; 01804 01805 if( x.IsNan() ) 01806 { 01807 if( err ) 01808 *err = err_improper_argument; 01809 01810 return x; 01811 } 01812 01813 result = x; 01814 01815 temp = 180; 01816 c += result.Mul(temp); 01817 01818 temp = 200; 01819 c += result.Div(temp); 01820 01821 if( err ) 01822 *err = c ? err_overflow : err_ok; 01823 01824 return result; 01825 } 01826 01827 01828 01829 01830 /* 01831 * 01832 * another functions 01833 * 01834 * 01835 */ 01836 01837 01838 /*! 01839 this function calculates the square root 01840 01841 Sqrt(9) = 3 01842 */ 01843 template<class ValueType> 01844 ValueType Sqrt (ValueType x, ErrorCode * err = 0) 01845 { 01846 if( x.IsNan() || x.IsSign() ) 01847 { 01848 if( err ) 01849 *err = err_improper_argument; 01850 01851 x.SetNan(); 01852 01853 return x; 01854 } 01855 01856 uint c = x.Sqrt(); 01857 01858 if( err ) 01859 *err = c ? err_overflow : err_ok; 01860 01861 return x; 01862 } 01863 01864 01865 01866 namespace auxiliaryfunctions 01867 { 01868 01869 template<class ValueType> 01870 bool RootCheckIndexSign(ValueType & x, const ValueType & index, ErrorCode * err) 01871 { 01872 if( index.IsSign() ) 01873 { 01874 // index cannot be negative 01875 if( err ) 01876 *err = err_improper_argument; 01877 01878 x.SetNan(); 01879 01880 return true; 01881 } 01882 01883 return false; 01884 } 01885 01886 01887 template<class ValueType> 01888 bool RootCheckIndexZero(ValueType & x, const ValueType & index, ErrorCode * err) 01889 { 01890 if( index.IsZero() ) 01891 { 01892 if( x.IsZero() ) 01893 { 01894 // there isn't root(0;0) - we assume it's not defined 01895 if( err ) 01896 *err = err_improper_argument; 01897 01898 x.SetNan(); 01899 01900 return true; 01901 } 01902 01903 // root(x;0) is 1 (if x!=0) 01904 x.SetOne(); 01905 01906 if( err ) 01907 *err = err_ok; 01908 01909 return true; 01910 } 01911 01912 return false; 01913 } 01914 01915 01916 template<class ValueType> 01917 bool RootCheckIndexOne(const ValueType & index, ErrorCode * err) 01918 { 01919 ValueType one; 01920 one.SetOne(); 01921 01922 if( index == one ) 01923 { 01924 //root(x;1) is x 01925 // we do it because if we used the PowFrac function 01926 // we would lose the precision 01927 if( err ) 01928 *err = err_ok; 01929 01930 return true; 01931 } 01932 01933 return false; 01934 } 01935 01936 01937 template<class ValueType> 01938 bool RootCheckIndexTwo(ValueType & x, const ValueType & index, ErrorCode * err) 01939 { 01940 if( index == 2 ) 01941 { 01942 x = Sqrt (x, err); 01943 01944 return true; 01945 } 01946 01947 return false; 01948 } 01949 01950 01951 template<class ValueType> 01952 bool RootCheckIndexFrac(ValueType & x, const ValueType & index, ErrorCode * err) 01953 { 01954 if( !index.IsInteger() ) 01955 { 01956 // index must be integer 01957 if( err ) 01958 *err = err_improper_argument; 01959 01960 x.SetNan(); 01961 01962 return true; 01963 } 01964 01965 return false; 01966 } 01967 01968 01969 template<class ValueType> 01970 bool RootCheckXZero(ValueType & x, ErrorCode * err) 01971 { 01972 if( x.IsZero() ) 01973 { 01974 // root(0;index) is zero (if index!=0) 01975 // RootCheckIndexZero() must be called beforehand 01976 x.SetZero(); 01977 01978 if( err ) 01979 *err = err_ok; 01980 01981 return true; 01982 } 01983 01984 return false; 01985 } 01986 01987 01988 template<class ValueType> 01989 bool RootCheckIndex(ValueType & x, const ValueType & index, ErrorCode * err, bool * change_sign) 01990 { 01991 *change_sign = false; 01992 01993 if( index.Mod2() ) 01994 { 01995 // index is odd (1,3,5...) 01996 if( x.IsSign() ) 01997 { 01998 *change_sign = true; 01999 x.Abs(); 02000 } 02001 } 02002 else 02003 { 02004 // index is even 02005 // x cannot be negative 02006 if( x.IsSign() ) 02007 { 02008 if( err ) 02009 *err = err_improper_argument; 02010 02011 x.SetNan(); 02012 02013 return true; 02014 } 02015 } 02016 02017 return false; 02018 } 02019 02020 02021 template<class ValueType> 02022 uint RootCorrectInteger(ValueType & old_x, ValueType & x, const ValueType & index) 02023 { 02024 if( !old_x.IsInteger() || x.IsInteger() || !index.exponent.IsSign() ) 02025 return 0; 02026 02027 // old_x is integer, 02028 // x is not integer, 02029 // index is relatively small (index.exponent<0 or index.exponent<=0) 02030 // (because we're using a special powering algorithm Big::PowUInt()) 02031 02032 uint c = 0; 02033 02034 ValueType temp(x); 02035 c += temp.Round(); 02036 02037 ValueType temp_round(temp); 02038 c += temp.PowUInt(index); 02039 02040 if( temp == old_x ) 02041 x = temp_round; 02042 02043 return (c==0)? 0 : 1; 02044 } 02045 02046 02047 02048 } // namespace auxiliaryfunctions 02049 02050 02051 02052 /*! 02053 indexth Root of x 02054 index must be integer and not negative <0;1;2;3....) 02055 02056 if index==0 the result is one 02057 if x==0 the result is zero and we assume root(0;0) is not defined 02058 02059 if index is even (2;4;6...) the result is x^(1/index) and x>0 02060 if index is odd (1;2;3;...) the result is either 02061 -(abs(x)^(1/index)) if x<0 or 02062 x^(1/index)) if x>0 02063 02064 (for index==1 the result is equal x) 02065 */ 02066 template<class ValueType> 02067 ValueType Root (ValueType x, const ValueType & index, ErrorCode * err = 0) 02068 { 02069 using namespace auxiliaryfunctions; 02070 02071 if( x.IsNan() || index.IsNan() ) 02072 { 02073 if( err ) 02074 *err = err_improper_argument; 02075 02076 x.SetNan(); 02077 02078 return x; 02079 } 02080 02081 if( RootCheckIndexSign(x, index, err) ) return x; 02082 if( RootCheckIndexZero(x, index, err) ) return x; 02083 if( RootCheckIndexOne ( index, err) ) return x; 02084 if( RootCheckIndexTwo (x, index, err) ) return x; 02085 if( RootCheckIndexFrac(x, index, err) ) return x; 02086 if( RootCheckXZero (x, err) ) return x; 02087 02088 // index integer and index!=0 02089 // x!=0 02090 02091 ValueType old_x(x); 02092 bool change_sign; 02093 02094 if( RootCheckIndex(x, index, err, &change_sign ) ) return x; 02095 02096 ValueType temp; 02097 uint c = 0; 02098 02099 // we're using the formula: root(x ; n) = exp( ln(x) / n ) 02100 c += temp.Ln(x); 02101 c += temp.Div(index); 02102 c += x.Exp(temp); 02103 02104 if( change_sign ) 02105 { 02106 // x is different from zero 02107 x.SetSign(); 02108 } 02109 02110 c += RootCorrectInteger(old_x, x, index); 02111 02112 if( err ) 02113 *err = c ? err_overflow : err_ok; 02114 02115 return x; 02116 } 02117 02118 02119 02120 /*! 02121 absolute value of x 02122 e.g. -2 = 2 02123 2 = 2 02124 */ 02125 template<class ValueType> 02126 ValueType Abs (const ValueType & x) 02127 { 02128 ValueType result( x ); 02129 result.Abs(); 02130 02131 return result; 02132 } 02133 02134 02135 /*! 02136 it returns the sign of the value 02137 e.g. -2 = -1 02138 0 = 0 02139 10 = 1 02140 */ 02141 template<class ValueType> 02142 ValueType Sgn (ValueType x) 02143 { 02144 x.Sgn(); 02145 02146 return x; 02147 } 02148 02149 02150 /*! 02151 the remainder from a division 02152 02153 e.g. 02154 mod( 12.6 ; 3) = 0.6 because 12.6 = 3*4 + 0.6 02155 mod(-12.6 ; 3) = -0.6 bacause -12.6 = 3*(-4) + (-0.6) 02156 mod( 12.6 ; -3) = 0.6 02157 mod(-12.6 ; -3) = -0.6 02158 */ 02159 template<class ValueType> 02160 ValueType Mod (ValueType a, const ValueType & b, ErrorCode * err = 0) 02161 { 02162 if( a.IsNan() || b.IsNan() ) 02163 { 02164 if( err ) 02165 *err = err_improper_argument; 02166 02167 a.SetNan(); 02168 02169 return a; 02170 } 02171 02172 uint c = a.Mod(b); 02173 02174 if( err ) 02175 *err = c ? err_overflow : err_ok; 02176 02177 return a; 02178 } 02179 02180 02181 02182 namespace auxiliaryfunctions 02183 { 02184 02185 /*! 02186 this function is used to store factorials in a given container 02187 'more' means how many values should be added at the end 02188 02189 e.g. 02190 std::vector<ValueType> fact; 02191 SetFactorialSequence(fact, 3); 02192 // now the container has three values: 1 1 2 02193 02194 SetFactorialSequence(fact, 2); 02195 // now the container has five values: 1 1 2 6 24 02196 */ 02197 template<class ValueType> 02198 void SetFactorialSequence (std::vector<ValueType> & fact, uint more = 20) 02199 { 02200 if( more == 0 ) 02201 more = 1; 02202 02203 uint start = static_cast<uint >(fact.size()); 02204 fact.resize(fact.size() + more); 02205 02206 if( start == 0 ) 02207 { 02208 fact[0] = 1; 02209 ++start; 02210 } 02211 02212 for(uint i=start ; i<fact.size() ; ++i) 02213 { 02214 fact[i] = fact[i-1]; 02215 fact[i].MulInt(i); 02216 } 02217 } 02218 02219 02220 /*! 02221 an auxiliary function used to calculate Bernoulli numbers 02222 02223 this function returns a sum: 02224 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)!)) 02225 02226 you should have sufficient factorials in cgamma.fact 02227 (cgamma.fact should have at least m items) 02228 02229 n_ should be equal 2 02230 */ 02231 template<class ValueType> 02232 ValueType SetBernoulliNumbersSum (CGamma<ValueType> & cgamma, const ValueType & n_, uint m, 02233 const volatile StopCalculating * stop = 0) 02234 { 02235 ValueType k_, temp, temp2, temp3, sum; 02236 02237 sum.SetZero(); 02238 02239 for(uint k=0 ; k<m ; ++k) // k<m means k<=m-1 02240 { 02241 if( stop && (k & 15)==0 ) // means: k % 16 == 0 02242 if( stop->WasStopSignal() ) 02243 return ValueType(); // NaN 02244 02245 if( k>1 && (k & 1) == 1 ) // for that k the Bernoulli number is zero 02246 continue; 02247 02248 k_ = k; 02249 02250 temp = n_; // n_ is equal 2 02251 temp.Pow(k_); 02252 // temp = 2^k 02253 02254 temp2 = cgamma.fact [m]; 02255 temp3 = cgamma.fact [k]; 02256 temp3.Mul(cgamma.fact [m-k]); 02257 temp2.Div(temp3); 02258 // temp2 = (m k) = m! / ( k! * (m-k)! ) 02259 02260 temp.Mul(temp2); 02261 temp.Mul(cgamma.bern [k]); 02262 02263 sum.Add(temp); 02264 // sum += 2^k * (m k) * B(k) 02265 02266 if( sum.IsNan() ) 02267 break; 02268 } 02269 02270 return sum; 02271 } 02272 02273 02274 /*! 02275 an auxiliary function used to calculate Bernoulli numbers 02276 start is >= 2 02277 02278 we use the recurrence formula: 02279 B(m) = 1 / (2*(1 - 2^m)) * sum(m) 02280 where sum(m) is calculated by SetBernoulliNumbersSum() 02281 */ 02282 template<class ValueType> 02283 bool SetBernoulliNumbersMore (CGamma<ValueType> & cgamma, uint start, const volatile StopCalculating * stop = 0) 02284 { 02285 ValueType denominator, temp, temp2, temp3, m_, sum, sum2, n_, k_; 02286 02287 const uint n = 2; 02288 n_ = n; 02289 02290 // start is >= 2 02291 for(uint m=start ; m<cgamma.bern .size() ; ++m) 02292 { 02293 if( (m & 1) == 1 ) 02294 { 02295 cgamma.bern [m].SetZero(); 02296 } 02297 else 02298 { 02299 m_ = m; 02300 02301 temp = n_; // n_ = 2 02302 temp.Pow(m_); 02303 // temp = 2^m 02304 02305 denominator.SetOne(); 02306 denominator.Sub(temp); 02307 if( denominator.exponent.AddOne() ) // it means: denominator.MulInt(2) 02308 denominator.SetNan(); 02309 02310 // denominator = 2 * (1 - 2^m) 02311 02312 cgamma.bern [m] = SetBernoulliNumbersSum (cgamma, n_, m, stop); 02313 02314 if( stop && stop->WasStopSignal() ) 02315 { 02316 cgamma.bern .resize(m); // valid numbers are in [0, m-1] 02317 return false; 02318 } 02319 02320 cgamma.bern [m].Div(denominator); 02321 } 02322 } 02323 02324 return true; 02325 } 02326 02327 02328 /*! 02329 this function is used to calculate Bernoulli numbers, 02330 returns false if there was a stop signal, 02331 'more' means how many values should be added at the end 02332 02333 e.g. 02334 typedef Big<1,2> MyBig; 02335 CGamma<MyBig> cgamma; 02336 SetBernoulliNumbers(cgamma, 3); 02337 // now we have three first Bernoulli numbers: 1 -0.5 0.16667 02338 02339 SetBernoulliNumbers(cgamma, 4); 02340 // now we have 7 Bernoulli numbers: 1 -0.5 0.16667 0 -0.0333 0 0.0238 02341 */ 02342 template<class ValueType> 02343 bool SetBernoulliNumbers (CGamma<ValueType> & cgamma, uint more = 20, const volatile StopCalculating * stop = 0) 02344 { 02345 if( more == 0 ) 02346 more = 1; 02347 02348 uint start = static_cast<uint >(cgamma.bern .size()); 02349 cgamma.bern .resize(cgamma.bern .size() + more); 02350 02351 if( start == 0 ) 02352 { 02353 cgamma.bern [0].SetOne(); 02354 ++start; 02355 } 02356 02357 if( cgamma.bern .size() == 1 ) 02358 return true; 02359 02360 if( start == 1 ) 02361 { 02362 cgamma.bern [1].Set05(); 02363 cgamma.bern [1].ChangeSign(); 02364 ++start; 02365 } 02366 02367 // we should have sufficient factorials in cgamma.fact 02368 if( cgamma.fact .size() < cgamma.bern .size() ) 02369 SetFactorialSequence (cgamma.fact , static_cast<uint>(cgamma.bern .size() - cgamma.fact .size())); 02370 02371 02372 return SetBernoulliNumbersMore (cgamma, start, stop); 02373 } 02374 02375 02376 /*! 02377 an auxiliary function used to calculate the Gamma() function 02378 02379 we calculate a sum: 02380 sum(n) = sum_{m=2} { B(m) / ( (m^2 - m) * n^(m-1) ) } = 1/(12*n) - 1/(360*n^3) + 1/(1260*n^5) + ... 02381 B(m) means a mth Bernoulli number 02382 the sum starts from m=2, we calculate as long as the value will not change after adding a next part 02383 */ 02384 template<class ValueType> 02385 ValueType GammaFactorialHighSum (const ValueType & n, CGamma<ValueType> & cgamma, ErrorCode & err, 02386 const volatile StopCalculating * stop) 02387 { 02388 ValueType temp, temp2, denominator, sum, oldsum; 02389 02390 sum.SetZero(); 02391 02392 for(uint m=2 ; m<TTMATH_ARITHMETIC_MAX_LOOP ; m+=2) 02393 { 02394 if( stop && (m & 3)==0 ) // (m & 3)==0 means: (m % 4)==0 02395 if( stop->WasStopSignal() ) 02396 { 02397 err = err_interrupt; 02398 return ValueType(); // NaN 02399 } 02400 02401 temp = (m-1); 02402 denominator = n; 02403 denominator.Pow(temp); 02404 // denominator = n ^ (m-1) 02405 02406 temp = m; 02407 temp2 = temp; 02408 temp.Mul(temp2); 02409 temp.Sub(temp2); 02410 // temp = m^2 - m 02411 02412 denominator.Mul(temp); 02413 // denominator = (m^2 - m) * n ^ (m-1) 02414 02415 if( m >= cgamma.bern .size() ) 02416 { 02417 if( !SetBernoulliNumbers (cgamma, m - cgamma.bern .size() + 1 + 3, stop) ) // 3 more than needed 02418 { 02419 // there was the stop signal 02420 err = err_interrupt; 02421 return ValueType(); // NaN 02422 } 02423 } 02424 02425 temp = cgamma.bern [m]; 02426 temp.Div(denominator); 02427 02428 oldsum = sum; 02429 sum.Add(temp); 02430 02431 if( sum.IsNan() || oldsum==sum ) 02432 break; 02433 } 02434 02435 return sum; 02436 } 02437 02438 02439 /*! 02440 an auxiliary function used to calculate the Gamma() function 02441 02442 we calculate a helper function GammaFactorialHigh() by using Stirling's series: 02443 n! = (n/e)^n * sqrt(2*pi*n) * exp( sum(n) ) 02444 where n is a real number (not only an integer) and is sufficient large (greater than TTMATH_GAMMA_BOUNDARY) 02445 and sum(n) is calculated by GammaFactorialHighSum() 02446 */ 02447 template<class ValueType> 02448 ValueType GammaFactorialHigh (const ValueType & n, CGamma<ValueType> & cgamma, ErrorCode & err, 02449 const volatile StopCalculating * stop) 02450 { 02451 ValueType temp, temp2, temp3, denominator, sum; 02452 02453 temp.Set2Pi(); 02454 temp.Mul(n); 02455 temp2 = Sqrt (temp); 02456 // temp2 = sqrt(2*pi*n) 02457 02458 temp = n; 02459 temp3.SetE(); 02460 temp.Div(temp3); 02461 temp.Pow(n); 02462 // temp = (n/e)^n 02463 02464 sum = GammaFactorialHighSum (n, cgamma, err, stop); 02465 temp3.Exp(sum); 02466 // temp3 = exp(sum) 02467 02468 temp.Mul(temp2); 02469 temp.Mul(temp3); 02470 02471 return temp; 02472 } 02473 02474 02475 /*! 02476 an auxiliary function used to calculate the Gamma() function 02477 02478 Gamma(x) = GammaFactorialHigh(x-1) 02479 */ 02480 template<class ValueType> 02481 ValueType GammaPlusHigh (ValueType n, CGamma<ValueType> & cgamma, ErrorCode & err, const volatile StopCalculating * stop) 02482 { 02483 ValueType one; 02484 02485 one.SetOne(); 02486 n.Sub(one); 02487 02488 return GammaFactorialHigh (n, cgamma, err, stop); 02489 } 02490 02491 02492 /*! 02493 an auxiliary function used to calculate the Gamma() function 02494 02495 we use this function when n is integer and a small value (from 0 to TTMATH_GAMMA_BOUNDARY] 02496 we use the formula: 02497 gamma(n) = (n-1)! = 1 * 2 * 3 * ... * (n-1) 02498 */ 02499 template<class ValueType> 02500 ValueType GammaPlusLowIntegerInt (uint n, CGamma<ValueType> & cgamma) 02501 { 02502 TTMATH_ASSERT( n > 0 ) 02503 02504 if( n - 1 < static_cast<uint>(cgamma.fact .size()) ) 02505 return cgamma.fact [n - 1]; 02506 02507 ValueType res; 02508 uint start = 2; 02509 02510 if( cgamma.fact .size() < 2 ) 02511 { 02512 res.SetOne(); 02513 } 02514 else 02515 { 02516 start = static_cast<uint >(cgamma.fact .size()); 02517 res = cgamma.fact [start-1]; 02518 } 02519 02520 for(uint i=start ; i<n ; ++i) 02521 res.MulInt(i); 02522 02523 return res; 02524 } 02525 02526 02527 /*! 02528 an auxiliary function used to calculate the Gamma() function 02529 02530 we use this function when n is integer and a small value (from 0 to TTMATH_GAMMA_BOUNDARY] 02531 */ 02532 template<class ValueType> 02533 ValueType GammaPlusLowInteger (const ValueType & n, CGamma<ValueType> & cgamma) 02534 { 02535 sint n_; 02536 02537 n.ToInt(n_); 02538 02539 return GammaPlusLowIntegerInt (n_, cgamma); 02540 } 02541 02542 02543 /*! 02544 an auxiliary function used to calculate the Gamma() function 02545 02546 we use this function when n is a small value (from 0 to TTMATH_GAMMA_BOUNDARY] 02547 we use a recurrence formula: 02548 gamma(z+1) = z * gamma(z) 02549 then: gamma(z) = gamma(z+1) / z 02550 02551 e.g. 02552 gamma(3.89) = gamma(2001.89) / ( 3.89 * 4.89 * 5.89 * ... * 1999.89 * 2000.89 ) 02553 */ 02554 template<class ValueType> 02555 ValueType GammaPlusLow (ValueType n, CGamma<ValueType> & cgamma, ErrorCode & err, const volatile StopCalculating * stop) 02556 { 02557 ValueType one, denominator, temp, boundary; 02558 02559 if( n.IsInteger() ) 02560 return GammaPlusLowInteger (n, cgamma); 02561 02562 one.SetOne(); 02563 denominator = n; 02564 boundary = TTMATH_GAMMA_BOUNDARY; 02565 02566 while( n < boundary ) 02567 { 02568 n.Add(one); 02569 denominator.Mul(n); 02570 } 02571 02572 n.Add(one); 02573 02574 // now n is sufficient big 02575 temp = GammaPlusHigh (n, cgamma, err, stop); 02576 temp.Div(denominator); 02577 02578 return temp; 02579 } 02580 02581 02582 /*! 02583 an auxiliary function used to calculate the Gamma() function 02584 */ 02585 template<class ValueType> 02586 ValueType GammaPlus (const ValueType & n, CGamma<ValueType> & cgamma, ErrorCode & err, const volatile StopCalculating * stop) 02587 { 02588 if( n > TTMATH_GAMMA_BOUNDARY ) 02589 return GammaPlusHigh (n, cgamma, err, stop); 02590 02591 return GammaPlusLow (n, cgamma, err, stop); 02592 } 02593 02594 02595 /*! 02596 an auxiliary function used to calculate the Gamma() function 02597 02598 this function is used when n is negative 02599 we use the reflection formula: 02600 gamma(1-z) * gamma(z) = pi / sin(pi*z) 02601 then: gamma(z) = pi / (sin(pi*z) * gamma(1-z)) 02602 02603 */ 02604 template<class ValueType> 02605 ValueType GammaMinus (const ValueType & n, CGamma<ValueType> & cgamma, ErrorCode & err, const volatile StopCalculating * stop) 02606 { 02607 ValueType pi, denominator, temp, temp2; 02608 02609 if( n.IsInteger() ) 02610 { 02611 // gamma function is not defined when n is negative and integer 02612 err = err_improper_argument; 02613 return temp; // NaN 02614 } 02615 02616 pi.SetPi(); 02617 02618 temp = pi; 02619 temp.Mul(n); 02620 temp2 = Sin (temp); 02621 // temp2 = sin(pi * n) 02622 02623 temp.SetOne(); 02624 temp.Sub(n); 02625 temp = GammaPlus (temp, cgamma, err, stop); 02626 // temp = gamma(1 - n) 02627 02628 temp.Mul(temp2); 02629 pi.Div(temp); 02630 02631 return pi; 02632 } 02633 02634 } // namespace auxiliaryfunctions 02635 02636 02637 02638 /*! 02639 this function calculates the Gamma function 02640 02641 it's multithread safe, you should create a CGamma<> object and use it whenever you call the Gamma() 02642 e.g. 02643 typedef Big<1,2> MyBig; 02644 MyBig x=234, y=345.53; 02645 CGamma<MyBig> cgamma; 02646 std::cout << Gamma(x, cgamma) << std::endl; 02647 std::cout << Gamma(y, cgamma) << std::endl; 02648 in the CGamma<> object the function stores some coefficients (factorials, Bernoulli numbers), 02649 and they will be reused in next calls to the function 02650 02651 each thread should have its own CGamma<> object, and you can use these objects with Factorial() function too 02652 */ 02653 template<class ValueType> 02654 ValueType Gamma (const ValueType & n, CGamma<ValueType> & cgamma, ErrorCode * err = 0, 02655 const volatile StopCalculating * stop = 0) 02656 { 02657 using namespace auxiliaryfunctions; 02658 02659 ValueType result; 02660 ErrorCode err_tmp; 02661 02662 if( n.IsNan() ) 02663 { 02664 if( err ) 02665 *err = err_improper_argument; 02666 02667 return n; 02668 } 02669 02670 if( cgamma.history .Get(n, result, err_tmp) ) 02671 { 02672 if( err ) 02673 *err = err_tmp; 02674 02675 return result; 02676 } 02677 02678 err_tmp = err_ok; 02679 02680 if( n.IsSign() ) 02681 { 02682 result = GammaMinus (n, cgamma, err_tmp, stop); 02683 } 02684 else 02685 if( n.IsZero() ) 02686 { 02687 err_tmp = err_improper_argument; 02688 result.SetNan(); 02689 } 02690 else 02691 { 02692 result = GammaPlus (n, cgamma, err_tmp, stop); 02693 } 02694 02695 if( result.IsNan() && err_tmp==err_ok ) 02696 err_tmp = err_overflow; 02697 02698 if( err ) 02699 *err = err_tmp; 02700 02701 if( stop && !stop->WasStopSignal() ) 02702 cgamma.history .Add(n, result, err_tmp); 02703 02704 return result; 02705 } 02706 02707 02708 /*! 02709 this function calculates the Gamma function 02710 02711 note: this function should be used only in a single-thread environment 02712 */ 02713 template<class ValueType> 02714 ValueType Gamma (const ValueType & n, ErrorCode * err = 0) 02715 { 02716 // warning: this static object is not thread safe 02717 static CGamma<ValueType> cgamma; 02718 02719 return Gamma (n, cgamma, err); 02720 } 02721 02722 02723 02724 namespace auxiliaryfunctions 02725 { 02726 02727 /*! 02728 an auxiliary function for calculating the factorial function 02729 02730 we use the formula: 02731 x! = gamma(x+1) 02732 */ 02733 template<class ValueType> 02734 ValueType Factorial2 (ValueType x, 02735 CGamma<ValueType> * cgamma = 0, 02736 ErrorCode * err = 0, 02737 const volatile StopCalculating * stop = 0) 02738 { 02739 ValueType result, one; 02740 02741 if( x.IsNan() || x.IsSign() || !x.IsInteger() ) 02742 { 02743 if( err ) 02744 *err = err_improper_argument; 02745 02746 x.SetNan(); 02747 02748 return x; 02749 } 02750 02751 one.SetOne(); 02752 x.Add(one); 02753 02754 if( cgamma ) 02755 return Gamma (x, *cgamma, err, stop); 02756 02757 return Gamma (x, err); 02758 } 02759 02760 } // namespace auxiliaryfunctions 02761 02762 02763 02764 /*! 02765 the factorial from given 'x' 02766 e.g. 02767 Factorial(4) = 4! = 1*2*3*4 02768 02769 it's multithread safe, you should create a CGamma<> object and use it whenever you call the Factorial() 02770 e.g. 02771 typedef Big<1,2> MyBig; 02772 MyBig x=234, y=54345; 02773 CGamma<MyBig> cgamma; 02774 std::cout << Factorial(x, cgamma) << std::endl; 02775 std::cout << Factorial(y, cgamma) << std::endl; 02776 in the CGamma<> object the function stores some coefficients (factorials, Bernoulli numbers), 02777 and they will be reused in next calls to the function 02778 02779 each thread should have its own CGamma<> object, and you can use these objects with Gamma() function too 02780 */ 02781 template<class ValueType> 02782 ValueType Factorial (const ValueType & x, CGamma<ValueType> & cgamma, ErrorCode * err = 0, 02783 const volatile StopCalculating * stop = 0) 02784 { 02785 return auxiliaryfunctions::Factorial2 (x, &cgamma, err, stop); 02786 } 02787 02788 02789 /*! 02790 the factorial from given 'x' 02791 e.g. 02792 Factorial(4) = 4! = 1*2*3*4 02793 02794 note: this function should be used only in a single-thread environment 02795 */ 02796 template<class ValueType> 02797 ValueType Factorial (const ValueType & x, ErrorCode * err = 0) 02798 { 02799 return auxiliaryfunctions::Factorial2 (x, (CGamma<ValueType> *)0, err, 0); 02800 } 02801 02802 02803 /*! 02804 this method prepares some coefficients: factorials and Bernoulli numbers 02805 stored in 'fact' and 'bern' objects 02806 02807 we're defining the method here because we're using Gamma() function which 02808 is not available in ttmathobjects.h 02809 02810 read the doc info in ttmathobjects.h file where CGamma<> struct is declared 02811 */ 02812 template<class ValueType> 02813 void CGamma<ValueType>::InitAll () 02814 { 02815 ValueType x = TTMATH_GAMMA_BOUNDARY + 1; 02816 02817 // history.Remove(x) removes only one object 02818 // we must be sure that there are not others objects with the key 'x' 02819 while( history.Remove(x) ) 02820 { 02821 } 02822 02823 // the simplest way to initialize is to call the Gamma function with (TTMATH_GAMMA_BOUNDARY + 1) 02824 // when x is larger then fewer coefficients we need 02825 Gamma (x, *this); 02826 } 02827 02828 02829 02830 } // namespace 02831 02832 02833 /*! 02834 this is for convenience for the user 02835 he can only use '#include <ttmath/ttmath.h>' 02836 */ 02837 #include "ttmathparser.h" 02838 02839 // Dec is not finished yet 02840 //#include "ttmathdec.h" 02841 02842 02843 02844 #ifdef _MSC_VER 02845 //warning C4127: conditional expression is constant 02846 #pragma warning( default: 4127 ) 02847 //warning C4702: unreachable code 02848 #pragma warning( default: 4702 ) 02849 //warning C4800: forcing value to bool 'true' or 'false' (performance warning) 02850 #pragma warning( default: 4800 ) 02851 #endif 02852 02853 #endif 02854
Generated on Tue Jul 12 2022 14:03:17 by
1.7.2