(************************************************

	   A LIBRARY OF TRIG FUNCTIONS

	  See TRIG.DOC for all comments.

*************************************************)

FUNCTION ARCCOS(X :REAL ) :REAL;
BEGIN
  -ARCTAN( X / SQRT(-X*X+1) ) +1.5708
END;

	(*******************************)

FUNCTION ARCSIN(X :REAL) :REAL;
BEGIN
  ARCTAN( X / SQRT(-X*X+1) )
END;

	(*******************************)

FUNCTION COSECANT(X :REAL) :REAL;
BEGIN
  COSECANT := 1/SIN(X)
END;

	(*******************************)

FUNCTION COSH(X :REAL) :REAL;
BEGIN
  COSH := (EXP(X) + EXP(-X)) /2
END;

	(*******************************)

FUNCTION LOG(B,X: REAL ) : REAL;
BEGIN
  LOG := LN(X) / LN(B)
END;

	(*******************************)

FUNCTION LOG10( X : REAL ) :REAL;
BEGIN
  LOG10 := LN(10) / LN(X)
END;

	(*******************************)

Function ARCTAN( x: real ): real;
{  Intrinsic function for calculating the ARC TANGENT  }
CONST
  a1 = 3.7092563;
  a2 = -7.10676;
  a3 = -0.26476862;
  b0 = 0.17465544;
  b1 = 6.762139;
  b2 = 3.3163354;
  b3 = 1.44863154;
VAR	i,k: real;
	signchg, adjust: boolean;
begin
  adjust := false;
  signchg := false;
  if x < 0.0 then
    begin
      signchg := true;
      x := -x
    end;
  if x > 1.0 then
    begin
      adjust := true;
      x := 1.0 / x
   end;
  i := x * x;
  k := x * (b0 + a1 / (b1 + i + a2 / (b2 + i + a3 / (b3 + i) ) ) );
  if adjust then k := halfpi - k;
  if signchg then k := - k;
  arctan := k
end;


Function LN( x: real ): real;
{ intrinsic function for calculating natural log }
const	c1 = 2.8853913;
	c3 = 0.9614706;
	c5 = 0.59897865;
	sqrtp5 = 0.7071068;
	loge2 = 0.6931472;
var	y: record
	     case char of
		'a': ( a: real );
		'b': ( b1, b2, b3, b4: 0..255 )
	   end;
    z, z2: real;
    i: integer;
begin
    If x = 1.0 then
      ln := 0.0
    Else
      begin
	y.a := x;
	i := (y.b1 + 128) mod 256 - 128;
	y.b1 := 0;
	z := (y.a - sqrtp5) / (y.a + sqrtp5);
	z2 := sqr( z );
	ln := (((c5*z2 + c3)*z2 + c1)*z - 0.5 + i) * loge2;
      end
end;


Function SQRT( x:real ): real;
{ Intrinsic function for square root }
VAR	j,i,k: real;
begin
  i := 0.0;
  If x < 0.0 then x := ABS( x );
  j := x / 2.0;
  If x <> 0.0 then
    Repeat
      k := i;
      i := j;
      j := (x / j + j )/2.0;
    Until (j = i) or (j = k);
  SQRT := j
End{of Square_Root};

	(*******************************)

FUNCTION Max(X, Y  :REAL):REAL;
BEGIN	IF X>Y THEN Max := X ELSE Max := Y
END;

	(*******************************)

FUNCTION POWER(X,Y :REAL) :REAL;
BEGIN
  POWER := EXP( LN(X) * Y)
END;

	(*******************************)

FUNCTION POWERI(X :REAL; N :INTEGER) :REAL;
  (* COMPUTES 'X' TO THE POWER 'N' RECURSIVELY.
	THE FUNCTION RETURNS 0 FOR X=0 *)
VAR T:REAL;
  BEGIN
  IF X=0.0 THEN
    POWERI:=0.0
  ELSE IF N=0 THEN
	 POWERI:=1.0
       ELSE IF N=1 THEN
	      POWERI:=X
	    ELSE IF N<0 THEN
		   BEGIN
		     T:=POWERI(1/X,(-N) DIV 2));
		     POWERI:=T*T*POWERI(1/X,(-N) MOD 2)
		   END
		 ELSE BEGIN
			T:=POWERI(X, N DIV 2);
			POWERI:=T*T*POWERI(X, N MOD 2)
		      END
