Missing SinCos() in Fortran

Language enhancement
1

In the ongoing quadpack discussion, this line of code was posted:

I have a general question about this common situation where both sin() and cos() are computed with the same arguments. This occurs fairly often in signal processing, Fourier analysis, and in linear algebra in the computation of Jacobi and Givens rotations. Some languages have a function, often called something like sincos(), that computes both values with a single call. This allows for a common argument reduction step to occur, and also there might be other common intermediates computed during the evaluation of the two values to improve efficiency. Fortran does not have such a function directly, but it does have the complex exponential, which of course returns the cos() in the real part and the sin() in the imaginary part of the return value. Is the use of the complex exponential in this way a practical replacement for the “missing” sincos() subroutine in fortran?

The alternative is to hope that the fortran compiler recognizes this optimization and automatically converts the two separate trig calls into the single sincos() equivalent in its mathematical library; is that a reasonable expectation?

2

The Cray compilers did this optimization starting in the late 1980s. (I asked the compiler group about it because the pattern kept coming up in some code I was working on, and they implemented it.)

A number of compilers do this optimization today, including gfortran and flang. Have to specify -O1 or higher. Interestingly gfortran calls cexp to do the work.

4 Likes
3

Wow. What timing. I spent more time explaining a cossin() function
a month ago used in an old code to someone only familiar with f90+
syntax than I did figuring out what it did (the comments were not
in the original. Comments used to be heavy and take up a lot of
space in the card box!). It looked something like this …

      program oldstyle
      complex cv, cossin
      real  trig(2), cosprt, sinprt
     0 equivalence 
     1( cv,trig(1) ),
     2( trig(1),cosprt), 
     3( trig(2),sinprt)
C     The "cossin()" function computes both sin(X) and cos(X) values with a single
C     call -- the complex value returned contains the cos() in the real part
C     and the sin() in the imaginary part of the return value.
      cossin(EXPRESSION)=cexp(cmplx(0.0,EXPRESSION))

      cv=cossin(3.0**2/4.0)
      write(*,*)'SIN part=',sinprt
      write(*,*)'COS part=',cosprt
      end

I think you are right and modern compilers typically optimize an adjacent
sin() and cos() call with the same argument into a call to an equivalent
of cossin().

So it would be interesting to time results, particularly since exp()
solves a lot more than just the special case of an angle of zero (0) so
it is potentially slower than the explicit separate calls (ie. exp()
will likely execute conditional code that will be evalulated because
the angle can be non-zero)..

The old method above was probably efficient because it used EQUIVALENCE and
statement functions. Any modern version of that that I can think of probably
adds considerably more overhead as well.

That being said, it would be interesting to time some cases.

If you can read the code you inherited some old code at some point or have been around for a while! :wink:

4

I would say that the modern fortran code to do this is relatively simple and straightforward by comparison.

cv = exp(cmplx(0.0,3.0**2/4.0))
write(*,*) 'SIN part=', cv%im
write(*,*) 'COS part=', cv%re

No equivalence is needed, no statement function is necessary, and the real and imaginary parts of the variable are individually available. Also, the expression works for any real kind, but the programmer must be careful with the cmplx() intrinsic because of its legacy.

The 0.0 real part of the exponent determines the magnitude, the imaginary part determines the angle. e^{(x+iy)} =e^x(cos(y)+i\ sin(y)). With x=0 the magnitude is 1. The complex exponential must test for this at runtime. In contrast, there is no equivalent magnitude test required for the real cos() and sin() expression. So in principle this is a little less efficient than a SinCos() subroutine would be, but I expect not in a significant way.

BTW, I’m assuming that is how the complex exponential is actually computed. There are other ways to compute the exponential, so this may all depend on the underlying math library.

edit: There is one more possible numerical issue when using the complex exp() function. When computing sin() and cos(), the results are in the range [-1,1]. The programmer never needs to check the result for numerical errors outside that range. I don’t think the fortran standard specifies this, but I think the IEEE standard does. In any case, if a fortran compiler did return a value outside that range, its users would complain. On the other hand, I don’t know if the complex exp() function is expected to satisfy that same requirement, even when the real part of the argument is exactly 0.0, despite the formal mathematical equivalence. If it is evaluated using an internal SinCos() function, then of course it would, but in general this might be relevant to the question of whether the complex exp() function can be used to replace the missing SinCos() function in fortran.

