mersenneforum.org Sin/cos algorithm question #2
 Register FAQ Search Today's Posts Mark Forums Read

 2020-01-11, 20:45 #1 Prime95 P90 years forever!     Aug 2002 Yeehaw, FL 740210 Posts Sin/cos algorithm question #2 Given cos(a) and sin(a), we want to compute cos(a*n) and sin(a*n) for n=2...10 My first implementation was a simple loop doing a complex multiply: Code:  base = step; for (i32 i = 1; i < MIDDLE; ++i) { u[i] = mul(u[i], base); // Apply twiddle base = mul(base, step); // Next twiddle } This site describes Chebyshev method Each cos / sin value can be computed with just one FMA instruction: Code:  steps[1] = step; u[1] = mul(u[1], steps[1]); // apply twiddle step1costimes2 = steps[1].cos * 2.0; steps[2].cos = step1xtimes2 * steps[1].cos - 1.0; steps[2].sin = step1xtimes2 * steps[1].sin; u[2] = mul(u[2], steps[2]); // apply twiddle for (i32 i = 3; i < MIDDLE; i++) { steps[i].cos = fma(step1costimes2, steps[i-1].cos, steps[i-2].cos); steps[i].sin = fma(step1costimes2, steps[i-1].sin, steps[i-2].sin); u[i] = mul(u[i], steps[i]); // apply twiddle } The multi-part question is: does Chebyshev's method have significantly worse accuracy issues? Testing indicates it does. can the accuracy issue be quantified or corrected?
 2020-01-11, 21:20 #2 Nick     Dec 2012 The Netherlands 2×3×52×11 Posts Unless I am missing the point here, the site you link to appears simply to be applying the standard formulae cos(A+B)=cos(A)cos(B)-sin(A)sin(B) and sin(A+B)=sin(A)cos(B)+cos(A)sin(B) and then iterating. So if the error for each trigonometric function is at most some epsilon then we can easily see how big the error gets on iteration. If you only need n up to 10, I would simply compute rational functions approximating sin(nx) and cos(nx) in the same way as for sin() and cos().
2020-01-11, 22:25   #3
Dr Sardonicus

Feb 2017
Nowhere

105348 Posts

Quote:
 Originally Posted by Prime95 Given cos(a) and sin(a), we want to compute cos(a*n) and sin(a*n) for n=2...10 My first implementation was a simple loop doing a complex multiply: Code:  base = step; for (i32 i = 1; i < MIDDLE; ++i) { u[i] = mul(u[i], base); // Apply twiddle base = mul(base, step); // Next twiddle } This site describes Chebyshev method Each cos / sin value can be computed with just one FMA instruction: The multi-part question is: does Chebyshev's method have significantly worse accuracy issues? Testing indicates it does. can the accuracy issue be quantified or corrected?
Ahh, the modified Chebyshev polynomials! As featured in this thread.

The use of n iterations would be bad for accuracy if n is large. However, the number of iterations can be reduced if n is factored. Since

Tn(2*cos(x)) = 2*cos(n*x), if n = a*b we have

2*cos(n*x) = Ta(Tb(2*cos(x)))

The composition could be extended to any number of factors. This might significantly reduce the number of iterations if n is the product of many small factors.

However, the cost of making this determination would probably make it useless in the present application.

I think using a polynomial or rational function approximation over a limited range is the way to go.

2020-01-12, 13:51   #4
R. Gerbicz

"Robert Gerbicz"
Oct 2005
Hungary

2×727 Posts

Quote:
 Originally Posted by Prime95 Given cos(a) and sin(a), we want to compute cos(a*n) and sin(a*n) for n=2...10 My first implementation was a simple loop doing a complex multiply:
Why not use formula for cos(x+y) and sin(x+y) to reduce the errors?
Get sin(10*n) from the previously computed sin(5*n),cos(5*n)
sin(9*n) from sin(4*n),cos(4*n) and sin(5*n),cos(5*n)
etc.

 2020-01-12, 16:16 #5 Dr Sardonicus     Feb 2017 Nowhere 444410 Posts Using the addition formulas might be numerically bad if k*a is at or near an integer multiple of pi/2, since either sin(k*a) or cos(k*a) will be very near 0, meaning there is nearly complete cancelation of the terms being added or subtracted in the addition formula, with a possible loss of sig figs.
2020-01-12, 16:48   #6
Prime95
P90 years forever!

Aug 2002
Yeehaw, FL

2×3,701 Posts

Quote:
 Originally Posted by R. Gerbicz Why not use formula for cos(x+y) and sin(x+y) to reduce the errors? Get sin(10*n) from the previously computed sin(5*n),cos(5*n) sin(9*n) from sin(4*n),cos(4*n) and sin(5*n),cos(5*n) etc.
That is one of my four implementations. The Chebyshev method generates the fewest F64 ops. Alas, it seems to have the worst accuracy.

Your suggested algorithm generates the most temporaries, which can be an issue on a GPU.

On AMD, with an GPU optimizer that sometimes serves as a random execution time generator, the straightforward compute of n+1 with a complex mul from cos(a*n) and sin(a*n) is fastest.

On nVidia, the fastest (other than the broken Chebyshev) computes the 3n and 5n values using 4*n and n as inputs (a complex mul and a complex mul by conjugate).

Last fiddled with by Prime95 on 2020-01-12 at 16:49

 2020-07-15, 06:09 #7 njuffa   Jul 2020 410 Posts I looked into a very similar question a couple of years ago and wrote up my conclusions in a question / answer pair on Stackoverflow: https://stackoverflow.com/questions/...lly-spaced-ang I compared four different rotation formulas to account for trade-offs between speed and accuracy.

 Similar Threads Thread Thread Starter Forum Replies Last Post neomacdev Computer Science & Computational Number Theory 1 2019-07-30 01:40 ET_ Software 20 2011-11-18 11:19 davieddy Miscellaneous Math 61 2011-07-06 20:13 geoff Twin Prime Search 5 2010-11-07 17:02 Angular Math 6 2002-09-26 00:13

All times are UTC. The time now is 13:41.

Sun Apr 11 13:41:34 UTC 2021 up 3 days, 8:22, 1 user, load averages: 1.72, 1.64, 1.55