|
Lines 838-851
Link Here
|
| 838 |
|
838 |
|
| 839 |
return aRet; |
839 |
return aRet; |
| 840 |
} |
840 |
} |
| 841 |
|
841 |
/** Calculates erf(x) by using Maclaurin expansion. Accurate up to |
| 842 |
|
842 |
x = 4. |
| 843 |
double Erf( double x ) |
843 |
*/ |
|
|
844 |
void ErfMaclaurin( double x, double& fErf ) |
| 844 |
{ |
845 |
{ |
| 845 |
if( x == 0.0 ) |
846 |
if( x == 0.0 ) |
| 846 |
return 0.0; |
847 |
{ |
|
|
848 |
fErf = 0.0; |
| 849 |
return; |
| 850 |
} |
| 847 |
|
851 |
|
| 848 |
double f, fZm, fN2, fn1, fOld, fT; |
852 |
double fZm, fN2, fn1, fOld, fT; |
| 849 |
sal_Bool bAdd = sal_True; |
853 |
sal_Bool bAdd = sal_True; |
| 850 |
sal_Int32 nMaxIter = 1000; |
854 |
sal_Int32 nMaxIter = 1000; |
| 851 |
|
855 |
|
|
Lines 857-869
Link Here
|
| 857 |
|
861 |
|
| 858 |
fT = x * fZm; |
862 |
fT = x * fZm; |
| 859 |
|
863 |
|
| 860 |
f = x - fT / fN2; |
864 |
fErf = x - fT / fN2; |
| 861 |
|
865 |
|
| 862 |
fOld = f * 0.9; |
866 |
fOld = fErf * 0.9; |
| 863 |
|
867 |
|
| 864 |
while( f != fOld && nMaxIter ) |
868 |
while( fErf != fOld && nMaxIter-- ) |
| 865 |
{ |
869 |
{ |
| 866 |
fOld = f; |
870 |
fOld = fErf; |
| 867 |
|
871 |
|
| 868 |
fN2 += 2.0; |
872 |
fN2 += 2.0; |
| 869 |
|
873 |
|
|
Lines 873-890
Link Here
|
| 873 |
fn1++; |
877 |
fn1++; |
| 874 |
|
878 |
|
| 875 |
if( bAdd ) |
879 |
if( bAdd ) |
| 876 |
f += fT / fN2; |
880 |
fErf += fT / fN2; |
| 877 |
else |
881 |
else |
| 878 |
f -= fT / fN2; |
882 |
fErf -= fT / fN2; |
| 879 |
|
883 |
|
| 880 |
bAdd = !bAdd; |
884 |
bAdd = !bAdd; |
|
|
885 |
} |
| 886 |
|
| 887 |
fErf *= 1.1283791670955125738961589031215452; |
| 888 |
} |
| 889 |
|
| 890 |
/** Calculates erf(x) by using Asymptotic expansion. Accurate for x > 4 |
| 891 |
but converges more quickly than octave/gnumeric. Calculated values |
| 892 |
are almost identical to Excel's. |
| 893 |
|
| 894 |
@author Kohei Yoshida <kohei@openoffice.org> |
| 895 |
|
| 896 |
@see #i19991# |
| 897 |
*/ |
| 898 |
void ErfAsymptotic( double x, double& fErf ) |
| 899 |
{ |
| 900 |
if( x == 0.0 ) |
| 901 |
{ |
| 902 |
fErf = 0.0; |
| 903 |
return; |
| 904 |
} |
| 905 |
|
| 906 |
long nMaxIter = 1000; |
| 907 |
|
| 908 |
double fFct2 = 1.0, fDenom = 1.0, fX = x, fn = 0.0; |
| 909 |
double fXCnt = 1.0; |
| 910 |
double f = 1.0 / fX; |
| 911 |
double fOld = f * 0.9; |
| 881 |
|
912 |
|
| 882 |
nMaxIter--; |
913 |
bool bAdd = false; |
|
|
914 |
double fDiff = 0.0, fDiffOld = 1000.0; |
| 915 |
while( f != fOld && nMaxIter-- ) |
| 916 |
{ |
| 917 |
fOld = f; |
| 918 |
|
| 919 |
fFct2 *= 2.0*++fn - 1.0; |
| 920 |
fDenom *= 2.0; |
| 921 |
fX *= x*x; |
| 922 |
fXCnt += 2.0; |
| 923 |
|
| 924 |
if ( bAdd ) |
| 925 |
f += fFct2 / ( fDenom * fX ); |
| 926 |
else |
| 927 |
f -= fFct2 / ( fDenom * fX ); |
| 928 |
|
| 929 |
bAdd = !bAdd; |
| 930 |
|
| 931 |
fDiff = fabs( f - fOld ); |
| 932 |
if ( fDiffOld < fDiff ) |
| 933 |
{ |
| 934 |
f = fOld; |
| 935 |
break; |
| 936 |
} |
| 937 |
fDiffOld = fDiff; |
| 883 |
} |
938 |
} |
| 884 |
|
939 |
|
| 885 |
return f * 1.128379167095512573896; // * 2/sqrt(PI) |
940 |
fErf = 1.0 - f*exp( -1.0*x*x )*0.5641895835477562869480794515607726; |
| 886 |
} |
941 |
} |
| 887 |
|
942 |
|
|
|
943 |
/** Parent erf function that calls two different algorithms based on value |
| 944 |
of x. It takes care of cases where x is negative ( erf is an odd function |
| 945 |
i.e. f(-x) = -f(x) ). |
| 946 |
|
| 947 |
@author Kohei Yoshida <kohei@openoffice.org> |
| 948 |
|
| 949 |
@see #i19991# |
| 950 |
*/ |
| 951 |
double Erf( double x ) |
| 952 |
{ |
| 953 |
bool bNegative = false; |
| 954 |
if ( x < 0.0 ) |
| 955 |
{ |
| 956 |
x = fabs( x ); |
| 957 |
bNegative = true; |
| 958 |
} |
| 959 |
|
| 960 |
double fErf; |
| 961 |
if ( x > 4.0 ) |
| 962 |
ErfAsymptotic( x, fErf ); |
| 963 |
else |
| 964 |
ErfMaclaurin( x, fErf ); |
| 965 |
|
| 966 |
if ( bNegative ) |
| 967 |
fErf *= -1.0; |
| 968 |
return fErf; |
| 969 |
} |
| 888 |
|
970 |
|
| 889 |
inline sal_Bool IsNum( sal_Unicode c ) |
971 |
inline sal_Bool IsNum( sal_Unicode c ) |
| 890 |
{ |
972 |
{ |