Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Special Functions #305

Open
LucaArgentiUCF opened this issue Jan 19, 2021 · 11 comments
Open

Special Functions #305

LucaArgentiUCF opened this issue Jan 19, 2021 · 11 comments
Labels
topic: mathematics linear algebra, sparse matrices, special functions, FFT, random numbers, statistics, ...

Comments

@LucaArgentiUCF
Copy link

There is a standard for the class of numerical functions that are needed in scientific computation. Traditionally, that was for a long time the "blue book" *Handbook of Mathematical Functions", by Abramowitz and Stegun. However, the American National Institute of Standard and Technology (NIST) has sing long ago taken the burden on its shoulder. After countless years of work and revision, the new
Digital Library of Mathematical Functions (DLMF)
has been published.
The DLMF is a the perfect starting point to broaden the standard library. All formulas are tested, some are implemented, and all rigorously referenced and documented. So there is no better blueprint.
Anywhere would be good to start, but I suggest

  • confluent hypergeometric functions
  • legendre and related functions
  • orthogonal polynomials
  • Coulomb functions
  • spherical-group coefficients (3j, 6j, 9j)
@Jim-215-Fisher
Copy link
Contributor

I like the idea. The special functions are needed in almost all areas. Currently, there are only gamma, log_gamma and bessel functions as standard intrinsic, we need more functions like beta, incomplete gamma, incomplete beta functions etc.

@arjenmarkus
Copy link
Member

arjenmarkus commented Jan 20, 2021 via email

@ivan-pi
Copy link
Member

ivan-pi commented Jan 20, 2021

Hello @LucaArgentiUCF ,

I think this issue might be a duplicate of #179 and #102. The former thread contains a link to a proposal with a specification for a special functions module. (The proposal to include such functions directly in the Fortran language standard was ultimately scrapped or rejected.)

If we could agree on an interface, that would be a good first step. We could also "borrow" the interfaces from the major vendor libraries:

I had a go at implementing the inverse error function last year using some preprocessing magic: #179 (comment)

I also started an implementation for a (non-special) cube root function: #214. Initially I thought it would be an easy one, but it turns out catching all the corner cases and guaranteeing full precision over the entire range of floating point variable range, while maintaining good performance is quite tricky. You definitely need a very good understanding of the mathematics, hardware, and numerical analysis if you want to produce high quality implementations.

Ultimately, I halted work on the erfinv function, because I lacked the knowledge (and time) to implement the tests recommed in post #179 (comment) to tabulate the units in last place (ULP). I would appreciate any help to push the erfinv further.

@arjenmarkus
Copy link
Member

arjenmarkus commented Jan 20, 2021 via email

@LucaArgentiUCF
Copy link
Author

I agree the issue overlaps with 179. However, I think that the Digital Library of Mathematical Functions (DLMF) is as authoritative and standard a reference as one can ever possibly get. My suggestion, therefore, would be to shape the special function libraries in terms of one module per DLMF chapter, and figure out which already existing libraries provide the functionalities.

Regarding TOMS, I do not quite get the statement that I have read in some threads that "it is not usable". Unless this sentiment meant to imply poor quality of the library (?), the copyright restrictions to CALGO are not particularly severe: for software past 2013, the rights are with the authors, so one can just ask them (I guess they would be happy to give the rights to use to a standard library they can boast about). For software prior to that date there is only the non-commercial-use clause. Maybe someone see that as a no-go. However, fortran users are for the most part from academia and research centers, so the library modules could include such a disclaimer. Or, if one has just one or two functions from TOMS in a module, one could try to convince ACM to issue a special perpetual no-cost license for those functions in the library. It may be seen by them as a way to promote CALGO, provided that their contribution is advertised in the description of the library.

Cheers,

Luca

@ivan-pi
Copy link
Member

ivan-pi commented Jan 20, 2021

For software prior to that date there is only the non-commercial-use clause. Maybe someone see that as a no-go. However, fortran users are for the most part from academia and research centers, so the library modules could include such a disclaimer.

I don't think this is fully the case. There are numerous engineering companies in the aeronautical, petrochemical, and other industries (optics, physics, astronomy, finance, defense, environment...) which are using Fortran. (I remember at a few Fortran monthly calls we had an employee from Boeing join in.) The projects listed on the fortran-lang page: https://fortran-lang.org/packages/scientific are likely only the tip of the Fortran iceberg.

The presence of the vendor libraries make me think there must be clients willing to pay for them. Personally, I work at the Technical University of Munich which is a high-ranking public research university, yet we don't have direct access to these libraries. (Access to the NAG library is available indirectly through the SuperMUC cluster funded by the Bavarian Academy of Sciences and Humanities. This resource is shared among many European universities, but even internally we need to submit a formal request to gain access.)

I have seen in multiple places of the Julia libraries they rolled their own versions of mathematical algorithms to avoid the ACM license. I lack experience with software licenses and their application in practice to judge whether this was really necessary or not.

@arjenmarkus
Copy link
Member

arjenmarkus commented Jan 20, 2021 via email

@ivan-pi
Copy link
Member

ivan-pi commented Jan 20, 2021

  • All these venerable codes use some means to determine machine parameters - Fortran now offers standard functions. Replace the old probing code by modern Fortran means?

I have seen this done in a few places previously. One example is here:
https://github.com/certik/fortran-utils/blob/master/src/legacy/amos/d1mach.f90

I can't remember which code it was, but I think I have seen the D1MACH function replaced by a static array somewhere.

Concerning SLATEC, I believe it is in the public domain. @jacobwilliams has refactored several components from SLATEC.

@ivan-pi
Copy link
Member

ivan-pi commented Jan 20, 2021

Matching the special functions introduced in C++17 is another reasonable goal:

https://en.cppreference.com/w/cpp/numeric/special_functions

@awvwgk awvwgk added the topic: mathematics linear algebra, sparse matrices, special functions, FFT, random numbers, statistics, ... label Feb 3, 2021
@ivan-pi
Copy link
Member

ivan-pi commented Oct 18, 2021

I found an old thread from @certik, which discusses how to implement the incomplete gamma function:

Fast and accurate double precision implementation of incomplete gamma function: https://scicomp.stackexchange.com/questions/3341/fast-and-accurate-double-precision-implementation-of-incomplete-gamma-function

The function is covered by Chapter 8 in DLMF.

@14NGiestas
Copy link
Member

I'm just moving and linking the old thread mentioned (#102) here, where it looks the discussion got further.

There are lots of implementations online.

As an example, here are the special functions that I needed in my projects over the past 10 years:

https://github.com/certik/fortran-utils/blob/b43bd24cd421509a5bc6d3b9c3eeae8ce856ed88/src/special.f90

https://github.com/certik/hfsolver/blob/b4c50c1979fb7e468b1852b144ba756f5a51788d/src/special_functions.f90

Being able to standardize on the API for all or most of them would be a huge deal. Other languages:

SciPy

https://docs.scipy.org/doc/scipy/reference/special.html

Matlab

https://www.mathworks.com/help/matlab/special-functions-1.html

Julia

Separate package: SpecialFunctions.jl

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
topic: mathematics linear algebra, sparse matrices, special functions, FFT, random numbers, statistics, ...
Projects
None yet
Development

No branches or pull requests

6 participants