Notebook

In [53]:

using BenchmarkTools
using SpecialFunctions

incbeta_I(x, a, b) = beta_inc(a, b, x)[1]

function ibetar_st(x::T, a, b;
        HUGE_VAL = T(Inf), FDeps = eps(T), rtol = eps(T), MAXITERS=200
    ) where T<:AbstractFloat
    a ≤ 0 && return HUGE_VAL
    if b ≤ 0
        x < 1 && return zero(T)
        x == 1 && return one(T)
        return HUGE_VAL
    end
    x > (a+1)/(a+b+2) && return 1 - ibetar_st(1 - x, b, a; HUGE_VAL, FDeps, rtol, MAXITERS)
    x ≤ FDeps && return zero(T)
    p1 = zero(T)
    q1 = one(T)
    p2 = exp(a*log(x) + b*log(1-x) + loggamma(a+b) - loggamma(a) - loggamma(b)) / a
    q2 = one(T)
    previous = HUGE_VAL
    k = 0
    while k < MAXITERS && !isapprox(p2, previous; rtol)
        previous = p2
        d = -(a+k)*(a+b+k)*x / ((a+2k)*(a+2k+1))
        p1 = p1*d + p2
        q1 = q1*d + q2
        k += 1
        d = k*(b-k)*x /((a+2k-1)*(a+2k))
        p2 = p2*d + p1
        q2 = q2*d + q1
        if q2 ≤ FDeps
            p2 = HUGE_VAL
        else
            p1 /= q2
            q1 /= q2
            p2 /= q2
            q2 = one(T)
        end
    end
    !isapprox(p2, previous; rtol) && error("ibetar did not converge.")
    p2
end

Out[53]:

ibetar_st (generic function with 1 method)

In [76]:

a, b = 100.0, 0.01
x = 0.1
@show incbeta_I(x, a, b)
@show ibetar_st(x, a, b);

incbeta_I(x, a, b) = 1.1675552034918687e-104
ibetar_st(x, a, b) = 1.1675552034918835e-104

In [77]:

@btime incbeta_I($x, $a, $b)
@btime ibetar_st($x, $a, $b);

  143.853 ns (0 allocations: 0 bytes)
  149.333 ns (0 allocations: 0 bytes)

In [78]:

@which SpecialFunctions._beta_inc(a, b, x)

Out[78]:

_beta_inc(a::Float64, b::Float64, x::Float64) in SpecialFunctions at D:\.julia\packages\SpecialFunctions\hefUc\src\beta_inc.jl:738

In [91]:

@code_warntype SpecialFunctions._beta_inc(a, b, x, 1-x)

MethodInstance for SpecialFunctions._beta_inc(::Float64, ::Float64, ::Float64, ::Float64)
  from _beta_inc(a::Float64, b::Float64, x::Float64, y::Float64) in SpecialFunctions at D:\.julia\packages\SpecialFunctions\hefUc\src\beta_inc.jl:738
Arguments
  #self#::Core.Const(SpecialFunctions._beta_inc)
  a::Float64
  b::Float64
  x::Float64
  y::Float64
Locals
  p::Float64
  q::Float64
  b0::Float64
  n::Int64
  x0::Float64
  y0::Float64
  a0::Float64
  ind::Bool
  lambda::Float64
  epps::Float64
  z::Float64
  @_17::Float64
