From: mcs@goblin.caltech.edu (Martin Shepherd)
Subject: Re: C vs. FORTRAN
Organization: California Institute of Technology
Date: Fri, 4 Sep 1992 00:53:37 GMT

----
In article 57339 of comp.lang.c, dwiseman@erim.org writes:
> I have a (two) simple programs, one in fortran, and one in C. I have compiled
> them both optimized, and timed the results.  The fortran version takes .1 s
> to run, and the C version takes 12.5s to run.  I know (hope) that there is a
> simple remedy for this.  The project that I am working on has to do many 
> sines and cosines, and I'd prefer to keep the entire program in C, and not
> mix the languages.  What can I do to improve the speed of the C version?
>----------------------------------------------------------------------------
>#include <stdio.h>
>#include <math.h>
>
>main (argc,argv)
>     int argc;
>     char **argv;
>{
>  double angle,delta,the_sin,the_cos;
>  int i;
>  
>  angle = 0.0;
>  delta = 0.01;
>  
>  for (i=0; i<1000000; i++)
>  {
>    the_sin = sin(angle);
>    the_cos = cos(angle);
>    angle += delta; 
>  }

[Fortran example deleted]

In the particular program cited, you can speed up the code considerably
by using trigonometric recurrence to remove the sin() and cos() function
calls from inside the loop.

Using the trig identities:

 cos(a+b) = cos(a)*cos(b) - sin(a)*sin(b)
 sin(a+b) = sin(a)*cos(b) + cos(a)*sin(a)

you can increment sin(angle) to sin(angle+delta) and cos(angle) to
cos(angle+delta) just by knowing the sin()'s and cos()'s of 'angle' and
'delta'. Thus you only have to call the trig functions four times,
just before the loop starts. I have coded this below.

---cut here-------------------------------------------------------------
#include <math.h>
#include <stdio.h>

#define START 0.0   /* The start 'angle' */
#define DELTA 0.01  /* The increment in 'angle' per iteration */

main (argc,argv)
     int argc;
     char **argv;
{
  double the_sin; /* The sin() of the last value of 'angle' */
  double the_cos; /* The cos() of the last value of 'angle' */
  double sin_inc; /* sin(delta) used to increment sin() and cos() */
  double cos_inc; /* cos(delta) used to increment sin() and cos() */
  double las_sin; /* Temporary intermediary in recurrence calculations */
  int i;
/*
 * Determine cos() and sin() of the start angle.
 */
  the_sin = sin(START);
  the_cos = cos(START);
/*
 * Precompute the sin() and cos() of the angle increment.
 */
  sin_inc = sin(DELTA);
  cos_inc = cos(DELTA);

  for (i=0; i<1000000; i++)
  {
/*
 * Use the current value of the_sin and the_cos here.
 */

/* Then compute the incremented sin() and cos() values here using
 * recurrence relations.
 */
    las_sin = the_sin;
    the_sin = the_sin*cos_inc + the_cos*sin_inc;
    the_cos = the_cos*cos_inc - las_sin*sin_inc;
  };
  return 0;
}

---cut here----------------------------------------------------------------

Given the number of times your loop repeats, precision errors will start to
creep in. If your loop really needs to have the above repeat
count and the precision errors prove to be a problem then you could always
add an outer loop to reinitialise the_sin and the_cos, using the library
sin() and cos() functions, every given number of iterations.
However you probably won't need to do this - I just inserted an printf()
in the program to check the value of 'the_sin' against the machine's
sin() function. After 10000 iterations the two values were still identical
within the 6 decimal places that printf() produced. The program as written
above takes 1 second of cpu time on a SUN IPX (using cc -O ).

mcs@phobos.caltech.edu  (Martin Shepherd)