Continued fraction expansion #2 for incomplete beta integral Cephes Math Library, Release 2.8: June, 2000 Copyright 1984, 1995, 2000 by Stephen L. Moshier *************************************************************************/
| 776 | Copyright 1984, 1995, 2000 by Stephen L. Moshier |
| 777 | *************************************************************************/ |
| 778 | static double incompletebetafe2(double a, |
| 779 | double b, |
| 780 | double x, |
| 781 | double big, |
| 782 | double biginv) |
| 783 | { |
| 784 | double result; |
| 785 | double xk; |
| 786 | double pk; |
| 787 | double pkm1; |
| 788 | double pkm2; |
| 789 | double qk; |
| 790 | double qkm1; |
| 791 | double qkm2; |
| 792 | double k1; |
| 793 | double k2; |
| 794 | double k3; |
| 795 | double k4; |
| 796 | double k5; |
| 797 | double k6; |
| 798 | double k7; |
| 799 | double k8; |
| 800 | double r; |
| 801 | double t; |
| 802 | double ans; |
| 803 | double z; |
| 804 | double thresh; |
| 805 | int n; |
| 806 | |
| 807 | k1 = a; |
| 808 | k2 = b-1.0; |
| 809 | k3 = a; |
| 810 | k4 = a+1.0; |
| 811 | k5 = 1.0; |
| 812 | k6 = a+b; |
| 813 | k7 = a+1.0; |
| 814 | k8 = a+2.0; |
| 815 | pkm2 = 0.0; |
| 816 | qkm2 = 1.0; |
| 817 | pkm1 = 1.0; |
| 818 | qkm1 = 1.0; |
| 819 | z = x/(1.0-x); |
| 820 | ans = 1.0; |
| 821 | r = 1.0; |
| 822 | n = 0; |
| 823 | thresh = 3.0*ap::machineepsilon; |
| 824 | do |
| 825 | { |
| 826 | xk = -z*k1*k2/(k3*k4); |
| 827 | pk = pkm1+pkm2*xk; |
| 828 | qk = qkm1+qkm2*xk; |
| 829 | pkm2 = pkm1; |
| 830 | pkm1 = pk; |
| 831 | qkm2 = qkm1; |
| 832 | qkm1 = qk; |
| 833 | xk = z*k5*k6/(k7*k8); |
| 834 | pk = pkm1+pkm2*xk; |
| 835 | qk = qkm1+qkm2*xk; |