search

fortran-lang / stdlib

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

Miscellaneous mathematical (non-special) functions #150

Open nshaffer opened 4 years ago

nshaffer commented 4 years ago

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:

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

certik commented 4 years ago

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

Some more worth considering:

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

epagone commented 4 years ago

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

jvdp1 commented 4 years ago

What about:

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

leonfoks commented 4 years ago

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

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

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

The CBRT/DCBRT function is also available.

nshaffer commented 4 years ago

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

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

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

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

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 https://github.com/fortran-lang/stdlib/issues/150#issuecomment-901406807, 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 commented 3 years ago

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

ivan-pi commented 3 years ago

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

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

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

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

@Ondřej Čertík @.***>, 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)

Op vr 27 aug. 2021 om 20:30 schreef Ondřej Čertík @.***

:

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

- https://gitlab.com/lfortran/lfortran/-/blob/06213849e4066eeed4dc0f15d6ceefbd5a5b64a1/src/runtime/impure/lfortran_intrinsic_math.f90

https://gitlab.com/lfortran/lfortran/-/blob/06213849e4066eeed4dc0f15d6ceefbd5a5b64a1/src/runtime/impure/lfortran_intrinsics.c

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:

- https://gitlab.com/lfortran/lfortran/-/blob/06213849e4066eeed4dc0f15d6ceefbd5a5b64a1/src/runtime/pure/lfortran_intrinsic_trig.f90#L59

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.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/fortran-lang/stdlib/issues/150#issuecomment-907393433, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAN6YRYBFRICDERXHMYEVGDT67KTLANCNFSM4KZEILLA . 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&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

ivan-pi commented 3 years ago

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

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

I have started looking into the relevant algorithms, using a book on the subject, The Mathematical Function Computation Handbook by Nelson Beebe. While that puts a lot of emphasis on getting results accurate in various floating-point formats, down to 1 or even 0.5 ulp, it is useful to get a feeling for what is involved. My first version will certainly not strive to achieve that kind of accuracy ;).

Op do 2 sep. 2021 om 01:00 schreef Ondřej Čertík @.***>:

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

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/fortran-lang/stdlib/issues/150#issuecomment-910866783, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAN6YR2XLW7ZKCZJ5Y6WZQ3T72WBFANCNFSM4KZEILLA . 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&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

arjenmarkus commented 3 years ago

Just a quick update: I implemented cosine and tangent in the straightforward, if not naïve, way. As expected the errors are a trifle larger for cosine than for sine (as I use the relation cos(x) = sin(pi/2-x)) and considerably larger for tan (calculated as sin(x)/cos(x)). I think the computation of the tangent function must be improved, although it is not a bound function in contrast to sine and cosine, so the large errors may be relatively small ;).

Op do 2 sep. 2021 om 08:34 schreef Arjen Markus @.***>:

I have started looking into the relevant algorithms, using a book on the subject, The Mathematical Function Computation Handbook by Nelson Beebe. While that puts a lot of emphasis on getting results accurate in various floating-point formats, down to 1 or even 0.5 ulp, it is useful to get a feeling for what is involved. My first version will certainly not strive to achieve that kind of accuracy ;).

Op do 2 sep. 2021 om 01:00 schreef Ondřej Čertík @.***

:

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

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/fortran-lang/stdlib/issues/150#issuecomment-910866783, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAN6YR2XLW7ZKCZJ5Y6WZQ3T72WBFANCNFSM4KZEILLA . 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&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

certik commented 3 years ago

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;

arjenmarkus commented 3 years ago

Right, I will concentrate on a fast implementation for now. BTW, a quick examination of the tan() results shows that the largest differences occur where the function itself is large. The relative error is not more than 1.0e-13.

Op zo 5 sep. 2021 om 05:39 schreef Ondřej Čertík @.***>:

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

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/fortran-lang/stdlib/issues/150#issuecomment-913080008, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAN6YR2PVUYSQKFT72K5ER3UALQ6ZANCNFSM4KZEILLA . 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&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

certik commented 3 years ago

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

Hi Ondrej,

I have completed the implementation of cos, tan, asin, acos and atan. See the attached code. I only tested them with the simple and rough test program, but it is looking quite good. (I based the inverse functions on code I found on Netlib). Next step: the exponentials, I'd say (I also have half an implementation of sqrt)

Op wo 8 sep. 2021 om 23:16 schreef Ondřej Čertík @.***>:

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:

- https://fortran-lang.discourse.group/t/what-are-some-applications-of-sin-x-for-x-1e16/1800

So far nobody gave me any application for this.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/fortran-lang/stdlib/issues/150#issuecomment-915577692, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAN6YR6JSIB5ILPKKWIJLSTUA7HDLANCNFSM4KZEILLA . 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&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

! test_intrinsic_trig.f90 -- ! Test program for the elementary trigonometric functions ! ! Note: rather naïve ! program test_intrinsic_trig use, intrinsic :: iso_fortran_env, only: sp => real32, dp => real64 use lfortran_intrinsic_trig, lf_sin => sin, lf_cos => cos, lf_tan => tan, lf_asin => asin, & lf_acos => acos, lf_atan => atan

implicit none

real(kind=dp)      :: dlfort, dnative, dr, ddiff, dneg_sin_diff,  dpos_sin_diff,  dabs_sin_diff,  &
                                                  dneg_cos_diff,  dpos_cos_diff,  dabs_cos_diff,  &
                                                  dneg_tan_diff,  dpos_tan_diff,  dabs_tan_diff,  &
                                                  dneg_asin_diff, dpos_asin_diff, dabs_asin_diff, &
                                                  dneg_acos_diff, dpos_acos_diff, dabs_acos_diff
real(kind=sp)      :: slfort, snative, sr, sdiff, sneg_sin_diff,  spos_sin_diff,  sabs_sin_diff,  &
                                                  sneg_cos_diff,  spos_cos_diff,  sabs_cos_diff,  &
                                                  sneg_tan_diff,  spos_tan_diff,  sabs_tan_diff,  &
                                                  sneg_asin_diff, spos_asin_diff, sabs_asin_diff, &
                                                  sneg_acos_diff, spos_acos_diff, sabs_acos_diff
real(kind=dp)      :: xd
real(kind=sp)      :: xs
integer            :: i
integer, parameter :: notest = 1000

real(kind=dp), parameter :: range = 20.0

dneg_sin_diff = 0.0_dp ; dpos_sin_diff = 0.0_dp
sneg_sin_diff = 0.0_sp ; spos_sin_diff = 0.0_sp
dneg_cos_diff = 0.0_dp ; dpos_cos_diff = 0.0_dp
sneg_cos_diff = 0.0_sp ; spos_cos_diff = 0.0_sp
dneg_tan_diff = 0.0_dp ; dpos_tan_diff = 0.0_dp
sneg_tan_diff = 0.0_sp ; spos_tan_diff = 0.0_sp

dabs_sin_diff = 0.0_dp ; dabs_sin_diff = 0.0_sp
dabs_cos_diff = 0.0_dp ; dabs_cos_diff = 0.0_sp
dabs_tan_diff = 0.0_dp ; dabs_tan_diff = 0.0_sp

do i = 1,notest
    call random_number( dr ); dr = 2.0_dp * range          * ( dr - 0.5_dp )
    call random_number( sr ); sr = 2.0_sp * real(range,sp) * ( sr - 0.5_sp )

    dlfort  = lf_sin( dr )
    slfort  = lf_sin( sr )

    dnative = sin( dr )
    snative = sin( sr )

    ddiff   = dlfort - dnative
    sdiff   = slfort - snative

    dneg_sin_diff = merge( ddiff, dneg_sin_diff, ddiff < dneg_sin_diff )
    sneg_sin_diff = merge( sdiff, sneg_sin_diff, sdiff < sneg_sin_diff )
    dpos_sin_diff = merge( ddiff, dpos_sin_diff, ddiff > dpos_sin_diff )
    spos_sin_diff = merge( sdiff, spos_sin_diff, sdiff > spos_sin_diff )
    dabs_sin_diff = dabs_sin_diff + abs(ddiff)
    sabs_sin_diff = sabs_sin_diff + abs(sdiff)

    dlfort  = lf_cos( dr )
    slfort  = lf_cos( sr )

    dnative = cos( dr )
    snative = cos( sr )

    ddiff   = dlfort - dnative
    sdiff   = slfort - snative

    dneg_cos_diff = merge( ddiff, dneg_cos_diff, ddiff < dneg_cos_diff )
    sneg_cos_diff = merge( sdiff, sneg_cos_diff, sdiff < sneg_cos_diff )
    dpos_cos_diff = merge( ddiff, dpos_cos_diff, ddiff > dpos_cos_diff )
    spos_cos_diff = merge( sdiff, spos_cos_diff, sdiff > spos_cos_diff )
    dabs_cos_diff = dabs_cos_diff + abs(ddiff)
    sabs_cos_diff = sabs_cos_diff + abs(sdiff)

    dlfort  = lf_tan( dr )
    slfort  = lf_tan( sr )

    dnative = tan( dr )
    snative = tan( sr )

    ddiff   = dlfort - dnative
    sdiff   = slfort - snative

    dneg_tan_diff = merge( ddiff, dneg_tan_diff, ddiff < dneg_tan_diff )
    sneg_tan_diff = merge( sdiff, sneg_tan_diff, sdiff < sneg_tan_diff )
    dpos_tan_diff = merge( ddiff, dpos_tan_diff, ddiff > dpos_tan_diff )
    spos_tan_diff = merge( sdiff, spos_tan_diff, sdiff > spos_tan_diff )
    dabs_tan_diff = dabs_tan_diff + abs(ddiff)
    sabs_tan_diff = sabs_tan_diff + abs(sdiff)
    write(88,'(10e15.6)') dlfort, dnative, ddiff, ddiff/dnative, slfort, snative, sdiff, sdiff/snative
