Xref: unixg.ubc.ca alt.fractals:2840 bit.listserv.frac-l:1253
Newsgroups: alt.fractals,bit.listserv.frac-l
Path: unixg.ubc.ca!ubc-cs!destroyer!caen!uakari.primate.wisc.edu!usenet.coe.montana.edu!news.u.washington.edu!milton.u.washington.edu!jesjones
From: jesjones@milton.u.washington.edu (Jesse Jones)
Subject: Fractal paper #2: Cordic.m
Message-ID: <1992Sep17.225337.11933@u.washington.edu>
Summary: Implementation of cordic algorithm in Modula-2.
Keywords: Fractal Cordic Integer Math Mandella
Sender: news@u.washington.edu (USENET News System)
Organization: University of Washington, Seattle
Date: Thu, 17 Sep 1992 22:53:37 GMT
Lines: 356
 
MPLEMENTATION MODULE Cordic;
 
(* This module is based on the cordic algorithm described in the 
   article "A Unified Algorithm for Elementary Functions" in 
   Computer Arithmetic edited by Earl Swartzlander Jr. *)
 
FROM Sane   IMPORT RoundDir, GetRound, SetRound, Rint, Scalb, INF;
IMPORT Functions;
 
 
CONST
   MaxBits = 32;
 
   PI     = 3.141592653589793238;   
   TwoPi  = 2.0*PI;
   HalfPi = PI/2.0;
   
   E = 2.718281828459045235;
 
TYPE
   Reply = RECORD
      x, y, z: REAL;
   END;
               
VAR
   Alpha: ARRAY [-1..1], [0..MaxBits] OF REAL;
   Scale: ARRAY [-1..1] OF REAL;
 
 
PROCEDURE InitLookUp;      
   VAR
      pow2: REAL;
      m, k: INTEGER;
      next: INTEGER;
