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.

Sign up for GitHub

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

Jump to bottom

Miscellaneous mathematical (non-special) functions #150

Open
nshaffer opened this issue Feb 21, 2020 · 36 comments
Open

Miscellaneous mathematical (non-special) functions #150

nshaffer opened this issue Feb 21, 2020 · 36 comments
Labels
topic: mathematics linear algebra, sparse matrices, special functions, FFT, random numbers, statistics, ...

Comments

@nshaffer
Copy link
Contributor

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
Copy link
Member

certik commented Feb 21, 2020

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
Copy link
Contributor Author

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
Copy link

epagone commented Feb 25, 2020

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

@jvdp1
Copy link
Member

jvdp1 commented Mar 1, 2020

What about:

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

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

@leonfoks
Copy link

leonfoks commented Mar 1, 2020

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
Copy link
Member

ivan-pi commented Jun 23, 2020

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

@ivan-pi
Copy link
Member

ivan-pi commented Jun 23, 2020

  • 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
Copy link
Contributor Author

nshaffer commented Jun 23, 2020 •
edited
Loading

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 ivan-pi mentioned this issue Jun 23, 2020
Open
@ivan-pi
Copy link
Member

ivan-pi commented Jun 23, 2020

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
Copy link

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
Copy link
Member

certik commented Aug 18, 2021 •
edited
Loading

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!

@nshaffer nshaffer mentioned this issue Aug 19, 2021
Closed
@arjenmarkus
Copy link
Member

arjenmarkus commented Aug 20, 2021 via email

@Beliavsky
Copy link

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

@ivan-pi
Copy link
Member

ivan-pi commented Aug 27, 2021

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

I don't think the rlibm library is common enough to be able to distribute it easily. The source files are located here. A module containing interfaces (as an fpm package) would certainly be a good start and could later be repurposed to stdlib if interest is high.

Thanks for the great find!

@certik
Copy link
Member

certik commented Aug 27, 2021

@arjenmarkus thanks for the question. LFortran currently has a runtime module that just interfaces a libc version of the math routines:

This module should satisfy most contributors to that thread saying to just use whatever is already done. Then we started implementing a pure Fortran version in a separate module:

Right now only double precision sin is implemented. We tested the accuracy, it seems as accurate as the gfortran's version of sin on an interval like (-10, 10). It is less accurate for large numbers due to errors in argument reduction. We do have an accurate C implementation of the argument reduction operation in: https://gitlab.com/lfortran/lfortran/-/blob/06213849e4066eeed4dc0f15d6ceefbd5a5b64a1/src/runtime/impure/lfortran_intrinsic_sin_c.c, but that's obviously not pure Fortran. So we have all the building blocks. As a start, we should implement highly accurate kernels of all functions, and use a pure Fortran reduction operation. I suspect some functions might be accurate enough. Some functions like sin might require special handling for large arguments, perhaps calling into the "impure" C version of the reduction, or reimplementing it in Fortran. If you want to help with any of this, please let me know! I would be happy to get you started. You don't need to know anything about LFortran --- this is pure Fortran, so you can just use gfortran to compile.

Right now thing are hardwired, but we will allow to select different versions from a command line.

@Beliavsky
Copy link

Beliavsky commented Aug 29, 2021 •
edited
Loading

It is common to compute the nth root of a real number. It would be nice if one could write x**(1/n), but that does not work. Instead one writes x**(1.0d0/n) or exp(log(x)/n). The latter does not look the formula in a scientific paper. Should stdlib add a function named something like nth_root(x,n) where x is real and n > 0? It returns x for n = 1, sqrt(x) for n = 2, and cubrt(x) for n = 3. For large n it uses exp(log(x)/n) but for example for n = 4 it can use sqrt(sqrt(x)) if that is faster. Maybe nth_root should be generalized to say power(x,i,j) which raises x to i/j.

@epagone
Copy link

epagone commented Aug 29, 2021

I like the idea of nth_root(x,n). I have been on the verge of proposing it myself a couple of times. Related discussion: #214.

@arjenmarkus
Copy link
Member

arjenmarkus commented Aug 30, 2021 via email

@ivan-pi
Copy link
Member

ivan-pi commented Aug 30, 2021

