ZX Sine

I spent the best part of my youth dwelling in front of my dear ZX Spectrum 48k. Imagine my surprise when I found out the best solution to a programming hack I tried to implement in January 2015 was to be found in the old ZX Spectrum ROM disassembly!

Timex Sinclair was a modified Spectrum, popular overseas.

The goal is to calculate sin(x) without using the libc math routine. This is commonly done using Taylor series. But as it turns out the ZX Spectrum method is both faster and more accurate!

As explained in Niels Moseley’s blogpost, the Spectrum uses Chebyshev polynomials to calculate sines. Without further ado, let me introduce some C macros!

#define T0(x) ( 1 )
#define T1(x) ( x )
#define T2(x) ( 2 * x*x - 1 )
#define T3(x) ( 4 * x*x*x - 3 * x )
#define T4(x) ( 8 * x*x*x*x - 8 * x*x + 1 )
#define T5(x) ( 16 * x*x*x*x*x - 20 * x*x*x + 5 * x )

#define C0 1.276278962f
#define C1 -.285261569f
#define C2 0.009118016f
#define C3 -.000136587f
#define C4 0.000001185f
#define C5 -.000000007f

#define P(z) ( C0 * T0(z) + C1 * T1(z) + C2 * T2(z) + C3 * T3(z) + C4 * T4(z) + C5 * T5(z) )

Here T0 – T5 are Chebyshev polynomials (first kind) and C0 – C5 are weights for the polynomial terms. In the Spectrum ROM C1 to C5 are given as half the values above and are multiplied by 2 inside P, but I have taken the liberty to multiply 2 directly into the C-coefficients here, like in the list of Chebyshev coefficients PDF down at the bottom.

Use the P macro like this:

float zx_sin(float x)
{
   double w = 4 * x;
   double z = 2 * w * w - 1;

   return P(z) * w;
}

In zx_sin() the argument x should be between -0.25 and 0.25 (percent of a full circle). In other words 0.25 of a full circle represents 90 degrees (π/2 radians).

Example

Calculate sin(π/8) radians.

/* 22.5 degrees */
zx_sin(22.5 / 360.0);

/* same as pi/8 radians */
zx_sin((M_PI / 8.0) / (2 * M_PI));

Let’s see how accurate the routine is:

Actual: 0.382683432
zx_sin: 0.382683426
Error:  0.000000006

Note. The error is further reduced if you use double precision floats as in/out arguments to zx_sin().

Symmetry

Because of the symmetry of the sine curve you can translate any sine value into the interval -0.25 to 0.25 (-π/2 to π/2 radians). The values in the range from π/2 to 3*π/2 radians form a horizontal mirror of the values from -π/2 to π/2 radians. See the image below:

Sine/Cosine symmetries (radians on X-axis)

This makes it possible to calculate all the sine values for a whole period. Note. The image also shows that cos(x) can be calculated as sin(x + π/2) radians.

Believe it or not, but even though it uses only single precision floats (doubles internally) – if x is within -0.25 to 0.25 – the zx_sin() routine above is correct to the eighth decimal!

According to the Chebyshev coefficients PDF (link at bottom) it is possible to directly get other trigonometrical functions, like cos() without shifting the phase, by using different C0-C5 coefficients.

See also

• Spectrum Basic implementation
• Tables of Chebyshev coefficients (page 145)