[J3] C TS 18661

[Bill Long]

> On Jan 30, 2019, at 1:50 PM, Steven G. Kargl <kargl at troutmask.apl.washington.edu> wrote:
> 
> 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)

I did try some other examples and the coss() function was called followed by extracting the real and imaginary parts of the answer. 

However, a more general issue is that writing sin(x) and cos(x) in the original source code makes the program much more readable - reflecting the intent of the algorithm.  Having the user try to “manually optimize” the code with opaque calls detracts from the code’s maintenance.  Better to write the code with clarity and let the compiler manage optimization. 

Cheers,
Bill

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

Bill Long                                                                       longb at cray.com
Principal Engineer, Fortran Technical Support &   voice:  651-605-9024
Bioinformatics Software Development                      fax:  651-605-9143
Cray Inc./ 2131 Lindau Lane/  Suite 1000/  Bloomington, MN  55425