Author 
Message 
Saruna #1 / 6

Alg for sinus?
Hi! I have been struggling with this for some time now  and I figured I'd better ask. Here is my problem: does anyone know an algorithm for calculating a sinus of a given angle x ? I know there is a standard function Sin(x), but I have to write my own implementation. The formula that I currently use is x x^3 x^5 x^7 sin(x) :=    +    .... etc until a member is < epsilon. 1! 3! 5! 7! I have had no problem implementing this, but the algorithm behaves very, very weirdly  if I calculate MySin(10) and TPSin(10, they are identical. However, as the angle increases, MySin(100) is pretty far off the TPSin(100), for instance. What is the problem? Any better algorithms??? Help!! Here is my source  try it with angle of 10, then 100... :( There are 2 implementations currently in there  OneKinfOf_Sinus is slower, MySinus can owerflow. {$N+} program Sinuses; const epsilon = 1E10; function OneKindOf_Sinus(x:integer):double; var ret, tSin: double; sign :1..1; expo, count : word; begin TSin := 0; expo := 1; sign := 1; repeat ret := 1.0; count := 0; repeat inc(count); ret := ret * (x / count); until count = expo; ret := ret * sign; Tsin := Tsin + ret; expo := expo + 2; sign := sign; until abs(ret) <= epsilon; OneKindOf_Sinus := TSin; end; function MySinus(Angle:integer):double; var tSin, A, B: double; CurrSum: word; Sign :1..1; begin CurrSum := 1; Sign := 1; A := Angle; tSin := A; repeat CurrSum := CurrSum + 2; B := (Angle / (CurrSum1)); A := A * (Angle / CurrSum); B := B * A; tSin := tSin + B * Sign; A := B; Sign := Sign * 1; until abs(B) < epsilon; {until CurrSum > 200;} MySinus := tSin; end; var x: integer; begin write('Enter an angle (I''l show ye'' a sinus if it):'); readln(x); x := x mod 360; writeln( 'My Sinus of ',x,' is ', MySinus(x)); writeln( 'Tp Sinus of ',x,' is ', Sin(x)); writeln( 'My Other_Sinus of ',x,' is ', OneKindOf_Sinus(x)); end.

Wed, 02 Dec 1998 03:00:00 GMT 


Chris Ryde #2 / 6