BEGIN
   pow2 := 1.0;
   FOR k := 0 TO MaxBits DO
      Alpha[1, k] := Functions.ArcTan(pow2);
      Alpha[0, k] := pow2;
      IF k = 0 THEN
         Alpha[-1, k] := 0.0;               (* Can't take ArcTanh of 1.0 *)
      ELSE
         Alpha[-1, k] := Functions.ArcTanh(pow2);
      END;
      pow2 := pow2/2.0;
   END;
   FOR m := -1 TO 1 DO
      k := ORD(m <> 1); next := 4;          
      Scale[m] := 1.0; 
      REPEAT
         Scale[m] := Scale[m]*Functions.Sqrt(1.0 + Scalb(-2*k, FLOAT(m)));
         IF (m = -1) AND (k = next) THEN
            next := 3*k + 1;
         ELSE
            INC(k);
         END;
      UNTIL k > NoBits;
      Scale[m] := 1.0/Scale[m];
   END;
END InitLookUp;
 
PROCEDURE DoCordic (x, y, z: REAL; m: INTEGER; reduceAngle: BOOLEAN): Reply;
   VAR 
      xn, yn, zn: REAL;
      start, k  : INTEGER;
      nextK     : INTEGER;
      subtract  : BOOLEAN;
      reply     : Reply;
BEGIN
   start := ORD(m <> 1);         (* start = 0 if m=1, else start = 1 *)
   k := start; nextK := 4;
   REPEAT
      IF reduceAngle THEN
         subtract := y > 0.0;
      ELSE
         subtract := z < 0.0;
      END;
      IF subtract THEN
         xn := x + Scalb(-k, FLOAT(m)*y);
         yn := y - Scalb(-k, x);
         zn := z + Alpha[m, k]
      ELSE
         xn := x - Scalb(-k, FLOAT(m)*y);
         yn := y + Scalb(-k, x);
         zn := z - Alpha[m, k]
      END;
      x := xn; y := yn; z := zn;
      IF (m = -1) AND (k = nextK) THEN
         nextK := 3*k + 1;
      ELSE
         INC(k);
      END;
   UNTIL k > NoBits;
   reply.x := xn; reply.y := yn; reply.z := zn;
   RETURN reply;
END DoCordic;
 
   (* x = multiple*interval + delta *)
PROCEDURE Reduce (x, interval: REAL; VAR multiple: INTEGER; VAR delta: REAL);
BEGIN
   x := ABS(x);
   IF x >= interval THEN
      multiple := TRUNC(x/interval);
      delta := x - FLOAT(multiple)*interval;
   ELSE
      multiple := 0;
      delta := x;
   END;
END Reduce;
 
(* ----------------------------------- Trig Functions -------------------------------------- *)
 
PROCEDURE Sin (x: REAL): REAL;               
   VAR 
      answer  : Reply;
      quadrant: INTEGER;
      angle   : REAL;
      result  : REAL;
BEGIN
   Reduce(x, HalfPi, quadrant, angle);
   answer := DoCordic(Scale[1], 0.0, angle, 1, FALSE);      
   CASE quadrant MOD 4 OF
      0: result := answer.y;   |      (* result =  Sin(angle) *)
      1: result := answer.x;   |      (* result =  Cos(angle) *)
      2: result := -answer.y;   |      (* result = -Sin(angle) *)
      3: result := -answer.x;   |      (* result = -Cos(angle) *)
   END;
   IF x < 0.0 THEN result := -result END;
   RETURN result; 
END Sin;
 
PROCEDURE Cos (x: REAL): REAL;
   VAR 
      answer  : Reply;
      quadrant: INTEGER;
      angle   : REAL;
      result  : REAL;
BEGIN
   Reduce(x, HalfPi, quadrant, angle);
   answer := DoCordic(Scale[1], 0.0, angle, 1, FALSE);      
   CASE quadrant MOD 4 OF
      0: result := answer.x;   |      (* result =  Cos(angle) *)
      1: result := -answer.y;   |      (* result = -Sin(angle) *)
      2: result := -answer.x;   |      (* result = -Cos(angle) *)
      3: result := answer.y;   |      (* result =  Sin(angle) *)
   END;
   RETURN result; 
END Cos;
 
PROCEDURE Tan (x: REAL): REAL;
   VAR 
      answer  : Reply;
      quadrant: INTEGER;
      angle   : REAL;
      sin, cos: REAL;
BEGIN
   Reduce(x, HalfPi, quadrant, angle);
   answer := DoCordic(Scale[1], 0.0, angle, 1, FALSE);      
   CASE quadrant MOD 4 OF
      0: sin := answer.y;
         cos := answer.x;   |      
      1: sin := answer.x;
         cos := -answer.y;   |      
      2: sin := -answer.y;
         cos := -answer.x;   |      
      3: sin := -answer.x;
         cos := answer.y;   |   
   END;
   IF x < 0.0 THEN sin := -sin END;
   RETURN sin/cos; 
END Tan;
 
PROCEDURE ArcTan (x: REAL): REAL;
   VAR result: Reply;
BEGIN
   result := DoCordic(1.0, x, 0.0, 1, TRUE);         
   RETURN result.z;         
END ArcTan;
 
(* -------------------------------- Hyperbolic Functions ----------------------------------- *)
 
PROCEDURE Sinh (x: REAL): REAL;               
   VAR 
      answer  : Reply;
      quadrant: INTEGER;
      angle   : REAL;
      result  : REAL;
BEGIN
   Reduce(x, Functions.Ln(2.0), quadrant, angle);
     answer := DoCordic(Scale[-1], 0.0, angle, -1, FALSE);   
   result := Scalb(quadrant-1, answer.x+answer.y - Scalb(-2*quadrant, answer.x-answer.y));
   IF x < 0.0 THEN result := - result END;
   RETURN result;
END Sinh;
 
PROCEDURE Cosh (x: REAL): REAL;
   VAR 
      answer  : Reply;
      quadrant: INTEGER;
      angle   : REAL;
      result  : REAL;
BEGIN
   Reduce(x, Functions.Ln(2.0), quadrant, angle);
     answer := DoCordic(Scale[-1], 0.0, angle, -1, FALSE);      
   result := Scalb(quadrant-1, answer.x+answer.y + Scalb(-2*quadrant, answer.x-answer.y));
   RETURN result;
END Cosh;
 
PROCEDURE Tanh (x: REAL): REAL;
   VAR 
      answer     : Reply;
      quadrant   : INTEGER;
      angle, temp: REAL;
      sinh, cosh : REAL;
BEGIN
   Reduce(x, Functions.Ln(2.0), quadrant, angle);
     answer := DoCordic(Scale[-1], 0.0, angle, -1, FALSE);      
   temp := Scalb(-2*quadrant, answer.x-answer.y);
   sinh := Scalb(quadrant-1, answer.x+answer.y - temp);
   IF x < 0.0 THEN sinh := -sinh END;
   cosh := Scalb(quadrant-1, answer.x+answer.y + temp);
   RETURN sinh/cosh;
END Tanh;
 
PROCEDURE ArcTanh (x: REAL): REAL;
   VAR 
      M, T  : REAL;
      E     : INTEGER;
      answer: Reply;
      result: REAL;
BEGIN
   IF x <= -1.0 THEN
      result := -INF();
   ELSIF x >= 1.0 THEN
      result := INF();
   ELSE
      E := 0; M := 1.0 - ABS(x);
      WHILE (M < 0.5) OR (M >= 1.0) DO
         IF M >= 1.0 THEN
            M := M/2.0;
            INC(E);
         ELSIF M < 0.5 THEN
            M := M*2.0;
            DEC(E);
         END;
      END;
      T := (2.0 - M - Scalb(E, M))/(2.0 + M - Scalb(E, M));
      answer := DoCordic(1.0, T, 0.0, -1, TRUE);      
      result := answer.z - Functions.Ln(2.0)*FLOAT(E)/2.0;
      IF x < 0.0 THEN result := -result END;
   END;
   RETURN result;
END ArcTanh;
 
(* -------------------------------- Logarithmic Functions ----------------------------------- *)
 
PROCEDURE Ln (x: REAL): REAL;                  
   VAR 
      frac  : REAL;
      exp   : INTEGER;
      answer: Reply;
      result: REAL;
BEGIN
   IF x <= 0.0 THEN
      result := -INF();
   ELSE
      exp := 0; frac := x;
      WHILE (frac < 0.5) OR (frac >= 1.0) DO
         IF frac >= 1.0 THEN
            frac := frac/2.0;
            INC(exp);
         ELSIF frac < 0.5 THEN
            frac := frac*2.0;
            DEC(exp);
         END;
      END;
      answer := DoCordic(1.0+frac, 1.0-frac, 0.0, -1, TRUE);      
      result := -2.0*answer.z + FLOAT(exp)*Functions.Ln(2.0);
   END;
   RETURN result;
END Ln;
 
PROCEDURE Exp (x: REAL): REAL;
   VAR 
      answer  : Reply;
      quadrant: INTEGER;
      delta   : REAL;
      result  : REAL;
BEGIN
   Reduce(x, Functions.Ln(2.0), quadrant, delta);
   answer := DoCordic(Scale[-1], 0.0, delta, -1, FALSE);      
   result := Scalb(quadrant, answer.y+answer.x);
   IF x < 0.0 THEN result := 1.0/result END;            (* can easily overflow... *)
   RETURN result;
END Exp;
 
(* ----------------------------------- Power Functions -------------------------------------- *)
 
PROCEDURE Sqrt (x: REAL): REAL;                  
   VAR 
      frac  : REAL;
      exp   : INTEGER;
      answer: Reply;
      result: REAL;
BEGIN
   IF x <= 0.0 THEN
      result := 0.0;
   ELSE
      exp := 0; frac := x;
      WHILE (frac < 0.5) OR (frac >= 1.0) DO
         IF frac >= 1.0 THEN
            frac := frac/2.0;
            INC(exp);
         ELSIF frac < 0.5 THEN
            frac := frac*2.0;
            DEC(exp);
         END;
      END;
      IF ODD(exp) THEN frac := frac/2.0 END;
        answer := DoCordic(frac+0.25, frac-0.25, 0.0, -1, TRUE);   
      IF ODD(exp) THEN
         result := Scalb((exp+1) DIV 2, answer.x*Scale[-1]);
      ELSE
         result := Scalb(exp DIV 2, answer.x*Scale[-1]);
      END;
   END;
   RETURN result;
END Sqrt;
 
   (* This doesn't use cordic math but is nonetheless very handy!
      The algorithm is from "SemiNumerical Algortihms" by Knuth. *)
PROCEDURE RaiseInt (x: REAL; n: INTEGER): REAL;
   VAR
      exp : INTEGER;
      Y, Z: REAL;
BEGIN
   exp := n; Y := 1.0; Z := x;
   LOOP
      IF ODD(exp) THEN
         exp := exp DIV 2;
         Y := Z*Y;
         IF exp = 0 THEN EXIT END;
      ELSE
         exp := exp DIV 2;
      END;
      Z := Z*Z;
   END;
   RETURN Y;
END RaiseInt;
 
 
BEGIN
   NoBits := 32;            (* default to 32 bits of precision *)
END Cordic.
 
  jesjones@milton.u.washington.edu