Continued fraction expansion #1 for incomplete beta integral Cephes Math Library, Release 2.8: June, 2000 Copyright 1984, 1995, 2000 by Stephen L. Moshier *************************************************************************/
| 662 | Copyright 1984, 1995, 2000 by Stephen L. Moshier |
| 663 | *************************************************************************/ |
| 664 | static double incompletebetafe(double a, |
| 665 | double b, |
| 666 | double x, |
| 667 | double big, |
| 668 | double biginv) |
| 669 | { |
| 670 | double result; |
| 671 | double xk; |
| 672 | double pk; |
| 673 | double pkm1; |
| 674 | double pkm2; |
| 675 | double qk; |
| 676 | double qkm1; |
| 677 | double qkm2; |
| 678 | double k1; |
| 679 | double k2; |
| 680 | double k3; |
| 681 | double k4; |
| 682 | double k5; |
| 683 | double k6; |
| 684 | double k7; |
| 685 | double k8; |
| 686 | double r; |
| 687 | double t; |
| 688 | double ans; |
| 689 | double thresh; |
| 690 | int n; |
| 691 | |
| 692 | k1 = a; |
| 693 | k2 = a+b; |
| 694 | k3 = a; |
| 695 | k4 = a+1.0; |
| 696 | k5 = 1.0; |
| 697 | k6 = b-1.0; |
| 698 | k7 = k4; |
| 699 | k8 = a+2.0; |
| 700 | pkm2 = 0.0; |
| 701 | qkm2 = 1.0; |
| 702 | pkm1 = 1.0; |
| 703 | qkm1 = 1.0; |
| 704 | ans = 1.0; |
| 705 | r = 1.0; |
| 706 | n = 0; |
| 707 | thresh = 3.0*ap::machineepsilon; |
| 708 | do |
| 709 | { |
| 710 | xk = -x*k1*k2/(k3*k4); |
| 711 | pk = pkm1+pkm2*xk; |
| 712 | qk = qkm1+qkm2*xk; |
| 713 | pkm2 = pkm1; |
| 714 | pkm1 = pk; |
| 715 | qkm2 = qkm1; |
| 716 | qkm1 = qk; |
| 717 | xk = x*k5*k6/(k7*k8); |
| 718 | pk = pkm1+pkm2*xk; |
| 719 | qk = qkm1+qkm2*xk; |
| 720 | pkm2 = pkm1; |
| 721 | pkm1 = pk; |