Miscellaneous mathematical (non-special) functions

[nshaffer]

I propose that we provide a module with miscellaneous mathematical functions, those which feel like they could/should be intrinsics, but are not for whatever reason. From my own experience, here are functions which I have felt the need to implement myself several times:

• expm1(x) and log1p(x) - I suspect that many compilers detect exp(x) - 1.0/log(1.0 + x) and call appropriate library routines, but I'd like to provide these explicitly.
• signum(x) - The intrinsic sign function is awkward to use when you just want to know the sign of a number.
• cbrt(x) - Again, I suspect many compilers detect x**(1.0/3.0) and call into their math library's cbrt. However, real exponents are ugly and error-prone and would like a more semantic replacement for this common case.
• sinc(x) - I am tired of writing my own implementation of this, and I suspect others are as well, even if it is trivial.
• phase(z) - Equivalent to atan2(z%im, z%re)

Let's brainstorm further functionality that you'd like to see.

[certik]

I agree, these would be all useful. Compilers should detect exp(x) - 1 and log(1 + x) and do the right thing, but it cannot hurt to have those functions explicitly.

Actually, sinc is not completely trivial. Here is the implementation that I use:

real(dp) elemental function sinc(x) result(r)
real(dp), intent(in) :: x
if (abs(x) < 1e-8_dp) then
    r = 1
else
    r = sin(x)/x
end if
end function

For high accuracy once must fine tune the cutoff, I ended up with 1e-8, but we would need to test this and ensure we are getting ~1e-15 accuracy for all x.

[nshaffer]

Some more worth considering:

• Functions to evaluate z*log(z) and log(z)/z which give the correct behavior near z=0. Not sure of good names.

Some with a distinctly statistical flavor, which maybe belong in the stats module:

• logit(x) for evaluating log(x)/log(1-x)
• logistic(x) for evaluating exp(x)/[1 + exp(x)] (the function whose inverse is logit)
• logsumexp(x) for evaluating log(sum(exp(x))) with x being an array.

[epagone]

Would the factorial and the binomial coefficient belong to this too?

[jvdp1]

What about:

• Cumulative product (cumprod)
• Cumulative sum (cumsum)

They would be probably in another module since they would not be elemental.

[leonfoks]

If a ”math” module is reserved for the elemental functions, Perhaps cumsum cumprod etc. Could go in a “numerical” module?

Since we are talking about summation, do we want to provide more robust implementations of sum and cumsum? The Kahan summation (and variants) spring to mind. I can reserve this discussion for the actual cumsum proposal too if need be. Related to #134 about multiple implementations of a concept.

[ivan-pi]

Do we have a name for this module already? Would stdlib_<experimental>_math be too broad?

[ivan-pi]

• expm1(x) and log1p(x) - I suspect that many compilers detect exp(x) - 1.0/log(1.0 + x) and call appropriate library routines, but I'd like to provide these explicitly.

These functions are available in the NSWC Library under the following names:

• REXP/DREXP for computing exp(x) - 1
• REXP1/DREXP1 for computing (exp(x) - 1)/x
• ALNREL/DLNREL for computing log(1 + a) for a > -1
• RLOG/DRLOG for computing x - 1 - log(x) for x > 0
• RLOG1/DRLOG1 for computing x - log(1 + x) for x > -1

The CBRT/DCBRT function is also available.

[nshaffer]

Yeah, NSWC has a lot of nice functionality, as long as you don't need complex numbers. But their implementations may still be useful reference material, especially for things like error tolerances for Taylor series approximations and the like.

As for the name, stdlib_experimental_math seems fine. Something more specific like mathutils also seems suitable. If we want to roll these in with special functions, that's fine too.

A good next step would be to write the markdown API for some obviously useful functions and implement the interfaces for them in a module (to be implemented later in submodules). I would take this up myself, but I have a new baby at home.

[ivan-pi]

A good next step would be to write the markdown API for some obviously useful functions and implement the interfaces for them in a module (to be implemented later in submodules). I would take this up myself, but I have a new baby at home.

I've already started with the cube root function, see #214. Congrats! 🎉

[Beliavsky]

A Hacker News thread discusses a post Speeding up atan2f by 50x, which provides C code. Are there Fortran programs where speeding up atan2 would be an important improvement?

[certik]

We started implementing pure Fortran versions of such functions as sin in LFortran: https://gitlab.com/lfortran/lfortran/-/blob/7e067df8f684f9762cd7e67c6a4bdd081e2971f9/src/runtime/pure/lfortran_intrinsic_trig.f90#L70

We plan to implement atan2 similarly as in the article @Beliavsky mentioned.

I suggested here to collaborate on such pure Fortran implementations: https://fortran-lang.discourse.group/t/fortran-runtime-math-library/755, unfortunately it was not received well. But we are doing exactly what I proposed there anyway and we welcome any collaborators!

[arjenmarkus]

Ondrej, do you have a list of functions that still need to be implemented? While I agree with the contributors to that thread, I can understand your point of view as well. I may be able to contribute. Op wo 18 aug. 2021 om 22:25 schreef Ondřej Čertík ***@***.***

…

: We started implementing pure Fortran versions of such functions as sin in LFortran: https://gitlab.com/lfortran/lfortran/-/blob/7e067df8f684f9762cd7e67c6a4bdd081e2971f9/src/runtime/pure/lfortran_intrinsic_trig.f90#L70 We plan to implement atan2 similarly as in the article @Beliavsky <https://github.com/Beliavsky> mentioned. I suggested here to collaborate on such pure Fortran implementations: https://fortran-lang.discourse.group/t/fortran-runtime-math-library/755, unfortunately it was not received well. But I am doing exactly what I proposed there anyway and I welcome any collaborators! — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#150 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAN6YR6PLCI5TBR7QNKZWM3T5QJKLANCNFSM4KZEILLA> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&utm_campaign=notification-email> .

[Beliavsky]

Should High Performance Correctly Rounded Math Libraries for 32-Bit Floating Point (in C) be callable from stdlib?