Unit MATHLIB;

(* Bibliotheque mathematique pour type real
  JD GAYRARD Fev. 94
 la bibliotheque est batie … partir des fonctions :
  ARCTAN, COS, EXP, LN, SIN, SQRT
  elle fournit les fonctions :
 ARCCOS, ARCSIN, ARCTAN2, LOG, TAN, PUISSANCE, SIGNE, MAX, MIN *)

(* Revision 1.0 de Jul. 95 pour :
- passage en double
- ajout de la fonction pwr_int, factorielle
- correction de puissance *)

(* revision 1.1 de Sep. 95 pour :
- passage en float
- correction de log (test valeur negative)
- ajout de ceiling, floor, Uran et Gran *)

(* revision 1.2 de Oct 95 pour :
- ajout de ten_to, module, deg_to_rad, rad_to_deg *)

{$G+}
{$N+}
{$E-}

interface

const author  = 'GAYRARD J-D';
      version = 'ver 1.2 - 10/95';

const PI_2      =  1.570796326794896619231322;   { pi / 2 }
      PI_3      =  1.047197551196597746154214;   { pi / 3 }
      PI_4      =  0.7853981633974483096156608;  { pi / 4 }
      SQRT_PI   =  1.772453850905516027298167;   { sqrt(pi) }
      SQRT_2PI  =  2.506628274631000502415765;   { sqrt(2.pi) }
      TWO_PI    =  6.283185307179586476925287;   { 2.pi }
      LN_PI     =  1.144729885849400174143427;   { ln(pi) }
      LOG_PI    =  0.4971498726941338543512683;  { log(pi) }
      LOG_E     =  0.4342944819032518276511289;  { log(e) }
      LN_10     =  2.302585092994045684017991;   { ln(10) }
      E         =  2.718281828459045235360287;   { exp(1) }
      ONE_RAD   = 57.295779513082320876798155;   { 1 rad in ø }
      ONE_DEG   =  0.017453292519943295769237;   { 1ø in rad }

