[J3] C TS 18661
Steven G. Kargl
kargl at troutmask.apl.washington.edu
Wed Jan 30 14:50:48 EST 2019
On Wed, Jan 30, 2019 at 06:03:25PM +0000, Bill Long wrote:
>
> > On Jan 30, 2019, at 11:29 AM, Steven G. Kargl via J3 <j3 at mailman.j3-fortran.org> wrote:
> >
> > In time harmonic analysis, one often computes exp(cmplx(0.,tau)).
> > Many people, including me, do cmplx(cos(tau),sin(tau)). Thus,
> > cos() and sin() are computed in pairs, and the argument reduction
> > for large arguments (tau > pi/4) is performed twice. With
> > sincos() argument reduction is done once and the number of function
> > calls is cut in half. As I happen to be the person who wrote the
> > sincos() implementation that appears in FreeBSD's libm, I have timing
> > codes
>
>
> Compilers can see that sin(x) and cos(x) with the same value of x are both computed, and replace it with an internal sincos(x) call that does exactly what you describe. This problem was solved long ago.
>
> Fro this example:
>
> subroutine sub (x, z)
> real(8),intent(in) :: x
> complex(8),intent(out) :: z
>
> z = cmplx(cos(x), sin(x), 8)
>
> end subroutine sub
>
>
> the internal dump of the compiled code is (Cray compiler, default options) is
>
> 1. subroutine sub( x, z )
> 5. 0[loc( z ),0].L = coss( 0[loc( x ),0].L )
> 7. return
>
> where coss is the internal name for what you call sincos.
>
> I assume other compilers do this obvious optimization as well.
>
Yes, that one is obvious. What happens if sin(x) and cos(x)
are displaced from each other:
c = cos(x)
few lines to a few dozens of lines of code. none changing x.
s = sin(x)
Does Cray merge something like this into sincos()?
I am somewhat amazed that there is pushback on giving Fortran
programmers direct access to an intrinsic subroutine with clear
benefits. I find it especially amazing consider J3 added Bessel
function routines, which are difficult to compute correctly. It
took only 17 years for someone (that would be me) to find a really
nasty bug in Sun Microsystem's fdlibm's jn.c [1]. That bugs
appeared in every libm that I could inspect, and is of course still
present in fdlibm from Netlib.
[1] https://svnweb.freebsd.org/base?view=revision&revision=215237
--
Steve
More information about the J3
mailing list