[WIP] port of Cephes implementation of Beta function (fixes #14256) #14349

CarloLucibello · 2015-12-10T14:07:13Z

port of https://github.com/scipy/scipy/blob/0cff7a5fe6226668729fc2551105692ce114c2b3/scipy/special/cephes/beta.c#L61-L133

Copyright issues (none actually) are discussed in #14256

see also #14165 for its relation with the binomial function

Comments appreciated @tkelman @stevengj @andreasnoack

TODO

extend test coverage
check if it works correctly for complex arguments
address type instability (help needed here):

julia> @code_warntype beta(3.,2.)
Variables:
  a::Float64
  b::Float64
  #s5::Bool
  #s8::Bool
  #s1::Int64
  y::Any
  s::Any
  #s6::Bool
  #s7::Bool
  #s2::Any
  #s3::Any
  yb::Any
  sb::Any
  #s4::Any
  ya::Any
  sa::Any
  ######fy#1781#7325#7333::Float64
  ######fy#1781#7325#7334::Float64

Body:
  begin  # /home/carlo/Git/julia/beta.jl, line 6:
      NewvarNode(symbol("#s5"))
      NewvarNode(symbol("#s8"))
      NewvarNode(symbol("#s1"))
      NewvarNode(:y)
      NewvarNode(:s)
      NewvarNode(symbol("#s6"))
      NewvarNode(symbol("#s7"))
      NewvarNode(symbol("#s2"))
      NewvarNode(symbol("#s3"))
      NewvarNode(:yb)
      NewvarNode(:sb)
      NewvarNode(symbol("#s4"))
      NewvarNode(:ya)
      NewvarNode(:sa)
      unless (Base.le_float)(a::Float64,0.0)::Bool goto 1 # floatfuncs.jl, line 27: # float.jl, line 309: # float.jl, line 268:
      ######fy#1781#7325#7333 = (Base.box)(Float64,(Base.sitofp)(Float64,0)) # float.jl, line 269:
      unless (Base.box)(Base.Bool,(Base.and_int)((Base.eq_float)((Base.box)(Base.Float64,(Base.trunc_llvm)(a::Float64)),a::Float64)::Bool,(Base.box)(Base.Bool,(Base.and_int)((Base.eq_float)((Base.box)(Base.Float64,(Base.sub_float)(a::Float64,a::Float64)),######fy#1781#7325#7333::Float64)::Bool,0 === (Base.box)(Int64,(Base.fptosi)(Int64,######fy#1781#7325#7333::Float64))::Bool)))) goto 0
      return (Main.beta_negint)(a::Float64,b::Float64)::Any
      0: 
      goto 1
      1:  # /home/carlo/Git/julia/beta.jl, line 7:
      unless (Base.le_float)(b::Float64,0.0)::Bool goto 3 # floatfuncs.jl, line 27: # float.jl, line 309: # float.jl, line 268:
      ######fy#1781#7325#7334 = (Base.box)(Float64,(Base.sitofp)(Float64,0)) # float.jl, line 269:
      unless (Base.box)(Base.Bool,(Base.and_int)((Base.eq_float)((Base.box)(Base.Float64,(Base.trunc_llvm)(b::Float64)),b::Float64)::Bool,(Base.box)(Base.Bool,(Base.and_int)((Base.eq_float)((Base.box)(Base.Float64,(Base.sub_float)(b::Float64,b::Float64)),######fy#1781#7325#7334::Float64)::Bool,0 === (Base.box)(Int64,(Base.fptosi)(Int64,######fy#1781#7325#7334::Float64))::Bool)))) goto 2
      return (Main.beta_negint)(b::Float64,a::Float64)::Any
      2: 
      goto 3
      3:  # /home/carlo/Git/julia/beta.jl, line 9:
      unless (Base.lt_float)((Base.box)(Base.Float64,(Base.abs_float)(a::Float64)),(Base.box)(Base.Float64,(Base.abs_float)(b::Float64)))::Bool goto 4 # /home/carlo/Git/julia/beta.jl, line 10:
      GenSym(0) = b::Float64
      GenSym(1) = a::Float64
      a = GenSym(0)
      b = GenSym(1)
      4:  # /home/carlo/Git/julia/beta.jl, line 13:
      unless (Base.lt_float)((Base.box)(Base.Float64,(Base.mul_float)(Main.ASYMP_FACTOR,(Base.box)(Base.Float64,(Base.abs_float)(b::Float64)))),(Base.box)(Base.Float64,(Base.abs_float)(a::Float64)))::Bool goto 5
      #s5 = (Base.lt_float)(Main.ASYMP_FACTOR,a::Float64)::Bool
      goto 6
      5: 
      #s5 = false
      6: 
      unless #s5::Bool goto 7 # /home/carlo/Git/julia/beta.jl, line 15:
      GenSym(2) = (Main.lbeta_asymp)(a::Float64,b::Float64)::Tuple{Any,Any}
      #s1 = 1
      GenSym(21) = (Base.getfield)(GenSym(2),1)::Any
      GenSym(22) = (Base.box)(Base.Int,(Base.add_int)(1,1))
      y = GenSym(21)
      #s1 = GenSym(22)
      GenSym(23) = (Base.getfield)(GenSym(2),2)::Any
      GenSym(24) = (Base.box)(Base.Int,(Base.add_int)(2,1))
      s = GenSym(23)
      #s1 = GenSym(24) # /home/carlo/Git/julia/beta.jl, line 16:
      return s * (Main.exp)(y)::Any::Any
      7:  # /home/carlo/Git/julia/beta.jl, line 19:
      y = (Base.box)(Base.Float64,(Base.add_float)(a::Float64,b::Float64)) # /home/carlo/Git/julia/beta.jl, line 20:
      #s6 = (Base.lt_float)(Main.MAXGAM,(Base.box)(Base.Float64,(Base.abs_float)(y::Float64)))::Bool
      unless #s6::Bool goto 8
      #s8 = #s6::Bool
      goto 10
      8: 
      #s7 = (Base.lt_float)(Main.MAXGAM,(Base.box)(Base.Float64,(Base.abs_float)(a::Float64)))::Bool
      unless #s7::Bool goto 9
      #s8 = #s7::Bool
      goto 10
      9: 
      #s8 = (Base.lt_float)(Main.MAXGAM,(Base.box)(Base.Float64,(Base.abs_float)(b::Float64)))::Bool
      10: 
      GenSym(5) = #s8::Bool
      unless GenSym(5) goto 11 # /home/carlo/Git/julia/beta.jl, line 21:
      GenSym(6) = (Main.lgamma_r)(y::Float64)::Any
      #s2 = (top(start))(GenSym(6))::Any
      GenSym(7) = (top(indexed_next))(GenSym(6),1,#s2)::Any
      y = (top(getfield))(GenSym(7),1)::Any
      #s2 = (top(getfield))(GenSym(7),2)::Any
      GenSym(8) = (top(indexed_next))(GenSym(6),2,#s2)::Any
      s = (top(getfield))(GenSym(8),1)::Any
      #s2 = (top(getfield))(GenSym(8),2)::Any # /home/carlo/Git/julia/beta.jl, line 22:
      GenSym(9) = (Main.lgamma_r)(b::Float64)::Any
      #s3 = (top(start))(GenSym(9))::Any
      GenSym(10) = (top(indexed_next))(GenSym(9),1,#s3)::Any
      yb = (top(getfield))(GenSym(10),1)::Any
      #s3 = (top(getfield))(GenSym(10),2)::Any
      GenSym(11) = (top(indexed_next))(GenSym(9),2,#s3)::Any
      sb = (top(getfield))(GenSym(11),1)::Any
      #s3 = (top(getfield))(GenSym(11),2)::Any # /home/carlo/Git/julia/beta.jl, line 23:
      y = yb - y::Any # /home/carlo/Git/julia/beta.jl, line 24:
      GenSym(12) = (Main.lgamma_r)(a::Float64)::Any
      #s4 = (top(start))(GenSym(12))::Any
      GenSym(13) = (top(indexed_next))(GenSym(12),1,#s4)::Any
      ya = (top(getfield))(GenSym(13),1)::Any
      #s4 = (top(getfield))(GenSym(13),2)::Any
      GenSym(14) = (top(indexed_next))(GenSym(12),2,#s4)::Any
      sa = (top(getfield))(GenSym(14),1)::Any
      #s4 = (top(getfield))(GenSym(14),2)::Any # /home/carlo/Git/julia/beta.jl, line 25:
      y = ya + y::Any # /home/carlo/Git/julia/beta.jl, line 29:
      return s * sa * sb * (Main.exp)(y)::Any::Any
      11:  # /home/carlo/Git/julia/beta.jl, line 32: # special/gamma.jl, line 3:
      GenSym(16) = (top(ccall))((top(tuple))(:tgamma,Base.Math.libm)::Tuple{Symbol,ASCIIString},Base.Math.Float64,(top(svec))(Base.Math.Float64)::SimpleVector,y::Float64,0)::Float64
      y = (Base.Math.nan_dom_err)(GenSym(16),y::Float64)::Float64 # /home/carlo/Git/julia/beta.jl, line 33: # special/gamma.jl, line 3:
      GenSym(18) = (top(ccall))((top(tuple))(:tgamma,Base.Math.libm)::Tuple{Symbol,ASCIIString},Base.Math.Float64,(top(svec))(Base.Math.Float64)::SimpleVector,a::Float64,0)::Float64
      a = (Base.Math.nan_dom_err)(GenSym(18),a::Float64)::Float64 # /home/carlo/Git/julia/beta.jl, line 34: # special/gamma.jl, line 3:
      GenSym(20) = (top(ccall))((top(tuple))(:tgamma,Base.Math.libm)::Tuple{Symbol,ASCIIString},Base.Math.Float64,(top(svec))(Base.Math.Float64)::SimpleVector,b::Float64,0)::Float64
      b = (Base.Math.nan_dom_err)(GenSym(20),b::Float64)::Float64 # /home/carlo/Git/julia/beta.jl, line 35:
      unless (Base.eq_float)(y::Float64,0.0)::Bool goto 12
      return Main.Inf
      12:  # /home/carlo/Git/julia/beta.jl, line 37:
      unless (Base.lt_float)((Base.box)(Base.Float64,(Base.abs_float)((Base.box)(Base.Float64,(Base.sub_float)((Base.box)(Base.Float64,(Base.abs_float)(b::Float64)),(Base.box)(Base.Float64,(Base.abs_float)(y::Float64)))))),(Base.box)(Base.Float64,(Base.abs_float)((Base.box)(Base.Float64,(Base.sub_float)((Base.box)(Base.Float64,(Base.abs_float)(a::Float64)),(Base.box)(Base.Float64,(Base.abs_float)(y::Float64)))))))::Bool goto 13 # /home/carlo/Git/julia/beta.jl, line 38:
      y = (Base.box)(Base.Float64,(Base.div_float)(b::Float64,y::Float64)) # /home/carlo/Git/julia/beta.jl, line 39:
      y = (Base.box)(Base.Float64,(Base.mul_float)(y::Float64,a::Float64))
      goto 14
      13:  # /home/carlo/Git/julia/beta.jl, line 41:
      y = (Base.box)(Base.Float64,(Base.div_float)(a::Float64,y::Float64)) # /home/carlo/Git/julia/beta.jl, line 42:
      y = (Base.box)(Base.Float64,(Base.mul_float)(y::Float64,b::Float64))
      14:  # /home/carlo/Git/julia/beta.jl, line 45:
      return y::Float64
  end::Any

tkelman · 2015-12-10T14:24:28Z

direct ports need to retain original copyright notices. either the header of this file should be modified to note the exception (along with license.md and contrib/add_license_to_files.jl) or we could make a separate file if we expect to port anything else from the same source.

KristofferC · 2015-12-10T14:30:53Z

base/special/gamma.jl

+
+ # Asymptotic expansion for  ln(|B(a, b)|) for a > ASYMP_FACTOR*max(|b|, 1).
+function lbeta_asymp(a::Number, b::Number)
+    r, s = lgammma_r(b)


CarloLucibello · 2015-12-10T14:40:08Z

sorry, I was working and testing a temporary file, should pass tests now

KristofferC · 2015-12-10T14:46:37Z

You still have 3 m in lgamma_r

CarloLucibello · 2015-12-10T14:49:12Z

direct ports need to retain original copyright notices. either the header of this file should be modified to note the exception (along with license.md and contrib/add_license_to_files.jl) or we could make a separate file if we expect to port anything else from the same source.

will do

stevengj · 2015-12-10T18:43:41Z

The basic problem is that the Cephes code is really only for double precision. You can't use the same constant factors MAXGAM and ASYMP_FACTOR in BigFloat precision.

stevengj · 2015-12-10T18:46:55Z

base/special/gamma.jl

+    real(a) <= 0.0 && isinteger(a) && return beta_negint(a, b)
+    real(b) <= 0.0 && isinteger(b) && return beta_negint(b, a)
+
+    if abs(a) < abs(b)


Note that abs is unnecessarily slow for complex numbers. I would tend to define a fastabs(x), where fastabs(x::Real) = abs(x) and fastabs(z::Complex) = abs(real(z)) + abs(imag(z)). This is within a factor of sqrt(2) of abs(z) but is much faster, and the sqrt(2) factor doesn't really matter for an application like this if you adjust the thresholds accordingly.

- use fastabs for complex numbers - add a test for asymptotic behaviour from issue JuliaLang#4301

CarloLucibello · 2016-01-19T12:55:29Z

todo:

make thresholds of asymptotic behaviours depend on input type. Not sure how to compute them...
check type stability

ararslan · 2019-03-18T23:39:25Z

beta and related functionality now lives in the SpecialFunctions package. I'll close this PR as it no longer applies to Base, but feel free to continue with these changes in the SpecialFunctions package.

CarloLucibello added 5 commits November 26, 2015 18:30

generalize binomial to real and complex arguments

ef1e530

Merge branch 'master' of https://github.com/CarloLucibello/julia

360e4d1

port of CEPHES implementation of Beta function

d0055a1

add some tests

67b47cf

some fix for complexes, still TODO though

088508c

KristofferC reviewed Dec 10, 2015
View reviewed changes

fixes

7fbfaff

fix lgammma

8601710

stevengj reviewed Dec 10, 2015
View reviewed changes

kshyatt added the maths Mathematical functions label Dec 10, 2015

simonbyrne mentioned this pull request Jan 10, 2016

beta function returns imprecise results #14620

Closed

CarloLucibello added 7 commits January 19, 2016 12:13

Merge branch 'master' of https://github.com/JuliaLang/julia

fb8e4ba

port of CEPHES implementation of Beta function

efce6cc

add some tests

d2de1b3

some fix for complexes, still TODO though

dab853d

fixes

0b93ab9

fix lgammma

8065f51

- Add credits to Cephes

5e2e09f

- use fastabs for complex numbers - add a test for asymptotic behaviour from issue JuliaLang#4301

rebase

cbea222

CarloLucibello added 2 commits January 24, 2016 12:50

fix typewarns for real types. Still to address complexes.

c5e346d

type stability for complex numbers

8bf6196

ararslan closed this Mar 18, 2019

Peiffap mentioned this pull request May 6, 2019

beta(-n,n) gives NaN, when n is an integer JuliaMath/SpecialFunctions.jl#149

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

[WIP] port of Cephes implementation of Beta function (fixes #14256) #14349

[WIP] port of Cephes implementation of Beta function (fixes #14256) #14349

CarloLucibello commented Dec 10, 2015

tkelman commented Dec 10, 2015

KristofferC Dec 10, 2015

CarloLucibello commented Dec 10, 2015

KristofferC commented Dec 10, 2015

CarloLucibello commented Dec 10, 2015

stevengj commented Dec 10, 2015

stevengj Dec 10, 2015

CarloLucibello commented Jan 19, 2016

ararslan commented Mar 18, 2019

[WIP] port of Cephes implementation of Beta function (fixes #14256) #14349

[WIP] port of Cephes implementation of Beta function (fixes #14256) #14349

Conversation

CarloLucibello commented Dec 10, 2015

tkelman commented Dec 10, 2015

KristofferC Dec 10, 2015

Choose a reason for hiding this comment

CarloLucibello commented Dec 10, 2015

KristofferC commented Dec 10, 2015

CarloLucibello commented Dec 10, 2015

stevengj commented Dec 10, 2015

stevengj Dec 10, 2015

Choose a reason for hiding this comment

CarloLucibello commented Jan 19, 2016

ararslan commented Mar 18, 2019