enddo

dabs_sin_diff = dabs_sin_diff / notest
sabs_sin_diff = sabs_sin_diff / notest
dabs_cos_diff = dabs_cos_diff / notest
sabs_cos_diff = sabs_cos_diff / notest
dabs_tan_diff = dabs_tan_diff / notest
sabs_tan_diff = sabs_tan_diff / notest

write( *, '(a,g0.2,a,g0.2,a)' ) 'Result for sin and cos in the range [',-range,',',range,']'
write( *, '(a)' )             'Maximum negative difference:'
write( *, '(a,g16.5)' )       '    sine (double precision):   ', dneg_sin_diff
write( *, '(a,g16.5)' )       '    cosine (double precision): ', dneg_cos_diff
write( *, '(a,g16.5)' )       '    tangent (double precision):', dneg_tan_diff
write( *, '(a,g16.5)' )       '    sine (single precision):   ', sneg_sin_diff
write( *, '(a,g16.5)' )       '    cosine (single precision): ', sneg_cos_diff
write( *, '(a,g16.5)' )       '    tangent (single precision):', sneg_tan_diff
write( *, '(a)' )             'Maximum positive difference:'
write( *, '(a,g16.5)' )       '    sine (double precision):   ', dpos_sin_diff
write( *, '(a,g16.5)' )       '    cosine (double precision): ', dpos_cos_diff
write( *, '(a,g16.5)' )       '    tangent (double precision):', dpos_tan_diff
write( *, '(a,g16.5)' )       '    sine (single precision):   ', spos_sin_diff
write( *, '(a,g16.5)' )       '    cosine (single precision): ', spos_cos_diff
write( *, '(a,g16.5)' )       '    tangent (single precision):', spos_tan_diff
write( *, '(a)' )             'Mean absolute difference:'
write( *, '(a,g16.5)' )       '    sine (double precision):   ', dabs_sin_diff
write( *, '(a,g16.5)' )       '    cosine (double precision): ', dabs_cos_diff
write( *, '(a,g16.5)' )       '    tangent (double precision):', dabs_tan_diff
write( *, '(a,g16.5)' )       '    sine (single precision):   ', sabs_sin_diff
write( *, '(a,g16.5)' )       '    cosine (single precision): ', sabs_cos_diff
write( *, '(a,g16.5)' )       '    tangent (single precision):', sabs_tan_diff

!
! Quick and dirty ...
!
write(*,*) 'Parameters ...'
xd = lf_asin(0.5_dp) ! Print the number of terms
xd = lf_atan(0.5_dp)

write( *, '(a)') 'asin: double precision'
do i = -10,10
    write(*,*) lf_asin(0.1_dp*i), asin(0.1_dp*i), lf_asin(0.1_dp*i) - asin(0.1_dp*i)
enddo
write( *, '(a)') 'asin: single precision'
do i = -10,10
    write(*,*) lf_asin(0.1_sp*i), asin(0.1_sp*i), lf_asin(0.1_sp*i) - asin(0.1_sp*i)
enddo

write( *, '(a)') 'asin: out of range argument'
write(*,*) lf_asin(1.1_dp), asin(1.1_dp), lf_asin(1.1_sp), asin(1.1_sp)

write( *, '(a)') 'acos: double precision'
do i = -10,10
    write(*,*) lf_acos(0.1_dp*i), acos(0.1_dp*i), lf_acos(0.1_dp*i) - acos(0.1_dp*i)
enddo
write( *, '(a)') 'acos: single precision'
do i = -10,10
    write(*,*) lf_acos(0.1_sp*i), acos(0.1_sp*i), lf_acos(0.1_sp*i) - acos(0.1_sp*i)
enddo
write( *, '(a)') 'atan: double precision'
do i = -10,10
    xd = tan(3.1415925_dp * i / 20.0_dp)
    write(*,*) xd, lf_atan(xd), atan(xd), lf_atan(xd) - atan(xd)
enddo
write( *, '(a)') 'atan: single precision'
do i = -10,10
    xs = tan(3.1415925_sp * i / 20.0_sp)
    write(*,*) xs, lf_atan(xs), atan(xs), lf_atan(xs) - atan(xs)
enddo

end program test_intrinsic_trig

module lfortran_intrinsic_trig use, intrinsic :: ieee_arithmetic implicit none private public :: sin, cos, tan, asin, acos, atan

integer, parameter :: sp = kind(1.0) integer, parameter :: dp = selected_real_kind(1+precision(1.0))

real(dp), parameter :: pi = 3.1415926535897932384626433832795_dp real(dp), parameter :: halfpi = pi / 2.0_dp

real(dp), parameter :: sqeps = sqrt( 6.0_dp * epsilon(1.0_dp) ) real(dp), parameter :: xbig = 1.0_dp / epsilon(1.0_dp)

interface sin module procedure ssin, dsin end interface

interface cos module procedure scos, dcos end interface

interface tan module procedure stan, dtan end interface

interface asin module procedure sasin, dasin end interface

interface acos module procedure sacos, dacos end interface

interface atan module procedure satan, datan end interface

contains

! sin --------------------

real(dp) function abs(x) result(r) real(dp), intent(in) :: x if (x >= 0) then r = x else r = -x end if end function

elemental integer function floor(x) result(r) real(dp), intent(in) :: x if (x >= 0) then r = x else r = x-1 end if end function

elemental real(dp) function modulo(x, y) result(r) real(dp), intent(in) :: x, y r = x-floor(x/y)*y end function

elemental real(dp) function min(x, y) result(r) real(dp), intent(in) :: x, y if (x < y) then r = x else r = y end if end function

elemental real(dp) function max(x, y) result(r) real(dp), intent(in) :: x, y if (x > y) then r = x else r = y end if end function

elemental real(dp) function dsin(x) result(r) real(dp), intent(in) :: x real(dp) :: y integer :: n y = modulo(x, 2*pi) y = min(y, pi - y) y = max(y, -pi - y) y = min(y, pi - y) r = kernel_dsin(y) end function

elemental real(dp) function dcos(x) result(r) real(dp), intent(in) :: x r = sin(halfpi-x) end function

elemental real(dp) function dtan(x) result(r) real(dp), intent(in) :: x r = sin(x)/sin(halfpi-x) end function

! Accurate on [-pi/2,pi/2] to about 1e-16 elemental real(dp) function kernel_dsin(x) result(res) real(dp), intent(in) :: x real(dp), parameter :: S1 = 0.9999999999999990771_dp real(dp), parameter :: S2 = -0.16666666666664811048_dp real(dp), parameter :: S3 = 8.333333333226519387e-3_dp real(dp), parameter :: S4 = -1.9841269813888534497e-4_dp real(dp), parameter :: S5 = 2.7557315514280769795e-6_dp real(dp), parameter :: S6 = -2.5051823583393710429e-8_dp real(dp), parameter :: S7 = 1.6046585911173017112e-10_dp real(dp), parameter :: S8 = -7.3572396558796051923e-13_dp real(dp) :: z z = xx res = x (S1+z(S2+z(S3+z(S4+z(S5+z(S6+z(S7+z*S8))))))) end function

elemental real(sp) function ssin(x) result(r) real(sp), intent(in) :: x real(sp) :: y real(dp) :: tmp, x2 x2 = x tmp = dsin(x2) r = tmp end function

elemental real(sp) function scos(x) result(r) real(sp), intent(in) :: x real(sp) :: y real(dp) :: tmp, x2 x2 = x tmp = dcos(x2) r = tmp end function

elemental real(sp) function stan(x) result(r) real(sp), intent(in) :: x real(sp) :: y real(dp) :: tmp, x2 x2 = x tmp = dtan(x2) r = tmp end function

real(dp) function dasin (x) real(dp), intent(in) :: x ! ! Transcribed from asin.f on netlib.org ! may 1980 edition. w. fullerton, c3, los alamos scientific lab. ! ! series for asin on the interval 0. to 5.00000d-01 ! with weighted error 1.60e-17 ! log weighted error 16.79 ! significant figures required 15.67 ! decimal places required 17.45 ! real(dp), dimension(20) :: asincs = [ & .10246391753227159e0_dp, & .054946487221245833e0_dp, & .004080630392544969e0_dp, & .000407890068546044e0_dp, & .000046985367432203e0_dp, & .000005880975813970e0_dp, & .000000777323124627e0_dp, & .000000106774233400e0_dp, & .000000015092399536e0_dp, & .000000002180972408e0_dp, & .000000000320759842e0_dp, & .000000000047855369e0_dp, & .000000000007225128e0_dp, & .000000000001101833e0_dp, & .000000000000169476e0_dp, & .000000000000026261e0_dp, & .000000000000004095e0_dp, & .000000000000000642e0_dp, & .000000000000000101e0_dp, & .000000000000000016e0_dp ]

integer, parameter :: nterms = size(asincs) real(dp) :: z, y

if ( x < -1.0_dp .or. x > 1.0_dp ) then dasin = ieee_value( x, ieee_quiet_nan ) else y = abs(x)

z = 0.
if ( y > sqeps ) z = y*y
if ( z <= 0.5_dp ) dasin = x*(1.0_dp + csevl (4.0_dp*z-1.0_dp, asincs, nterms))
if ( z >  0.5_dp ) dasin = halfpi - sqrt(1.0_dp-z)*(1.0_dp + csevl (3.0_dp-4.0_dp*z, asincs, nterms))
if ( x /= 0.0_dp ) dasin = sign (dasin, x)

endif end function dasin

real(sp) function sasin (x) real(sp), intent(in) :: x

real(dp) :: tmp, x2 x2 = x tmp = dasin(x2) sasin = tmp end function sasin

real(dp) function dacos (x) real(dp), intent(in) :: x

dacos = halfpi - dasin(x) end function dacos

real(sp) function sacos (x) real(sp), intent(in) :: x

real(dp) :: tmp, x2 x2 = x tmp = dacos(x2) sacos = tmp end function sacos

real(dp) function datan (x) real(dp), intent(in) :: x

! ! Transcribed from asin.f on netlib.org ! jan 1978 edition. w. fullerton, c3, los alamos scientific lab. ! ! series for atan on the interval 0. to 4.00000d-02 ! with weighted error 1.00e-17 ! log weighted error 17.00 ! significant figures required 16.38 ! decimal places required 17.48 ! real(dp), dimension(9), parameter :: atancs = [ & .48690110349241406e0_dp, & -.006510831636717464e0_dp, & .000038345828265245e0_dp, & -.000000268722128762e0_dp, & .000000002050093098e0_dp, & -.000000000016450717e0_dp, & .000000000000136509e0_dp, & -.000000000000001160e0_dp, & .000000000000000010e0_dp ]

! tanp8(n) = tan(n*pi/8.) real(dp), dimension(3), parameter :: tanp8 = [ & .414213562373095048e+0_dp, & 1.0e0_dp, & 2.41421356237309504e+0_dp ]

! xbndn = tan((2n-1)pi/16.0) real(dp), parameter :: xbnd1 = +.198912367379658006e+0_dp real(dp), parameter :: xbnd2 = +.668178637919298919e+0_dp real(dp), parameter :: xbnd3 = +1.49660576266548901e+0_dp real(dp), parameter :: xbnd4 = +5.02733949212584810e+0_dp

! conpi8(n) + pi8(n) = n*pi/8.0 real(dp), dimension(4), parameter :: conpi8 = [ & 0.375e0_dp, & 0.75e0_dp, & 1.125e0_dp, & 1.5e0_dp ] real(dp), dimension(4), parameter :: pi8 = [ & +.176990816987241548e-1_dp, & +.353981633974483096e-1_dp, & +.530972450961724644e-1_dp, & 0.0707963267948966192e0_dp ]

integer, parameter :: nterms = size(atancs)

integer :: n real(dp) :: y, t

y = abs(x) if ( y <= xbnd1 ) then datan = x if ( y > sqeps ) datan = x(0.75_dp+csevl(50.0_dpy**2-1.0_dp, atancs, nterms)) else if ( y <= xbnd4 ) then n = 1 if ( y > xbnd2 ) n = 2 if ( y > xbnd3 ) n = 3

    t     = (y - tanp8(n)) / (1.0_dp + y*tanp8(n))
    datan = sign (conpi8(n) + (pi8(n) + t*(0.75_dp + csevl(50.0_dp*t**2-1.0_dp, atancs, nterms)) ), x)
else
    datan = conpi8(4) + pi8(4)
    if ( y < xbig ) datan = conpi8(4) + (pi8(4) - (0.75_dp + csevl (50.0_dp/y**2-1.0_dp, atancs, nterms))/y )
    datan = sign (datan, x)
endif

endif end function datan

real(sp) function satan (x) real(sp), intent(in) :: x

real(dp) :: tmp, x2 x2 = x tmp = datan(x2) satan = tmp end function satan

! ! Auxiliary functions for asin ! integer function inits (os, nos, eta) ! ! april 1977 version. w. fullerton, c3, los alamos scientific lab. ! ! initialize the orthogonal series so that inits is the number of terms ! needed to insure the error is no larger than eta. ordinarily, eta ! will be chosen to be one-tenth machine precision. ! ! input arguments -- ! os array of nos coefficients in an orthogonal series. ! nos number of coefficients in os. ! eta requested accuracy of series. ! real(dp), dimension(:), intent(in) :: os integer, intent(in) :: nos real(dp), intent(in) :: eta

real(dp)                           :: err
integer                            :: i

err = 0.
do i = nos,1,-1
  err = err + abs(os(i))
  if (err > eta) exit
enddo

inits = i

end function inits

real(dp) function csevl (x, cs, n) ! ! april 1977 version. w. fullerton, c3, los alamos scientific lab. ! ! evaluate the n-term chebyshev series cs at x. adapted from ! r. broucke, algorithm 446, c.a.c.m., 16, 254 (1973). also see fox ! and parker, chebyshev polys in numerical analysis, oxford press, p.56. ! ! x value at which the series is to be evaluated. ! cs array of n terms of a chebyshev series. in eval- ! uating cs, only half the first coef is summed. ! n number of terms in array cs. ! ! Note: in this version sanity checking has bee nremoved, as the ! function is only used within this module, hence we know ! the arguments are within range ! real(dp), intent(in) :: x real(dp), dimension(:), intent(in) :: cs integer, intent(in) :: n

integer :: i real(dp) :: b0, b1, b2, twox

b1 = 0.0_dp b0 = 0.0_dp twox = 2.0_dpx do i = n,1,-1 b2 = b1 b1 = b0 b0 = twoxb1 - b2 + cs(i) enddo

csevl = 0.5_dp * (b0-b2) end function csevl

end module

arjenmarkus commented 3 years ago

I have an update for the trig routines (now elemental) and a new set of routines for the exponential functions and their counterparts.

Note: because of restrictions on the functions you can use in specification expressions, I still need to fill in a few literal values, because now log() functions are used. These values mark the bounds of a few regimes and from a casual glance I conclude that they are not really critical, so that machine-independent values can be used (the bounds are actually dependent on epsilon(), tiny() and huge()).

Ondrej, okay to put them here or do you prefer another means?

Op do 9 sep. 2021 om 21:39 schreef Arjen Markus @.***>:

Hi Ondrej,

I have completed the implementation of cos, tan, asin, acos and atan. See the attached code. I only tested them with the simple and rough test program, but it is looking quite good. (I based the inverse functions on code I found on Netlib). Next step: the exponentials, I'd say (I also have half an implementation of sqrt)

Op wo 8 sep. 2021 om 23:16 schreef Ondřej Čertík @.***

:

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:

- https://fortran-lang.discourse.group/t/what-are-some-applications-of-sin-x-for-x-1e16/1800

So far nobody gave me any application for this.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/fortran-lang/stdlib/issues/150#issuecomment-915577692, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAN6YR6JSIB5ILPKKWIJLSTUA7HDLANCNFSM4KZEILLA . 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&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

! TODO: exp, log ! ! dlog not allowed in parameter definitions ! ! sqrt will have to be removed as well ! ! Do we want to use literal numbers instead? ! ! Carefully define the common parameters as module parameters, others locally

module lfortran_intrinsic_exp use, intrinsic :: ieee_arithmetic implicit none private public :: exp, cosh, sinh, tanh, acosh, asinh, atanh, log, log10

integer, parameter :: sp = kind(1.0) integer, parameter :: dp = selected_real_kind(1+precision(1.0))

real(dp), parameter :: xbig = 1.0_dp / epsilon(1.0_dp) real(dp), parameter :: sqeps = sqrt (3.0_dpepsilon(1.0_dp)) ! TODO: Make them local? real(dp), parameter :: xmax = -0.5_dplog(epsilon(1.0_dp)) ! TODO: literal? real(dp), parameter :: aln2 = 0.69314718055994530942e0_dp real(dp), parameter :: dxrel = sqrt (2.0_dp * epsilon(1.0_dp))

interface exp module procedure sexp, dexp end interface

interface sinh module procedure ssinh, dsinh end interface

interface cosh module procedure scosh, dcosh end interface

interface tanh module procedure stanh, dtanh end interface

interface asinh module procedure sasinh, dasinh end interface

interface acosh module procedure sacosh, dacosh end interface

interface atanh module procedure satanh, datanh end interface

interface log module procedure slog, dlog end interface

interface log10 module procedure slog10, dlog10 end interface

contains

elemental real(dp) function dexp(x) real(dp), intent(in) :: x ! ! may 1980 edition. w. fullerton c3, los alamos scientific lab. ! ! series for exp on the interval -1.00000d+00 to 1.00000d+00 ! with weighted error 4.81e-18 ! log weighted error 17.32 ! significant figures required 15.95 ! decimal places required 17.77 ! real(dp), parameter, dimension(8) :: expcs = [ & .086656949331498571e0_dp, & .000938494869299839e0_dp, & .000006776039709981e0_dp, & .000000036693120039e0_dp, & .000000000158959053e0_dp, & .000000000000573859e0_dp, & .000000000000001775e0_dp, & .000000000000000004e0_dp ] integer, parameter :: nterms = size(expcs) real(dp), parameter, dimension(17) :: twon16 = [ & 0.e0_dp, & .44273782427413840e-1_dp, & .90507732665257659e-1_dp, & .13878863475669165e0_dp, & .18920711500272107e0_dp, & .24185781207348405e0_dp, & .29683955465100967e0_dp, & .35425554693689273e0_dp, & .41421356237309505e0_dp, & .47682614593949931e0_dp, & .54221082540794082e0_dp, & .61049033194925431e0_dp, & .68179283050742909e0_dp, & .75625216037329948e0_dp, & .83400808640934246e0_dp, & .91520656139714729e0_dp, & 1.0e0_dp ] ! aln216 is 16.0/alog(2.) - 23.0 real(dp), parameter :: aln216 = .083120654223414518e0_dp real(dp), parameter :: xmin = log(tiny(1.0_dp)) + 0.01_dp ! TODO: literal? real(dp), parameter :: xmax = log(huge(1.0_dp)) - 0.001_dp ! TODO: literal? real(dp) :: xint, y, f integer :: n, n16, ndx

if (x >= xmin) then if (x > xmax) then dexp = ieee_value( dexp, ieee_positive_inf ) else xint = aint(x) y = x - xint

    y    = 23.0_dp*y + x*aln216
    n    = y
    f    = y - real(n,dp)
    n    = 23.0_dp*xint + real(n,dp)
    n16  = n/16
    if (n < 0) n16 = n16 - 1
    ndx = n - 16*n16 + 1

    dexp = 1.0_dp + (twon16(ndx) + f*(1.0+twon16(ndx)) * csevl (f, expcs, nterms) )
    dexp = set_exponent( dexp, n16+1)
endif

else dexp = 0.0_dp endif end function dexp

elemental real(dp) function dsinh(x) result(r) real(dp), intent(in) :: x real(dp) :: y if ( abs(x) > xbig ) then r = 0.5_dp exp(x) else y = exp(x) r = 0.5_dp (y - 1.0_dp/y) endif end function dsinh

elemental real(dp) function dcosh(x) result(r) real(dp), intent(in) :: x real(dp) :: y if ( abs(x) > xbig ) then r = 0.5_sp exp(x) else y = exp(x) r = 0.5_dp (y + 1.0_dp/y) endif end function dcosh

elemental real(dp) function dtanh (x) real(dp), intent(in) :: x ! ! april 1977 edition. w. fullerton, c3, los alamos scientific lab. ! ! series for tanh on the interval 0. to 1.00000d+00 ! with weighted error 9.88e-18 ! log weighted error 17.01 ! significant figures required 16.25 ! decimal places required 17.62 ! real(dp), dimension(17), parameter :: tanhcs = [ & -.25828756643634710e0_dp, & -.11836106330053497e0_dp, & .009869442648006398e0_dp, & -.000835798662344582e0_dp, & .000070904321198943e0_dp, & -.000006016424318120e0_dp, & .000000510524190800e0_dp, & -.000000043320729077e0_dp, & .000000003675999055e0_dp, & -.000000000311928496e0_dp, & .000000000026468828e0_dp, & -.000000000002246023e0_dp, & .000000000000190587e0_dp, & -.000000000000016172e0_dp, & .000000000000001372e0_dp, & -.000000000000000116e0_dp, & .000000000000000009e0_dp ] integer, parameter :: nterms = size(tanhcs) real(dp) :: y

y = abs(x) if (y <= 1.0_dp ) then dtanh = x if (y > sqeps) dtanh = x(1. + csevl (2.0_dpx*x-1.0_dp, tanhcs, nterms)) elseif (y <= xmax) then y = exp(y) dtanh = sign ((y-1.0_dp/y)/(y+1.0_dp/y), x) else dtanh = sign (1.0_dp, x) endif end function dtanh

elemental real(dp) function dasinh (x) real(dp), intent(in) :: x ! ! april 1977 edition. w. fullerton, c3, los alamos scientific lab. ! ! series for asnh on the interval 0. to 1.00000d+00 ! with weighted error 2.19e-17 ! log weighted error 16.66 ! significant figures required 15.60 ! decimal places required 17.31 ! real(dp), parameter, dimension(20) :: asnhcs = [ & -.12820039911738186e0_dp, & -.058811761189951768e0_dp, & .004727465432212481e0_dp, & -.000493836316265361e0_dp, & .000058506207058557e0_dp, & -.000007466998328931e0_dp, & .000001001169358355e0_dp, & -.000000139035438587e0_dp, & .000000019823169483e0_dp, & -.000000002884746841e0_dp, & .000000000426729654e0_dp, & -.000000000063976084e0_dp, & .000000000009699168e0_dp, & -.000000000001484427e0_dp, & .000000000000229037e0_dp, & -.000000000000035588e0_dp, & .000000000000005563e0_dp, & -.000000000000000874e0_dp, & .000000000000000138e0_dp, & -.000000000000000021e0_dp ] integer, parameter :: nterms = size(asnhcs) real(dp), parameter :: sqeps = sqrt(epsilon(1.0_dp)) ! Local definitions real(dp), parameter :: xmax = 1.0_dp/sqeps real(dp) :: y

y = abs(x) if (y <= 1.0) then dasinh = x if (y > sqeps) dasinh = x(1.0_dp + csevl (2.0_dpx2-1.0_dp, asnhcs,nterms)) else if (y < xmax) then dasinh = dlog (y + sqrt(y2+1.0_dp)) else dasinh = aln2 + dlog(y) endif dasinh = sign (dasinh, x) endif end function dasinh

elemental real(dp) function dacosh (x) real(dp), intent(in) :: x ! ! april 1977 edition. w. fullerton, c3, los alamos scientific lab. !

if (x < 1.0_dp) then dacosh = ieee_value( x, ieee_quiet_nan ) else if (x < xmax) then dacosh = dlog (x + sqrt(x**2-1.0_dp)) else dacosh = aln2 + dlog(x) endif endif end function dacosh

elemental real(dp) function datanh (x) real(dp), intent(in) :: x ! ! april 1977 edition. w. fullerton, c3, los alamos scientific lab. ! ! series for atnh on the interval 0. to 2.50000d-01 ! with weighted error 6.70e-18 ! log weighted error 17.17 ! significant figures required 16.01 ! decimal places required 17.76 ! real(dp), dimension(15), parameter :: atnhcs = [ & .094395102393195492e0_dp, & .049198437055786159e0_dp, & .002102593522455432e0_dp, & .000107355444977611e0_dp, & .000005978267249293e0_dp, & .000000350506203088e0_dp, & .000000021263743437e0_dp, & .000000001321694535e0_dp, & .000000000083658755e0_dp, & .000000000005370503e0_dp, & .000000000000348665e0_dp, & .000000000000022845e0_dp, & .000000000000001508e0_dp, & .000000000000000100e0_dp, & .000000000000000006e0_dp ] integer, parameter :: nterms = size(atnhcs) real(dp) :: y

y = abs(x) if (y >= 1.0_dp) then datanh = ieee_value( datanh, ieee_quiet_nan ) else datanh = x if (y > sqeps .and. y <= 0.5_dp) then datanh = x(1.0_dp + csevl (8.0_dpxx-1.0_dp, atnhcs, nterms)) elseif (y.gt.0.5) then datanh = 0.5_dpdlog((1.0_dp+x)/(1.0_dp-x)) endif endif end function datanh