Body::Tuple{Float64, Float64}
1 ──        Core.NewvarNode(:(b0))
│           Core.NewvarNode(:(n))
│           Core.NewvarNode(:(x0))
│           Core.NewvarNode(:(y0))
│           Core.NewvarNode(:(a0))
│           Core.NewvarNode(:(ind))
│           Core.NewvarNode(:(lambda))
│           Core.NewvarNode(:(epps))
│           Core.NewvarNode(:(z))
│           (p = 0.0)
│           (q = 0.0)
│    %12  = (a < 0.0)::Bool
└───        goto #3 if not %12
2 ──        goto #4
3 ── %15  = (b < 0.0)::Bool
└───        goto #5 if not %15
4 ┄─ %17  = Core.tuple(a, b)::Tuple{Float64, Float64}
│    %18  = SpecialFunctions.DomainError(%17, "a or b is negative")::Core.PartialStruct(DomainError, Any[Tuple{Float64, Float64}, String])
│           SpecialFunctions.throw(%18)
└───        Core.Const(:(return %19))
5 ┄─ %21  = (a == 0.0)::Bool
└───        goto #8 if not %21
6 ── %23  = (b == 0.0)::Bool
└───        goto #8 if not %23
7 ── %25  = Core.tuple(a, b)::Tuple{Float64, Float64}
│    %26  = SpecialFunctions.DomainError(%25, "a and b are 0.0")::Core.PartialStruct(DomainError, Any[Tuple{Float64, Float64}, String])
│           SpecialFunctions.throw(%26)
└───        Core.Const(:(return %27))
8 ┄─ %29  = (x < 0.0)::Bool
└───        goto #10 if not %29
9 ──        goto #11
10 ─ %32  = (x > 1.0)::Bool
└───        goto #12 if not %32
11 ┄ %34  = SpecialFunctions.DomainError(x, "x < 0 or x > 1")::Core.PartialStruct(DomainError, Any[Float64, String])
│           SpecialFunctions.throw(%34)
└───        Core.Const(:(return %35))
12 ┄ %37  = (y < 0.0)::Bool
└───        goto #14 if not %37
13 ─        goto #15
14 ─ %40  = (y > 1.0)::Bool
└───        goto #16 if not %40
15 ┄ %42  = SpecialFunctions.DomainError(y, "y < 0 or y > 1")::Core.PartialStruct(DomainError, Any[Float64, String])
│           SpecialFunctions.throw(%42)
└───        Core.Const(:(return %43))
16 ┄ %45  = (x + y)::Float64
│           (z = %45 - 1.0)
│    %47  = SpecialFunctions.abs(z)::Float64
│    %48  = SpecialFunctions.eps()::Core.Const(2.220446049250313e-16)
│    %49  = (3.0 * %48)::Core.Const(6.661338147750939e-16)
│    %50  = (%47 > %49)::Bool
└───        goto #18 if not %50
17 ─ %52  = Core.tuple(x, y)::Tuple{Float64, Float64}
│    %53  = SpecialFunctions.DomainError(%52, "x + y != 1.0")::Core.PartialStruct(DomainError, Any[Tuple{Float64, Float64}, String])
│           SpecialFunctions.throw(%53)
└───        Core.Const(:(return %54))
18 ┄ %56  = (x == 0.0)::Bool
└───        goto #20 if not %56
19 ─ %58  = Core.tuple(0.0, 1.0)::Core.Const((0.0, 1.0))
└───        return %58
20 ─ %60  = (y == 0.0)::Bool
└───        goto #22 if not %60
21 ─ %62  = Core.tuple(1.0, 0.0)::Core.Const((1.0, 0.0))
└───        return %62
22 ─ %64  = (a == 0.0)::Bool
└───        goto #24 if not %64
23 ─ %66  = Core.tuple(1.0, 0.0)::Core.Const((1.0, 0.0))
└───        return %66
24 ─ %68  = (b == 0.0)::Bool
└───        goto #26 if not %68
25 ─ %70  = Core.tuple(0.0, 1.0)::Core.Const((0.0, 1.0))
└───        return %70
26 ─ %72  = SpecialFunctions.eps()::Core.Const(2.220446049250313e-16)
│           (epps = SpecialFunctions.max(%72, 1.0e-15))
│    %74  = SpecialFunctions.max(a, b)::Float64
│    %75  = (0.001 * epps::Core.Const(1.0e-15))::Core.Const(1.0e-18)
│    %76  = (%74 < %75)::Bool
└───        goto #28 if not %76
27 ─ %78  = (a + b)::Float64
│    %79  = (b / %78)::Float64
│    %80  = (a + b)::Float64
│    %81  = (a / %80)::Float64
│    %82  = Core.tuple(%79, %81)::Tuple{Float64, Float64}
└───        return %82
28 ─        (ind = false)
│           (a0 = a)
│           (b0 = b)
│           (x0 = x)
│           (y0 = y)
│    %89  = SpecialFunctions.min(a0, b0)::Float64
│    %90  = (%89 > 1.0)::Bool
└───        goto #61 if not %90
29 ─ %92  = (a > b)::Bool
└───        goto #31 if not %92
30 ─ %94  = (a + b)::Float64
│    %95  = (%94 * y)::Float64
│           (@_17 = %95 - b)
└───        goto #32
31 ─ %98  = (a + b)::Float64
│    %99  = (%98 * x)::Float64
└───        (@_17 = a - %99)
32 ┄        (lambda = @_17)
│    %102 = (lambda < 0.0)::Bool
└───        goto #34 if not %102
33 ─        (ind = true)
│           (a0 = b)
│           (b0 = a)
│           (x0 = y)
│           (y0 = x)
└───        (lambda = SpecialFunctions.abs(lambda))
34 ┄ %110 = (b0 < 40.0)::Bool
└───        goto #37 if not %110
35 ─ %112 = (b0 * x0)::Float64
│    %113 = (%112 <= 0.7)::Bool
└───        goto #37 if not %113
36 ─        (p = SpecialFunctions.beta_inc_power_series(a0, b0, x0, epps::Core.Const(1.0e-15)))
│           (q = 1.0 - p)
└───        goto #58
37 ┄ %118 = (b0 < 40.0)::Bool
└───        goto #46 if not %118
38 ─        (n = SpecialFunctions.trunc(SpecialFunctions.Int, b0))
│           (b0 = b0 - n)
│    %122 = (b0 == 0.0)::Bool
└───        goto #40 if not %122
39 ─        (n = n - 1)
└───        (b0 = 1.0)
40 ┄        (p = SpecialFunctions.beta_inc_diff(b0, a0, y0, x0, n, epps::Core.Const(1.0e-15)))
│    %127 = (x0 <= 0.7)::Bool
└───        goto #42 if not %127
41 ─ %129 = p::Float64
│    %130 = SpecialFunctions.beta_inc_power_series(a0, b0, x0, epps::Core.Const(1.0e-15))::Float64
│           (p = %129 + %130)
│           (q = 1.0 - p)
└───        goto #45
42 ─ %134 = (a0 <= 15.0)::Bool
└───        goto #44 if not %134
43 ─        (n = 20)
│    %137 = p::Float64
│    %138 = SpecialFunctions.beta_inc_diff(a0, b0, x0, y0, n::Core.Const(20), epps::Core.Const(1.0e-15))::Float64
│           (p = %137 + %138)
└───        (a0 = a0 + n::Core.Const(20))
44 ┄ %141 = a0::Float64
│    %142 = b0::Float64
│    %143 = x0::Float64
│    %144 = y0::Float64
│    %145 = p::Float64
│    %146 = SpecialFunctions.eps()::Core.Const(2.220446049250313e-16)
│    %147 = (15.0 * %146)::Core.Const(3.3306690738754696e-15)
│           (p = SpecialFunctions.beta_inc_asymptotic_asymmetric(%141, %142, %143, %144, %145, %147))
└───        (q = 1.0 - p)
45 ┄        goto #58
46 ─ %151 = (a0 > b0)::Bool
└───        goto #53 if not %151
47 ─ %153 = (b0 <= 100.0)::Bool
└───        goto #49 if not %153
48 ─        goto #50
49 ─ %156 = lambda::Float64
│    %157 = (0.03 * b0)::Float64
│    %158 = (%156 > %157)::Bool
└───        goto #51 if not %158
50 ┄ %160 = a0::Float64
│    %161 = b0::Float64
│    %162 = x0::Float64
│    %163 = y0::Float64
│    %164 = lambda::Float64
│    %165 = SpecialFunctions.eps()::Core.Const(2.220446049250313e-16)
│    %166 = (15.0 * %165)::Core.Const(3.3306690738754696e-15)
│           (p = SpecialFunctions.beta_inc_cont_fraction(%160, %161, %162, %163, %164, %166))
│           (q = 1.0 - p)
└───        goto #52
51 ─ %170 = a0::Float64
│    %171 = b0::Float64
│    %172 = lambda::Float64
│    %173 = SpecialFunctions.eps()::Core.Const(2.220446049250313e-16)
│    %174 = (100.0 * %173)::Core.Const(2.220446049250313e-14)
│           (p = SpecialFunctions.beta_inc_asymptotic_symmetric(%170, %171, %172, %174))
└───        (q = 1.0 - p)
52 ┄        goto #58
53 ─ %178 = (a0 <= 100.0)::Bool
└───        goto #55 if not %178
54 ─        goto #56
55 ─ %181 = lambda::Float64
│    %182 = (0.03 * a0)::Float64
│    %183 = (%181 > %182)::Bool
└───        goto #57 if not %183
56 ┄ %185 = a0::Float64
│    %186 = b0::Float64
│    %187 = x0::Float64
│    %188 = y0::Float64
│    %189 = lambda::Float64
│    %190 = SpecialFunctions.eps()::Core.Const(2.220446049250313e-16)
│    %191 = (15.0 * %190)::Core.Const(3.3306690738754696e-15)
│           (p = SpecialFunctions.beta_inc_cont_fraction(%185, %186, %187, %188, %189, %191))
│           (q = 1.0 - p)
└───        goto #58
57 ─ %195 = a0::Float64
│    %196 = b0::Float64
│    %197 = lambda::Float64
│    %198 = SpecialFunctions.eps()::Core.Const(2.220446049250313e-16)
│    %199 = (100.0 * %198)::Core.Const(2.220446049250313e-14)
│           (p = SpecialFunctions.beta_inc_asymptotic_symmetric(%195, %196, %197, %199))
└───        (q = 1.0 - p)
58 ┄        goto #60 if not ind
59 ─ %203 = Core.tuple(q, p)::Tuple{Float64, Float64}
└───        return %203
60 ─ %205 = Core.tuple(p, q)::Tuple{Float64, Float64}
└───        return %205
61 ─ %207 = (x > 0.5)::Bool
└───        goto #63 if not %207
62 ─        (ind = true)
│           (a0 = b)
│           (b0 = a)
│           (y0 = x)
└───        (x0 = y)
63 ┄ %214 = b0::Float64
│    %215 = epps::Core.Const(1.0e-15)
│    %216 = (epps::Core.Const(1.0e-15) * a0)::Float64
│    %217 = SpecialFunctions.min(%215, %216)::Float64
│    %218 = (%214 < %217)::Bool
└───        goto #65 if not %218
64 ─        (p = SpecialFunctions.beta_inc_power_series2(a0, b0, x0, epps::Core.Const(1.0e-15)))
│           (q = 1.0 - p)
└───        goto #89
65 ─ %223 = a0::Float64
│    %224 = epps::Core.Const(1.0e-15)
│    %225 = (epps::Core.Const(1.0e-15) * b0)::Float64
│    %226 = SpecialFunctions.min(%224, %225)::Float64
│    %227 = (%223 < %226)::Bool
└───        goto #68 if not %227
66 ─ %229 = (b0 * x0)::Float64
│    %230 = (%229 <= 1.0)::Bool
└───        goto #68 if not %230
67 ─        (q = SpecialFunctions.beta_inc_power_series1(a0, b0, x0, epps::Core.Const(1.0e-15)))
│           (p = 1.0 - q)
└───        goto #89
68 ┄ %235 = SpecialFunctions.max(a0, b0)::Float64
│    %236 = (%235 > 1.0)::Bool
└───        goto #82 if not %236
69 ─ %238 = (b0 <= 1.0)::Bool
└───        goto #71 if not %238
70 ─        (p = SpecialFunctions.beta_inc_power_series(a0, b0, x0, epps::Core.Const(1.0e-15)))
│           (q = 1.0 - p)
└───        goto #81
71 ─ %243 = (x0 >= 0.3)::Bool
└───        goto #73 if not %243
72 ─        (q = SpecialFunctions.beta_inc_power_series(b0, a0, y0, epps::Core.Const(1.0e-15)))
│           (p = 1.0 - q)
└───        goto #81
73 ─ %248 = (x0 >= 0.1)::Bool
└───        goto #78 if not %248
74 ─ %250 = (b0 > 15.0)::Bool
└───        goto #76 if not %250
75 ─ %252 = b0::Float64
│    %253 = a0::Float64
│    %254 = y0::Float64
│    %255 = x0::Float64
│    %256 = q::Core.Const(0.0)
│    %257 = SpecialFunctions.eps()::Core.Const(2.220446049250313e-16)
│    %258 = (15.0 * %257)::Core.Const(3.3306690738754696e-15)
│           (q = SpecialFunctions.beta_inc_asymptotic_asymmetric(%252, %253, %254, %255, %256, %258))
│           (p = 1.0 - q)
└───        goto #77
76 ─        (n = 20)
│           (q = SpecialFunctions.beta_inc_diff(b0, a0, y0, x0, n::Core.Const(20), epps::Core.Const(1.0e-15)))
│           (b0 = b0 + n::Core.Const(20))
│    %265 = b0::Float64
│    %266 = a0::Float64
│    %267 = y0::Float64
│    %268 = x0::Float64
│    %269 = q::Float64
│    %270 = SpecialFunctions.eps()::Core.Const(2.220446049250313e-16)
│    %271 = (15.0 * %270)::Core.Const(3.3306690738754696e-15)
│           (q = SpecialFunctions.beta_inc_asymptotic_asymmetric(%265, %266, %267, %268, %269, %271))
└───        (p = 1.0 - q)
77 ┄        goto #81
78 ─ %275 = (x0 * b0)::Float64
│    %276 = (%275 ^ a0)::Float64
│    %277 = (%276 <= 0.7)::Bool
└───        goto #80 if not %277
79 ─        (p = SpecialFunctions.beta_inc_power_series(a0, b0, x0, epps::Core.Const(1.0e-15)))
│           (q = 1.0 - p)
└───        goto #81
80 ─        (n = 20)
│           (q = SpecialFunctions.beta_inc_diff(b0, a0, y0, x0, n::Core.Const(20), epps::Core.Const(1.0e-15)))
│           (b0 = b0 + n::Core.Const(20))
│    %285 = b0::Float64
│    %286 = a0::Float64
│    %287 = y0::Float64
│    %288 = x0::Float64
│    %289 = q::Float64
│    %290 = SpecialFunctions.eps()::Core.Const(2.220446049250313e-16)
│    %291 = (15.0 * %290)::Core.Const(3.3306690738754696e-15)
│           (q = SpecialFunctions.beta_inc_asymptotic_asymmetric(%285, %286, %287, %288, %289, %291))
└───        (p = 1.0 - q)
81 ┄        goto #89
82 ─ %295 = a0::Float64
│    %296 = SpecialFunctions.min(0.2, b0)::Float64
│    %297 = (%295 >= %296)::Bool
└───        goto #84 if not %297
83 ─        (p = SpecialFunctions.beta_inc_power_series(a0, b0, x0, epps::Core.Const(1.0e-15)))
│           (q = 1.0 - p)
└───        goto #89
84 ─ %302 = (x0 ^ a0)::Float64
│    %303 = (%302 <= 0.9)::Bool
└───        goto #86 if not %303
85 ─        (p = SpecialFunctions.beta_inc_power_series(a0, b0, x0, epps::Core.Const(1.0e-15)))
│           (q = 1.0 - p)
└───        goto #89
86 ─ %308 = (x0 >= 0.3)::Bool
└───        goto #88 if not %308
87 ─        (q = SpecialFunctions.beta_inc_power_series(b0, a0, y0, epps::Core.Const(1.0e-15)))
│           (p = 1.0 - q)
└───        goto #89
88 ─        (n = 20)
│           (q = SpecialFunctions.beta_inc_diff(b0, a0, y0, x0, n::Core.Const(20), epps::Core.Const(1.0e-15)))
│           (b0 = b0 + n::Core.Const(20))
│    %316 = b0::Float64
│    %317 = a0::Float64
│    %318 = y0::Float64
│    %319 = x0::Float64
│    %320 = q::Float64
│    %321 = SpecialFunctions.eps()::Core.Const(2.220446049250313e-16)
│    %322 = (15.0 * %321)::Core.Const(3.3306690738754696e-15)
│           (q = SpecialFunctions.beta_inc_asymptotic_asymmetric(%316, %317, %318, %319, %320, %322))
└───        (p = 1.0 - q)
89 ┄        goto #91 if not ind
90 ─ %326 = Core.tuple(q, p)::Tuple{Float64, Float64}
└───        return %326
91 ─ %328 = Core.tuple(p, q)::Tuple{Float64, Float64}
└───        return %328

In [79]:

epps = max(eps(), 1.0e-15)

Out[79]:

1.0e-15

In [81]:

@btime SpecialFunctions.beta_inc_power_series($a, $b, $x, $epps)

  130.464 ns (0 allocations: 0 bytes)

Out[81]:

1.1675552034918687e-104

In [83]:

y = 1 - x
lambda = a > b ? (a + b)*y - b : a - (a + b)*x

Out[83]:

89.999

In [89]:

SpecialFunctions.beta_inc_asymptotic_symmetric(b, a, lambda, 100.0*eps())

AssertionError: a >= 15.0

Stacktrace:
 [1] beta_inc_asymptotic_symmetric(a::Float64, b::Float64, lambda::Float64, epps::Float64)
   @ SpecialFunctions D:\.julia\packages\SpecialFunctions\hefUc\src\beta_inc.jl:263
 [2] top-level scope
   @ In[89]:1
 [3] eval
   @ .\boot.jl:368 [inlined]
 [4] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base .\loading.jl:1428

In [ ]: