Changeset 1023
Legend:
 Unmodified
 Added
 Removed

trunk/src/real.cpp
r1022 r1023 679 679 } 680 680 681 static real fast_exp(real const &x) 682 { 683 /* This fast exp method is tuned to work on the [0..1] range and 681 real log10(real const &x) 682 { 683 return log(x) * real::R_LOG10E; 684 } 685 686 static real fast_exp_sub(real const &x, real const &y) 687 { 688 /* This fast exp method is tuned to work on the [1..1] range and 684 689 * no effort whatsoever was made to improve convergence outside this 685 * domain of validity. */ 686 real ret = real::R_1, fact = real::R_1, xn = x; 690 * domain of validity. The argument y is used for cases where we 691 * don't want the leading 1 in the Taylor series. */ 692 real ret = real::R_1  y, fact = real::R_1, xn = x; 687 693 688 694 for (int i = 1; ; i++) … … 722 728 int e0 = x / real::R_LN2; 723 729 real x0 = x  (real)e0 * real::R_LN2; 724 real x1 = fast_exp (x0);730 real x1 = fast_exp_sub(x0, real::R_0); 725 731 x1.m_signexp += e0; 726 732 return x1; … … 732 738 int e0 = x; 733 739 real x0 = x  (real)e0; 734 real x1 = fast_exp (x0 * real::R_LN2);740 real x1 = fast_exp_sub(x0 * real::R_LN2, real::R_0); 735 741 x1.m_signexp += e0; 736 742 return x1; 743 } 744 745 real sinh(real const &x) 746 { 747 /* We cannot always use (exp(x)exp(x))/2 because we'll lose 748 * accuracy near zero. We only use this identity for x>0.5. If 749 * x<=0.5, we compute exp(x)1 and exp(x)1 instead. */ 750 bool near_zero = (fabs(x) < real::R_1 >> 1); 751 real x1 = near_zero ? fast_exp_sub(x, real::R_1) : exp(x); 752 real x2 = near_zero ? fast_exp_sub(x, real::R_1) : exp(x); 753 return (x1  x2) >> 1; 754 } 755 756 real tanh(real const &x) 757 { 758 /* See sinh() for the strategy here */ 759 bool near_zero = (fabs(x) < real::R_1 >> 1); 760 real x1 = near_zero ? fast_exp_sub(x, real::R_1) : exp(x); 761 real x2 = near_zero ? fast_exp_sub(x, real::R_1) : exp(x); 762 real x3 = near_zero ? x1 + x2 + real::R_2 : x1 + x2; 763 return (x1  x2) / x3; 764 } 765 766 real cosh(real const &x) 767 { 768 /* No need to worry about accuracy here; maybe the last bit is slightly 769 * off, but that's about it. */ 770 return (exp(x) + exp(x)) >> 1; 737 771 } 738 772 
trunk/src/real.h
r1022 r1023 73 73 74 74 /* Hyperbolic functions */ 75 /* FIXME: sinh cosh tanh */ 75 friend real sinh(real const &x); 76 friend real cosh(real const &x); 77 friend real tanh(real const &x); 76 78 77 79 /* Exponential and logarithmic functions */ … … 80 82 friend real log(real const &x); 81 83 friend real log2(real const &x); 82 /* FIXME: frexp ldexp log10 modf */ 84 friend real log10(real const &x); 85 /* FIXME: frexp ldexp modf */ 83 86 84 87 /* Power functions */
Note: See TracChangeset
for help on using the changeset viewer.