elemental real(dp) function dlog (x) real(dp), intent(in) :: x ! ! june 1977 edition. w. fullerton, c3, los alamos scientific lab. ! ! series for aln on the interval 0. to 3.46021d-03 ! with weighted error 1.50e-16 ! log weighted error 15.82 ! significant figures required 15.65 ! decimal places required 16.21 ! real(dp), parameter, dimension(6) :: alncs = [ & 1.3347199877973882e0_dp, & .000693756283284112e0_dp, & .000000429340390204e0_dp, & .000000000289338477e0_dp, & .000000000000205125e0_dp, & .000000000000000150e0_dp ] integer, parameter :: nterms = size(alncs) real(dp), parameter, dimension(4) :: center = [ & 1.0_dp, & 1.25_dp, & 1.50_dp, & 1.75_dp ] real(dp), parameter, dimension(5) :: alncen = [ & 0.0e0_dp, & +.223143551314209755e+0_dp, & +.405465108108164381e+0_dp, & +.559615787935422686e+0_dp, & +.693147180559945309e+0_dp ] ! aln2 = alog(2.0) - 0.625 real(dp), parameter :: aln2 = 0.068147180559945309e0_dp

integer :: n, ntrval real(dp) :: xn, y, t, t2

  if ( x <= 0.0_dp ) then
      dlog = ieee_value( dlog, ieee_quiet_nan )
  else
      n  = exponent(x)
      y  = set_exponent(x,0)
      xn = n - 1
      y  = 2.0_dp*y
      ntrval = 4.0_dp*y - 2.5_dp
      if (ntrval == 5) t = ((y-1.0_dp)-1.0_dp) / (y+2.0_dp)
      if (ntrval < 5)  t = (y-center(ntrval))/(y+center(ntrval))
      t2 = t**2

      dlog = 0.625_dp*xn + (aln2*xn + alncen(ntrval) + 2.0_dp*t + t*t2*csevl(578.0_dp*t2-1.0_dp, alncs, nterms) )
  endif

end function dlog

elemental real(dp) function dlog10 (x) real(dp), intent(in) :: x ! ! april 1977 version. w. fullerton, c3, los alamos scientific lab. ! real(dp), parameter :: aloge = 0.43429448190325182765e0_dp

dlog10 = aloge * dlog(x) end function dlog10

elemental real(sp) function sexp(x) result(r) real(sp), intent(in) :: x real(dp) :: y, tmp y = x tmp = dexp(y) r = tmp end function sexp

elemental real(sp) function ssinh(x) result(r) real(sp), intent(in) :: x real(dp) :: y, tmp y = x tmp = dsinh(y) r = tmp end function ssinh

elemental real(sp) function scosh(x) result(r) real(sp), intent(in) :: x real(dp) :: y, tmp y = x tmp = dcosh(y) r = tmp end function scosh

elemental real(sp) function stanh(x) result(r) real(sp), intent(in) :: x real(dp) :: y, tmp y = x tmp = dtanh(y) r = tmp end function stanh

elemental real(sp) function sasinh(x) result(r) real(sp), intent(in) :: x real(dp) :: y, tmp y = x tmp = dasinh(y) r = tmp end function sasinh

elemental real(sp) function sacosh(x) result(r) real(sp), intent(in) :: x real(dp) :: y, tmp y = x tmp = dacosh(y) r = tmp end function sacosh

elemental real(sp) function satanh(x) result(r) real(sp), intent(in) :: x real(dp) :: y, tmp y = x tmp = datanh(y) r = tmp end function satanh

elemental real(sp) function slog(x) result(r) real(sp), intent(in) :: x real(dp) :: y, tmp y = x tmp = dlog(y) r = tmp end function slog

elemental real(sp) function slog10(x) result(r) real(sp), intent(in) :: x real(dp) :: y, tmp y = x tmp = dlog10(y) r = tmp end function slog10

pure real(dp) function csevl (x, cs, n) ! ! april 1977 version. w. fullerton, c3, los alamos scientific lab. ! ! evaluate the n-term chebyshev series cs at x. adapted from ! r. broucke, algorithm 446, c.a.c.m., 16, 254 (1973). also see fox ! and parker, chebyshev polys in numerical analysis, oxford press, p.56. ! ! x value at which the series is to be evaluated. ! cs array of n terms of a chebyshev series. in eval- ! uating cs, only half the first coef is summed. ! n number of terms in array cs. ! ! Note: in this version sanity checking has bee nremoved, as the ! function is only used within this module, hence we know ! the arguments are within range ! real(dp), intent(in) :: x real(dp), dimension(:), intent(in) :: cs integer, intent(in) :: n

integer :: i real(dp) :: b0, b1, b2, twox

b1 = 0.0_dp b0 = 0.0_dp twox = 2.0_dpx do i = n,1,-1 b2 = b1 b1 = b0 b0 = twoxb1 - b2 + cs(i) enddo

csevl = 0.5_dp * (b0-b2) end function csevl

end module

! test_intrinsic_trig.f90 -- ! Test program for the elementary trigonometric functions ! ! Note: rather naïve ! program test_intrinsic_trig use, intrinsic :: iso_fortran_env, only: sp => real32, dp => real64 use lfortran_intrinsic_trig, lf_sin => sin, lf_cos => cos, lf_tan => tan, lf_asin => asin, & lf_acos => acos, lf_atan => atan

implicit none

real(kind=dp)      :: dlfort, dnative, dr, ddiff, dneg_sin_diff,  dpos_sin_diff,  dabs_sin_diff,  &
                                                  dneg_cos_diff,  dpos_cos_diff,  dabs_cos_diff,  &
                                                  dneg_tan_diff,  dpos_tan_diff,  dabs_tan_diff,  &
                                                  dneg_asin_diff, dpos_asin_diff, dabs_asin_diff, &
                                                  dneg_acos_diff, dpos_acos_diff, dabs_acos_diff
real(kind=sp)      :: slfort, snative, sr, sdiff, sneg_sin_diff,  spos_sin_diff,  sabs_sin_diff,  &
                                                  sneg_cos_diff,  spos_cos_diff,  sabs_cos_diff,  &
                                                  sneg_tan_diff,  spos_tan_diff,  sabs_tan_diff,  &
                                                  sneg_asin_diff, spos_asin_diff, sabs_asin_diff, &
                                                  sneg_acos_diff, spos_acos_diff, sabs_acos_diff
real(kind=dp)      :: xd
real(kind=sp)      :: xs
integer            :: i
integer, parameter :: notest = 1000

real(kind=dp), parameter :: range = 20.0

dneg_sin_diff = 0.0_dp ; dpos_sin_diff = 0.0_dp
sneg_sin_diff = 0.0_sp ; spos_sin_diff = 0.0_sp
dneg_cos_diff = 0.0_dp ; dpos_cos_diff = 0.0_dp
sneg_cos_diff = 0.0_sp ; spos_cos_diff = 0.0_sp
dneg_tan_diff = 0.0_dp ; dpos_tan_diff = 0.0_dp
sneg_tan_diff = 0.0_sp ; spos_tan_diff = 0.0_sp

dabs_sin_diff = 0.0_dp ; dabs_sin_diff = 0.0_sp
dabs_cos_diff = 0.0_dp ; dabs_cos_diff = 0.0_sp
dabs_tan_diff = 0.0_dp ; dabs_tan_diff = 0.0_sp

do i = 1,notest
    call random_number( dr ); dr = 2.0_dp * range          * ( dr - 0.5_dp )
    call random_number( sr ); sr = 2.0_sp * real(range,sp) * ( sr - 0.5_sp )

    dlfort  = lf_sin( dr )
    slfort  = lf_sin( sr )

    dnative = sin( dr )
    snative = sin( sr )

    ddiff   = dlfort - dnative
    sdiff   = slfort - snative

    dneg_sin_diff = merge( ddiff, dneg_sin_diff, ddiff < dneg_sin_diff )
    sneg_sin_diff = merge( sdiff, sneg_sin_diff, sdiff < sneg_sin_diff )
    dpos_sin_diff = merge( ddiff, dpos_sin_diff, ddiff > dpos_sin_diff )
    spos_sin_diff = merge( sdiff, spos_sin_diff, sdiff > spos_sin_diff )
    dabs_sin_diff = dabs_sin_diff + abs(ddiff)
    sabs_sin_diff = sabs_sin_diff + abs(sdiff)

    dlfort  = lf_cos( dr )
    slfort  = lf_cos( sr )

    dnative = cos( dr )
    snative = cos( sr )

    ddiff   = dlfort - dnative
    sdiff   = slfort - snative

    dneg_cos_diff = merge( ddiff, dneg_cos_diff, ddiff < dneg_cos_diff )
    sneg_cos_diff = merge( sdiff, sneg_cos_diff, sdiff < sneg_cos_diff )
    dpos_cos_diff = merge( ddiff, dpos_cos_diff, ddiff > dpos_cos_diff )
    spos_cos_diff = merge( sdiff, spos_cos_diff, sdiff > spos_cos_diff )
    dabs_cos_diff = dabs_cos_diff + abs(ddiff)
    sabs_cos_diff = sabs_cos_diff + abs(sdiff)

    dlfort  = lf_tan( dr )
    slfort  = lf_tan( sr )

    dnative = tan( dr )
    snative = tan( sr )

    ddiff   = dlfort - dnative
    sdiff   = slfort - snative

    dneg_tan_diff = merge( ddiff, dneg_tan_diff, ddiff < dneg_tan_diff )
    sneg_tan_diff = merge( sdiff, sneg_tan_diff, sdiff < sneg_tan_diff )
    dpos_tan_diff = merge( ddiff, dpos_tan_diff, ddiff > dpos_tan_diff )
    spos_tan_diff = merge( sdiff, spos_tan_diff, sdiff > spos_tan_diff )
    dabs_tan_diff = dabs_tan_diff + abs(ddiff)
    sabs_tan_diff = sabs_tan_diff + abs(sdiff)
    write(88,'(10e15.6)') dlfort, dnative, ddiff, ddiff/dnative, slfort, snative, sdiff, sdiff/snative
