home *** CD-ROM | disk | FTP | other *** search
Modula Implementation | 1986-03-09 | 3.6 KB | 159 lines |
-
- IMPLEMENTATION MODULE MathLibS;
- (* *)
- (* Math routines that work with the 8087
- contributed by Steve Eckhart *)
- (* *)
- CONST
- pi = 3.14159265358979;
- twopi = pi*2.0;
- piover2 = pi*0.5;
- piover4 = pi*0.25;
- MaxErr = 1.0E-14;
- PROCEDURE sin(x : REAL) : REAL;
- VAR
- sx, sign : REAL;
- BEGIN
- sign := 1.0;
- IF x<0.0 THEN
- sign := -1.0;
- x := -x;
- END;
- WHILE x>pi DO
- x := x-twopi;
- END;
- IF x<0.0 THEN
- sign := -sign;
- x := -x;
- END;
- IF x>piover2 THEN
- x := pi-x;
- END;
- IF x<=piover4 THEN
- sx := SinExp(x);
- ELSE
- sx := CosExp(piover2-x);
- END;
- RETURN sign*sx;
- END sin;
- (*series expansion for sin(x) *)
- PROCEDURE SinExp(s : REAL) : REAL;
- VAR
- term, sum, iter : REAL;
- BEGIN
- term := s;
- sum := s;
- iter := 1.0;
- WHILE ABS(term)>MaxErr DO
- iter := iter+1.0;
- term := -term*s*s/((2.0*iter-2.0)*(2.0*iter-1.0));
- sum := sum+term;
- END;
- RETURN sum;
- END SinExp;
- PROCEDURE cos(x : REAL) : REAL;
- VAR
- cx, sign : REAL;
- BEGIN
- sign := 1.0;
- IF x<0.0 THEN
- x := -x;
- END;
- WHILE x>pi DO
- x := x-twopi;
- END;
- IF x<0.0 THEN
- x := -x;
- END;
- IF x>piover2 THEN
- x := pi-x;
- sign := -1.0;
- END;
- IF x<=piover4 THEN
- cx := CosExp(x);
- ELSE
- cx := SinExp(piover2-x);
- END;
- RETURN sign*cx;
- END cos;
- (*series expansion for cos(x) *)
- PROCEDURE CosExp(c : REAL) : REAL;
- VAR
- term, sum, iter : REAL;
- BEGIN
- term := 1.0;
- sum := 1.0;
- iter := 1.0;
- WHILE ABS(term)>MaxErr DO
- term := -term*c*c/(2.0*iter*(2.0*iter-1.0));
- sum := sum+term;
- iter := iter+1.0;
- END;
- RETURN sum;
- END CosExp;
- PROCEDURE atan(x : REAL) : REAL;
- VAR
- i : CARDINAL;
- Sign, Reduc, x2, Term, Sum : REAL;
- BEGIN
- IF x<0.0 THEN
- Sign := -1.0;
- x := -x;
- ELSE
- Sign := 1.0;
- END;
- IF x<0.4142 THEN
- Reduc := 0.0;
- ELSE
- x := (x-1.0)/(x+1.0);
- IF x<0.4142 THEN
- Reduc := 1.0;
- ELSE
- x := (x-1.0)/(x+1.0);
- Reduc := 2.0;
- END;
- END;
- IF ABS(x)<MaxErr THEN
- RETURN Sign*Reduc*piover4;
- END;
- x2 := x*x;
- i := 1;
- Term := 1.0;
- Sum := 1.0;
- WHILE ABS(Term)>MaxErr DO
- INC(i);
- Term := -x2*Term;
- Sum := Sum+Term/FLOAT(2*i-1);
- END;
- Sum := Sign*(x*Sum+Reduc*piover4);
- RETURN Sum;
- END atan;
- PROCEDURE sqrt(x : REAL) : REAL;
- VAR
- exp, guess, newguess : REAL;
- BEGIN
- IF x<=0.0 THEN
- RETURN 0.0;
- END;
- exp := 1.0;
- WHILE x>=100.0 DO
- x := x*0.01;
- exp := exp*10.0;
- END;
- WHILE x<1.0 DO
- x := x*100.0;
- exp := exp*0.1;
- END;
- IF x=1.0 THEN
- RETURN exp;
- END;
- newguess := 4.0;
- REPEAT
- guess := newguess;
- newguess := (x/guess+guess)*0.5;
- UNTIL ABS(newguess-guess)<MaxErr;
- RETURN newguess*exp;
- END sqrt;
- BEGIN
- END MathLibS.