type float = double; { a modifier suivant l'utilisation }

(* utilisable avec tout types de reel et avec controle du domaine
de definition des fonctions *)

function tan(x : float): float;
function arcsin(x : float): float;
function arccos(x : float): float;
function arctan2(x, y : float): float;          { retourne arctan (y/x) }
function log(x : float): float;
function y_to_x( y, x : float): float;          { retourne y^x}
function signe(x, y : float): float;            { retourne x avec le signe de y }
function max(x, y : float): float;
function min(x, y : float): float;
function pwr_int(x : float; n : integer): float;   { retourne x^n }
function ten_to(x : float): float;                 { retourne 10^x }
function fac(n : integer): float;                  { retourne x! }
function Uran: float;                              { uniform law }
function Gran: float;                              { gaussian law }
function ceiling(x : float): float;
function floor(x : float): float;
function module(x, y : float): float;
function deg_to_rad(x : float): float;
function rad_to_deg(x : float): float;

implementation

function Uran: float;
(* loi uniforme *)
begin
uran := random
end;

function Gran: float;
(* loi gaussienne *)
var k : integer;
    sum : float;
begin
sum := 0.0;
for k := 0 to 16 do sum := sum + random;
gran := sum / 16.0
end;

function signe(x, y : float): float;
(* retourne x avec le signe de y *)
begin
if x > 0.0 then if y > 0.0 then signe := x
                           else signe := -x
           else if y < 0.0 then signe := x
                           else signe := -x
end;

function min(x, y : float): float;
(* retourne le plus petit *)
begin
if x > y then min := y
         else min := x
end;

function max(x, y : float): float;
(* retourne le plus grand *)
begin
if x < y then max := y
         else max := x
end;

function ceiling(x : float): float;
{ return the nearest integer value above x }
begin
if x <> int(x) then ceiling := int(x) + 1
               else ceiling := x
end;

function floor(x : float): float;
{ return the nearest integer value below x }
begin
if x <> int(x) then floor := int(x) - 1
               else floor := x
end;

function module(x, y : float): float;
{ retourne sqrt( x.x + y.y), distance du point (x,y) a l'origine (0,0) }
begin
module := sqrt(x * x + y * y)
end;

function tan(x : float): float;
(* retourne la tangente de x (en radian) *)
var cosx : float;
begin
cosx := cos(x);
if cosx = 0.0
   then begin
        writeln('******* Fonction tan ********');
        writeln('********* OVERFLOW **********');
        halt
        end
   else tan := sin(x) / cosx
end;

function arcsin(x : float): float;
(* retourne l'arcsin de x, x compris entre -1 et 1 *)
{                           ________
   arcsin(x) = arctan( x / V 1 - x.x ) }
begin
if (x > 1.0) or (x < -1.0)
   then begin
        writeln('****** Fonction arcsin ******');
        writeln('********* OVERFLOW **********');
        halt
        end
   else if x = 0.0
           then arcsin := 0.0
           else if x = 1.0
                   then arcsin := pi_2
                   else if x = -1.0
                           then arcsin := - pi_2
                           else arcsin := arctan(x / sqrt( 1.0 - x * x))
end;

function arccos(x : float): float;
(* retourne l'arccos de x, x compris entre -1 et 1 *)
{                       ________
   arcsin(x) = arctan( V 1 - x.x / x ) }
var y : float;
begin
if (x > 1.0) or (x < -1.0)
   then begin
        writeln('****** Fonction arccos ******');
        writeln('********* OVERFLOW **********');
        halt
        end
   else if x = 0.0
           then arccos := pi_2
           else if x = 1.0
                then arccos := 0.0
                else if x = -1.0
                        then arccos := pi
                        else begin
                             y := arctan(sqrt( 1.0 - x * x) / x);
                             if x > 0.0
                                then arccos := y
                                else arccos := y + pi;
                             end
end;

function arctan2(x, y : float): float;
{ retourne l'arctan de y/x }
begin
if x = 0.0
   then arctan2 := signe(pi_2, y)
   else if x > 0.0
        then arctan2 := arctan(y/x)
        else arctan2 := arctan(y/x) + signe(pi,y)
end;

function y_to_x (y, x : float): float;
(* retourne y^x, y positif par la methode e^x.ln(y) *)
begin
if y >= 0 then y_to_x := exp( x * ln(y))
          else begin
               writeln('****** Fonction puissance ******');
               writeln('****** NEGATIVE ARGUMENT *******');
               halt
               end
end;

function ten_to(x : float): float;
begin
ten_to := exp(x * LN_10)
end;

function log(x : float): float;
(* retourne de logarithme decimal de x, x positif
utilise la methode log10(x) = ln(x)/ln(10) *)
begin
if x >= 0 then log := log_E * ln(x)
          else begin
               writeln('********* Fonction log *********');
               writeln('****** NEGATIVE ARGUMENT *******');
               halt
               end
end;

function pwr_int(x : float; n : integer) : float;
{ retourne x^n, n entier, utilise la methode multiplicative }
var       temp : float;
          i : integer;
begin
if n = 0 then pwr_int := 1.0
         else
   if (x = 0.0) or (n = 1) then pwr_int := x
                           else
             begin
             temp := 1.0;
             for i := 1 to abs(n) do temp := temp * x;
             if n > 0 then pwr_int := temp
                      else pwr_int := 1.0 / temp
             end
end;

function fac(n : integer): float;
(* returne n! , n > 0 *)
var temp : float;
    i : integer;
begin
if n <= 0 then begin
               writeln('********* Fonction fac *********');
               writeln('****** NEGATIVE ARGUMENT *******');
               halt
               end
          else begin
               temp := 1.0;
               for i := 2 to n do temp := temp * i;
               fac := temp
               end
end;

function deg_to_rad(x : float): float;
{ conversion degres vers radians }
begin
deg_to_rad := one_deg * x
end;

function rad_to_deg(x : float): float;
{conversion radians vers degres }
begin
rad_to_deg := one_rad * x
end;

begin
randomize
end.