brainwagon

I can never remember these formulas, so I wrote this program. I’m putting it here so I won’t lose it, and so others may benefit.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* $Id$
 *
 * Written by Mark VandeWettering.
 *
 * Any copyright would be pretty stupid.  Use this code as you see
 * fit, for any purpose, commercial or non-commercial, with or
 * without attribution.
 *
 * I wrote this little test program mostly as documentation for some
 * interesting recurrence relationships which I keep forgetting and
 * have to dig up from time to time.
 *
 * To generate a continuous tone, you basically are repeatedly
 * evaluating sin(n * ang), for n = 0... however many you need (or,
 * alternatively, cos).  Writing the loop in this naive way is pretty
 * horrible though, it takes many cycles to evaluate a sine or cosine.
 * Another way to do so is to compute the complex number m = cos(i *
 * ang) (where i is the sqrt(-1)) and then maintain an accumulator which
 * is repeatedly multiplied by m.  Unfortunately, eventually numerical
 * roundoff makes the resulting vector diverge.
 *
 * Luckily, there is a third option: to use a recurrence relationship
 * for cos(n * ang) and sin(n * ang).
 *
 * From section 5.4 of Numerical Recipes: The Art of Scientific Computing:
 *
 * cos(n*ang) = 2 cos(ang) cos((n-1) * ang) - cos((n-2)*ang)
 * sin(n*ang) = 2 cos(ang) sin((n-1) * ang) - sin((n-2)*ang)
 *
 * Somewhat amazingly, these recurrence relationships are _incredibly_
 * stable.  The program below implements both of them (generating a
 * sin and a cosine in each iteration) using only two multiplies and
 * two adds per iteration.  On the 3.4Ghz Intel machine I have, this
 * entire program runs in about 11 seconds, and the vectors traced
 * out by the program are still unit vectors to the precision of
 * a normal printf.  You can try replacing the "doubles", with "floats",
 * and you will see the precision isn't as good, but might be adequate
 * for your needs.  It should be noted, however, that the program runs
 * in exactly the same time in either case (doubles and floats are both
 * single cycle operations) with my machine/compiler.
 */

/* This test program runs for a billion cycles, which at the sample rate
 * specified is just about a full day.  The idea here is to see if the
 * generated sines and * cosines continue to trace out unit vectors. */

#define NCYCLES         1000000000
#define SAMPLE_RATE     11025
#define FREQ            (800.0)

main()
{
    double ang = 2.0 * M_PI * FREQ / SAMPLE_RATE ;
    double c0, c1, c2 ;
    double s0, s1, s2 ;
    double cm =  2.0 * cos(ang) ;
    int i ;

    c0 = cos(-2.0 * ang) ;
    c1 = cos(-1.0 * ang) ;

    s0 = sin(-2.0 * ang) ;
    s1 = sin(-1.0 * ang) ;

    for (i=0; i<NCYCLES; i++) {
            c2 = cm * c1 - c0 ; /* c2 is cos(i * ang) ; */
            s2 = cm * s1 - s0 ; /* s2 is sin(i * ang) ; */
            c0 = c1 ;
            c1 = c2 ;
            s0 = s1 ;
            s1 = s2 ;
    }

    for (i=0; i<100; i++) {
            c2 = cm * c1 - c0 ; /* c2 is cos(i * ang) ; */
            s2 = cm * s1 - s0 ; /* s2 is sin(i * ang) ; */
            printf("%f %f %f\n", s2, c2, s2*s2+c2*c2) ;
            c0 = c1 ;
            c1 = c2 ;
            s0 = s1 ;
            s1 = s2 ;
    }
}

/* $Log$ */

Addendum: Jim Blinn has a chapter on circle drawing in his book, JIm Blinn’s Corner which is probably relevent. It’s been some time since I read it (and I can’t find it on my shelf at the moment), but I believe it also points you to Minksy’s circle drawing algorithm that appeared in HAKMEM. It’s clever, but doesn’t really draw circles, and is no faster than the method above. What’s really odd is that it uses less storage for temporaries, which makes it look as if it were buggy.

Addendum²: Well, I coded that program up, just for fun. It turns out that you can use it as well. Gosper explains how to compute the desired value of epsilon in the Minsky’s algorithm, and unlike the code above, it is stable when evaluated using floating point arithmetic. Unfortunately, it generates a fairly notable ellipse, which doesn’t change over repeated evaluation, but isn’t very round either. As the speed the points go around the curve decrease, the curve becomes more circular.

Still, pretty interesting.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/*
 * According to Gosper, Item 151 in HAKMEM,
 * 1-epsilon^2/2 = cos(theta), where theta
 * is the desired angular increment.
 * This means that...
 * 1 - cos(theta) = epsilon^2 / 2
 * epsilon = sqrt(2.0 * (1-cos(theta))) ;
 */

#define SAMPLE_RATE 11025
#define FREQ        800.0

main()
{
    float x = 1.0, y = 0.0 ;
    float ang = 2.0 * M_PI * FREQ / SAMPLE_RATE ;
    int i ;
    float epsilon = sqrt(2.0 * (1-cos(ang))) ;

    for (i=0; i<1000000000; i++) {
        x -= epsilon * y ;
        y += epsilon * x ;
    }
    for (i=0; i<1000; i++) {
        x -= epsilon * y ;
        y += epsilon * x ;
        printf("%f %f %f\n", x, y, sqrt(x*x+y*y)) ;
    }
}

Addendum³: I’m not sure what I really said about the “complex multiply” method really is true. Using the builtin gcc complex support, it runs in about the same time (even though I think it should have twice as many multiplies, perhaps it pipelines better?) and works just as well with respect to roundoff.

For completeness:

#include <stdio.h>
#include <stdlib.h>
#include <complex.h>
#include <math.h>

#define SAMPLE_RATE 11025
#define FREQ        800.0

main()
{
        int i ;
        complex m = cexp(2.0i * M_PI * FREQ / SAMPLE_RATE) ;
        complex acc = 1 ;

        for (i=0; i<1000000000; i++)
                acc *= m ;
        for (i=0; i<1000; i++) {
                printf("%f %f %f\n", creal(acc), cimag(acc), cabs(acc)) ;
                acc *= m ;
        }
}