From 2208822aff34b15b5d9e5d93e470e922d0de54f6 Mon Sep 17 00:00:00 2001 From: Johanni Brea Date: Fri, 24 Jan 2025 10:31:36 +0100 Subject: [PATCH 1/6] add owent --- Project.toml | 2 + src/SpecialFunctions.jl | 5 ++ src/owent.jl | 139 ++++++++++++++++++++++++++++++++++++++++ test/owent.jl | 104 ++++++++++++++++++++++++++++++ test/runtests.jl | 3 +- 5 files changed, 252 insertions(+), 1 deletion(-) create mode 100644 src/owent.jl create mode 100644 test/owent.jl diff --git a/Project.toml b/Project.toml index 45862755..e735732f 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,7 @@ version = "2.5.0" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" IrrationalConstants = "92d709cd-6900-40b7-9082-c6be49f344b6" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688" OpenLibm_jll = "05823500-19ac-5b8b-9628-191a04bc5112" OpenSpecFun_jll = "efe28fd5-8261-553b-a9e1-b2916fc3738e" @@ -19,6 +20,7 @@ SpecialFunctionsChainRulesCoreExt = "ChainRulesCore" ChainRulesCore = "0.9.44, 0.10, 1" ChainRulesTestUtils = "0.6.8, 0.7, 1" IrrationalConstants = "0.1, 0.2" +LinearAlgebra = "1.11.0" LogExpFunctions = "0.3.2" OpenLibm_jll = "0.7, 0.8" OpenSpecFun_jll = "0.5" diff --git a/src/SpecialFunctions.jl b/src/SpecialFunctions.jl index 0d296cc0..74fa88ea 100644 --- a/src/SpecialFunctions.jl +++ b/src/SpecialFunctions.jl @@ -1,5 +1,7 @@ module SpecialFunctions +import LinearAlgebra: dot + using IrrationalConstants: twoπ, halfπ, @@ -7,6 +9,7 @@ using IrrationalConstants: sqrt2π, invπ, inv2π, + inv4π, invsqrt2, invsqrt2π, logtwo, @@ -58,6 +61,7 @@ export logerfcx, faddeeva, eta, + owent, # Gamma functions gamma, @@ -103,6 +107,7 @@ include("gamma.jl") include("gamma_inc.jl") include("betanc.jl") include("beta_inc.jl") +include("owent.jl") if !isdefined(Base, :get_extension) include("../ext/SpecialFunctionsChainRulesCoreExt.jl") end diff --git a/src/owent.jl b/src/owent.jl new file mode 100644 index 00000000..943e28b1 --- /dev/null +++ b/src/owent.jl @@ -0,0 +1,139 @@ +# Owen's T Function +# Written by Andy Gough; August 2021 (see https://github.com/JuliaStats/StatsFuns.jl/issues/99#issuecomment-1124581689) +# Edited by Johanni Brea to make type stable; January 2025 +# Rev 1.09 +# MIT License +# +# dependencies +# IrrationalConstants +# SpecialFunctions +# LinearAlgebra +# +# HISTORY +# In the past 20 or so years, most implementations of Owen's T function have followed the algorithms given in "Fast and accurate Calculation of Owen's +# T-Function", by M. Patefield and D. Tandy, Journal of Statistical Software, 5 (5), 1 - 25 (2000) +# +# Six algorithms were given, and which is was used depends on the values of (h,a) +# +# T1: first m terms of series expansion of Owen (1956) +# T2: approximates 1/(1+x^2) by power series expansion up to order 2m +# T3: approximates 1/(1+x^2) by chebyshev polynomials of degree 2m in x +# T4: new expression for zi from T2 +# T5: Gauss 2m-point quadrature; 30 figures accuracy with m=48 (p. 18) +# T6: For when a is very close to 1, use formula derived from T(h,1) = 1/2 Φ(h)[1-Φ(h)] +# +# They developed code for these algorithms on a DEC VAX 750. The VAX 750 came out in 1980 and had a processor clock speed of 3.125 MHz. +# +# The reason for 6 algorithms was to speed up the function when possible, with T1 being faster than T2, T2 faster than T3, etc. +# +# THIS FUNCTION +# A native Julia implementation, based on the equations in the paper. The FORTRAN source code was not analyzed, translated, or used. This is a new +# implementation that takes advantages of Julia's unique capabilities (and those of modern computers). +# +# T1 through T4 are not implemented. Instead, if a < 0.999999, T5 is used to calculate Owen's T (using 48 point Gauss-Legendre quadrature) +# For 0.999999 < a < 1.0, T6 is implemented. +# +# REFERENCES +# [1] "Fast and accurate Calculation of Owen's T-Function", by M. Patefield and D. Tandy, Journal of Statistical Software, 5 (5), 1 - 25 (2000) +# [2] "Tables for Computing Bivariate Normal Probabilities", by Donald P. Owen, The Annals of Mathematical Statistics, Vol. 27, No. 4 (Dec 1956), pp. 1075-1090 +# +# Partial Derivatives (FYI) +# D[owent[x,a],x] = -exp(-0.5*x^2)*erf(a*x/sqrt2)/(2*sqrt2π) +# D[owent[x,a],a] = exp(-0.5*(1+a^2)*(x^2))/((1+a^2)*2π) +# +@doc raw""" +owent(h,a) : Returns the value of Owen's T function for (h,a) + +Owen's T function: +```math +T(h,a) = \frac{1}{2\pi } \int_{0}^{a} \frac{e^{-\frac{1}{2}h^2(1+x^2)}}{1+x^2}dx\quad(-\infty < h,a < +\infty) +``` + +For *h* and *a* > 0, *T(h,a)* gives the volume of the uncorrelated bivariate normal distribution with zero mean and unit variance +over the area from *y = ax* and *y = 0* and to the right of *x = h*. + +EXAMPLE: +``` +julia> owent(0.0625, 0.025) +0.003970281304296922 +``` + +Worst case accuracy is about 2e-16. +""" +function owent(h::T, a::T) where {T <: Real} + + invsqrt2_T = T(invsqrt2) + inv2π_T = T(inv2π) + + #********************* + # shortcut evaluations + #********************* + + if h < 0 + return owent(abs(h),a) + end + + if h == 0 + return atan(a)*inv2π_T + end + + if a < 0 + return -owent(h,abs(a)) + end + + if a == 0 + return zero(a) + end + + if a == 1 + return T(0.125)*erfc(-h*invsqrt2_T)*erfc(h*invsqrt2_T) + end + + if a == Inf + return T(0.25)*erfc(sqrt(h^2)*invsqrt2_T) + end + + # below reduces the range from -inf < h,a < +inf to h ≥ 0, 0 ≤ a ≤ 1 + if a > 1 + return T(0.25)*(erfc(-h*invsqrt2_T) + erfc(-a*h*invsqrt2_T)) - T(0.25)*erfc(-h*invsqrt2_T)*erfc(-a*h*invsqrt2_T) - owent(a*h,one(a)/a) + end + + # calculate Owen's T + + if a ≤ T(0.999999) + + + x = (-T(0.9987710072524261), -T(0.9935301722663508), -T(0.9841245837228269), -T(0.9705915925462473), -T(0.9529877031604309), -T(0.9313866907065543), -T(0.9058791367155696) + , -T(0.8765720202742479), -T(0.8435882616243935), -T(0.8070662040294426), -T(0.7671590325157404), -T(0.7240341309238146), -T(0.6778723796326639), -T(0.6288673967765136) + , -T(0.5772247260839727), -T(0.5231609747222331), -T(0.4669029047509584), -T(0.4086864819907167), -T(0.34875588629216075), -T(0.28736248735545555), -T(0.22476379039468905) + , -T(0.1612223560688917), -T(0.0970046992094627), -T(0.03238017096286937), T(0.03238017096286937), T(0.0970046992094627), T(0.1612223560688917), T(0.22476379039468905) + , T(0.28736248735545555), T(0.34875588629216075), T(0.4086864819907167), T(0.4669029047509584), T(0.5231609747222331), T(0.5772247260839727), T(0.6288673967765136) + , T(0.6778723796326639), T(0.7240341309238146), T(0.7671590325157404), T(0.8070662040294426), T(0.8435882616243935), T(0.8765720202742479), T(0.9058791367155696) + , T(0.9313866907065543), T(0.9529877031604309), T(0.9705915925462473), T(0.9841245837228269), T(0.9935301722663508), T(0.9987710072524261)) + + w = (T(0.0031533460523059122), T(0.0073275539012762885), T(0.011477234579234613), T(0.015579315722943824), T(0.01961616045735561), T(0.023570760839324363) + , T(0.027426509708356944), T(0.031167227832798003), T(0.03477722256477052), T(0.038241351065830737), T(0.04154508294346467), T(0.0446745608566943), T(0.04761665849249045) + , T(0.05035903555385445), T(0.05289018948519363), T(0.055199503699984116), T(0.05727729210040322), T(0.05911483969839564), T(0.06070443916589387), T(0.06203942315989268) + , T(0.06311419228625402), T(0.06392423858464813), T(0.06446616443594998), T(0.06473769681268386), T(0.06473769681268386), T(0.06446616443594998), T(0.06392423858464813) + , T(0.06311419228625402), T(0.06203942315989268), T(0.06070443916589387), T(0.05911483969839564), T(0.05727729210040322), T(0.055199503699984116) + , T(0.05289018948519363), T(0.05035903555385445), T(0.04761665849249045), T(0.0446745608566943), T(0.04154508294346467), T(0.038241351065830737), T(0.03477722256477052) + , T(0.031167227832798003), T(0.027426509708356944), T(0.023570760839324363), T(0.01961616045735561), T(0.015579315722943824), T(0.011477234579234613), T(0.0073275539012762885) + , T(0.0031533460523059122)) + + towen = dot(w, t2.(h,a,x)) + + return towen + else + # a > 0.999999, T6 from paper (quadrature using QuadGK would also work, but be slower) + + j = T(0.5)*erfc(-h*invsqrt2_T) + k = atan((one(a)-a)/(one(a)+a)) + towen = T(0.5)*j*(one(h)-j)-inv2π_T*k*exp((-T(0.5)*(one(a)-a)*h^2)/k) + + return towen + end +end + +t2(h::T, a, x) where T = T(inv4π)*a*exp(-T(0.5)*(h^2)*(one(h)+(a*x)^2))/(one(h)+(a*x)^2) + +owent(h::Real, a::Real) = owent(promote(h,a)...) diff --git a/test/owent.jl b/test/owent.jl new file mode 100644 index 00000000..eb1f2a1d --- /dev/null +++ b/test/owent.jl @@ -0,0 +1,104 @@ +# Owen's T function tests + +using Test +using IrrationalConstants +using LinearAlgebra +using SpecialFunctions + +# test values for accurate and precise calculation +hvec = [0.0625, 6.5, 7.0, 4.78125, 2.0, 1.0, 0.0625, 1, 1, 1, 1, 0.5, 0.5, 0.5, 0.5, 0.25, 0.25, 0.25, 0.25, 0.125, 0.125, 0.125, 0.125, 0.0078125 +, 0.0078125, 0.0078125, 0.0078125, 0.0078125, 0.0078125, 0.0625, 0.5, 0.9, 2.5, 7.33, 0.6, 1.6, 2.33, 2.33] +avec = [0.25, 0.4375, 0.96875, 0.0625, 0.5, 0.9999975, 0.999999125, 0.5, 1, 2, 3, 0.5, 1, 2, 3, 0.5, 1, 2, 3, 0.5, 1, 2, 3, 0.5, 1, 2, 3, 10, 100 +, 0.999999999999999, 0.999999999999999, 0.999999999999999, 0.999999999999999, 0.999999999999999, 0.999999999999999, 0.999999999999999, 0.999999999999999 +, 0.99999] +cvec = [big"0.0389119302347013668966224771378", big"2.00057730485083154100907167685e-11", big"6.399062719389853083219914429e-13" +, big"1.06329748046874638058307112826e-7", big"0.00862507798552150713113488319155", big"0.0667418089782285927715589822405" +, big"0.1246894855262192" +, big"0.04306469112078537", big"0.06674188216570097", big"0.0784681869930841", big"0.0792995047488726", big"0.06448860284750375", big"0.1066710629614485" +, big"0.1415806036539784", big"0.1510840430760184", big"0.07134663382271778", big"0.1201285306350883", big"0.1666128410939293", big"0.1847501847929859" +, big"0.07317273327500386", big"0.1237630544953746", big"0.1737438887583106", big"0.1951190307092811", big"0.07378938035365545" +, big"0.1249951430754052", big"0.1761984774738108", big"0.1987772386442824", big"0.2340886964802671", big"0.2479460829231492" +, big"0.1246895548850743676554299881345328280176736760739893903915691894" +, big"0.1066710629614484543187382775527753945264849005582264731161129477" +, big"0.0750909978020473015080760056431386286348318447478899039422181015" +, big"0.0030855526911589942124216949767707430484322201889568086810922629" +, big"5.7538182971139187466647478665179637676531179007295252996453e-14", big"0.0995191725772188724714794470740785702586033387786949658229016920" +, big"0.0258981646643923680014142514442989928165349517076730515952020227" +, big"0.0049025023268168675126146823752680242063832053551244071400100690" +, big"0.0049024988349089450612896251009169062698683918433614542387524648"] + +# test values for type stability checking +ht = [0.0625, 0.0, 0.5, 0.5, 0.5] +at = [0.025, 0.5, 0.0, 1.0, +Inf] + +#= +# Interactive tests, displayin error and error magnitude +mvbck = "\033[1A\033[58G" +mva = "\033[72G" +mve = "\033[90G" + +warmup = owent(0.0625,0.025) +warmup = owent(0.0625,0.999999125) + +for i in 1:size(hvec,1) + if i == 1 + println("\t\tExecution Time","\033[58G","h",mva,"a",mve,"error\t\t","log10 error") + end + h = hvec[i] + a = avec[i] + c = cvec[i] + print(i,"\t") + @time tx = owent(h,a) + err = Float64(tx-c) + logerr = Float64(round(log10(abs(tx-c)),sigdigits=3)) + println(mvbck,h,mva,a,mve,err,"\t",logerr) +end +=# + +@testset "Owen's T" begin + + # check that error for calc is within specification + @testset "Owen's T value checks" begin + + for i in 1:size(hvec,1) + h = hvec[i] # h test value + a = avec[i] # a test value + c = cvec[i] # the "correct" answer + t = owent(h,a) + + err = round(log10(abs(t-c)),digits=0) + + @test err ≤ -16.0 + end + + end + + @testset "Owen's T Shortcut Evaluations" begin + @test owent(-0.0625,0.025) == owent(0.0625,0.025) # if h<0 equals owent(abs(h),a) + @test owent(0.0,0.025) == atan(0.025)*inv2π # if h=0, t = atan(a)*inv2π + @test owent(0.0625,0.0) == 0.0 # when a=0, owent=0 + @test owent(0.0625,-0.025) == -owent(0.0625,0.025) # if a<0, t = -owent(h,abs(a)) + @test owent(0.0625,1.0) == 0.125*erfc(-0.0625*invsqrt2)*erfc(0.0625*invsqrt2) # if a=1, t = (formula) + @test owent(0.0625,+Inf) == 0.25*erfc(sqrt(0.0625^2)*invsqrt2) # if a=∞, t = (formula) + @test owent(0.0625,-Inf) == -0.25*erfc(sqrt(0.0625^2)*invsqrt2) # should also work, due to a<0 condition + end + + @testset "Owen's T type stability" begin + + for T1 in [Float16, Float32, Float64, BigFloat] + for T2 in [BigFloat, Float64, Float32, Float16] + for i in 1:size(ht,1) + h=T1(ht[i]) + a=T2(at[i]) + (p1, p2) = promote(h,a) + t=owent(h,a) + @test typeof(t) == typeof(p1) + #println("T1: ",T1," ",typeof(p1),"\tT2: ",T2," ",typeof(p2),"\t",i,"\th ",h,"\ta ",a,"\tt ",t,"\t",typeof(t)," ",typeof(t)==typeof(p1)) + end + end + end + + end # test type stability + +end +# Owen's T Function Tests diff --git a/test/runtests.jl b/test/runtests.jl index 69dcafe7..197b9825 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -33,7 +33,8 @@ tests = [ "logabsgamma", "sincosint", "other_tests", - "chainrules" + "chainrules", + "owent" ] const testdir = dirname(@__FILE__) From 4c5659fbae5334bf2830e9631577d9ade150c04e Mon Sep 17 00:00:00 2001 From: Johanni Brea Date: Fri, 24 Jan 2025 10:59:00 +0100 Subject: [PATCH 2/6] rm LinearAlgebra dependency --- Project.toml | 2 -- src/SpecialFunctions.jl | 2 -- src/owent.jl | 5 ++++- 3 files changed, 4 insertions(+), 5 deletions(-) diff --git a/Project.toml b/Project.toml index e735732f..45862755 100644 --- a/Project.toml +++ b/Project.toml @@ -5,7 +5,6 @@ version = "2.5.0" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" IrrationalConstants = "92d709cd-6900-40b7-9082-c6be49f344b6" -LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688" OpenLibm_jll = "05823500-19ac-5b8b-9628-191a04bc5112" OpenSpecFun_jll = "efe28fd5-8261-553b-a9e1-b2916fc3738e" @@ -20,7 +19,6 @@ SpecialFunctionsChainRulesCoreExt = "ChainRulesCore" ChainRulesCore = "0.9.44, 0.10, 1" ChainRulesTestUtils = "0.6.8, 0.7, 1" IrrationalConstants = "0.1, 0.2" -LinearAlgebra = "1.11.0" LogExpFunctions = "0.3.2" OpenLibm_jll = "0.7, 0.8" OpenSpecFun_jll = "0.5" diff --git a/src/SpecialFunctions.jl b/src/SpecialFunctions.jl index 74fa88ea..10927697 100644 --- a/src/SpecialFunctions.jl +++ b/src/SpecialFunctions.jl @@ -1,7 +1,5 @@ module SpecialFunctions -import LinearAlgebra: dot - using IrrationalConstants: twoπ, halfπ, diff --git a/src/owent.jl b/src/owent.jl index 943e28b1..90c7e682 100644 --- a/src/owent.jl +++ b/src/owent.jl @@ -120,7 +120,10 @@ function owent(h::T, a::T) where {T <: Real} , T(0.031167227832798003), T(0.027426509708356944), T(0.023570760839324363), T(0.01961616045735561), T(0.015579315722943824), T(0.011477234579234613), T(0.0073275539012762885) , T(0.0031533460523059122)) - towen = dot(w, t2.(h,a,x)) + towen = zero(a) + @inbounds for i in eachindex(w) + towen += w[i] * t2(h,a,x[i]) + end return towen else From b6c96e72966a45c2c54913f39fce16fdb321e328 Mon Sep 17 00:00:00 2001 From: Johanni Brea Date: Fri, 24 Jan 2025 11:05:02 +0100 Subject: [PATCH 3/6] fix --- test/owent.jl | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/test/owent.jl b/test/owent.jl index eb1f2a1d..4a9347f0 100644 --- a/test/owent.jl +++ b/test/owent.jl @@ -1,9 +1,6 @@ # Owen's T function tests -using Test -using IrrationalConstants -using LinearAlgebra -using SpecialFunctions +using IrrationalConstants: inv2π, invsqrt2 # test values for accurate and precise calculation hvec = [0.0625, 6.5, 7.0, 4.78125, 2.0, 1.0, 0.0625, 1, 1, 1, 1, 0.5, 0.5, 0.5, 0.5, 0.25, 0.25, 0.25, 0.25, 0.125, 0.125, 0.125, 0.125, 0.0078125 From db56d5d59afd7916f06e2ab8c0a1ca5f48742963 Mon Sep 17 00:00:00 2001 From: Johanni Brea Date: Mon, 27 Jan 2025 12:02:53 +0100 Subject: [PATCH 4/6] use sum --- src/owent.jl | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/owent.jl b/src/owent.jl index 90c7e682..d94ad07f 100644 --- a/src/owent.jl +++ b/src/owent.jl @@ -120,12 +120,7 @@ function owent(h::T, a::T) where {T <: Real} , T(0.031167227832798003), T(0.027426509708356944), T(0.023570760839324363), T(0.01961616045735561), T(0.015579315722943824), T(0.011477234579234613), T(0.0073275539012762885) , T(0.0031533460523059122)) - towen = zero(a) - @inbounds for i in eachindex(w) - towen += w[i] * t2(h,a,x[i]) - end - - return towen + return sum(w .* t2.(h, a, x)) else # a > 0.999999, T6 from paper (quadrature using QuadGK would also work, but be slower) From 6bd5170135516810ce1a71073a26ee4041f99ce5 Mon Sep 17 00:00:00 2001 From: Johanni Brea Date: Mon, 27 Jan 2025 12:03:08 +0100 Subject: [PATCH 5/6] improve docs --- src/owent.jl | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/src/owent.jl b/src/owent.jl index d94ad07f..81346f07 100644 --- a/src/owent.jl +++ b/src/owent.jl @@ -1,13 +1,12 @@ # Owen's T Function # Written by Andy Gough; August 2021 (see https://github.com/JuliaStats/StatsFuns.jl/issues/99#issuecomment-1124581689) -# Edited by Johanni Brea to make type stable; January 2025 +# Edited by Johanni Brea; January 2025 # Rev 1.09 # MIT License # # dependencies # IrrationalConstants # SpecialFunctions -# LinearAlgebra # # HISTORY # In the past 20 or so years, most implementations of Owen's T function have followed the algorithms given in "Fast and accurate Calculation of Owen's @@ -33,32 +32,31 @@ # T1 through T4 are not implemented. Instead, if a < 0.999999, T5 is used to calculate Owen's T (using 48 point Gauss-Legendre quadrature) # For 0.999999 < a < 1.0, T6 is implemented. # -# REFERENCES -# [1] "Fast and accurate Calculation of Owen's T-Function", by M. Patefield and D. Tandy, Journal of Statistical Software, 5 (5), 1 - 25 (2000) -# [2] "Tables for Computing Bivariate Normal Probabilities", by Donald P. Owen, The Annals of Mathematical Statistics, Vol. 27, No. 4 (Dec 1956), pp. 1075-1090 -# # Partial Derivatives (FYI) # D[owent[x,a],x] = -exp(-0.5*x^2)*erf(a*x/sqrt2)/(2*sqrt2π) # D[owent[x,a],a] = exp(-0.5*(1+a^2)*(x^2))/((1+a^2)*2π) # @doc raw""" -owent(h,a) : Returns the value of Owen's T function for (h,a) + owent(h, a) -Owen's T function: +Returns the value of Owen's T function ```math T(h,a) = \frac{1}{2\pi } \int_{0}^{a} \frac{e^{-\frac{1}{2}h^2(1+x^2)}}{1+x^2}dx\quad(-\infty < h,a < +\infty) ``` -For *h* and *a* > 0, *T(h,a)* gives the volume of the uncorrelated bivariate normal distribution with zero mean and unit variance -over the area from *y = ax* and *y = 0* and to the right of *x = h*. +For *h* and *a* > 0, *T(h,a)* gives the volume of the uncorrelated bivariate normal distribution with zero mean and unit variance over the area from *y = ax* and *y = 0* and to the right of *x = h*. -EXAMPLE: +## Example ``` julia> owent(0.0625, 0.025) 0.003970281304296922 ``` -Worst case accuracy is about 2e-16. +## References +"Fast and accurate Calculation of Owen's T-Function", by M. Patefield and D. Tandy, Journal of Statistical Software, 5 (5), 1 - 25 (2000) + +"Tables for Computing Bivariate Normal Probabilities", by Donald P. Owen, The Annals of Mathematical Statistics, Vol. 27, No. 4 (Dec 1956), pp. 1075-1090 +# """ function owent(h::T, a::T) where {T <: Real} From e4b203be5877ff0c2d5b3f455397afeec7c7bc42 Mon Sep 17 00:00:00 2001 From: Johanni Brea Date: Mon, 27 Jan 2025 14:43:31 +0100 Subject: [PATCH 6/6] Float16, Float32, Float64, BigFloat --- src/owent.jl | 36 +++++++++++++++++------------------- 1 file changed, 17 insertions(+), 19 deletions(-) diff --git a/src/owent.jl b/src/owent.jl index 81346f07..f048e45d 100644 --- a/src/owent.jl +++ b/src/owent.jl @@ -99,25 +99,7 @@ function owent(h::T, a::T) where {T <: Real} # calculate Owen's T if a ≤ T(0.999999) - - - x = (-T(0.9987710072524261), -T(0.9935301722663508), -T(0.9841245837228269), -T(0.9705915925462473), -T(0.9529877031604309), -T(0.9313866907065543), -T(0.9058791367155696) - , -T(0.8765720202742479), -T(0.8435882616243935), -T(0.8070662040294426), -T(0.7671590325157404), -T(0.7240341309238146), -T(0.6778723796326639), -T(0.6288673967765136) - , -T(0.5772247260839727), -T(0.5231609747222331), -T(0.4669029047509584), -T(0.4086864819907167), -T(0.34875588629216075), -T(0.28736248735545555), -T(0.22476379039468905) - , -T(0.1612223560688917), -T(0.0970046992094627), -T(0.03238017096286937), T(0.03238017096286937), T(0.0970046992094627), T(0.1612223560688917), T(0.22476379039468905) - , T(0.28736248735545555), T(0.34875588629216075), T(0.4086864819907167), T(0.4669029047509584), T(0.5231609747222331), T(0.5772247260839727), T(0.6288673967765136) - , T(0.6778723796326639), T(0.7240341309238146), T(0.7671590325157404), T(0.8070662040294426), T(0.8435882616243935), T(0.8765720202742479), T(0.9058791367155696) - , T(0.9313866907065543), T(0.9529877031604309), T(0.9705915925462473), T(0.9841245837228269), T(0.9935301722663508), T(0.9987710072524261)) - - w = (T(0.0031533460523059122), T(0.0073275539012762885), T(0.011477234579234613), T(0.015579315722943824), T(0.01961616045735561), T(0.023570760839324363) - , T(0.027426509708356944), T(0.031167227832798003), T(0.03477722256477052), T(0.038241351065830737), T(0.04154508294346467), T(0.0446745608566943), T(0.04761665849249045) - , T(0.05035903555385445), T(0.05289018948519363), T(0.055199503699984116), T(0.05727729210040322), T(0.05911483969839564), T(0.06070443916589387), T(0.06203942315989268) - , T(0.06311419228625402), T(0.06392423858464813), T(0.06446616443594998), T(0.06473769681268386), T(0.06473769681268386), T(0.06446616443594998), T(0.06392423858464813) - , T(0.06311419228625402), T(0.06203942315989268), T(0.06070443916589387), T(0.05911483969839564), T(0.05727729210040322), T(0.055199503699984116) - , T(0.05289018948519363), T(0.05035903555385445), T(0.04761665849249045), T(0.0446745608566943), T(0.04154508294346467), T(0.038241351065830737), T(0.03477722256477052) - , T(0.031167227832798003), T(0.027426509708356944), T(0.023570760839324363), T(0.01961616045735561), T(0.015579315722943824), T(0.011477234579234613), T(0.0073275539012762885) - , T(0.0031533460523059122)) - + x, w = gauss_legendre(T) return sum(w .* t2.(h, a, x)) else # a > 0.999999, T6 from paper (quadrature using QuadGK would also work, but be slower) @@ -133,3 +115,19 @@ end t2(h::T, a, x) where T = T(inv4π)*a*exp(-T(0.5)*(h^2)*(one(h)+(a*x)^2))/(one(h)+(a*x)^2) owent(h::Real, a::Real) = owent(promote(h,a)...) + +# 48-point Gauss-Legendre quadrature (Arblib.hypgeom_legendre_p_ui_root!) +gauss_legendre(::Type{Float64}) = + (0.9987710072524261, 0.9935301722663508, 0.9841245837228269, 0.9705915925462473, 0.9529877031604309, 0.9313866907065543, 0.9058791367155696, 0.8765720202742479, 0.8435882616243935, 0.8070662040294426, 0.7671590325157404, 0.7240341309238146, 0.6778723796326639, 0.6288673967765136, 0.5772247260839727, 0.523160974722233, 0.4669029047509584, 0.4086864819907167, 0.34875588629216075, 0.28736248735545555, 0.22476379039468905, 0.1612223560688917, 0.0970046992094627, 0.03238017096286936, -0.03238017096286936, -0.0970046992094627, -0.1612223560688917, -0.22476379039468905, -0.28736248735545555, -0.34875588629216075, -0.4086864819907167, -0.4669029047509584, -0.523160974722233, -0.5772247260839727, -0.6288673967765136, -0.6778723796326639, -0.7240341309238146, -0.7671590325157404, -0.8070662040294426, -0.8435882616243935, -0.8765720202742479, -0.9058791367155696, -0.9313866907065543, -0.9529877031604309, -0.9705915925462473, -0.9841245837228269, -0.9935301722663508, -0.9987710072524261), + (0.0031533460523058385, 0.0073275539012762625, 0.01147723457923454, 0.015579315722943849, 0.01961616045735553, 0.02357076083932438, 0.027426509708356948, 0.03116722783279809, 0.03477722256477044, 0.03824135106583071, 0.04154508294346475, 0.04467456085669428, 0.04761665849249048, 0.05035903555385447, 0.05289018948519367, 0.055199503699984165, 0.057277292100403214, 0.059114839698395635, 0.06070443916589388, 0.062039423159892665, 0.06311419228625402, 0.06392423858464819, 0.06446616443595009, 0.06473769681268392, 0.06473769681268392, 0.06446616443595009, 0.06392423858464819, 0.06311419228625402, 0.062039423159892665, 0.06070443916589388, 0.059114839698395635, 0.057277292100403214, 0.055199503699984165, 0.05289018948519367, 0.05035903555385447, 0.04761665849249048, 0.04467456085669428, 0.04154508294346475, 0.03824135106583071, 0.03477722256477044, 0.03116722783279809, 0.027426509708356948, 0.02357076083932438, 0.01961616045735553, 0.015579315722943849, 0.01147723457923454, 0.0073275539012762625, 0.0031533460523058385) +# 24-point Gauss-Legendre quadrature (Arblib.hypgeom_legendre_p_ui_root!) +gauss_legendre(::Type{Float32}) = + (0.9951872f0, 0.9747286f0, 0.93827456f0, 0.88641554f0, 0.82000196f0, 0.74012417f0, 0.64809364f0, 0.5454215f0, 0.43379351f0, 0.31504267f0, 0.19111887f0, 0.064056896f0, -0.064056896f0, -0.19111887f0, -0.31504267f0, -0.43379351f0, -0.5454215f0, -0.64809364f0, -0.74012417f0, -0.82000196f0, -0.88641554f0, -0.93827456f0, -0.9747286f0, -0.9951872f0), + (0.01234123f0, 0.02853139f0, 0.044277437f0, 0.059298586f0, 0.07334648f0, 0.086190164f0, 0.097618654f0, 0.10744427f0, 0.115505666f0, 0.12167047f0, 0.12583746f0, 0.1279382f0, 0.1279382f0, 0.12583746f0, 0.12167047f0, 0.115505666f0, 0.10744427f0, 0.097618654f0, 0.086190164f0, 0.07334648f0, 0.059298586f0, 0.044277437f0, 0.02853139f0, 0.01234123f0) +gauss_legendre(::Type{Float16}) = + (Float16(0.9814), Float16(0.9043), Float16(0.77), Float16(0.5874), Float16(0.368), Float16(0.1252), Float16(-0.1252), Float16(-0.368), Float16(-0.5874), Float16(-0.77), Float16(-0.9043), Float16(-0.9814)), + (Float16(0.04718), Float16(0.10693), Float16(0.16), Float16(0.2031), Float16(0.2335), Float16(0.2491), Float16(0.2491), Float16(0.2335), Float16(0.2031), Float16(0.16), Float16(0.10693), Float16(0.04718)) +# 48-point Gauss-Legendre quadrature (Arblib.hypgeom_legendre_p_ui_root!) +gauss_legendre(::Type{BigFloat}) = + (big"0.9987710072524261186005414915631136400889376502767210386129404813754588436074878", big"0.9935301722663507575479287508490741183566147495946719296171518380987546182067713", big"0.9841245837228268577445836000265988305892392234173847299576501679855297780009794", big"0.9705915925462472504614119838006600573024339116308837060283723521653233091284874", big"0.9529877031604308607229606660257183432085413318239187368639476034939458705853333", big"0.9313866907065543331141743801016012677199970856189504298706048642530730422171056", big"0.9058791367155696728220748356710117883122621998274108453524854254710168231209838", big"0.8765720202742478859056935548050967545616485337299619927478757518746727101403824", big"0.8435882616243935307110898445196560498708870117375524015149131998988410546898503", big"0.8070662040294426270825530430245384459730130294604153865758629418121821540044232", big"0.7671590325157403392538554375229690536226423308482073722351285886640508368078524", big"0.7240341309238146546744822334936652465850928122807223627293663025733514606200864", big"0.6778723796326639052118512806759090588499546790260486130710406429754946468798164", big"0.6288673967765136239951649330699946520249089997901617709817329945195319139770715", big"0.5772247260839727038178092385404787728539972861401955280523973994277369963343583", big"0.5231609747222330336782258691375085262891876218118841075802295472194144547473473", big"0.4669029047509584045449288616507985092368121042585169441818691951347943934426029", big"0.4086864819907167299162254958146332864599228429948880647711509833256205384841253", big"0.3487558862921607381598179372704079161343096499683925760321229677812815940686729", big"0.2873624873554555767358864613167976878515583058010397789085000321689998442687597", big"0.224763790394689061224865440174692277438561804041654806164742641045181941897513", big"0.161222356068891718056437390783497694774374379741895117703242637556516342099581", big"0.09700469920946269893005395585362452015273622930093698643058076594480403626262214", big"0.03238017096286936203332224315213444204596280236151809242500322001737781920338223", big"-0.03238017096286936203332224315213444204596280236151809242500322001737781920338223", big"-0.09700469920946269893005395585362452015273622930093698643058076594480403626262214", big"-0.161222356068891718056437390783497694774374379741895117703242637556516342099581", big"-0.224763790394689061224865440174692277438561804041654806164742641045181941897513", big"-0.2873624873554555767358864613167976878515583058010397789085000321689998442687597", big"-0.3487558862921607381598179372704079161343096499683925760321229677812815940686729", big"-0.4086864819907167299162254958146332864599228429948880647711509833256205384841253", big"-0.4669029047509584045449288616507985092368121042585169441818691951347943934426029", big"-0.5231609747222330336782258691375085262891876218118841075802295472194144547473473", big"-0.5772247260839727038178092385404787728539972861401955280523973994277369963343583", big"-0.6288673967765136239951649330699946520249089997901617709817329945195319139770715", big"-0.6778723796326639052118512806759090588499546790260486130710406429754946468798164", big"-0.7240341309238146546744822334936652465850928122807223627293663025733514606200864", big"-0.7671590325157403392538554375229690536226423308482073722351285886640508368078524", big"-0.8070662040294426270825530430245384459730130294604153865758629418121821540044232", big"-0.8435882616243935307110898445196560498708870117375524015149131998988410546898503", big"-0.8765720202742478859056935548050967545616485337299619927478757518746727101403824", big"-0.9058791367155696728220748356710117883122621998274108453524854254710168231209838", big"-0.9313866907065543331141743801016012677199970856189504298706048642530730422171056", big"-0.9529877031604308607229606660257183432085413318239187368639476034939458705853333", big"-0.9705915925462472504614119838006600573024339116308837060283723521653233091284874", big"-0.9841245837228268577445836000265988305892392234173847299576501679855297780009794", big"-0.9935301722663507575479287508490741183566147495946719296171518380987546182067713", big"-0.9987710072524261186005414915631136400889376502767210386129404813754588436074878"), + (big"0.003153346052305838632677311543891487578283938831693622295209493250319586438316842", big"0.007327553901276262102383979621786550058707902559201353274881829548806980072502799", big"0.01147723457923453948959266760909162808642050630874764065376681674103503658508731", big"0.01557931572294384872817695583446031397637626899155246951309343105269243335619984", big"0.01961616045735552781446071965221270969581303773413223918112083050740924629812146", big"0.02357076083932437914051930137844923022172973852218859873423906486456506379639118", big"0.0274265097083569482000738362625058204511841551616509759972809374993765019410236", big"0.03116722783279808890206575684635441945428534148356953550954371886143141262424302", big"0.03477722256477043889254858596380241059728139690706809871800663617967672335903626", big"0.03824135106583070631721725652371561786382396835498228892925819103405053922410909", big"0.04154508294346474921405882236106479775347282603403806308273482122272582562965843", big"0.04467456085669428041944858712585039498846278686250200843292144633919149051230188", big"0.04761665849249047482590662347892983015799806674344968539676989627880988507905503", big"0.05035903555385447495780761908786560603299409302590633069379205724693441466024811", big"0.05289018948519366709550505626469891466172648563310918638649123384829276249063296", big"0.05519950369998416286820349519163543900445092560756100054805625793058523675145725", big"0.05727729210040321570515023468470057624152712300411207753884993747681745421856388", big"0.05911483969839563574647481743351991065965560255705499855629113348583514270048057", big"0.06070443916589388005296923202782047788526086425647775511151144466063789427975123", big"0.06203942315989266390419778413759851830638339966509146156903781450273903590161649", big"0.06311419228625402565712602275023331812741364337110079121114724790803811921086588", big"0.06392423858464818662390620182551540891897408498264299989087420749955378258611148", big"0.06446616443595008220650419365770506572569192445553030876055845653739235337295456", big"0.06473769681268392250302493873659155355208191894663651001456309552308307891126462", big"0.06473769681268392250302493873659155355208191894663651001456309552308307891126462", big"0.06446616443595008220650419365770506572569192445553030876055845653739235337295456", big"0.06392423858464818662390620182551540891897408498264299989087420749955378258611148", big"0.06311419228625402565712602275023331812741364337110079121114724790803811921086588", big"0.06203942315989266390419778413759851830638339966509146156903781450273903590161649", big"0.06070443916589388005296923202782047788526086425647775511151144466063789427975123", big"0.05911483969839563574647481743351991065965560255705499855629113348583514270048057", big"0.05727729210040321570515023468470057624152712300411207753884993747681745421856388", big"0.05519950369998416286820349519163543900445092560756100054805625793058523675145725", big"0.05289018948519366709550505626469891466172648563310918638649123384829276249063296", big"0.05035903555385447495780761908786560603299409302590633069379205724693441466024811", big"0.04761665849249047482590662347892983015799806674344968539676989627880988507905503", big"0.04467456085669428041944858712585039498846278686250200843292144633919149051230188", big"0.04154508294346474921405882236106479775347282603403806308273482122272582562965843", big"0.03824135106583070631721725652371561786382396835498228892925819103405053922410909", big"0.03477722256477043889254858596380241059728139690706809871800663617967672335903626", big"0.03116722783279808890206575684635441945428534148356953550954371886143141262424302", big"0.0274265097083569482000738362625058204511841551616509759972809374993765019410236", big"0.02357076083932437914051930137844923022172973852218859873423906486456506379639118", big"0.01961616045735552781446071965221270969581303773413223918112083050740924629812146", big"0.01557931572294384872817695583446031397637626899155246951309343105269243335619984", big"0.01147723457923453948959266760909162808642050630874764065376681674103503658508731", big"0.007327553901276262102383979621786550058707902559201353274881829548806980072502799", big"0.003153346052305838632677311543891487578283938831693622295209493250319586438316842")