Unit MATH387;

(* Bibliotheque mathematique pour utilisation du coprocesseur flottant
  JD GAYRARD Sept. 95
 la bibliotheque est batie … partir des fonctions du coprocesseur
 du type 386, elle fournit les fonctions suivantes:
      fsin, fcos, ftan, farctan, farctan2,
      farcsin, farccos, fmod, mod_2PI,
      ten_to, y_to_x, fexp, fln, flog...

Aucune verification du domaine de definition des fonctions n'est faite,
pas plus qu'un controle de la validite des operandes. Il est conseille
d'utiliser cette bibliotheque pour les types single et double exclusivement *)

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

{ table opcode du 387 non comprise par turbo pascal V7 }
{ FSIN    : D9 FE
  FCOS    : D9 FF
  FSINCOS : D9 FB
  FPREM1  : D9 F5 }

interface

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

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

(* use only with 86387, 86486 or pentium for type single, double andt extended,
no check of definition domain of the function or range (FPU limitation).
The f prefix avoids function redefinition of system runtime library *)

function fsin(x : float): float;
function fcos(x : float): float;
procedure sincos(x : float; var sinus, cosinus : float);
function ftan(x : float): float;
function farcsin(x : float): float;
function farccos(x : float): float;
function farctan(x : float): float;
function farctan2(x, y : float): float;   { return arctan(y/x) }
function fmod(x , y : float): float;      { return x mod y }
function fmod_2PI(x : float): float;
function fln(x : float): float;
function flog(x : float): float;
function fexp(x : float): float;
function ften_to(x : float) : float;      { return 10**x }
function fy_to_x(y , x : float): float;   { return y**x }
function module(x, y : float): float;

implementation

function fsin(x : float): float; assembler;
{if x < pi.2^62, then C2 is set to 0 and ST = sin(x)
                 else C2 is set to 1 and ST = x }
{no check range validity is performed in this function}
{$IFDEF fpc}
{$ASMMODE intel}
{$ENDIF}
asm
   FLD x           { load x }
   DB $D9, $FE     { opcode for FSIN }
end;

function fcos(x : float): float; assembler;
{ if x < pi.2^62, then C2 is set to 0 and ST = sin(x)
                  else C2 is set to 1 and ST = x }
{no range validity check is performed in this function}
{$IFDEF fpc}
{$ASMMODE intel}
{$ENDIF}
asm
   FLD x         { load angle }
   DB $D9, $FF   { opcode for FCOS }
end;

procedure sincos(x : float; var sinus, cosinus : float);
(* retourne sinus et cosinus(x), utilisable uniquement
 avec 80387, 80468 et pentium *)
var lcos, lsin : float;
begin
  {$IFDEF fpc}
  {$ASMMODE intel}
  {$ENDIF}
     asm                 { ST(0)     ST(1) }
          FLD x          {  x          -   }
          DB $D9, $FB    { cos(x)   sin(x) }
          FSTP lcos      { sin(x)      -   }
          FSTP lsin      {  -          -   }
     end;
cosinus := lcos;
sinus := lsin
end;

function ftan(x : float): float; assembler;
{ if x < pi.2^62, then C2 is set to 0 and ST = 1 and ST(1) = tan(x)
                  else C2 is set to 1 and ST = x }
{no range validity check is performed in this function}
{$IFDEF fpc}
{$ASMMODE intel}
{$ENDIF}
asm              { ST(0)    ST(1) }
   FLD x         {  x        -    }
   FPTAN         {  1      tan(x) }
   FSTP ST(0)    { tan(x)    -    }
end;

function farcsin(x : float): float; assembler;
(* retourne l'arcsin de x *)
{  methode :                ________
   arcsin(x) = arctan( x / V 1 - x.x ) }
{no range validity check is performed in this function |x| > 1 }
{$IFDEF fpc}
{$ASMMODE intel}
{$ENDIF}
asm                 { ST(0)     ST(1)     ST(2) }
   FLD X            {  x         -         -    }
   FLD ST(0)        {  x         x         -    }
   FMUL ST(0), ST   { x.x        x         -    }
   FLD1             {  1         x.x       x    }
   FSUBRP ST(1), ST { 1 - xý      x        -    }
   FSQRT            { sqrt(1-xý)  x        -    }
   FPATAN           { arcsin(x)   -        -    }
end;

function farccos(x : float): float; assembler;
{ retourne arccos(x)
   methode :            ________
   arcsin(x) = arctan( V 1 - x.x / x ) }
{ pas de controle de domaine de definition |x| > 1 }
{$IFDEF fpc}
{$ASMMODE intel}
{$ENDIF}
asm                 { ST(0)     ST(1)     ST(2) }
   FLD X            {  x         -         -    }
   FLD ST(0)        {  x         x         -    }
   FMUL ST(0), ST   { x.x        x         -    }
   FLD1             {  1         x.x       x    }
   FSUBRP ST(1), ST { 1 - xý      x        -    }
   FSQRT            { sqrt(1-xý)  x        -    }
   FXCH             {   x         z        -    }
   FPATAN           { arccos(x)   -        -    }
end;

function farctan(x : float): float; assembler;
{$IFDEF fpc}
{$ASMMODE intel}
{$ENDIF}
asm              { ST(0)    ST(1) }
   FLD x         {  x         -   }
   FLD1          {  1         x   }
   FPATAN        { atan(x/1)  -   }
end;

function farctan2(x, y : float): float; assembler;
{ retoune arctan (y / x) }
{$IFDEF fpc}
{$ASMMODE intel}
{$ENDIF}
asm              { ST(0)    ST(1) }
   FLD y         {  y         -   }
   FLD x         {  x         y   }
   FPATAN        { atan(y/x)  -   }
end;

function fmod(x, y : float): float; assembler;
(* retourne x mod y *)
{$IFDEF fpc}
{$ASMMODE intel}
{$ENDIF}
asm                 { ST(0)    ST(1) }
   FLD Y            {  y        -    }
   FLD X            {  x        y    }
@repeat_mod:
   FPREM            { x mod y   y    }
   FSTSW AX
   SAHF
   JP @repeat_mod
   FSTP ST(1)       { x mod y   -    }
end;

function fmod_2PI( x : float): float; assembler;
(* retourne x mod 2.pi *)
{$IFDEF fpc}
{$ASMMODE intel}
{$ENDIF}
asm                   { ST(0)    ST(1) }
   FLDPI              {  pi        -   }
   FADD ST, ST        { 2.pi       -   }
   FLD x              {  x        2.pi }
@unit_circle:
   FPREM              { x mod 2pi  2pi }
   FSTSW AX
   SAHF
   JP @unit_circle
   FSTP ST(1)         { x mod 2pi   -  }
end;

function fln(x : float): float; assembler;
{ retourne le logarithme naturel de x, utilise
 la methode loge(x) = loge(2).log2(x)}
{ pas de verification du domaine de definition (x < 0) }
{$IFDEF fpc}
{$ASMMODE intel}
{$ENDIF}
asm             {  ST(0)          ST(1)  }
   FLDLN2       { ln(2)            -     }
   FLD X        {   x             ln(2)  }
   FYL2X        { ln(2).log2(x)    -     }
end;

function flog(x : float): float; assembler;
{ retourne le logarithme base 10 de x, utilise
 la methode log10(x) = log10(2).log2(x) }
{ pas de verification du domaine de definition (x < 0) }
{$IFDEF fpc}
{$ASMMODE intel}
{$ENDIF}
asm             {  ST(0)          ST(1)  }
   FLDLG2       { log10(2)         -     }
   FLD X        {   x           log10(2) }
   FYL2X        {log2(x).log10(2)  -     }
end;

function fexp(x : float): float; assembler;
{ retourne e^x, par la methode e^x = 2^(y.log2(e)) }
{ 2^z = 2^f.2^i with f = frac(z) and i = int(z) }
{ 2^f is computed with F2XM1, 2^i with FSCALE }
const round_down : word = $177F;
var control_ww : word;
{$IFDEF fpc}
{$ASMMODE intel}
{$ENDIF}
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          -   }
   FSTP ST(1)         { 10^x       -          -   }
end;

function ften_to(x : float): float; assembler;
{ retourne 10^x, par la methode 10^x = 2^(y.log2(10)) }
{ 2^z = 2^f.2^i with f = frac(z) and i = int(z)       }
{ 2^f is computed with F2XM1, 2^i with FSCALE }
const round_down : word = $177F;
var control_ww : word;
{$IFDEF fpc}
{$ASMMODE intel}
{$ENDIF}
asm                   { ST(0)     ST(1)     ST(2) }
   FLD X              {  x         -         -    }
   FLDL2T             { log2(10)   x         -    }
   FMULP ST(1), ST    { x.log2(10) -         -    }
   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          -   }
   FSTP ST(1)         { 10^x       -          -   }
end;

function fy_to_x(y, x : float): float; assembler;
{ retourne y^x, par la methode y^x = 2^(y.log2(y))    }
{no range validity check is performed in this function (y > 0) }
{ 2^z = 2^f.2^i with f = frac(z) and i = int(z)  }
{ 2^f is computed with F2XM1, 2^i with FSCALE }
const round_down : word = $177F;
var   control_ww : word;
{$IFDEF fpc}
{$ASMMODE intel}
{$ENDIF}
asm                   { ST(0)     ST(1)     ST(2) }
   FLD Y              {  y         -         -    }
   FLD X              {  x         y         -    }
   FYL2X              { x.log2(y)  -         -    }
   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          -   }
   FSTP ST(1)         { y^x        -          -   }
end;

function module(x, y : float): float; assembler;
(* retourne le module du complexe (x,y) *)
{$IFDEF fpc}
{$ASMMODE intel}
{$ENDIF}
asm                  { ST(0)     ST(1) }
   FLD Y             {  y         -    }
   FMUL ST(0), ST    { y.y        -    }
   FLD X             {  x         y.y  }
   FMUL ST(0), ST    { x.x        y.y  }
   FADDP ST(1), ST   { d.d        -    }
   FSQRT             {  d         -    }
end;

end.