enddo

dabs_sin_diff = dabs_sin_diff / notest
sabs_sin_diff = sabs_sin_diff / notest
dabs_cos_diff = dabs_cos_diff / notest
sabs_cos_diff = sabs_cos_diff / notest
dabs_tan_diff = dabs_tan_diff / notest
sabs_tan_diff = sabs_tan_diff / notest

write( *, '(a,g0.2,a,g0.2,a)' ) 'Result for sin and cos in the range [',-range,',',range,']'
write( *, '(a)' )             'Maximum negative difference:'
write( *, '(a,g16.5)' )       '    sine (double precision):   ', dneg_sin_diff
write( *, '(a,g16.5)' )       '    cosine (double precision): ', dneg_cos_diff
write( *, '(a,g16.5)' )       '    tangent (double precision):', dneg_tan_diff
write( *, '(a,g16.5)' )       '    sine (single precision):   ', sneg_sin_diff
write( *, '(a,g16.5)' )       '    cosine (single precision): ', sneg_cos_diff
write( *, '(a,g16.5)' )       '    tangent (single precision):', sneg_tan_diff
write( *, '(a)' )             'Maximum positive difference:'
write( *, '(a,g16.5)' )       '    sine (double precision):   ', dpos_sin_diff
write( *, '(a,g16.5)' )       '    cosine (double precision): ', dpos_cos_diff
write( *, '(a,g16.5)' )       '    tangent (double precision):', dpos_tan_diff
write( *, '(a,g16.5)' )       '    sine (single precision):   ', spos_sin_diff
write( *, '(a,g16.5)' )       '    cosine (single precision): ', spos_cos_diff
write( *, '(a,g16.5)' )       '    tangent (single precision):', spos_tan_diff
write( *, '(a)' )             'Mean absolute difference:'
write( *, '(a,g16.5)' )       '    sine (double precision):   ', dabs_sin_diff
write( *, '(a,g16.5)' )       '    cosine (double precision): ', dabs_cos_diff
write( *, '(a,g16.5)' )       '    tangent (double precision):', dabs_tan_diff
write( *, '(a,g16.5)' )       '    sine (single precision):   ', sabs_sin_diff
write( *, '(a,g16.5)' )       '    cosine (single precision): ', sabs_cos_diff
write( *, '(a,g16.5)' )       '    tangent (single precision):', sabs_tan_diff

!
! Quick and dirty ...
!
write(*,*) 'Parameters ...'
xd = lf_asin(0.5_dp) ! Print the number of terms
xd = lf_atan(0.5_dp)

write( *, '(a)') 'asin: double precision'
do i = -10,10
    write(*,*) lf_asin(0.1_dp*i), asin(0.1_dp*i), lf_asin(0.1_dp*i) - asin(0.1_dp*i)
enddo
write( *, '(a)') 'asin: single precision'
do i = -10,10
    write(*,*) lf_asin(0.1_sp*i), asin(0.1_sp*i), lf_asin(0.1_sp*i) - asin(0.1_sp*i)
enddo

write( *, '(a)') 'asin: out of range argument'
write(*,*) lf_asin(1.1_dp), asin(1.1_dp), lf_asin(1.1_sp), asin(1.1_sp)

write( *, '(a)') 'acos: double precision'
do i = -10,10
    write(*,*) lf_acos(0.1_dp*i), acos(0.1_dp*i), lf_acos(0.1_dp*i) - acos(0.1_dp*i)
enddo
write( *, '(a)') 'acos: single precision'
do i = -10,10
    write(*,*) lf_acos(0.1_sp*i), acos(0.1_sp*i), lf_acos(0.1_sp*i) - acos(0.1_sp*i)
enddo
write( *, '(a)') 'atan: double precision'
do i = -10,10
    xd = tan(3.1415925_dp * i / 20.0_dp)
    write(*,*) xd, lf_atan(xd), atan(xd), lf_atan(xd) - atan(xd)
enddo
write( *, '(a)') 'atan: single precision'
do i = -10,10
    xs = tan(3.1415925_sp * i / 20.0_sp)
    write(*,*) xs, lf_atan(xs), atan(xs), lf_atan(xs) - atan(xs)
enddo

end program test_intrinsic_trig

! test_intrinsic_exp.f90 -- ! Test program for the elementary trigonometric functions ! ! Note: rather naïve ! program test_intrinsic_exp use, intrinsic :: iso_fortran_env, only: sp => real32, dp => real64 use lfortran_intrinsic_exp, lf_sinh => sinh, lf_cosh => cosh, lf_tanh => tanh, lf_asinh => asinh, & lf_acosh => acosh, lf_atanh => atanh, lf_exp => exp

implicit none

real(kind=dp)      :: dlfort, dnative, dr, ddiff, dneg_sinh_diff,  dpos_sinh_diff,  dabs_sinh_diff,  &
                                                  dneg_cosh_diff,  dpos_cosh_diff,  dabs_cosh_diff,  &
                                                  dneg_tanh_diff,  dpos_tanh_diff,  dabs_tanh_diff,  &
                                                  dneg_asinh_diff, dpos_asinh_diff, dabs_asinh_diff, &
                                                  dneg_acosh_diff, dpos_acosh_diff, dabs_acosh_diff
real(kind=sp)      :: slfort, snative, sr, sdiff, sneg_sinh_diff,  spos_sinh_diff,  sabs_sinh_diff,  &
                                                  sneg_cosh_diff,  spos_cosh_diff,  sabs_cosh_diff,  &
                                                  sneg_tanh_diff,  spos_tanh_diff,  sabs_tanh_diff,  &
                                                  sneg_asinh_diff, spos_asinh_diff, sabs_asinh_diff, &
                                                  sneg_acosh_diff, spos_acosh_diff, sabs_acosh_diff
real(kind=dp)      :: xd
real(kind=sp)      :: xs
integer            :: i
integer, parameter :: notest = 1000

real(kind=dp), parameter :: range = 20.0

dneg_sinh_diff = 0.0_dp ; dpos_sinh_diff = 0.0_dp
sneg_sinh_diff = 0.0_sp ; spos_sinh_diff = 0.0_sp
dneg_cosh_diff = 0.0_dp ; dpos_cosh_diff = 0.0_dp
sneg_cosh_diff = 0.0_sp ; spos_cosh_diff = 0.0_sp
dneg_tanh_diff = 0.0_dp ; dpos_tanh_diff = 0.0_dp
sneg_tanh_diff = 0.0_sp ; spos_tanh_diff = 0.0_sp

dabs_sinh_diff = 0.0_dp ; dabs_sinh_diff = 0.0_sp
dabs_cosh_diff = 0.0_dp ; dabs_cosh_diff = 0.0_sp
dabs_tanh_diff = 0.0_dp ; dabs_tanh_diff = 0.0_sp