5

Just as a side remark: in that case the original formulas were written with complex exponents of course (this is wave amplitude and corresponding transformations). I intentionally coded everything as two real arrays to understand what was happening & debug. Beside from that, since the numerical integration is used afterwords, any speed-ups from these two lines would be completely negligible. But yes, if a complex exponent (for me double precision is necessary) is calculated with the same speed as a trigonometric function, then I will keep it in mind. Thanks for suggestion!

6

There is a long history of this in fortran. Although f77 (and earlier) had single precision complex types, it did not have a double precision complex type. If double precision complex was supported by a compiler of that era, it was as an extension, not as part of the standard. So many programmers of that era would code complex algorithms using just real arithmetic because they needed double precision or because their code might possibly need to be converted in the future to double precision. Also, compilers would typically optimize the real expressions better than they would the complex expressions. This was due to, among other reasons, the stride between the real (or imaginary) elements in a complex array was two numerical units rather than one, and also the chicken-or-egg problem that if programmers didn’t use complex variables, there was little incentive for vendors to spend time optimizing complex operations. This convention was further reinforced by popular algorithms, such as complex FFT, that took separate arrays for the real and imaginary arguments. That single-precision-only legacy is also why the generic cmplx() intrinsic looks so out-of-place relative to other intrinsics in the modern language; it was originally defined to always return a single precision complex result, regardless of its argument type, because that was the only standard complex type there was. The dcmplx() function was an extension, not a standard intrinsic. So now, to be careful, it is usually recommended to always include the kind argument when using that generic intrinsic so it is clear what the programmer intended (e.g. cmplx(x,y,kind(x))). With that convention, you can write code that is portable w.r.t. the real/complex working precision kind value.

1 Like
7

A few years ago when I wrote this subroutine I found that depending of the compiler in use a call to this routine was either marginally slower that two sequential calls to sin and cos, and potentially up to 2 times faster.

…