It is common to compute the nth root of a real number. It would be nice if one could write x**(1/n), but that does not work. Instead one writes x**(1.0d0/n) or exp(log(x)/n). The latter does not look the formula in a scientific paper. Should stdlib add a function named something like nth_root(x,n) where x is real and n > 0? It returns x for n = 1, sqrt(x) for n = 2, and cubrt(x) for n = 3. For large n it uses exp(log(x)/n) but for example for n = 4 it can use sqrt(sqrt(x)) if that is faster. Maybe nth_root should be generalized to say power(x,i,j) which raises x to i/j.

I think a nthroot function would be in scope of stdlib. MATLAB has one limited to real roots. If you read the discussion in #214 it becomes partially clear why. Since there can be multiple complex roots it becomes difficult to settle on a convention, i.e. which complex root should the function return.

Concerning implementation, I would generally avoid trying to outsmart the compiler. For the expression x**(1./3.) some compilers already call an internal cbrt function. Perhaps x**(1.0/4.0) already generates two square root calls?

@certik
Copy link
Member

certik commented Sep 1, 2021

yes, I am interested to help out. It is
an exercise in practical mathematics :). (FYI: I have experience with
implementing special functions in a different language)

Perfect! You can simply try to implement some of the trigonometric functions in pure Fortran, and compare accuracy and speed with the default ones. Here is how we have done the sin function in pure Fortran: https://gitlab.com/lfortran/lfortran/-/blob/92e11002b133d909f8429bfe75d57ec554f6e557/src/runtime/pure/lfortran_intrinsic_trig.f90#L59

@arjenmarkus
Copy link
Member

arjenmarkus commented Sep 2, 2021 via email

@arjenmarkus
Copy link
Member

arjenmarkus commented Sep 3, 2021 via email

@certik
Copy link
Member

certik commented Sep 5, 2021

The sin implementation that I linked to is 1e-16 accurate based on our measurements. Ultimately there are two approaches, so we should have two versions;

  • The best possible accuracy across the whole range of single or double precision numbers; as performing as one can get withing this constrain
  • The fastest speed possible, as accurate as one can get without sacrificing speed

@arjenmarkus
Copy link
Member

arjenmarkus commented Sep 8, 2021 via email

@certik
Copy link
Member

certik commented Sep 8, 2021

