{
"cells": [
{
"cell_type": "code",
"execution_count": 53,
"id": "ff8267c2-386a-4ba0-b4a3-ac5c49b1a9e4",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"ibetar_st (generic function with 1 method)"
]
},
"execution_count": 53,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using BenchmarkTools\n",
"using SpecialFunctions\n",
"\n",
"incbeta_I(x, a, b) = beta_inc(a, b, x)[1]\n",
"\n",
"function ibetar_st(x::T, a, b;\n",
" HUGE_VAL = T(Inf), FDeps = eps(T), rtol = eps(T), MAXITERS=200\n",
" ) where T<:AbstractFloat\n",
" a ≤ 0 && return HUGE_VAL\n",
" if b ≤ 0\n",
" x < 1 && return zero(T)\n",
" x == 1 && return one(T)\n",
" return HUGE_VAL\n",
" end\n",
" x > (a+1)/(a+b+2) && return 1 - ibetar_st(1 - x, b, a; HUGE_VAL, FDeps, rtol, MAXITERS)\n",
" x ≤ FDeps && return zero(T)\n",
" p1 = zero(T)\n",
" q1 = one(T)\n",
" p2 = exp(a*log(x) + b*log(1-x) + loggamma(a+b) - loggamma(a) - loggamma(b)) / a\n",
" q2 = one(T)\n",
" previous = HUGE_VAL\n",
" k = 0\n",
" while k < MAXITERS && !isapprox(p2, previous; rtol)\n",
" previous = p2\n",
" d = -(a+k)*(a+b+k)*x / ((a+2k)*(a+2k+1))\n",
" p1 = p1*d + p2\n",
" q1 = q1*d + q2\n",
" k += 1\n",
" d = k*(b-k)*x /((a+2k-1)*(a+2k))\n",
" p2 = p2*d + p1\n",
" q2 = q2*d + q1\n",
" if q2 ≤ FDeps\n",
" p2 = HUGE_VAL\n",
" else\n",
" p1 /= q2\n",
" q1 /= q2\n",
" p2 /= q2\n",
" q2 = one(T)\n",
" end\n",
" end\n",
" !isapprox(p2, previous; rtol) && error(\"ibetar did not converge.\")\n",
" p2\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 76,
"id": "5d6154b4-f1f6-4ab3-9ede-34266506abe8",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"incbeta_I(x, a, b) = 1.1675552034918687e-104\n",
"ibetar_st(x, a, b) = 1.1675552034918835e-104\n"
]
}
],
"source": [
"a, b = 100.0, 0.01\n",
"x = 0.1\n",
"@show incbeta_I(x, a, b)\n",
"@show ibetar_st(x, a, b);"
]
},
{
"cell_type": "code",
"execution_count": 77,
"id": "39cf8786-7a2b-494f-a4e0-94125a3983ff",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 143.853 ns (0 allocations: 0 bytes)\n",
" 149.333 ns (0 allocations: 0 bytes)\n"
]
}
],
"source": [
"@btime incbeta_I($x, $a, $b)\n",
"@btime ibetar_st($x, $a, $b);"
]
},
{
"cell_type": "code",
"execution_count": 78,
"id": "389d70f6-3437-4249-8ddd-ed19c49456fa",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"_beta_inc(a::Float64, b::Float64, x::Float64) in SpecialFunctions at D:\\.julia\\packages\\SpecialFunctions\\hefUc\\src\\beta_inc.jl:738"
],
"text/plain": [
"_beta_inc(a::Float64, b::Float64, x::Float64) in SpecialFunctions at D:\\.julia\\packages\\SpecialFunctions\\hefUc\\src\\beta_inc.jl:738"
]
},
"execution_count": 78,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@which SpecialFunctions._beta_inc(a, b, x)"
]
},