pure elemental subroutine sincosD(x, s, c)
!--------------------------------------------------------------
! SUBROUTINE: sincos
!
! DESCRIPTION:
!   Computes both the sine and cosine of an angle simultaneously. 
!
! INTERFACE:
!   x (Input)  : The angle in radians, 64-bit real.
!   s (Output) : The calculated sine value, 64-bit real.
!   c (Output) : The calculated cosine value, 64-bit real.
!
! METHODOLOGY:
!   1. Range Reduction: Uses high-precision Cody-Waite 
!      constants (P1, P2, P3) to reduce the input angle x 
!      modulo pi/2, preventing catastrophic cancellation error
!      for large values of x.
!   2. Approximation: Computes polynomial expansions for sine
!      and cosine using Maclaurin Series.
!   3. Reconstruction: Uses a lookup table (LUT) based on the 
!      quadrant to correct the signs and swap sine/cosine 
!      functions as required.
!--------------------------------------------------------------
use iso_fortran_env, only: real64, int32, int64
implicit none
integer, parameter :: wp = real64
real(wp), intent(in) :: x
real(wp), intent(out) :: s, c
real(wp) :: x_reduced, x2, s_poly, c_poly
integer(int64) :: quadrant
integer(int32) :: q_idx
! Constants (sufficient for real128)
real(wp), parameter :: PI_4 = 0.78539816339744830961566084581987572_wp
real(wp), parameter :: INV_PI_2 = 0.63661977236758134307553505349005745_wp
! High-Precision Cody-Waite Split Constants (for real64)
real(wp), parameter :: P1 = 1.570796326734125614166259765625_wp
real(wp), parameter :: P2 = 6.0771005063039659766010730103343e-11_wp
real(wp), parameter :: P3 = 2.0222662487959506315414243419361e-21_wp
! Cosine Factorial Reciprocals
real(wp), parameter :: C12 = 1.0_wp / 479001600.0_wp
real(wp), parameter :: C10 = -1.0_wp / 3628800.0_wp
real(wp), parameter :: C8  = 1.0_wp / 40320.0_wp
real(wp), parameter :: C6  = -1.0_wp / 720.0_wp
real(wp), parameter :: C4  = 1.0_wp / 24.0_wp
real(wp), parameter :: C2  = -0.5_wp
! Sine Factorial Reciprocals
real(wp), parameter :: S13 = 1.0_wp / 6227020800.0_wp
real(wp), parameter :: S11 = -1.0_wp / 39916800.0_wp
real(wp), parameter :: S9  = 1.0_wp / 362880.0_wp
real(wp), parameter :: S7  = -1.0_wp / 5040.0_wp
real(wp), parameter :: S5  = 1.0_wp / 120.0_wp
real(wp), parameter :: S3  = -1.0_wp / 6.0_wp
! Lookup Tables
real(wp), parameter :: cos_sign_lut(0:3) = [1.0_wp,-1.0_wp,-1.0_wp, 1.0_wp]
real(wp), parameter :: sin_sign_lut(0:3) = [1.0_wp, 1.0_wp,-1.0_wp,-1.0_wp]
real(wp), parameter :: cos_weight(0:3)   = [1.0_wp, 0.0_wp, 1.0_wp, 0.0_wp]
real(wp), parameter :: sin_weight(0:3)   = [0.0_wp, 1.0_wp, 0.0_wp, 1.0_wp]
    if (x == 0.0_wp) then
        s = x        
        c = 1.0_wp
    else
        ! Range reduction
        if (abs(x) <= PI_4) then
            x_reduced = x
            q_idx = 0
        else
            quadrant = int(x * INV_PI_2 + sign(0.5_wp, x), kind=int64)
            x_reduced = ((x - real(quadrant, kind=wp) * P1) &
                            - real(quadrant, kind=wp) * P2) &
                            - real(quadrant, kind=wp) * P3
            ! Zero-based masking
            q_idx = int(iand(quadrant, 3_int64), kind=int32)
        end if
        x2 = x_reduced * x_reduced
        ! Cosine
        c_poly = C10 + x2 * C12
        c_poly = C8  + x2 * c_poly
        c_poly = C6  + x2 * c_poly
        c_poly = C4  + x2 * c_poly
        c_poly = C2  + x2 * c_poly
        c_poly = 1.0_wp + x2 * c_poly
        ! Sine
        s_poly = S11 + x2 * S13
        s_poly = S9  + x2 * s_poly
        s_poly = S7  + x2 * s_poly
        s_poly = S5  + x2 * s_poly
        s_poly = S3  + x2 * s_poly
        s_poly = x_reduced * (1.0_wp + x2 * s_poly)
        ! Reconstruction
        c = (c_poly * cos_weight(q_idx) + s_poly * sin_weight(q_idx)) * cos_sign_lut(q_idx)
        s = (s_poly * cos_weight(q_idx) + c_poly * sin_weight(q_idx)) * sin_sign_lut(q_idx)
    end if
end subroutine sincosD

…

1 Like
8

Nice. Using six different compilers and just calling the C sincos procedure versus using the call to exp with a zero angle and complex values; and calling sin() and cos() explicitly I found a huge range of performances.

The fastest by far was with the AOCC flang compiler and ifx, which vectorized the simple expression so (given enough threads/cores) running up to thirty times faster than some of the other methods. Calling an external optimized version from a math library might get some good results too.

Like in anything regarding timing using a simple case instead of the actual application there were a lot of little gotchas so I did not post that here but it left me wishing there was a concerted effort to provide opensource versions of such routines for optimization by the community. So far the focus of the stdlib project has been more to add to what is available, not to duplicate the existing intrinsics but it would be interesting to provide sin() and cos() as well.

Some of the flaws as you hint at to using EXP was it was inaccurate with large values whereas the C sincos() and sin() and cos() were much better in that regard.

