Unit HYPER387;

(* Bibliotheque mathematique des fonctions hyperboliques *)
(* JD GAYRARD oct 95 *)
(* utilisables uniquement avec un 86387, 86486 et pentium,
pour type single, double et extended, sans controle de domaine de
definition (de la fonction) ou d'utilisation (limitation du FPU).
le prefixe f est pour eviter la redefinition *)

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

interface

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

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

(* fonctions trigonometriques directes *)
function fch(x : float): float;
function fsh(x : float): float;
function fth(x : float): float;

(* fonctions trigonometriques inverses *)
function farg_ch(x : float): float;
function farg_sh(x : float): float;
function farg_th(x : float): float;

implementation

(* fonctions trigonometriques directes *)

function fch(x : float): float; assembler;
(* retourne le cosinus hyperbolique de l'argument *)
{ ch(x) = [exp(x) + exp(-x)] / 2
methode : z = exp(x), ch(x) = 1/2 (z + 1/z)
          z = 2^y, y = x.log2(e),
          z = 2^f.2^i, f = frac(y), i = int(y)
{ 2^f is computed with F2XM1, 2^i with FSCALE }
const round_down : word = $177F;
      one_half : float = 0.5;
var   control_ww : word;
asm                   { ST(0)     ST(1)     ST(2) }
   FLD X              {  x         -         -    }
   FLDL2E             { log2(e)    x         -    }
   FMULP ST(1), ST    { x.log2(e)  -         -    }
   FSTCW control_ww
   FLDCW round_down
   FLD ST(0)          {  z         z          -   }
   FRNDINT            { int(z)     z          -   }
   FLDCW control_ww
   FXCH               {  z         i          -   }
   FSUB ST, ST(1)     {  f         i          -   }
   F2XM1              { 2^f-1      i          -   }
   FLD1               {  1        2^f-1       i   }
   FADDP ST(1), ST    { 2^f        i          -   }
   FSCALE             { 2^f.2^i    i          -   }
   FST ST(1)          { e^x       e^x         -   }
   FLD1               {  1         z          z   }
   FDIVRP ST(1), ST   { 1/z        z          -   }
   FADDP ST(1), ST    { z+1/z      -          -   }
   FLD one_half       { 0.5       z+1/z       -   }
   FMULP ST(1), ST    { ch(x)      -          -   }
end;

function fsh(x : float): float; assembler;
(* retourne le sinus hyperbolique de l'argument *)
{ sh(x) = [exp(x) - exp(-x)] / 2
methode : z = exp(x), ch(x) = 1/2 (z - 1/z)
          z = 2^y, y = x.log2(e),
          z = 2^f.2^i, f = frac(y), i = int(y)
{ 2^f is computed with F2XM1, 2^i with FSCALE }
const round_down : word = $177F;
      one_half : float = 0.5;
var   control_ww : word;
asm                   { ST(0)     ST(1)     ST(2) }
   FLD X              {  x         -         -    }
   FLDL2E             { log2(e)    x         -    }
   FMULP ST(1), ST    { x.log2(e)  -         -    }
   FSTCW control_ww
   FLDCW round_down
   FLD ST(0)          {  y         y          -   }
   FRNDINT            { int(y)     y          -   }
   FLDCW control_ww
   FXCH               {  y         i          -   }
   FSUB ST, ST(1)     {  f         i          -   }
   F2XM1              { 2^f-1      i          -   }
   FLD1               {  1        2^f-1       i   }
   FADDP ST(1), ST    { 2^f        i          -   }
   FSCALE             { 2^f.2^i    i          -   }
   FST ST(1)          { e^x       e^x         -   }
   FLD1               {  1         z          z   }
   FDIVRP ST(1), ST   { 1/z        z          -   }
   FSUBP ST(1), ST    { z-1/z      -          -   }
   FLD one_half       { 0.5      z-1/z)       -   }
   FMULP ST(1), ST    { sh(x)      -          -   }
end;

function fth(x : float): float; assembler;
(* retourne la tangente hyperbolique de l'argument *)
(* th(x) = sh(x) / ch(x) *)
{ th(x) = [exp(x) - exp(x)] / [exp(x) + exp(-x)]
methode : z = exp(x), ch(x) = (z - 1/z) / (z + 1/z)
          z = 2^y, y = x.log2(e),
          z = 2^f.2^i, f = frac(y), i = int(y)
{ 2^f is computed with F2XM1, 2^i with FSCALE }
const round_down : word = $177F;
      one_half : float = 0.5;
var   control_ww : word;
asm                   { ST(0)     ST(1)     ST(2) }
   FLD X              {  x         -         -    }
   FLDL2E             { log2(e)    x         -    }
   FMULP ST(1), ST    { x.log2(e)  -         -    }
   FSTCW control_ww
   FLDCW round_down
   FLD ST(0)          {  z         z          -   }
   FRNDINT            { int(z)     z          -   }
   FLDCW control_ww
   FXCH               {  z         i          -   }
   FSUB ST, ST(1)     {  f         i          -   }
   F2XM1              { 2^f-1      i          -   }
   FLD1               {  1        2^f-1       i   }
   FADDP ST(1), ST    { 2^f        i          -   }
   FSCALE             { 2^f.2^i    i          -   }
   FST ST(1)          { e^x       e^x         -   }
   FLD1               {  1         z          z   }
   FDIV  ST, ST(1)    { 1/z        z          z   }
   FSUB  ST(2), ST    { 1/z        z        z-1/z }
   FADDP ST(1), ST    { z+1/z    z-1/z        -   }
   FDIVP ST(1), ST    { th(x)      -          -   }
end;

(* fonctions trigonometriques inverses *)

function farg_ch(x : float): float; assembler;
(* retourne l'arc cosinus hyperbolique de l'argument *)
(*                       ________          *)
(* arg ch(x) = ln ( x + V x.x - 1 )  x >=1 *)
asm                 {  ST(0)         ST(1)          ST(2)  }
   FLDLN2           {  ln(2)          -              -     }
   FLD X            {   x            ln(2)           -     }
   FLD ST(0)        {   x             x             ln(2)  }
   FMUL ST(0), ST   {   x.x           x             ln(2)  }
   FLD1             {   1             x.x            x     }
   FSUBP ST(1), ST  { x.x - 1         x             ln(2)  }
   FSQRT            { sqrt(x2-1)      x             ln(2)  }
   FADDP ST(1), ST  { x + z          ln(2)           -     }
   FYL2X            { arg_ch(x)       -              -     }
end;

function farg_sh(x : float): float; assembler;
(* retourne l'arc sinus hyperbolique de l'argument *)
(*                       _________   *)
(* arg sh(x) = ln ( x + V x.x + 1 )  *)
asm                 {  ST(0)         ST(1)          ST(2)  }
   FLDLN2           {  ln(2)          -              -     }
   FLD X            {   x            ln(2)           -     }
   FLD ST(0)        {   x             x             ln(2)  }
   FMUL ST(0), ST   {   x.x           x             ln(2)  }
   FLD1             {   1             x.x            x     }
   FADDP ST(1), ST  { x.x + 1         x             ln(2)  }
   FSQRT            { sqrt(x.x+1)     x             ln(2)  }
   FADDP ST(1), ST  { x + z          ln(2)           -     }
   FYL2X            { arg_sh(x)       -              -     }
end;

function farg_th(x : float): float; assembler;
(* retourne l'arc tangente hyperbolique de l'argument *)
(* arg th(x) = 1/2 ln [ (1 + x) / (1 - x) *)
asm                 {  ST(0)         ST(1)          ST(2)  }
   FLDLN2           {  ln(2)          -              -     }
   FLD X            {   x            ln(2)           -     }
   FLD ST(0)        {   x             x             ln(2)  }
   FLD1             {   1             x              x     }
   FADDP ST(1), ST  { 1 + x           x             ln(2)  }
   FXCH             {   x            1 + x          ln(2)  }
   FLD1             {   1             x             1 + x  }
   FSUBRP ST(1), ST { 1 - x          1 + x          ln(2)  }
   FDIVP ST(1), ST  { 1+x/1-x        ln(2)           -     }
   FSQRT            {                ln(2)           -     }
   FYL2X            { ln(z)           -              -     }
end;

end.