search

fortran-lang / stdlib

Fortran Standard Library
https://stdlib.fortran-lang.org
MIT License
1.03k stars 163 forks source link

Special Functions #305

Open LucaArgentiUCF opened 3 years ago

LucaArgentiUCF commented 3 years ago

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

Jim-215-Fisher commented 3 years ago

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 commented 3 years ago

Definitely. There are quite a number of libraries out there that claim to calculate these functions. But a good reference is necessary to ensure the quality. In addition to the NIST library, I have a digital copy of the book "The Mathematical Function Computation Handbook" by Nelson Beebe, He pays a lot of attention to the accuracy and to the corner cases. It could serve as a second opinion :).

Op di 19 jan. 2021 om 21:44 schreef Jing notifications@github.com:

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.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/fortran-lang/stdlib/issues/305#issuecomment-763128046, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAN6YR4FENMILBTKIDQHO7DS2XVLVANCNFSM4WJL6ZSQ .

ivan-pi commented 3 years ago

Hello @LucaArgentiUCF ,

I think this issue might be a duplicate of https://github.com/fortran-lang/stdlib/issues/179 and https://github.com/fortran-lang/stdlib/issues/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: https://github.com/fortran-lang/stdlib/issues/179#issuecomment-647178168

I also started an implementation for a (non-special) cube root function: https://github.com/fortran-lang/stdlib/issues/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 https://github.com/fortran-lang/stdlib/issues/179#issuecomment-647697547 to tabulate the units in last place (ULP). I would appreciate any help to push the erfinv further.

arjenmarkus commented 3 years ago

The SLATEC library on the Netlib site seems to offer quite a few routines for special functions as well.

Things I am wondering about:

Mostly thinking out loud :).

Op wo 20 jan. 2021 om 12:09 schreef Ivan Pribec notifications@github.com:

Hello @LucaArgentiUCF https://github.com/LucaArgentiUCF ,

I think this issue might be a duplicate of #179 https://github.com/fortran-lang/stdlib/issues/179 and #102 https://github.com/fortran-lang/stdlib/issues/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) https://github.com/fortran-lang/stdlib/issues/179#issuecomment-647178168

I also started an implementation for a (non-special) cube root function:

214 https://github.com/fortran-lang/stdlib/issues/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 can be very tricky. You definitely need a very good understand of both the mathematics, hardware, and numerical analysis if you want to produce high quality implementations.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/fortran-lang/stdlib/issues/305#issuecomment-763528813, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAN6YRYRPP7GQCGXERJBYMLS222W7ANCNFSM4WJL6ZSQ .

LucaArgentiUCF commented 3 years ago

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 commented 3 years ago

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 commented 3 years ago

The Netlib site is a trifle vague, On their FAQ page:

So, unless code is associated with a TOMS algorithm, it is likely to come with a very permissive license.

Op wo 20 jan. 2021 om 15:23 schreef Ivan Pribec notifications@github.com:

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.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/fortran-lang/stdlib/issues/305#issuecomment-763638421, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAN6YR34HWOOZMQS2PVEBD3S23RO3ANCNFSM4WJL6ZSQ .

ivan-pi commented 3 years ago
  • 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 commented 3 years ago

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

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

ivan-pi commented 2 years ago

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 commented 2 years ago

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