Alg for sinus?
Quote:
> I have had no problem implementing this, but the algorithm behaves very, very > weirdly  if I calculate MySin(10) and TPSin(10, they are identical. However, > as the angle increases, MySin(100) is pretty far off the TPSin(100), for > instance. What is the problem? Any better algorithms??? Help!!
I think you'll find it's an accuracy problem. As the size of the angle increases the series requires more terms in order to create an accurate estimation of the actual sin. Fortunately this problem is easily overcome as sin(x) has a period of 2*pi, so you need to do x := x MOD twoPi at the start of each of your functions to make x into a value in the interval [0, 2*pi). (Ok, so you can't use MOD on reals, but I'm sure you can find a way around that  it's not too hard.) Enjoy. ChrisR:

Fri, 04 Dec 1998 03:00:00 GMT 


Leo Sarasu #3 / 6

Alg for sinus?
The method you are using is called the Taylor expansion series. As you probably know, this method finds an aproximation (as precise as you want) of the function based on a known value of this function. In this case you are implicitly using 0 as the reference value, for which sin(0) = 0. If you were using a different value x0, for which you knew the result, say sin(x0)=y0, then you should use the following series: sin (x) = 1*y0 + (xx0)*cos(x0) + (xx0)^2*(y0)/2 + (xx0)^3*(cos(x0))/6 + .... + (xx0)^n*(nthDerivative of sin(x0))/n! (Notice that cos(x0) = sqrt (1y0^2)) This well know method converges faster the closer x is to x0. Therefore, for a value like 100, it's no surprise that your Sin procedure converges slowly, since you are taking x0=0. Since you want to use your procedure for any value of x, it's useless to find a given x0 which will give a fast convergence. There are other simpler tricks instead: You use the fact that sin (x+2*Pi) = sin (x). This means that for a given x you can calculate the sin (x mod 2*Pi) and the result is the same. Therefore, the maximum value you will be inputting to the series will be 2*Pi, or about 6,28. That's already an improvement. A finer trick is by using sin (x+Pi) =  sin (x). This gives a maximum input value for the series of Pi. I'm sure you know how to implement these checks at the beginning of your routine. Otherwise, ask me. Cheers, Leo

Fri, 04 Dec 1998 03:00:00 GMT 


PZAN #4 / 6

Alg for sinus?
Quote:
> Hi! > I have been struggling with this for some time now  and I figured I'd better > ask. Here is my problem: does anyone know an algorithm for calculating > a sinus of a given angle x ? I know there is a standard function Sin(x), but I > have to write my own implementation. The formula that I currently use > is x x^3 x^5 x^7 > sin(x) :=    +    .... etc until a member is < epsilon. > 1! 3! 5! 7! > I have had no problem implementing this, but the algorithm behaves very, very > weirdly  if I calculate MySin(10) and TPSin(10, they are identical. However, > as the angle increases, MySin(100) is pretty far off the TPSin(100), for > instance. What is the problem? Any better algorithms??? Help!! > Here is my source  try it with angle of 10, then 100... :( > There are 2 implementations currently in there  OneKinfOf_Sinus is slower, > MySinus can owerflow. > {$N+} > program Sinuses; > const epsilon = 1E10; > function OneKindOf_Sinus(x:integer):double; > var ret, tSin: double; > sign :1..1; > expo, count : word; > begin > TSin := 0; > expo := 1; > sign := 1; > repeat > ret := 1.0; > count := 0; > repeat > inc(count); > ret := ret * (x / count); > until count = expo; > ret := ret * sign; > Tsin := Tsin + ret; > expo := expo + 2; > sign := sign; > until abs(ret) <= epsilon; > OneKindOf_Sinus := TSin; > end; > function MySinus(Angle:integer):double; > var tSin, A, B: double; > CurrSum: word; > Sign :1..1; > begin > CurrSum := 1; > Sign := 1; > A := Angle; > tSin := A; > repeat > CurrSum := CurrSum + 2; > B := (Angle / (CurrSum1)); > A := A * (Angle / CurrSum); > B := B * A; > tSin := tSin + B * Sign; > A := B; > Sign := Sign * 1; > until abs(B) < epsilon; > {until CurrSum > 200;} > MySinus := tSin; > end; > var x: integer; > begin > write('Enter an angle (I''l show ye'' a sinus if it):'); > readln(x); > x := x mod 360; > writeln( 'My Sinus of ',x,' is ', MySinus(x)); > writeln( 'Tp Sinus of ',x,' is ', Sin(x)); > writeln( 'My Other_Sinus of ',x,' is ', OneKindOf_Sinus(x)); > end.Some one else mentioned that the greater the absolute value of x, the
longer it takes to converge and the more sensitive it is to fixed length math errors. They are correct. If my memory servers me, there is an op code for the '387, 486 and beyond processors which bring an arbitrary argument into the range of pi to pi. This is done in extended reals to minimize errors. You might also want to declare your variables to be extended to improve accuracy. If you have a '387 or 486 or better, execution time will be about the same. Good luck P Zank

Sat, 05 Dec 1998 03:00:00 GMT 


Dr John Stockton NPL #5 / 6

Alg for sinus?
Quote: >>... >... >If my memory servers me, there is an op >code for the '387, 486 and beyond processors which bring an arbitrary >argument into the range of pi to pi. This is done in extended reals to >minimize errors.
There is FPREM, and FPREM1, according to information received; I'm not clear on the definition, but FPREM tests to be a sort of Push(Pop mod Pop)  so it could be used to reduce into the range 2pi..0 / 0..+2pi. Anyone got more information? it would of course be a useful function in AC measurement work, even though it can be written in standard Pascal. BPW didn't recognise FPREM1  maybe it's exactly what one wants? It's in the TP6 asm book, as '387 only (pre486 book).. Offshore news still takes up to a week to get here; please Email me a copy of nonUK nonpascal nonTV followups ! Posting also slow ! Regret that the system puts unzoned (UK civil) time on messages. 
National Physical Laboratory, Teddington, Middlesex, TW11 0LW, UK Direct Telephone +44 181943 6087, Nearby Fax +44 181943 7138

Mon, 07 Dec 1998 03:00:00 GMT 


Voti #6 / 6

Alg for sinus?
Quote: >have to write my own implementation. The formula that I currently use >is x x^3 x^5 x^7 > sin(x) :=    +    .... etc until a member is < epsilon. > 1! 3! 5! 7! >I have had no problem implementing this, but the algorithm behaves very, very >weirdly  if I calculate MySin(10) and TPSin(10, they are identical. However, >as the angle increases, MySin(100) is pretty far off the TPSin(100), for >instance. What is the problem? Any better algorithms??? Help!!
I think the broblem is that when the argument x increases the intermediate calculations are likely to introduce relatively big rounding errors (e.g. 700000001.15700000000.78=0.37 only if your numbers have 11 or more digits of accuracy, otherwise it will be rounded to 0), which are accumulated and destroy the final result. You should try to reduce any argument x to an equivalent x0 between 0 and pi, and compute the series for this x0. Votis  Votis Kokavessis Mathematics teacher & pascal programmer Thessaloniki, Greece


Mon, 07 Dec 1998 03:00:00 GMT 