do i = 1,notest
    call random_number( dr ); dr = 2.0_dp * range          * ( dr - 0.5_dp )
    call random_number( sr ); sr = 2.0_sp * real(range,sp) * ( sr - 0.5_sp )

    dlfort  = lf_sinh( dr )
    slfort  = lf_sinh( sr )

    dnative = sinh( dr )
    snative = sinh( sr )

    ddiff   = dlfort - dnative
    sdiff   = slfort - snative

    dneg_sinh_diff = merge( ddiff, dneg_sinh_diff, ddiff < dneg_sinh_diff )
    sneg_sinh_diff = merge( sdiff, sneg_sinh_diff, sdiff < sneg_sinh_diff )
    dpos_sinh_diff = merge( ddiff, dpos_sinh_diff, ddiff > dpos_sinh_diff )
    spos_sinh_diff = merge( sdiff, spos_sinh_diff, sdiff > spos_sinh_diff )
    dabs_sinh_diff = dabs_sinh_diff + abs(ddiff)
    sabs_sinh_diff = sabs_sinh_diff + abs(sdiff)

    dlfort  = lf_cosh( dr )
    slfort  = lf_cosh( sr )

    dnative = cosh( dr )
    snative = cosh( sr )

    ddiff   = dlfort - dnative
    sdiff   = slfort - snative

    dneg_cosh_diff = merge( ddiff, dneg_cosh_diff, ddiff < dneg_cosh_diff )
    sneg_cosh_diff = merge( sdiff, sneg_cosh_diff, sdiff < sneg_cosh_diff )
    dpos_cosh_diff = merge( ddiff, dpos_cosh_diff, ddiff > dpos_cosh_diff )
    spos_cosh_diff = merge( sdiff, spos_cosh_diff, sdiff > spos_cosh_diff )
    dabs_cosh_diff = dabs_cosh_diff + abs(ddiff)
    sabs_cosh_diff = sabs_cosh_diff + abs(sdiff)

    dlfort  = lf_tanh( dr )
    slfort  = lf_tanh( sr )

    dnative = tanh( dr )
    snative = tanh( sr )

    ddiff   = dlfort - dnative
    sdiff   = slfort - snative

    dneg_tanh_diff = merge( ddiff, dneg_tanh_diff, ddiff < dneg_tanh_diff )
    sneg_tanh_diff = merge( sdiff, sneg_tanh_diff, sdiff < sneg_tanh_diff )
    dpos_tanh_diff = merge( ddiff, dpos_tanh_diff, ddiff > dpos_tanh_diff )
    spos_tanh_diff = merge( sdiff, spos_tanh_diff, sdiff > spos_tanh_diff )
    dabs_tanh_diff = dabs_tanh_diff + abs(ddiff)
    sabs_tanh_diff = sabs_tanh_diff + abs(sdiff)
    write(88,'(10e15.6)') dlfort, dnative, ddiff, ddiff/dnative, slfort, snative, sdiff, sdiff/snative
enddo

dabs_sinh_diff = dabs_sinh_diff / notest
sabs_sinh_diff = sabs_sinh_diff / notest
dabs_cosh_diff = dabs_cosh_diff / notest
sabs_cosh_diff = sabs_cosh_diff / notest
dabs_tanh_diff = dabs_tanh_diff / notest
sabs_tanh_diff = sabs_tanh_diff / notest

write( *, '(a,g0.2,a,g0.2,a)' ) 'Result for sinh and cosh in the range [',-range,',',range,']'
write( *, '(a)' )             'Maximum negative difference:'
write( *, '(a,g16.5)' )       '    sinh (double precision): ', dneg_sinh_diff
write( *, '(a,g16.5)' )       '    cosh (double precision): ', dneg_cosh_diff
write( *, '(a,g16.5)' )       '    tanh (double precision): ', dneg_tanh_diff
write( *, '(a,g16.5)' )       '    sinh (single precision): ', sneg_sinh_diff
write( *, '(a,g16.5)' )       '    cosh (single precision): ', sneg_cosh_diff
write( *, '(a,g16.5)' )       '    tanh (single precision): ', sneg_tanh_diff
write( *, '(a)' )             'Maximum positive difference:'
write( *, '(a,g16.5)' )       '    sinh (double precision): ', dpos_sinh_diff
write( *, '(a,g16.5)' )       '    cosh (double precision): ', dpos_cosh_diff
write( *, '(a,g16.5)' )       '    tanh (double precision): ', dpos_tanh_diff
write( *, '(a,g16.5)' )       '    sinh (single precision): ', spos_sinh_diff
write( *, '(a,g16.5)' )       '    cosh (single precision): ', spos_cosh_diff
write( *, '(a,g16.5)' )       '    tanh (single precision): ', spos_tanh_diff
write( *, '(a)' )             'Mean absolute difference:'
write( *, '(a,g16.5)' )       '    sinh (double precision): ', dabs_sinh_diff
write( *, '(a,g16.5)' )       '    cosh (double precision): ', dabs_cosh_diff
write( *, '(a,g16.5)' )       '    tanh (double precision): ', dabs_tanh_diff
write( *, '(a,g16.5)' )       '    sinh (single precision): ', sabs_sinh_diff
write( *, '(a,g16.5)' )       '    cosh (single precision): ', sabs_cosh_diff
write( *, '(a,g16.5)' )       '    tanh (single precision): ', sabs_tanh_diff

!
! Quick and dirty ...
!
write(*,*) 'Parameters ...'
xd = lf_asinh(0.5_dp) ! Print the number of terms
xd = lf_atanh(0.5_dp)

write( *, '(a)') 'asinh: double precision'
do i = -10,10
    write(*,*) lf_asinh(0.1_dp*i*abs(i)), asinh(0.1_dp*i*abs(i)), lf_asinh(0.1_dp*i*abs(i)) - asinh(0.1_dp*i*abs(i))
enddo
write( *, '(a)') 'asinh: single precision'
do i = -10,10
    write(*,*) lf_asinh(0.1_sp*i*abs(i)), asinh(0.1_sp*i*abs(i)), lf_asinh(0.1_sp*i*abs(i)) - asinh(0.1_sp*i*abs(i))
enddo

write( *, '(a)') 'acosh: out of range argument'
write(*,*) lf_acosh(0.9_dp), acosh(0.9_dp), lf_acosh(0.9_sp), acosh(0.9_sp)

write( *, '(a)') 'acosh: double precision'
do i = 0,10
    write(*,*) lf_acosh(1.0_dp+0.1_dp*i*abs(i)), acosh(1.0_dp+0.1_dp*i*abs(i)), &
        lf_acosh(1.0_dp+0.1_dp*i*abs(i)) - acosh(1.0_dp+0.1_dp*i*abs(i))
enddo
write( *, '(a)') 'acosh: single precision'
do i = 0,10
    write(*,*) lf_acosh(1.0_sp+0.1_sp*i*abs(i)), acosh(1.0_sp+0.1_sp*i*abs(i)), &
        lf_acosh(1.0_sp+0.1_sp*i*abs(i)) - acosh(1.0_sp+0.1_sp*i*abs(i))
enddo
write( *, '(a)') 'atanh: double precision'
do i = -10,10
    xd = tanh(3.1415925_dp * i / 20.0_dp)
    write(*,*) xd, lf_atanh(xd), atanh(xd), lf_atanh(xd) - atanh(xd)
enddo
write( *, '(a)') 'atanh: single precision'
do i = -10,10
    xs = tanh(3.1415925_sp * i / 20.0_sp)
    write(*,*) xs, lf_atanh(xs), atanh(xs), lf_atanh(xs) - atanh(xs)
enddo

!
! Exponential function
!
write( *, '(a)') 'exp: double precision'
do i = -20,20
    xd = i
    write(*,*) xd, lf_exp(xd), exp(xd), lf_exp(xd)-exp(xd)
enddo
write( *, '(a)') 'exp: single precision'
do i = -20,20
    xs = i
    write(*,*) xs, lf_exp(xs), exp(xs), lf_exp(xs)-exp(xs)
enddo

end program test_intrinsic_exp

module lfortran_intrinsic_trig use, intrinsic :: ieee_arithmetic implicit none private public :: sin, cos, tan, asin, acos, atan

integer, parameter :: sp = kind(1.0) integer, parameter :: dp = selected_real_kind(1+precision(1.0))

real(dp), parameter :: pi = 3.1415926535897932384626433832795_dp real(dp), parameter :: halfpi = pi / 2.0_dp

real(dp), parameter :: sqeps = sqrt( 6.0_dp * epsilon(1.0_dp) ) real(dp), parameter :: xbig = 1.0_dp / epsilon(1.0_dp)

interface sin module procedure ssin, dsin end interface

interface cos module procedure scos, dcos end interface

interface tan module procedure stan, dtan end interface

interface asin module procedure sasin, dasin end interface

interface acos module procedure sacos, dacos end interface

interface atan module procedure satan, datan end interface

contains

! sin --------------------

elemental real(dp) function abs(x) result(r) real(dp), intent(in) :: x if (x >= 0) then r = x else r = -x end if end function

elemental integer function floor(x) result(r) real(dp), intent(in) :: x if (x >= 0) then r = x else r = x-1 end if end function

elemental real(dp) function modulo(x, y) result(r) real(dp), intent(in) :: x, y r = x-floor(x/y)*y end function

elemental real(dp) function min(x, y) result(r) real(dp), intent(in) :: x, y if (x < y) then r = x else r = y end if end function

elemental real(dp) function max(x, y) result(r) real(dp), intent(in) :: x, y if (x > y) then r = x else r = y end if end function