END(* POWERI *);

	(*******************************)

(* Another version of the Integer Power function *)
FUNCTION POWERI(x :REAL; y :INTEGER) :REAL;
VAR
  w,z	:REAL;
  i	:INTEGER;
BEGIN
w := x;
z := 1;
i := y;
WHILE i > 0 DO
  BEGIN (*   z * (w ** i) = x ** y   *)
    IF ODD(i) THEN z := z * w;
    i := i DIV 2;
    w := SQR(w)
  END;
Poweri := z  (*   z = x ** y	*)
END.

	(*******************************)

PROCEDURE SEEDRAND;
(* INITIAL VALUES FOR SEED1 AND SEED2 MAY BE
   INPUT HERE  *)
BEGIN
   SEED1 := 10946;
   SEED2 := 17711;
END;

FUNCTION RANDOM : INTEGER;
(* Implement a Fibonacci series Random number generator.
   Written for PASCAL/Z By Raymond E. Penley, September 1979
   RANDOM will return numbers from 0 to 32767
   Call RANDOM with the following convention:
	 Range		 Use
	  0 - 32	RANDOM DIV 1000
	  0 - 327	RANDOM DIV 100
	  0 - 32767	RANDOM

GLOBAL
	SEED1, SEED2 : INTEGER
*)
CONST
	HALFINT = 16383; (* 1/2 OF MAXINT *)
VAR
	h1, h2, hadded : INTEGER;

BEGIN
  h1 := SEED1 DIV 2;
  h2 := SEED2 DIV 2;
  IF (h1+h2) >= HALFINT THEN
    hadded := h1 + h2 - HALFINT
  ELSE
    hadded := h1 + h2;
  SEED1 := SEED2;
  SEED2 := hadded * 2;(* Restore from previous DIVision *)
  RANDOM := SEED2
END(* RANDOM *);

	(*******************************)

FUNCTION SECANT( X :REAL ) :REAL;
BEGIN
  SECANT := 1/COS(X)
END;


FUNCTION SINH( X :REAL ) :REAL;
BEGIN
  SINH := (EXP(X) - EXP(-X))/2
END;

	(*******************************)

FUNCTION SQUAREROOT (VALUE : REAL ) : REAL;
CONST	EPSILON = 0.000001; (* 1E-6 *)
VAR	ROOT : REAL;
BEGIN
  IF VALUE = 0 THEN
    SQUAREROOT := 0
  ELSE BEGIN
	 ROOT:=1;
	 REPEAT
	   ROOT:=(VALUE/ROOT+ROOT)/2
	 UNTIL ABS(VALUE/SQR(ROOT)-1) < EPSILON;
	 SQUAREROOT := ROOT
       END
END(*** SQUAREROOT ***);

	(*******************************)

FUNCTION TAN(X :REAL) :REAL;
BEGIN
  TAN := SIN(X)/COS(X)
END;


FUNCTION TANH(X :REAL) :REAL;
BEGIN
  TANH := SINH(X)/COSH(X)
END;

	(*******************************)

FUNCTION COMPOUNDINT(PRINCIPAL, RATE:REAL;
		     PERIODNBR:INTEGER):REAL;
BEGIN
   COMPOUNDINT:=PRINCIPAL * POWERI(1.0+RATE, PERIODNBR)
END(* COMPOUNDINT *);

	(*******************************)

FUNCTION GCD ( NUM, DEN : INTEGER ) :INTEGER;
VAR	REMAINDER : INTEGER;
BEGIN
  WHILE DEN <> 0 DO
  BEGIN
    REMAINDER := NUM MOD DEN;
    NUM := DEN;
    DEN := REMAINDER
  END;(* WHILE *)
GCD := NUM
END; (*  GCD  *)

	(*******************************)

PROCEDURE LOWTERM ( VAR NUM, DEN : INTEGER );
(* MUST BE USED WITH FUNCTION GCD *)
VAR	DIVISOR : INTEGER;
BEGIN
  DIVISOR := GCD(NUM, DEN);
  IF DIVISOR > 1 THEN
    BEGIN
      NUM := NUM DIV DIVISOR;
      DEN := DEN DIV DIVISOR
    END
END; (*  LOWTERM  *)

	(*******************************)
   NUM := NUM DIV DIVISOR;
      DEN := DEN DIV DIVISOR
    END
END; (*  LOWTERM  *)

	(**********