sincos.pas

embedded magazine source code for the year 1992

{ Listing 1 -- The Naive Approach }

Function Sin(x: real): real;

Const Eps = 1.0e-8;
Var n: integer;
    s, Last, Sign: real;
Begin
   n := 1;
   Last := x;
   s := Last;
   Sign := 1.0;
   Repeat
      Sign := - Sign;
      n := n + 2;
      Last := (Last * x * x)/(n * (n-1));
      s := s + Sign * Last;
   Until Last < Eps;
   Sin := s;
END;


{ Listing 2 -- Horner's Method }

{ Find the Sine of an Angle <= 45 }
Function _Sin(x: real): real;
Const s1 = 1.0/2.0/3.0;
      s2 = 1.0/4.0/5.0;
      s3 = 1.0/6.0/7.0;
      s4 = 1.0/8.0/9.0;
Var z: real;
Begin
   z := x*x;
   _Sin := ((((s4*z-1.0)*s3*z+1.0)*s2*z-1.0)*s1*z+1.0)*x;
End;

{ Find the Cosine of an Angle <= 45 }
Function _Cos(x: real): real;
Const c1 = 1.0/1.0/2.0;
      c2 = 1.0/3.0/4.0;
      c3 = 1.0/5.0/6.0;
      c4 = 1.0/7.0/8.0;
Var z: real;
Begin
   z := x*x;
   _Cos := (((c4*z-1.0)*c3*z+1.0)*c2*z-1.0)*c1*z+1.0;
End;


{ Listing 3 -- A Naive Range Reduction }

Function sin(x) is
   x = x mod 360;
   if x > 180 then
      return - sin_1(x - 180)
   else
      return sin_1(x)
   endif;
end;

Function sin_1(x) is
   if x > 90 then
      return cos_2(x - 90)
   else
      return sin_2(x)
   endif;
end;

Function sin_2(x) is
   if x > 45 then
      return _cos(90 - x)
   else
      return _sin(x)
   endif;
end;

Function cos_2(x) is
   if x > 45 then
      return _sin(90 - x)
   else
      return _cos(x)
   endif;
end;


{ Listing 4 -- Sign Function ... The Right Way }

Function Sin(x: real): real;
Var n: integer;
Begin
   n := Round(x/Pi_Over_Two);
   x := x - n * Pi_Over_Two;
   n := n Mod 4;
   if n < 0 Then
      Inc(n, 4);
   Case n Of
      0: Sin :=   _Sin(x);
      1: Sin :=   _Cos(x);
      2: Sin := - _Sin(x);
      3: Sin := - _Cos(x);
   End;
End;


Function Cos(x: real): real;
Begin
   Cos := Sin(x + Pi_Over_Two);
End;


{ Listing 5 -- Fast Sine/Cosine Using +/- 15 Degree Range }

Procedure SinCos(x: Real; Var s, c: Real);
Var n: Integer;
   s1, c1: Real;
Begin
   n := Round(x/Pi_Over_Six);
   x := x - n * Pi_Over_Six;
   n := n Mod 12;
   if n < 0 Then
      Inc(n, 12);
   s1 := _Sin(x);
   c1 := _Cos(x);
   Case n Of
      0: Begin
            s :=  s1;
            c :=  c1;
         End;
      1: Begin
            s :=   Cos_30 * s1 + Sin_30 * c1;
            c :=  -Sin_30 * s1 + Cos_30 * c1;
         End;
      2: Begin
            s :=   Sin_30 * s1 + Cos_30 * c1;
            c :=  -Cos_30 * s1 + Sin_30 * c1;
         End;
      3: Begin
            s :=   c1;
            c :=  -s1;
         End;
      4: Begin
            s :=  -Sin_30 * s1 + Cos_30 * c1;
                   c :=  -Cos_30 * s1 - Sin_30 * c1;
         End;
      5: Begin
            s :=  -Cos_30 * s1 + Sin_30 * c1;
            c :=  -Sin_30 * s1 - Cos_30 * c1;
         End;
      6: Begin
            s :=  -s1;
            c :=  -c1;
         End;
      7: Begin
            s :=  -Cos_30 * s1 - Sin_30 * c1;
            c :=   Sin_30 * s1 - Cos_30 * c1;
         End;
      8: Begin
            s :=  -Sin_30 * s1 - Cos_30 * c1;
            c :=   Cos_30 * s1 - Sin_30 * c1;
         End;
      9: Begin
            s :=  -c1;
            c :=   s1;
         End;
      10:Begin
            s :=  Sin_30 * s1 - Cos_30 * c1;
            c :=  Cos_30 * s1 + Sin_30 * c1;
         End;
      11:Begin
            s :=  Cos_30 * s1 - Sin_30 * c1;
            c :=  Sin_30 * s1 + Cos_30 * c1;
         End;
   End;
End;


{ Listing 6 -- Integer Sine/Cosine }

function sin(y) is
begin
   y = y << 1;
   z = (y * y) << 2;
   return (-y*(z*(z*(z*s4-s3)+s2)-s1) << 1;
end;

function cos(y) is
begin
   y = y << 1;
   z = (y * y) << 2;
   return c0 - z*(z*(x*c3-c2)+c1);
end;