B4J Code Snippet Gamma and Beta functions for computing Cumulative Distribution Functions

There are numerous implementations of Gamma and Beta functions in many computer languages
This is my implementation for B4X as applied to statistics.

http://en.wikipedia.org/wiki/Gamma_function
http://www.stat.tamu.edu/~jnewton/604/chap3.pdf
http://www.aip.de/groups/soe/local/numres/bookfpdf/f6-1.pdf

Sub AppStart (Args() As String)
    Dim markTime As Long = DateTime.Now
    Log("2-tailed z = 3.291   p=" & (NumberFormat2((1-cdfZ(3.291)), 1,4,4,False)))
    Log("right-tail chiSq(14) = 36.1   p=" & NumberFormat2((1-cdfC(36.1, 14)),1,4,4,False))
    Log("right-tail F(3,5) = 33.2025   p=" & NumberFormat2((1-cdfF(33.2025, 3, 5)), 1,4,4,False))
    Log("2-tailed Pearson R(14) = 0.742    p=" & NumberFormat2((1-cdfR(0.742, 14)), 1,4,4,False))
    Log("2-tailed t-Value(13) = 4.221   p=" & NumberFormat2((1-cdfT(4.221, 13)), 1,4,4,False))
    Log((DateTime.now - markTime))   '~15 milliseconds in release mode
End Sub

Sub cdfZ(x As Double) As Double
    If x>=0 Then Return (gammaNormal(.5, x * x / 2)) Else Return (1 - cdfZ(-x))
End Sub

Sub cdfC(x As Double, v As Double) As Double
    If x>=0 Then Return (gammaNormal(.5 * v, x / 2)) Else Return (1 - cdfC(-x, v))
End Sub

Sub cdfT(x As Double, v As Double) As Double
    If x>=0 Then Return (1 - betaNormal(v / 2, .5, v / (v + x*x))) Else Return (1 - cdfT(-x, v))
End Sub

Sub cdfF(x As Double, v1 As Double, v2 As Double) As Double
    If x>=0 Then Return (1 - betaNormal(v2 / 2, v1 / 2, v2 / (v2 + v1*x))) Else Return (1 - cdfF(-x, v1, v2))
End Sub

Sub cdfR(x As Double, v As Double) As Double
    If x>=0 Then Return cdfT(x * Sqrt(v)/Sqrt(1.0 - x*x), v) Else Return cdfR(-x, v)
End Sub

Sub gammaNormal(a As Double, x As Double) As Double
    Dim sum, delta, aStep, result As Double
    result = 0.0
    If (x = 0.0) Then
        result = .5
    Else
        aStep = a: delta = 1.0 / a
        sum = delta
        Dim n As Int = 0
        Do While (Abs(delta) > Abs(sum) * .0000001)
            n = n + 1
            If n>100 Then Return 0.0
            aStep = aStep + 1
            delta = delta * x / aStep
            sum = sum + delta
        Loop
        result = sum * Power(cE, -x + a * Logarithm(x,cE) - logGamma(a))
    End If
    Return result
End Sub

Sub betaNormal(a As Double, b As Double, x As Double) As Double
    If (x < 0.0 Or x > 1.0) Then Return 0.0
    Dim result As Double
    If (x = 0.0 Or x = 1.0) Then
        result = 0.0
    Else
        result = Power(cE, logGamma(a + b) - logGamma(a) - logGamma(b) + a * Logarithm(x, cE) + b * Logarithm(1.0 - x, cE))
    End If
    If (x < (a + 1.0)/(a + b + 2.0)) Then
        result = result * continuousFractions(a, b, x) / a
    Else
        result =  1.0 - result * continuousFractions(b, a, 1.0-x) / b
    End If
    Return result
End Sub

Sub continuousFractions(a As Double, b As Double, x As Double) As Double
    Dim m, m2 As Int
    Dim aa, c, d, delta, result As Double
    c = 1.0
    d = 1.0 / Abs(1.0 - (a +  b) * x / (a + 1))
    result = d
    m = 0: m2 = 0
    Do While Abs(delta-1.0) > .0000001
        m = m + 1
        If m>100 Then Return 0.0
        m2 = m2 + 2
        aa = m * (b - m) * x / ((a + m2) * (a - 1 + m2))
        c = (1.0 + aa/c)
        d = 1.0 / (1.0 + aa*d)
        result = result * c * d
    'alternating two different fractional terms
        aa = -(a + m) * (a + b + m) * x / ((a + m2) * (a + 1 + m2))
        c = (1.0 + aa/c)
        d = 1.0 / (1.0 + aa*d)
        delta = c * d
        result = result * delta
    Loop
    Return result
End Sub

Sub logGamma(z As Double) As Double
    Dim  zStep, denominator, series, Lanczos() As Double   'Lanzos coefficients g=5 n=7
    Lanczos = Array As Double( 1.000000000190015, _
        76.18009172947146, -86.50532032941677, 24.01409824083091, _
        -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5)
    denominator = (z + 5.5) - (z + 0.5) * Logarithm((z + 5.5), cE)
    series = Lanczos(0)
    zStep = z
    For j = 1 To 6
        zStep = zStep + 1
        series = series + Lanczos(j) / zStep
    Next
    Return Logarithm(Sqrt(2*cPI) * series / z, cE) - denominator
End Sub