This is a nice example I will add to my little timing test. The test is more useful in exposing the range of optimizations at the machine code level, what various compiler switches can do, and exploring all the troubles with simple timing tests (I almost always try to profile and time the real application instead of an external timing because there are so many variations possible between a test problem and a real application) than being a guide on calling sin() and cos() in compbination but all that being said I am curious what the results will be as I was surprised by some of the ones I got from what on the surface was a simple test program.

This looks like it should have the accuracy of the intrinsics. Nice example. If anything particularly interesting happens adding it to my timing test I will post the timings and the test program (and a “caveat emptor” about the perils of using little test timing programs as well, as they can be great learning experiences to use, but often teach a lot of unexpected issues more than the initial goal!).

Looks like a keeper though. Thanks for posting it.

9

This looks good. It does accurate argument reduction, which is shared for both the cos and sin evaluation, and it precomputes the divisions so that only multiplies and adds are used in the code. I’m unfamiliar with some of these details. First, the argument reduction is only down to pi/2 (I think); is there any advantage to reducing it down to pi/4? Or maybe that is what the cos_weight and sin_weight arrays do? An alternative approach for the argument reduction is to shift the binary exponent of the argument (with scale()), and then do recursive squaring to compute the result (think of the complex exp() equivalent); would that approach be more or less accurate? The polynomials are evaluated manually rather than in a loop; does the compiler compute those statements optimally, or would the loop code be better? There is no test to ensure the result values are <=1.0; is it alright to always avoid that test? Finally, the intrinsic trig functions that end in D take degree arguments, so in modern codes perhaps the subroutine name should be changed to avoid confusion.

10

When the Cray CFT compiler was taught to do the optimization, there were already pre-existing COSS routines (both scalar and vector) in the math library. Thus the common code for argument reduction and such was only done once. According to this old reference, CEXP called COSS internally (see page A-4):

1 Like
11

I’m curious about the argument reduction step. The argument is limited to 2^24, and the floating point fraction had 48 bits, so I’m thinking that it did not do anything special to preserve precision for large arguments. For those with more Cray experience, does that seem correct?

12

I can’t answer your question directly. But I do know the math library ($ARLIB on COS, libm on Unicos) was written in CAL (Cray Assembly Language) and pretty highly optimized for speed. A few years later (1990-ish?) the fellow who supported it went through all the intrinsics in an attempt to improve their accuracy to within (IIRC) 1 ULP or less. Some intrinsics ran a bit slower, and some he found speedups. Since the answers could differ a bit, it was released as libmv2. Folks could then use the old libm for reproducability, or the new one for best accuracy.

I used to have a copy of an internal paper he wrote on the results. No idea if I still have it somewhere…

13

The range reduction in the subroutine I posted is from -pi/4 to pi/4 i.e. spanning pi/2. This minimises the number of terms required in the Maclaurin series.

Yes there should probably be check to ensure that the magnitude of returned s and c values are <= 1, clamping where necessary.

There are various ways to perform the reconstruction. Look up tables are not strictly required and the same functionality can be achieved for example using SELECT CASE.

As there are likely to be very few cases where x is exactly equal to zero, this specific test can be removed, which may be better for vectorisation when dealing with large arrays?

I suspect a compiler would unroll any do-loops to arrive at the “manual” polynomial evaluation as the code I posted. Originally I had loops, calculating the factorial reciprocals on the fly and only exiting when convergence was achieved, which is definitely not the right way to do it.

The fast Cody-Waite range reduction starts to breaks down above |X| = 10^5 so another path is required for very large arguments.

…

1 Like
14

I think it would be a useful feature in fortran if all intrinsics had this kind of option. Sometimes the programmer wants speed, even it it costs a bit or two of accuracy, and sometimes the programmer wants accuracy, even if it costs a little more to compute. It would be nice if the programmer had that option. As the language is now, and has been for 70+ years, if the programmer wants speed then he needs to write his own code, and if the programmer wants accuracy, then he needs to write his own code, the language specification neither requires nor suggests either. IEEE arithmetic does make some recommendations in this regard, but I don’t think they are strict requirements. IEEE support in fortran (which is triggered by the USE statement) does begin to give the programmer some control. That might be a reasonable model for the more general feature.