Special functions to Base

[lindahua]

@andreasnoack implements a great collection of special functions for Distributions.jl. Some of these functions are probably useful outside of the distribution package.

Source file: https://github.com/JuliaStats/Distributions.jl/blob/12d5aea1d9bea01e4883c120dc870dcf8912ecf9/src/specialfuns.jl

We have a discussion in issue JuliaStats/Distributions.jl#148. Per @StefanKarpinski's suggestion, I open an issue here to invite broader attention.

Basically, we want to discuss: which of those functions should be moved to the Julia base.

[andreasnoack]

@JeffBezanson Follow up on you question in JuliaStats/Distributions.jl#148 (comment). The source is located at http://www.ualberta.ca/CNS/RESEARCH/Software/NumericalNSWC/site.html but it was develop in the US Navy and made public available. The last version is from the beginning of the 90s.

We could just call that code except for the function that initiated my translation, i.e. the incomplete beta function. It is written in single precision, but they had some special computers for naval surface warfare computations in the beginning of the 90s so they could have 14 digits precision in single precision. Therefore the incomplete beta function would need a translation to double precision anyway and my choice was Julia over Fortran.

Also, if you look at the double precision constants they are given with up to 40 digits so by having the code in Julian we can maybe use that extra precision by calling the function with BigFloats.

[simonbyrne]

Also, as @andreasnoack pointed out, there are some accuracy issues with the current beta and lbeta in Base:

julia> lbeta(1e10,1e-10)
23.025848388671875

julia> float64(lbeta(big(1e10),big(1e-10)))
23.02585092758015

julia> Distributions.dbetln(1e10,1e-10)
23.02585092758015

[pao]

@andreasnoack the NSWC library is one we've discussed before. We've been unable to resolve the unclear license of the code. It's probably not eligible for US copyright as the work of public servants, but that's not clearly stated. There's a longer discussion around here somewhere.

[pao]

Oh, found it--you were in that discussion. Carry on. JuliaMath/openlibm#20

[stevengj]

I feel like there should be a JuliaLang/SpecialFunctions.jl package that can be used as a repository for new special-function implementations (or ports), focusing on pure-Julia implementations (as opposed to C or Fortran code, which we can continue to put into openspecfun). Distributions.jl seems like the wrong place.

@nolta has an unregistered https://github.com/nolta/SpecialFunctions.jl package that might be a starting point. (I feel like the https://github.com/nolta/Elliptic.jl package should be folded into a larger special-functions package.)

I have a feeling that some of the more esoteric functions in Base, like the polygamma functions, might be better off moved to such a package too.

I've also been playing with implementing exponential integrals (#7089) and polylogarithms, and it would be helpful to have such a package to put them into since I'm not sure they belong in Base either.

[tkelman]

We should maybe create a new separate org for generic mathematically oriented packages. JuliaLang's CI resources are shared between base and all the packages here, so we should start migrating some of the packages into other organizations. Otherwise agree it's best not to tie these to Base's development cycle.

[johnmyleswhite]

The StatsFuns package is a step in this direction, but I agree that we should just have an actual special functions package.

[stevengj]

@johnmyleswhite, StatsFun is good for functions that are clearly limited to statistics, like the various CDF functions, but for more "generic" transcendental functions like log1psq it would be nice to move them to a separate package.

[simonbyrne]

+1 for a pure Julia SpecialFunctions.jl. Happy to help out, though probably don't have time to drive it.

[JeffBezanson]

Superseded by #8536 --- we should have a SpecialFunctions package.

[simonbyrne]

There is still the question of what should be moved out to a package (or stdlib, or whatever): the currently exported functions from Base.Math are:

• trigonmetric:
  • libm: sin, cos, tan, asin, acos, atan,atan2
  • Julia: sec, csc, cot, asec, acsc, acot,
• hyperbolic:
  • libm: sinh, cosh, tanh, asinh, acosh, atanh,
  • Julia: sech, csch, coth, asech, acsch, acoth,
• sinc:
  • Julia: sinc, cosc
• trig pi multiple:
  • Julia: sinpi, cospi,
• trig degree:
  • Julia: cosd, cotd, cscd, secd, sind, tand, acosd, acotd, acscd, asecd, asind, atand
• degree-radian conversions:
  • Julia: rad2deg, deg2rad
• log:
  • libm: log, log2, log10, log1p
• exp:
  • libm: exp, exp2, expm1
  • Julia: exp10
• roots:
  • intrinsic:sqrt
  • libm: cbrt
• error function:
  • libm: erf, erfc
  • Faddeeva (openspecfun): erfcx, erfi, dawson
  • Julia: erfinv, erfcinv
• gamma:
  • libm: gamma, lgamma
  • Julia: lfact, beta, lbeta, digamma, trigamma,polygamma, invdigamma, eta, zeta
• min/max:
  • Julia: max, min, minmax, clamp, clamp!
• floating point specials:
  • Julia: exponent, significand, ldexp, frexp,
• Airy functions:
  • AMOS (openspecfun): airy, airyai, airyprime, airyaiprime, airybi, airybiprime, airyx
• Bessel functions:
  • AMOS (openspecfun): various besselX and hankelX (there's a lot of these)
• misc:
  • libm: modf, ^
  • Julia: hypot, @evalpoly
  • rem_pio2 (openspecfun): mod2pi

[JeffBezanson]

Categories I'd move to a package(s): error function, gamma, airy, bessel

[simonbyrne]

What about ones that also have mpfr (e.g. bessely0 and erf)? Should we move those to the module?

[JeffBezanson]

Yes those can move to the package as well.

[ararslan]

In reference to what should be moved, see also JuliaMath/SpecialFunctions.jl#12.