mersenneforum.org  

Go Back   mersenneforum.org > Extra Stuff > Programming

Reply
 
Thread Tools
Old 2020-01-11, 20:45   #1
Prime95
P90 years forever!
 
Prime95's Avatar
 
Aug 2002
Yeehaw, FL

740210 Posts
Default 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?
Prime95 is offline   Reply With Quote
Old 2020-01-11, 21:20   #2
Nick
 
Nick's Avatar
 
Dec 2012
The Netherlands

2×3×52×11 Posts
Default

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().
Nick is online now   Reply With Quote
Old 2020-01-11, 22:25   #3
Dr Sardonicus
 
Dr Sardonicus's Avatar
 
Feb 2017
Nowhere

105348 Posts
Default

Quote:
Originally Posted by Prime95 View Post
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:

<snip code>

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.
Dr Sardonicus is online now   Reply With Quote
Old 2020-01-12, 13:51   #4
R. Gerbicz
 
R. Gerbicz's Avatar
 
"Robert Gerbicz"
Oct 2005
Hungary

2×727 Posts
Default

Quote:
Originally Posted by Prime95 View Post
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.
R. Gerbicz is offline   Reply With Quote
Old 2020-01-12, 16:16   #5
Dr Sardonicus
 
Dr Sardonicus's Avatar
 
Feb 2017
Nowhere

444410 Posts
Default

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.
Dr Sardonicus is online now   Reply With Quote
Old 2020-01-12, 16:48   #6
Prime95
P90 years forever!
 
Prime95's Avatar
 
Aug 2002
Yeehaw, FL

2×3,701 Posts
Default

Quote:
Originally Posted by R. Gerbicz View Post
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
Prime95 is offline   Reply With Quote
Old 2020-07-15, 06:09   #7
njuffa
 
Jul 2020

410 Posts
Default

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.
njuffa is offline   Reply With Quote
Reply

Thread Tools


Similar Threads
Thread Thread Starter Forum Replies Last Post
Karatsuba algorithm relation to Matrix Multiplication (Strassen Algorithm)? neomacdev Computer Science & Computational Number Theory 1 2019-07-30 01:40
Doubt about P-1 algorithm ET_ Software 20 2011-11-18 11:19
TF algorithm(s) davieddy Miscellaneous Math 61 2011-07-06 20:13
Sieve algorithm geoff Twin Prime Search 5 2010-11-07 17:02
Prime95's FFT Algorithm 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

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2021, Jelsoft Enterprises Ltd.

This forum has received and complied with 0 (zero) government requests for information.

Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation.
A copy of the license is included in the FAQ.