elemental real(dp) function dsin(x) result(r) real(dp), intent(in) :: x real(dp) :: y integer :: n y = modulo(x, 2*pi) y = min(y, pi - y) y = max(y, -pi - y) y = min(y, pi - y) r = kernel_dsin(y) end function

elemental real(dp) function dcos(x) result(r) real(dp), intent(in) :: x r = sin(halfpi-x) end function

elemental real(dp) function dtan(x) result(r) real(dp), intent(in) :: x r = sin(x)/sin(halfpi-x) end function

! Accurate on [-pi/2,pi/2] to about 1e-16 elemental real(dp) function kernel_dsin(x) result(res) real(dp), intent(in) :: x real(dp), parameter :: S1 = 0.9999999999999990771_dp real(dp), parameter :: S2 = -0.16666666666664811048_dp real(dp), parameter :: S3 = 8.333333333226519387e-3_dp real(dp), parameter :: S4 = -1.9841269813888534497e-4_dp real(dp), parameter :: S5 = 2.7557315514280769795e-6_dp real(dp), parameter :: S6 = -2.5051823583393710429e-8_dp real(dp), parameter :: S7 = 1.6046585911173017112e-10_dp real(dp), parameter :: S8 = -7.3572396558796051923e-13_dp real(dp) :: z z = xx res = x (S1+z(S2+z(S3+z(S4+z(S5+z(S6+z(S7+z*S8))))))) end function

elemental real(sp) function ssin(x) result(r) real(sp), intent(in) :: x real(sp) :: y real(dp) :: tmp, x2 x2 = x tmp = dsin(x2) r = tmp end function

elemental real(sp) function scos(x) result(r) real(sp), intent(in) :: x real(sp) :: y real(dp) :: tmp, x2 x2 = x tmp = dcos(x2) r = tmp end function

elemental real(sp) function stan(x) result(r) real(sp), intent(in) :: x real(sp) :: y real(dp) :: tmp, x2 x2 = x tmp = dtan(x2) r = tmp end function

elemental real(dp) function dasin (x) real(dp), intent(in) :: x ! ! Transcribed from asin.f on netlib.org ! may 1980 edition. w. fullerton, c3, los alamos scientific lab. ! ! series for asin on the interval 0. to 5.00000d-01 ! with weighted error 1.60e-17 ! log weighted error 16.79 ! significant figures required 15.67 ! decimal places required 17.45 ! real(dp), dimension(20), parameter :: asincs = [ & .10246391753227159e0_dp, & .054946487221245833e0_dp, & .004080630392544969e0_dp, & .000407890068546044e0_dp, & .000046985367432203e0_dp, & .000005880975813970e0_dp, & .000000777323124627e0_dp, & .000000106774233400e0_dp, & .000000015092399536e0_dp, & .000000002180972408e0_dp, & .000000000320759842e0_dp, & .000000000047855369e0_dp, & .000000000007225128e0_dp, & .000000000001101833e0_dp, & .000000000000169476e0_dp, & .000000000000026261e0_dp, & .000000000000004095e0_dp, & .000000000000000642e0_dp, & .000000000000000101e0_dp, & .000000000000000016e0_dp ]

integer, parameter :: nterms = size(asincs) real(dp) :: z, y

if ( x < -1.0_dp .or. x > 1.0_dp ) then dasin = ieee_value( x, ieee_quiet_nan ) else y = abs(x)

z = 0.
if ( y > sqeps ) z = y*y
if ( z <= 0.5_dp ) dasin = x*(1.0_dp + csevl (4.0_dp*z-1.0_dp, asincs, nterms))
if ( z >  0.5_dp ) dasin = halfpi - sqrt(1.0_dp-z)*(1.0_dp + csevl (3.0_dp-4.0_dp*z, asincs, nterms))
if ( x /= 0.0_dp ) dasin = sign (dasin, x)

endif end function dasin

elemental real(sp) function sasin (x) real(sp), intent(in) :: x

real(dp) :: tmp, x2 x2 = x tmp = dasin(x2) sasin = tmp end function sasin

elemental real(dp) function dacos (x) real(dp), intent(in) :: x

dacos = halfpi - dasin(x) end function dacos

elemental real(sp) function sacos (x) real(sp), intent(in) :: x

real(dp) :: tmp, x2 x2 = x tmp = dacos(x2) sacos = tmp end function sacos

elemental real(dp) function datan (x) real(dp), intent(in) :: x

! ! Transcribed from asin.f on netlib.org ! jan 1978 edition. w. fullerton, c3, los alamos scientific lab. ! ! series for atan on the interval 0. to 4.00000d-02 ! with weighted error 1.00e-17 ! log weighted error 17.00 ! significant figures required 16.38 ! decimal places required 17.48 ! real(dp), dimension(9), parameter :: atancs = [ & .48690110349241406e0_dp, & -.006510831636717464e0_dp, & .000038345828265245e0_dp, & -.000000268722128762e0_dp, & .000000002050093098e0_dp, & -.000000000016450717e0_dp, & .000000000000136509e0_dp, & -.000000000000001160e0_dp, & .000000000000000010e0_dp ]

! tanp8(n) = tan(n*pi/8.) real(dp), dimension(3), parameter :: tanp8 = [ & .414213562373095048e+0_dp, & 1.0e0_dp, & 2.41421356237309504e+0_dp ]

! xbndn = tan((2n-1)pi/16.0) real(dp), parameter :: xbnd1 = +.198912367379658006e+0_dp real(dp), parameter :: xbnd2 = +.668178637919298919e+0_dp real(dp), parameter :: xbnd3 = +1.49660576266548901e+0_dp real(dp), parameter :: xbnd4 = +5.02733949212584810e+0_dp

! conpi8(n) + pi8(n) = n*pi/8.0 real(dp), dimension(4), parameter :: conpi8 = [ & 0.375e0_dp, & 0.75e0_dp, & 1.125e0_dp, & 1.5e0_dp ] real(dp), dimension(4), parameter :: pi8 = [ & +.176990816987241548e-1_dp, & +.353981633974483096e-1_dp, & +.530972450961724644e-1_dp, & 0.0707963267948966192e0_dp ]

integer, parameter :: nterms = size(atancs)

integer :: n real(dp) :: y, t

y = abs(x) if ( y <= xbnd1 ) then datan = x if ( y > sqeps ) datan = x(0.75_dp+csevl(50.0_dpy**2-1.0_dp, atancs, nterms)) else if ( y <= xbnd4 ) then n = 1 if ( y > xbnd2 ) n = 2 if ( y > xbnd3 ) n = 3

    t     = (y - tanp8(n)) / (1.0_dp + y*tanp8(n))
    datan = sign (conpi8(n) + (pi8(n) + t*(0.75_dp + csevl(50.0_dp*t**2-1.0_dp, atancs, nterms)) ), x)
else
    datan = conpi8(4) + pi8(4)
    if ( y < xbig ) datan = conpi8(4) + (pi8(4) - (0.75_dp + csevl (50.0_dp/y**2-1.0_dp, atancs, nterms))/y )
    datan = sign (datan, x)
endif

endif end function datan

elemental real(sp) function satan (x) real(sp), intent(in) :: x

real(dp) :: tmp, x2 x2 = x tmp = datan(x2) satan = tmp end function satan

pure real(dp) function csevl (x, cs, n) ! ! april 1977 version. w. fullerton, c3, los alamos scientific lab. ! ! evaluate the n-term chebyshev series cs at x. adapted from ! r. broucke, algorithm 446, c.a.c.m., 16, 254 (1973). also see fox ! and parker, chebyshev polys in numerical analysis, oxford press, p.56. ! ! x value at which the series is to be evaluated. ! cs array of n terms of a chebyshev series. in eval- ! uating cs, only half the first coef is summed. ! n number of terms in array cs. ! ! Note: in this version sanity checking has bee nremoved, as the ! function is only used within this module, hence we know ! the arguments are within range ! real(dp), intent(in) :: x real(dp), dimension(:), intent(in) :: cs integer, intent(in) :: n

integer :: i real(dp) :: b0, b1, b2, twox

b1 = 0.0_dp b0 = 0.0_dp twox = 2.0_dpx do i = n,1,-1 b2 = b1 b1 = b0 b0 = twoxb1 - b2 + cs(i) enddo

csevl = 0.5_dp * (b0-b2) end function csevl

end module

certik commented 3 years ago

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

Yes, that is fine with me - it is closer to LFortran than if I owned the repository :).

Op wo 15 sep. 2021 om 06:52 schreef Ondřej Čertík @.***

:

Awesome, thanks @arjenmarkus https://github.com/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 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 ?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/fortran-lang/stdlib/issues/150#issuecomment-919697476, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAN6YRY2XK6GSODDZYT3P6LUCAQ7VANCNFSM4KZEILLA . 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&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

Beliavsky commented 2 years ago

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

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

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

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

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

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 commented 1 year ago

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.