How to compute sine wave with accuracy over the time

Question

Use case is to generate a sine wave for digital synthesis, so, we need to compute all values of sin(d t) where:

t is an integer number, representing the sample number. This is variable. Range is from 0 to 158,760,000 for one hour sound of CD quality.

d is double, representing the delta of the angle. This is constant. And the range is: greater than 0 , less than pi.

Goal is to achieve high accuracy with traditional int and double data types. Performance is not important.

Naive implementation is:

double next()
{
    t++;
    return sin( ((double) t) * (d) );
}

But, the problem is when t increases, accuracy gets reduced because big numbers provided to "sin" function.

An improved version is the following:

double next()
{
    d_sum += d;
    if (d_sum >= (M_PI*2)) d_sum -= (M_PI*2);

    return sin(d_sum);
}

Here, I make sure to provide numbers in range from 0 to 2*pi to the "sin" function.

But, now, the problem is when d is small, there are many small additions which decreases the accuracy every time.

The question here is how to improve the accuracy.

---

Appendix 1

"accuracy gets reduced because big numbers provided to "sin" function":

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

#define TEST      (300000006.7846112)
#define TEST_MOD  (0.0463259891528704262050786960234519968548937998410258872449766)
#define SIN_TEST  (0.0463094209176730795999323058165987662490610492247070175523420)

int main()
{
    double a = sin(TEST);
    double b = sin(TEST_MOD);

    printf("a=%0.20f \n" , a);
    printf("diff=%0.20f \n" , a - SIN_TEST);
    printf("b=%0.20f \n" , b);
    printf("diff=%0.20f \n" , b - SIN_TEST);
    return 0;
}

Output:

a=0.04630944601888796475
diff=0.00000002510121488442
b=0.04630942091767308033
diff=0.00000000000000000000

Answer 1

You can try an approach that is used is some implementations of fast Fourier transformation. Values of trigonometric function are calculated based on previous values and delta.

Sin(A + d) = Sin(A) * Cos(d) + Cos(A) * Sin(d)

Here we have to store and update cosine value too and store constant (for given delta) factors Cos(d) and Sin(d).

Now about precision: cosine(d) for small d is very close to 1, so there is risk of precision loss (there are only few significant digits in numbers like 0.99999987). To overcome this issue, we can store constant factors as

dc = Cos(d) - 1 =  - 2 * Sin(d/2)^2
ds = Sin(d) 

using another formulas to update current value
(here sa = Sin(A) for current value, ca = Cos(A) for current value)

ts = sa //remember last values
tc = ca
sa = sa * dc + ca * ds
ca = ca * dc - ts * ds
sa = sa + ts
ca = ca + tc

P.S. Some FFT implementations periodically (every K steps) renew sa and ca values through trig. functions to avoid error accumulation.

Example result. Calculations in doubles.

d=0.000125
800000000 iterations
finish angle 100000 radians

                             cos               sin
described method       -0.99936080743598  0.03574879796994 
Cos,Sin(100000)         -0.99936080743821  0.03574879797202
windows Calc           -0.9993608074382124518911354141448 
                            0.03574879797201650931647050069581           

Answer 2

sin(x) = sin(x + 2N∙π), so the problem can be boiled down to accurately finding a small number which is equal to a large number x modulo 2π.

For example, –1.61059759 ≅ 256 mod 2π, and you can calculate sin(-1.61059759) with more precision than sin(256)

So let's choose some integer number to work with, 256. First find small numbers which are equal to powers of 256, modulo 2π:

// to be calculated once for a given frequency
// approximate hard-coded numbers for d = 1 below:
double modB = -1.61059759;  // = 256  mod (2π / d)
double modC =  2.37724612;  // = 256² mod (2π / d)
double modD = -0.89396887;  // = 256³ mod (2π / d)

and then split your index as a number in base 256:

// split into a base 256 representation
int a = i         & 0xff;
int b = (i >> 8)  & 0xff;
int c = (i >> 16) & 0xff;
int d = (i >> 24) & 0xff;

You can now find a much smaller number x which is equal to i modulo 2π/d

// use our smaller constants instead of the powers of 256
double x = a + modB * b + modC * c + modD * d;
double the_answer = sin(d * x);

For different values of d you'll have to calculate different values modB, modC and modD, which are equal to those powers of 256, but modulo (2π / d). You could use a high precision library for these couple of calculations.