Alg for sinus? 
Author Message
 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 = 1E-10;

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 / (CurrSum-1));
 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  
 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  
 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 + (x-x0)*cos(x0) + (x-x0)^2*(-y0)/2 +
(x-x0)^3*(-cos(x0))/6 + .... + (x-x0)^n*(nth-Derivative of sin(x0))/n!

(Notice that cos(x0) = sqrt (1-y0^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  
 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 = 1E-10;

> 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 / (CurrSum-1));
>  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  
 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 (pre-486 book)..

Offshore news still takes up to a week to get here;  please E-mail me
a copy of non-UK non-pascal non-TV 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 181-943 6087, Nearby Fax +44 181-943 7138  



Mon, 07 Dec 1998 03:00:00 GMT  
 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.15-700000000.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  
 
 [ 6 post ] 

 Relevant Pages 
 

 
Powered by phpBB® Forum Software