On a related note, I have figured out how to create fast argument reduction, 1e-16 accurate up to x=1e16 using the double double approach (emulating quadruple precision using doubles, but it's very simple and no branching, thus should vectorize well). Beyond 1e16 I have to implement the general (slow) algorithm, so that we can have everything in pure Fortran. I've been struggling to figure out applications in that range, since the distance between the closest floating point numbers around x=1e16 is exactly 2.0, and it only gets worse for higher x. So that's roughly 3 numbers per period of sin. I can't figure out any application where I would want such a sparse grid. It seems the floating point around 1e16 are just not very usable for numerical calculations. I asked here with the same question:

So far nobody gave me any application for this.

@arjenmarkus
Copy link
Member

arjenmarkus commented Sep 9, 2021 via email

@arjenmarkus
Copy link
Member

arjenmarkus commented Sep 13, 2021 via email

@certik
Copy link
Member

certik commented Sep 15, 2021 •
edited
Loading

Awesome, thanks @arjenmarkus ! I tried to edit your comments to fold the long code, but GitHub does not allow that.

We can either put these into LFortran, or we can start a new separate repository. Starting a separate repository would be better, as we can then more easily maintain it as a community, have multiple versions (with different choices of speed vs accuracy), infrastructure to do accuracy plots, benchmarks, etc. Then LFortran (and hopefully other compilers too!) can just copy the version that they like into their runtime library.

I was hoping these would be ideally maintained by fortran-lang itself, but given the (hopefully only temporary) opposition at https://fortran-lang.discourse.group/t/fortran-runtime-math-library/755, we can maybe start this repository separately, and then move it under fortran-lang later.

Do you want me to start a repository under my name at github.com/certik ?

@arjenmarkus
Copy link
Member

arjenmarkus commented Sep 15, 2021 via email

@awvwgk awvwgk added the topic: mathematics linear algebra, sparse matrices, special functions, FFT, random numbers, statistics, ... label Sep 18, 2021
@Beliavsky
Copy link

Should there be a function (subroutine?) to compute the sine and cosine of the same argument together?

In the context of generating normal variates using Box-Muller, Ron Shephard said on Fortran Discourse

"Many math libraries have a subroutine that computes both sine and cosine values of the argument. In addition to the range reduction, some of the polynomial evaluation steps can be shared, making the total effort less than the sum of the effort for the two separate evaluations."

@urbanjost
Copy link

It is indeed very common for expressions requiring a sine of a value to also require a cosine. Using a subroutine that calculates both would require a separate call, and IMO generally be less readable though. "X=SIN(T)+3.0*COS(T)" IMO opinion is clearer
than

call trig(t,s,c)
x=s+3.0*c

If tables and interpolation are being used I would guess the computational efficiency would not change much; but if polynomial expansion or something else was used it might be more significant. Given that trig. functions are often highly optimized I think i would want to see some numbers. I am wondering what languages have the routine and if anyone has programs using it. I can imagine some HPC apps where it probably would add up; but personally when profiling codes I am involved with I cannot remember trig functions showing up as a major player in any of them although many of them call trig. functions throughout.
I could imagine that being different in other languages.
Anyone have any hard numbers?

@urbanjost
Copy link

Seems to usually be called "sincos"; I was toying with the idea that Fortran sort of already had this with complex
numbers since e**(iT)=cos(T)+isin(T) for a limited range of T and did not get at all what I expected, so need to look at this later; but
I was hoping I would make a fast sincos() in just a few lines using complex values, and instead found it to be slower and in some ranges got odd results.

program main
!   e ^(i · T) = cos(T) + i · sin(T).
use,intrinsic :: iso_fortran_env, only : int8, int16, int32, int64, real32, real64, real128
implicit none
complex,parameter   :: i=sqrt((-1,0))
real(kind=real64)   :: time_begin, time_end
complex             :: sc
real                :: T
! "e" is the base of the natural logarithm system; named in honor of Euler, but known as Napier's constant.
real,parameter :: e = 2.71828182845904523536028747135266249775724709369995_real64
real,parameter :: pi = 3.141592653589793238462643383279502884197169399375105820974944592307_real64
real,parameter :: step=epsilon(0.0)*4

write(*,*)'estimated passes=',2*pi/step

T=0.0
call CPU_TIME(time_begin)
do while (T <= 2*pi)
   sc=e**(i*T)
   T=T+step
enddo
write(*,*)T,sc
call CPU_TIME (time_end)
print *, 'Elapsed time is ', time_end - time_begin

T=0.0
call CPU_TIME(time_begin)
do while (T <= 2*pi)
   sc=cmplx(cos(T),sin(T))
   T=T+step
enddo
write(*,*)T,sc
call CPU_TIME (time_end)
print *, 'Elapsed time is ', time_end - time_begin

end program main

@Romendakil
Copy link

Jjust a stupid question, why do you define the complex unit this way:

implicit none
complex,parameter   :: i=sqrt((-1,0))

and not instead

complex, parameter :: i = cmplx (0, 1, kind=<your kind>)

?

@urbanjost
Copy link

No reason it would be any better, I like yours. I was just dashing something together to see if Fortran basically already had a "sincos" function because it supports complex variables, which a lot of languages do not; and how it would perform. It surprised me by being faster than a sin and cos call over some ranges but overall being slower and such. Wondering if anyone else has comments on going down that path, etc.

@certik
Copy link
Member

certik commented Jun 13, 2022 •
edited
Loading

I think for sin/cos in particular, compiles are usually clever enough to internally convert code like:

a = sin(x)
b = cos(x)

to something like:

call sincos(x, a, b)

We'll definitely do this and similar optimizations in LFortran in the future.

But one still has to have a fast / vectorizable sincos implementation that the compiler would call. So adding it to stdlib would be nice.

@Beliavsky
Copy link

Beliavsky commented Jul 31, 2023 •
edited
Loading

Rational powers appear fairly often in programs, and a recent preprint Generalising the Fast Reciprocal Square Root Algorithm discusses how to compute them quickly. They could be a candidate for stdlib.

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

No branches or pull requests