Using the GNU C Standard Library Function sincos() in Fortran

David S. Smith, © 2008

Introduction

Most floating point units on modern CPUs have a built-in instruction called sincos that returns both the sine and the cosine of the argument on the stack. In scientific computations, one rarely uses sin() or cos() without the also calling the other, so the utility of such a function is obvious. Unfortunately, a large fraction of scientific programs are still written in Fortran, and most Fortran implementations don't have sincos() in their libraries. One workaround is to call the C standard library on systems with GNU C installed.

Procedure

Create a C source file named "sincos.c" with the following content:

#define _GNU_SOURCE
#include <math.h>

void sincos_ (double *x, double *sin, double *cos)
{
   sincos(*x, sin, cos);
}

Compile the source file with the command

% gcc -O3 -c sincos.c  

The "-c" switch tells the compiler to compile only and not try to link this file with the main program yet. The sincos() function is a GNU extension to the C standard library, so you must have the GNU C compiler installed (or another implementation of C with this function included).

In your Fortran source, freely call the sincos() subroutine. (Note that it is a subroutine, not a function, like sin() and cos().) Pass the argument of the trig functions as the first argument to sincos(), and pass two 8-byte reals (real*8 or double precision) as the second and third arguments to hold the returned sin and cos values, respectively. For example, here is a program to print the sin and cos of pi/2:

      program sincos_ex1
      implicit none
      real*8 x, PI, s, c
      parameter (PI = 3.14159265358979323846)
      x = 0.5 * PI
      call sincos(x, s, c)
      print *, x, s, c
      end 

The final step is to add in the sincos.o object file when you compile the Fortran source so that the compiled machine instructions in sincos.o can be found by the Fortran compiler. For example,

% g77 sincos.o sincos_ex1.f -o sincos_ex1
% ./sincos_ex1
  1.57079637  1. -4.371139E-08
%

Benchmarks

To benchmark the speed difference, I ran two sample programs that crudely integrated sin and cos from 0 to 100 pi on my AMD Opteron 270-based machine:

Using sincos()

      program sincos_test
      implicit none
      real*8 x, s, c, PI, dx, is, ic
      parameter (PI = 3.14159265358979323846)
      integer i, j, N
      parameter (N = 1000000)
      dx = PI / dble(N);
      x = 0.5 * dx;
      is = 0.0
      ic = 0.0
      do i = 1, 100
        do j = 1, N - 1
           call sincos(x, s, c)
           is = is + s * dx
           ic = ic + c * dx
           x = x + dx
        end do
      end do
      print *,  x, is, ic
      end 

% gfortran -O3 sincos_test.f sincos.o -o sincos_test 
% time ./sincos_test
   314.158961642532       4.747585443014451E-008 -3.052872859691430E-004

real    0m7.889s
user    0m7.820s
sys     0m0.004s
%

The source using sin() and cos()

      program sincos_test2
      implicit none
      real*8 x, s, c, PI, dx, is, ic
      parameter (PI = 3.14159265358979323846)
      integer i, j, N
      parameter (N = 1000000)
      dx = PI / dble(N);
      x = 0.5 * dx;
      is = 0.0
      ic = 0.0
      do i = 1, 100
        do j = 1, N - 1
           s = sin(x)
           c = cos(x)
           is = is + s * dx
           ic = ic + c * dx
           x = x + dx
        end do
      end do
      print *,  x, is, ic
      end 

% gfortran -O3 sincos_test2.f sincos.o -o sincos_test2
% time ./sincos_test2
   314.158961642532       4.747585443013708E-008 -3.052872859689210E-004

real    0m15.413s
user    0m15.381s
sys     0m0.008s
%

Using sincos() is roughly twice as fast as calling sin() and cos() separately. For codes that spend a majority of time calculating sines and cosines, this optimization is critical.