[J3] C TS 18661

[Steven G. Kargl]

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