
{
"cell_type": "code",
"execution_count": 91,
"id": "2b17e006-c5d0-41ab-8872-db50f4c15d70",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"MethodInstance for SpecialFunctions._beta_inc(::Float64, ::Float64, ::Float64, ::Float64)\n",
" from _beta_inc(a::Float64, b::Float64, x::Float64, y::Float64) in SpecialFunctions at D:\\.julia\\packages\\SpecialFunctions\\hefUc\\src\\beta_inc.jl:738\n",
"Arguments\n",
" #self#\u001b[36m::Core.Const(SpecialFunctions._beta_inc)\u001b[39m\n",
" a\u001b[36m::Float64\u001b[39m\n",
" b\u001b[36m::Float64\u001b[39m\n",
" x\u001b[36m::Float64\u001b[39m\n",
" y\u001b[36m::Float64\u001b[39m\n",
"Locals\n",
" p\u001b[36m::Float64\u001b[39m\n",
" q\u001b[36m::Float64\u001b[39m\n",
" b0\u001b[36m::Float64\u001b[39m\n",
" n\u001b[36m::Int64\u001b[39m\n",
" x0\u001b[36m::Float64\u001b[39m\n",
" y0\u001b[36m::Float64\u001b[39m\n",
" a0\u001b[36m::Float64\u001b[39m\n",
" ind\u001b[36m::Bool\u001b[39m\n",
" lambda\u001b[36m::Float64\u001b[39m\n",
" epps\u001b[36m::Float64\u001b[39m\n",
" z\u001b[36m::Float64\u001b[39m\n",
" @_17\u001b[36m::Float64\u001b[39m\n",
"Body\u001b[36m::Tuple{Float64, Float64}\u001b[39m\n",
"\u001b[90m1 ──\u001b[39m Core.NewvarNode(:(b0))\n",
"\u001b[90m│ \u001b[39m Core.NewvarNode(:(n))\n",
"\u001b[90m│ \u001b[39m Core.NewvarNode(:(x0))\n",
"\u001b[90m│ \u001b[39m Core.NewvarNode(:(y0))\n",
"\u001b[90m│ \u001b[39m Core.NewvarNode(:(a0))\n",
"\u001b[90m│ \u001b[39m Core.NewvarNode(:(ind))\n",
"\u001b[90m│ \u001b[39m Core.NewvarNode(:(lambda))\n",
"\u001b[90m│ \u001b[39m Core.NewvarNode(:(epps))\n",
"\u001b[90m│ \u001b[39m Core.NewvarNode(:(z))\n",
"\u001b[90m│ \u001b[39m (p = 0.0)\n",
"\u001b[90m│ \u001b[39m (q = 0.0)\n",
"\u001b[90m│ \u001b[39m %12 = (a < 0.0)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #3 if not %12\n",
"\u001b[90m2 ──\u001b[39m goto #4\n",
"\u001b[90m3 ──\u001b[39m %15 = (b < 0.0)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #5 if not %15\n",
"\u001b[90m4 ┄─\u001b[39m %17 = Core.tuple(a, b)\u001b[36m::Tuple{Float64, Float64}\u001b[39m\n",
"\u001b[90m│ \u001b[39m %18 = SpecialFunctions.DomainError(%17, \"a or b is negative\")\u001b[36m::Core.PartialStruct(DomainError, Any[Tuple{Float64, Float64}, String])\u001b[39m\n",
"\u001b[90m│ \u001b[39m SpecialFunctions.throw(%18)\n",
"\u001b[90m└───\u001b[39m Core.Const(:(return %19))\n",
"\u001b[90m5 ┄─\u001b[39m %21 = (a == 0.0)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #8 if not %21\n",
"\u001b[90m6 ──\u001b[39m %23 = (b == 0.0)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #8 if not %23\n",
"\u001b[90m7 ──\u001b[39m %25 = Core.tuple(a, b)\u001b[36m::Tuple{Float64, Float64}\u001b[39m\n",
"\u001b[90m│ \u001b[39m %26 = SpecialFunctions.DomainError(%25, \"a and b are 0.0\")\u001b[36m::Core.PartialStruct(DomainError, Any[Tuple{Float64, Float64}, String])\u001b[39m\n",
"\u001b[90m│ \u001b[39m SpecialFunctions.throw(%26)\n",
"\u001b[90m└───\u001b[39m Core.Const(:(return %27))\n",
"\u001b[90m8 ┄─\u001b[39m %29 = (x < 0.0)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #10 if not %29\n",
"\u001b[90m9 ──\u001b[39m goto #11\n",
"\u001b[90m10 ─\u001b[39m %32 = (x > 1.0)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #12 if not %32\n",
"\u001b[90m11 ┄\u001b[39m %34 = SpecialFunctions.DomainError(x, \"x < 0 or x > 1\")\u001b[36m::Core.PartialStruct(DomainError, Any[Float64, String])\u001b[39m\n",
"\u001b[90m│ \u001b[39m SpecialFunctions.throw(%34)\n",
"\u001b[90m└───\u001b[39m Core.Const(:(return %35))\n",
"\u001b[90m12 ┄\u001b[39m %37 = (y < 0.0)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #14 if not %37\n",
"\u001b[90m13 ─\u001b[39m goto #15\n",
"\u001b[90m14 ─\u001b[39m %40 = (y > 1.0)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #16 if not %40\n",
"\u001b[90m15 ┄\u001b[39m %42 = SpecialFunctions.DomainError(y, \"y < 0 or y > 1\")\u001b[36m::Core.PartialStruct(DomainError, Any[Float64, String])\u001b[39m\n",
"\u001b[90m│ \u001b[39m SpecialFunctions.throw(%42)\n",
"\u001b[90m└───\u001b[39m Core.Const(:(return %43))\n",
"\u001b[90m16 ┄\u001b[39m %45 = (x + y)\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m (z = %45 - 1.0)\n",
"\u001b[90m│ \u001b[39m %47 = SpecialFunctions.abs(z)\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %48 = SpecialFunctions.eps()\u001b[36m::Core.Const(2.220446049250313e-16)\u001b[39m\n",
"\u001b[90m│ \u001b[39m %49 = (3.0 * %48)\u001b[36m::Core.Const(6.661338147750939e-16)\u001b[39m\n",
"\u001b[90m│ \u001b[39m %50 = (%47 > %49)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #18 if not %50\n",
"\u001b[90m17 ─\u001b[39m %52 = Core.tuple(x, y)\u001b[36m::Tuple{Float64, Float64}\u001b[39m\n",
"\u001b[90m│ \u001b[39m %53 = SpecialFunctions.DomainError(%52, \"x + y != 1.0\")\u001b[36m::Core.PartialStruct(DomainError, Any[Tuple{Float64, Float64}, String])\u001b[39m\n",
"\u001b[90m│ \u001b[39m SpecialFunctions.throw(%53)\n",
"\u001b[90m└───\u001b[39m Core.Const(:(return %54))\n",
"\u001b[90m18 ┄\u001b[39m %56 = (x == 0.0)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #20 if not %56\n",
"\u001b[90m19 ─\u001b[39m %58 = Core.tuple(0.0, 1.0)\u001b[36m::Core.Const((0.0, 1.0))\u001b[39m\n",
"\u001b[90m└───\u001b[39m return %58\n",
"\u001b[90m20 ─\u001b[39m %60 = (y == 0.0)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #22 if not %60\n",
"\u001b[90m21 ─\u001b[39m %62 = Core.tuple(1.0, 0.0)\u001b[36m::Core.Const((1.0, 0.0))\u001b[39m\n",
"\u001b[90m└───\u001b[39m return %62\n",
"\u001b[90m22 ─\u001b[39m %64 = (a == 0.0)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #24 if not %64\n",
"\u001b[90m23 ─\u001b[39m %66 = Core.tuple(1.0, 0.0)\u001b[36m::Core.Const((1.0, 0.0))\u001b[39m\n",
"\u001b[90m└───\u001b[39m return %66\n",
"\u001b[90m24 ─\u001b[39m %68 = (b == 0.0)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #26 if not %68\n",
"\u001b[90m25 ─\u001b[39m %70 = Core.tuple(0.0, 1.0)\u001b[36m::Core.Const((0.0, 1.0))\u001b[39m\n",
"\u001b[90m└───\u001b[39m return %70\n",
"\u001b[90m26 ─\u001b[39m %72 = SpecialFunctions.eps()\u001b[36m::Core.Const(2.220446049250313e-16)\u001b[39m\n",
"\u001b[90m│ \u001b[39m (epps = SpecialFunctions.max(%72, 1.0e-15))\n",
"\u001b[90m│ \u001b[39m %74 = SpecialFunctions.max(a, b)\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %75 = (0.001 * epps::Core.Const(1.0e-15))\u001b[36m::Core.Const(1.0e-18)\u001b[39m\n",
"\u001b[90m│ \u001b[39m %76 = (%74 < %75)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #28 if not %76\n",
"\u001b[90m27 ─\u001b[39m %78 = (a + b)\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %79 = (b / %78)\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %80 = (a + b)\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %81 = (a / %80)\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %82 = Core.tuple(%79, %81)\u001b[36m::Tuple{Float64, Float64}\u001b[39m\n",
"\u001b[90m└───\u001b[39m return %82\n",
"\u001b[90m28 ─\u001b[39m (ind = false)\n",
"\u001b[90m│ \u001b[39m (a0 = a)\n",
"\u001b[90m│ \u001b[39m (b0 = b)\n",
"\u001b[90m│ \u001b[39m (x0 = x)\n",
"\u001b[90m│ \u001b[39m (y0 = y)\n",
"\u001b[90m│ \u001b[39m %89 = SpecialFunctions.min(a0, b0)\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %90 = (%89 > 1.0)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #61 if not %90\n",
"\u001b[90m29 ─\u001b[39m %92 = (a > b)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #31 if not %92\n",
"\u001b[90m30 ─\u001b[39m %94 = (a + b)\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %95 = (%94 * y)\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m (@_17 = %95 - b)\n",
"\u001b[90m└───\u001b[39m goto #32\n",
"\u001b[90m31 ─\u001b[39m %98 = (a + b)\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %99 = (%98 * x)\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m└───\u001b[39m (@_17 = a - %99)\n",
"\u001b[90m32 ┄\u001b[39m (lambda = @_17)\n",
"\u001b[90m│ \u001b[39m %102 = (lambda < 0.0)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #34 if not %102\n",
"\u001b[90m33 ─\u001b[39m (ind = true)\n",
"\u001b[90m│ \u001b[39m (a0 = b)\n",
"\u001b[90m│ \u001b[39m (b0 = a)\n",
"\u001b[90m│ \u001b[39m (x0 = y)\n",
"\u001b[90m│ \u001b[39m (y0 = x)\n",
"\u001b[90m└───\u001b[39m (lambda = SpecialFunctions.abs(lambda))\n",
"\u001b[90m34 ┄\u001b[39m %110 = (b0 < 40.0)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #37 if not %110\n",
"\u001b[90m35 ─\u001b[39m %112 = (b0 * x0)\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %113 = (%112 <= 0.7)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #37 if not %113\n",
"\u001b[90m36 ─\u001b[39m (p = SpecialFunctions.beta_inc_power_series(a0, b0, x0, epps::Core.Const(1.0e-15)))\n",
"\u001b[90m│ \u001b[39m (q = 1.0 - p)\n",
"\u001b[90m└───\u001b[39m goto #58\n",
"\u001b[90m37 ┄\u001b[39m %118 = (b0 < 40.0)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #46 if not %118\n",
"\u001b[90m38 ─\u001b[39m (n = SpecialFunctions.trunc(SpecialFunctions.Int, b0))\n",
"\u001b[90m│ \u001b[39m (b0 = b0 - n)\n",
"\u001b[90m│ \u001b[39m %122 = (b0 == 0.0)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #40 if not %122\n",
"\u001b[90m39 ─\u001b[39m (n = n - 1)\n",
"\u001b[90m└───\u001b[39m (b0 = 1.0)\n",
"\u001b[90m40 ┄\u001b[39m (p = SpecialFunctions.beta_inc_diff(b0, a0, y0, x0, n, epps::Core.Const(1.0e-15)))\n",
"\u001b[90m│ \u001b[39m %127 = (x0 <= 0.7)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #42 if not %127\n",
"\u001b[90m41 ─\u001b[39m %129 = p\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %130 = SpecialFunctions.beta_inc_power_series(a0, b0, x0, epps::Core.Const(1.0e-15))\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m (p = %129 + %130)\n",
"\u001b[90m│ \u001b[39m (q = 1.0 - p)\n",
"\u001b[90m└───\u001b[39m goto #45\n",
"\u001b[90m42 ─\u001b[39m %134 = (a0 <= 15.0)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #44 if not %134\n",
"\u001b[90m43 ─\u001b[39m (n = 20)\n",
"\u001b[90m│ \u001b[39m %137 = p\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %138 = SpecialFunctions.beta_inc_diff(a0, b0, x0, y0, n::Core.Const(20), epps::Core.Const(1.0e-15))\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m (p = %137 + %138)\n",
"\u001b[90m└───\u001b[39m (a0 = a0 + n::Core.Const(20))\n",
"\u001b[90m44 ┄\u001b[39m %141 = a0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %142 = b0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %143 = x0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %144 = y0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %145 = p\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %146 = SpecialFunctions.eps()\u001b[36m::Core.Const(2.220446049250313e-16)\u001b[39m\n",
"\u001b[90m│ \u001b[39m %147 = (15.0 * %146)\u001b[36m::Core.Const(3.3306690738754696e-15)\u001b[39m\n",
"\u001b[90m│ \u001b[39m (p = SpecialFunctions.beta_inc_asymptotic_asymmetric(%141, %142, %143, %144, %145, %147))\n",
"\u001b[90m└───\u001b[39m (q = 1.0 - p)\n",
"\u001b[90m45 ┄\u001b[39m goto #58\n",
"\u001b[90m46 ─\u001b[39m %151 = (a0 > b0)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #53 if not %151\n",
"\u001b[90m47 ─\u001b[39m %153 = (b0 <= 100.0)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #49 if not %153\n",
"\u001b[90m48 ─\u001b[39m goto #50\n",
"\u001b[90m49 ─\u001b[39m %156 = lambda\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %157 = (0.03 * b0)\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %158 = (%156 > %157)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #51 if not %158\n",
"\u001b[90m50 ┄\u001b[39m %160 = a0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %161 = b0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %162 = x0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %163 = y0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %164 = lambda\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %165 = SpecialFunctions.eps()\u001b[36m::Core.Const(2.220446049250313e-16)\u001b[39m\n",
"\u001b[90m│ \u001b[39m %166 = (15.0 * %165)\u001b[36m::Core.Const(3.3306690738754696e-15)\u001b[39m\n",
"\u001b[90m│ \u001b[39m (p = SpecialFunctions.beta_inc_cont_fraction(%160, %161, %162, %163, %164, %166))\n",
"\u001b[90m│ \u001b[39m (q = 1.0 - p)\n",
"\u001b[90m└───\u001b[39m goto #52\n",
"\u001b[90m51 ─\u001b[39m %170 = a0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %171 = b0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %172 = lambda\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %173 = SpecialFunctions.eps()\u001b[36m::Core.Const(2.220446049250313e-16)\u001b[39m\n",
"\u001b[90m│ \u001b[39m %174 = (100.0 * %173)\u001b[36m::Core.Const(2.220446049250313e-14)\u001b[39m\n",
"\u001b[90m│ \u001b[39m (p = SpecialFunctions.beta_inc_asymptotic_symmetric(%170, %171, %172, %174))\n",
"\u001b[90m└───\u001b[39m (q = 1.0 - p)\n",
"\u001b[90m52 ┄\u001b[39m goto #58\n",
"\u001b[90m53 ─\u001b[39m %178 = (a0 <= 100.0)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #55 if not %178\n",
"\u001b[90m54 ─\u001b[39m goto #56\n",
"\u001b[90m55 ─\u001b[39m %181 = lambda\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %182 = (0.03 * a0)\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %183 = (%181 > %182)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #57 if not %183\n",
"\u001b[90m56 ┄\u001b[39m %185 = a0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %186 = b0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %187 = x0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %188 = y0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %189 = lambda\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %190 = SpecialFunctions.eps()\u001b[36m::Core.Const(2.220446049250313e-16)\u001b[39m\n",
"\u001b[90m│ \u001b[39m %191 = (15.0 * %190)\u001b[36m::Core.Const(3.3306690738754696e-15)\u001b[39m\n",
"\u001b[90m│ \u001b[39m (p = SpecialFunctions.beta_inc_cont_fraction(%185, %186, %187, %188, %189, %191))\n",
"\u001b[90m│ \u001b[39m (q = 1.0 - p)\n",
"\u001b[90m└───\u001b[39m goto #58\n",
"\u001b[90m57 ─\u001b[39m %195 = a0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %196 = b0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %197 = lambda\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %198 = SpecialFunctions.eps()\u001b[36m::Core.Const(2.220446049250313e-16)\u001b[39m\n",
"\u001b[90m│ \u001b[39m %199 = (100.0 * %198)\u001b[36m::Core.Const(2.220446049250313e-14)\u001b[39m\n",
"\u001b[90m│ \u001b[39m (p = SpecialFunctions.beta_inc_asymptotic_symmetric(%195, %196, %197, %199))\n",
"\u001b[90m└───\u001b[39m (q = 1.0 - p)\n",
"\u001b[90m58 ┄\u001b[39m goto #60 if not ind\n",
"\u001b[90m59 ─\u001b[39m %203 = Core.tuple(q, p)\u001b[36m::Tuple{Float64, Float64}\u001b[39m\n",
"\u001b[90m└───\u001b[39m return %203\n",
"\u001b[90m60 ─\u001b[39m %205 = Core.tuple(p, q)\u001b[36m::Tuple{Float64, Float64}\u001b[39m\n",
"\u001b[90m└───\u001b[39m return %205\n",
"\u001b[90m61 ─\u001b[39m %207 = (x > 0.5)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #63 if not %207\n",
"\u001b[90m62 ─\u001b[39m (ind = true)\n",
"\u001b[90m│ \u001b[39m (a0 = b)\n",
"\u001b[90m│ \u001b[39m (b0 = a)\n",
"\u001b[90m│ \u001b[39m (y0 = x)\n",
"\u001b[90m└───\u001b[39m (x0 = y)\n",
"\u001b[90m63 ┄\u001b[39m %214 = b0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %215 = epps\u001b[36m::Core.Const(1.0e-15)\u001b[39m\n",
"\u001b[90m│ \u001b[39m %216 = (epps::Core.Const(1.0e-15) * a0)\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %217 = SpecialFunctions.min(%215, %216)\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %218 = (%214 < %217)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #65 if not %218\n",
"\u001b[90m64 ─\u001b[39m (p = SpecialFunctions.beta_inc_power_series2(a0, b0, x0, epps::Core.Const(1.0e-15)))\n",
"\u001b[90m│ \u001b[39m (q = 1.0 - p)\n",
"\u001b[90m└───\u001b[39m goto #89\n",
"\u001b[90m65 ─\u001b[39m %223 = a0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %224 = epps\u001b[36m::Core.Const(1.0e-15)\u001b[39m\n",
"\u001b[90m│ \u001b[39m %225 = (epps::Core.Const(1.0e-15) * b0)\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %226 = SpecialFunctions.min(%224, %225)\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %227 = (%223 < %226)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #68 if not %227\n",
"\u001b[90m66 ─\u001b[39m %229 = (b0 * x0)\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %230 = (%229 <= 1.0)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #68 if not %230\n",
"\u001b[90m67 ─\u001b[39m (q = SpecialFunctions.beta_inc_power_series1(a0, b0, x0, epps::Core.Const(1.0e-15)))\n",
"\u001b[90m│ \u001b[39m (p = 1.0 - q)\n",
"\u001b[90m└───\u001b[39m goto #89\n",
"\u001b[90m68 ┄\u001b[39m %235 = SpecialFunctions.max(a0, b0)\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %236 = (%235 > 1.0)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #82 if not %236\n",
"\u001b[90m69 ─\u001b[39m %238 = (b0 <= 1.0)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #71 if not %238\n",
"\u001b[90m70 ─\u001b[39m (p = SpecialFunctions.beta_inc_power_series(a0, b0, x0, epps::Core.Const(1.0e-15)))\n",
"\u001b[90m│ \u001b[39m (q = 1.0 - p)\n",
"\u001b[90m└───\u001b[39m goto #81\n",
"\u001b[90m71 ─\u001b[39m %243 = (x0 >= 0.3)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #73 if not %243\n",
"\u001b[90m72 ─\u001b[39m (q = SpecialFunctions.beta_inc_power_series(b0, a0, y0, epps::Core.Const(1.0e-15)))\n",
"\u001b[90m│ \u001b[39m (p = 1.0 - q)\n",
"\u001b[90m└───\u001b[39m goto #81\n",
"\u001b[90m73 ─\u001b[39m %248 = (x0 >= 0.1)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #78 if not %248\n",
"\u001b[90m74 ─\u001b[39m %250 = (b0 > 15.0)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #76 if not %250\n",
"\u001b[90m75 ─\u001b[39m %252 = b0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %253 = a0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %254 = y0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %255 = x0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %256 = q\u001b[36m::Core.Const(0.0)\u001b[39m\n",
"\u001b[90m│ \u001b[39m %257 = SpecialFunctions.eps()\u001b[36m::Core.Const(2.220446049250313e-16)\u001b[39m\n",
"\u001b[90m│ \u001b[39m %258 = (15.0 * %257)\u001b[36m::Core.Const(3.3306690738754696e-15)\u001b[39m\n",
"\u001b[90m│ \u001b[39m (q = SpecialFunctions.beta_inc_asymptotic_asymmetric(%252, %253, %254, %255, %256, %258))\n",
"\u001b[90m│ \u001b[39m (p = 1.0 - q)\n",
"\u001b[90m└───\u001b[39m goto #77\n",
"\u001b[90m76 ─\u001b[39m (n = 20)\n",
"\u001b[90m│ \u001b[39m (q = SpecialFunctions.beta_inc_diff(b0, a0, y0, x0, n::Core.Const(20), epps::Core.Const(1.0e-15)))\n",
"\u001b[90m│ \u001b[39m (b0 = b0 + n::Core.Const(20))\n",
"\u001b[90m│ \u001b[39m %265 = b0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %266 = a0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %267 = y0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %268 = x0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %269 = q\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %270 = SpecialFunctions.eps()\u001b[36m::Core.Const(2.220446049250313e-16)\u001b[39m\n",
"\u001b[90m│ \u001b[39m %271 = (15.0 * %270)\u001b[36m::Core.Const(3.3306690738754696e-15)\u001b[39m\n",
"\u001b[90m│ \u001b[39m (q = SpecialFunctions.beta_inc_asymptotic_asymmetric(%265, %266, %267, %268, %269, %271))\n",
"\u001b[90m└───\u001b[39m (p = 1.0 - q)\n",
"\u001b[90m77 ┄\u001b[39m goto #81\n",
"\u001b[90m78 ─\u001b[39m %275 = (x0 * b0)\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %276 = (%275 ^ a0)\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %277 = (%276 <= 0.7)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #80 if not %277\n",
"\u001b[90m79 ─\u001b[39m (p = SpecialFunctions.beta_inc_power_series(a0, b0, x0, epps::Core.Const(1.0e-15)))\n",
"\u001b[90m│ \u001b[39m (q = 1.0 - p)\n",
"\u001b[90m└───\u001b[39m goto #81\n",
"\u001b[90m80 ─\u001b[39m (n = 20)\n",
"\u001b[90m│ \u001b[39m (q = SpecialFunctions.beta_inc_diff(b0, a0, y0, x0, n::Core.Const(20), epps::Core.Const(1.0e-15)))\n",
"\u001b[90m│ \u001b[39m (b0 = b0 + n::Core.Const(20))\n",
"\u001b[90m│ \u001b[39m %285 = b0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %286 = a0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %287 = y0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %288 = x0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %289 = q\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %290 = SpecialFunctions.eps()\u001b[36m::Core.Const(2.220446049250313e-16)\u001b[39m\n",
"\u001b[90m│ \u001b[39m %291 = (15.0 * %290)\u001b[36m::Core.Const(3.3306690738754696e-15)\u001b[39m\n",
"\u001b[90m│ \u001b[39m (q = SpecialFunctions.beta_inc_asymptotic_asymmetric(%285, %286, %287, %288, %289, %291))\n",
"\u001b[90m└───\u001b[39m (p = 1.0 - q)\n",
"\u001b[90m81 ┄\u001b[39m goto #89\n",
"\u001b[90m82 ─\u001b[39m %295 = a0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %296 = SpecialFunctions.min(0.2, b0)\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %297 = (%295 >= %296)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #84 if not %297\n",
"\u001b[90m83 ─\u001b[39m (p = SpecialFunctions.beta_inc_power_series(a0, b0, x0, epps::Core.Const(1.0e-15)))\n",
"\u001b[90m│ \u001b[39m (q = 1.0 - p)\n",
"\u001b[90m└───\u001b[39m goto #89\n",
"\u001b[90m84 ─\u001b[39m %302 = (x0 ^ a0)\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %303 = (%302 <= 0.9)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #86 if not %303\n",
"\u001b[90m85 ─\u001b[39m (p = SpecialFunctions.beta_inc_power_series(a0, b0, x0, epps::Core.Const(1.0e-15)))\n",
"\u001b[90m│ \u001b[39m (q = 1.0 - p)\n",
"\u001b[90m└───\u001b[39m goto #89\n",
"\u001b[90m86 ─\u001b[39m %308 = (x0 >= 0.3)\u001b[36m::Bool\u001b[39m\n",
"\u001b[90m└───\u001b[39m goto #88 if not %308\n",
"\u001b[90m87 ─\u001b[39m (q = SpecialFunctions.beta_inc_power_series(b0, a0, y0, epps::Core.Const(1.0e-15)))\n",
"\u001b[90m│ \u001b[39m (p = 1.0 - q)\n",
"\u001b[90m└───\u001b[39m goto #89\n",
"\u001b[90m88 ─\u001b[39m (n = 20)\n",
"\u001b[90m│ \u001b[39m (q = SpecialFunctions.beta_inc_diff(b0, a0, y0, x0, n::Core.Const(20), epps::Core.Const(1.0e-15)))\n",
"\u001b[90m│ \u001b[39m (b0 = b0 + n::Core.Const(20))\n",
"\u001b[90m│ \u001b[39m %316 = b0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %317 = a0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %318 = y0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %319 = x0\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %320 = q\u001b[36m::Float64\u001b[39m\n",
"\u001b[90m│ \u001b[39m %321 = SpecialFunctions.eps()\u001b[36m::Core.Const(2.220446049250313e-16)\u001b[39m\n",
"\u001b[90m│ \u001b[39m %322 = (15.0 * %321)\u001b[36m::Core.Const(3.3306690738754696e-15)\u001b[39m\n",
"\u001b[90m│ \u001b[39m (q = SpecialFunctions.beta_inc_asymptotic_asymmetric(%316, %317, %318, %319, %320, %322))\n",
"\u001b[90m└───\u001b[39m (p = 1.0 - q)\n",
"\u001b[90m89 ┄\u001b[39m goto #91 if not ind\n",
"\u001b[90m90 ─\u001b[39m %326 = Core.tuple(q, p)\u001b[36m::Tuple{Float64, Float64}\u001b[39m\n",
"\u001b[90m└───\u001b[39m return %326\n",
"\u001b[90m91 ─\u001b[39m %328 = Core.tuple(p, q)\u001b[36m::Tuple{Float64, Float64}\u001b[39m\n",
"\u001b[90m└───\u001b[39m return %328\n",
"\n"
]
}
],
"source": [
"@code_warntype SpecialFunctions._beta_inc(a, b, x, 1-x)"
]
},