Power series for incomplete beta integral. Use when b*x is small and x not too close to 1. Cephes Math Library, Release 2.8: June, 2000 Copyright 1984, 1995, 2000 by Stephen L. Moshier *************************************************************************/
| 892 | Copyright 1984, 1995, 2000 by Stephen L. Moshier |
| 893 | *************************************************************************/ |
| 894 | static double incompletebetaps(double a, double b, double x, double maxgam) |
| 895 | { |
| 896 | double result; |
| 897 | double s; |
| 898 | double t; |
| 899 | double u; |
| 900 | double v; |
| 901 | double n; |
| 902 | double t1; |
| 903 | double z; |
| 904 | double ai; |
| 905 | double sg; |
| 906 | |
| 907 | ai = 1.0/a; |
| 908 | u = (1.0-b)*x; |
| 909 | v = u/(a+1.0); |
| 910 | t1 = v; |
| 911 | t = u; |
| 912 | n = 2.0; |
| 913 | s = 0.0; |
| 914 | z = ap::machineepsilon*ai; |
| 915 | while(ap::fp_greater(fabs(v),z)) |
| 916 | { |
| 917 | u = (n-b)*x/n; |
| 918 | t = t*u; |
| 919 | v = t/(a+n); |
| 920 | s = s+v; |
| 921 | n = n+1.0; |
| 922 | } |
| 923 | s = s+t1; |
| 924 | s = s+ai; |
| 925 | u = a*log(x); |
| 926 | if( ap::fp_less(a+b,maxgam)&&ap::fp_less(fabs(u),log(ap::maxrealnumber)) ) |
| 927 | { |
| 928 | t = gamma(a+b)/(gamma(a)*gamma(b)); |
| 929 | s = s*t*pow(x, a); |
| 930 | } |
| 931 | else |
| 932 | { |
| 933 | t = lngamma(a+b, sg)-lngamma(a, sg)-lngamma(b, sg)+u+log(s); |
| 934 | if( ap::fp_less(t,log(ap::minrealnumber)) ) |
| 935 | { |
| 936 | s = 0.0; |
| 937 | } |
| 938 | else |
| 939 | { |
| 940 | s = exp(t); |
| 941 | } |
| 942 | } |
| 943 | result = s; |
| 944 | return result; |
| 945 | } |
| 946 | |
| 947 | |
| 948 |
no test coverage detected