
{
"cell_type": "code",
"execution_count": 79,
"id": "713939dc-f203-4a2b-b017-b96e95a8afe3",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.0e-15"
]
},
"execution_count": 79,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"epps = max(eps(), 1.0e-15)"
]
},
{
"cell_type": "code",
"execution_count": 81,
"id": "a5c58dfd-9aa1-4c72-bd38-3fd2772ba9e6",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 130.464 ns (0 allocations: 0 bytes)\n"
]
},
{
"data": {
"text/plain": [
"1.1675552034918687e-104"
]
},
"execution_count": 81,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@btime SpecialFunctions.beta_inc_power_series($a, $b, $x, $epps)"
]
},
{
"cell_type": "code",
"execution_count": 83,
"id": "859b1068-59f5-4d8d-a880-66e55910102e",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"89.999"
]
},
"execution_count": 83,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"y = 1 - x\n",
"lambda = a > b ? (a + b)*y - b : a - (a + b)*x"
]
},
{
"cell_type": "code",
"execution_count": 89,
"id": "ddad4115-0517-4dd0-a1a6-3302240fa3c3",
"metadata": {},
"outputs": [
{
"ename": "LoadError",
"evalue": "AssertionError: a >= 15.0",
"output_type": "error",
"traceback": [
"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"
]
}
],
"source": [
"SpecialFunctions.beta_inc_asymptotic_symmetric(b, a, lambda, 100.0*eps())"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "326da83d-8dc3-4ea3-a0b4-fe5d6e405ed1",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"jupytext": {
"formats": "ipynb,jl:hydrogen"
},
"kernelspec": {
"display_name": "Julia 1.8.1",
"language": "julia",
"name": "julia-1.8"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.8.1"
}
},
"nbformat": 4,
"nbformat_minor": 5
}