{
    $Id$
    This file is part of the Free Pascal run time library.
    Copyright (c) 2000 by Jonas Maebe and other members of the
    Free Pascal development team

    Implementation of mathamatical Routines (only for real)

    See the file COPYING.FPC, included in this distribution,
    for details about the copyright.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

 **********************************************************************}


{****************************************************************************
                       EXTENDED data type routines
 ****************************************************************************}

    function pi : double;[internconst:in_pi];
      begin
        pi := 3.14159265358979320;
      end;

    function abs(d : extended) : extended;[internproc:in_abs_extended];
    function sqr(d : extended) : extended;[internproc:in_sqr_extended];
    function sqrt(d : extended) : extended;[internproc:in_sqrt_extended];

    function arctan(d : extended) : extended;[internconst:in_arctan_extended];
      begin
        runerror(207);
      end;

    function ln(d : extended) : extended;[internconst:in_ln_extended];
      begin
        runerror(207);
      end;

    function sin(d : extended) : extended;[internconst: in_sin_extended];
      begin
        runerror(207);
      end;

    function cos(d : extended) : extended;[internconst:in_cos_extended];
      begin
        runerror(207);
      end;

    function exp(d : extended) : extended;assembler;[internconst:in_const_exp];
      begin
        runerror(207);
      end;


    function frac(d : extended) : extended;assembler;[internconst:in_const_frac];
      begin
        runerror(207);
      end;


    function int(d : extended) : extended;assembler;[internconst:in_const_int];
      begin
        runerror(207);
      end;


    function trunc(d : extended) : longint;assembler;[internconst:in_const_trunc];
      { input: d in fr1      }
      { output: result in r3 }
      assembler;
      var
        temp:
          record
            case byte of
              0: (l1,l2: longint);
              1: (d: double);
          end;
      asm
        fctiwz   fr1,fr1
        stfd     fr1,temp.d
        lwz      r3,temp.l2
      end ['r3','fr1'];


    function round(d : extended) : longint;assembler;[internconst:in_const_round];
      { input: d in fr1      }
      { output: result in r3 }
      assembler;
      var
        temp:
          record
            case byte of
              0: (l1,l2: longint);
              1: (d: double);
          end;
      asm
        fctiw    fr1,fr1
        stfd     fr1,temp.d
        lwz      r3,temp.l2
      end ['r3','fr1'];


   function power(bas,expo : extended) : extended;
     begin
        if bas=0 then
          begin
            if expo<>0 then
              power:=0.0
            else
              HandleError(207);
          end
        else if expo=0 then
         power:=1
        else
        { bas < 0 is not allowed }
         if bas<0 then
          handleerror(207)
         else
          power:=exp(ln(bas)*expo);
     end;


{****************************************************************************
                       Longint data type routines
 ****************************************************************************}

   function power(bas,expo : longint) : longint;
     begin
        if bas=0 then
          begin
            if expo<>0 then
              power:=0
            else
              HandleError(207);
          end
        else if expo=0 then
         power:=1
        else
         begin
           if bas<0 then
            begin
              if odd(expo) then
               power:=-round(exp(ln(-bas)*expo))
              else
               power:=round(exp(ln(-bas)*expo));
            end
           else
            power:=round(exp(ln(bas)*expo));
         end;
     end;

{****************************************************************************
                    Helper routines to support old TP styled reals
 ****************************************************************************}

    { warning: the following converts a little-endian TP-style real }
    { to a big-endian double. So don't byte-swap the TP real!       }
    function real2double(r : real48) : double;

      var
         res : array[0..7] of byte;
         exponent : word;

      begin
         { copy mantissa }
         res[6]:=0;
         res[5]:=r[1] shl 5;
         res[4]:=(r[1] shr 3) or (r[2] shl 5);
         res[3]:=(r[2] shr 3) or (r[3] shl 5);
         res[2]:=(r[3] shr 3) or (r[4] shl 5);
         res[1]:=(r[4] shr 3) or (r[5] and $7f) shl 5;
         res[0]:=(r[5] and $7f) shr 3;

         { copy exponent }
         { correct exponent: }
         exponent:=(word(r[0])+(1023-129));
         res[1]:=res[1] or ((exponent and $f) shl 4);
         res[0]:=exponent shr 4;

         { set sign }
         res[0]:=res[0] or (r[5] and $80);
         real2double:=double(res);
      end;


{****************************************************************************
                         Int to real helpers
 ****************************************************************************}

function fpc_int64_to_real(i: int64): double; compilerproc;
assembler;
{ input: high(i) in r3, low(i) in r4 }
{ output: double(i) in f0            }
var
  temp:
    record
      case byte of
        0: (l1,l2: cardinal);
        1: (d: double);
    end;
asm
           lis    r0,$4330
           stw    r0,temp.l1
           xoris  r3,r3,$8000
           stw    r3,temp.l2
           lis    r3,longint_to_real_helper@ha
           lfd    fr1,longint_to_real_helper@l(r3)
           lfd    fr0,temp.d
           stw    r4,temp.l2
           fsub   fr0,fr0,fr1
           lis    r4,cardinal_to_real_helper@ha
           lfd    fr1,cardinal_to_real_helper@l(r4)
           lis    r3,int_to_real_factor@ha
           lfd    fr3,temp
           lfd    fr2,int_to_real_factor@l(r3)
           fsub   fr3,fr3,fr1
           fmadd  fr1,fr0,fr2,fr3
end ['r0','r3','r4','fr0','fr1','fr2','fr3'];


function fpc_qword_to_real(q: qword): double; compilerproc;
assembler;
{ input: high(q) in r3, low(q) in r4 }
{ output: double(q) in f0            }
var
  temp:
    record
      case byte of
        0: (l1,l2: cardinal);
        1: (d: double);
    end;
asm
           lis    r0,$4330
           stw    r0,temp.l1
           stw    r3,temp.l2
           lfd    fr0,temp.d
           lis    r3,cardinal_to_real_helper@ha
           lfd    fr1,cardinal_to_real_helper@l(r3)
           stw    r4,temp.l2
           fsub   fr0,fr0,fr1
           lfd    fr3,temp
           lis    r3,int_to_real_factor@ha
           lfd    fr2,int_to_real_factor@l(r3)
           fsub   fr3,fr3,fr1
           fmadd  fr1,fr0,fr2,fr3
end ['r0','r3','fr0','fr1','fr2','fr3'];



{
  $Log$
  Revision 1.3  2001-12-02 16:19:45  jonas
    * fpu results are returned in fr1, not fr0

  Revision 1.2  2001/10/30 17:18:14  jonas
    * fixed fpc_int64_to_double and fpc_int64_to_double (fpc_int64_to_double
      is now mostly tested and should work fine, fpc_qword_to_double should
      work too since it's almost the same)

  Revision 1.1  2001/10/28 14:09:13  jonas
    + initial implementation, lots of things still missing


}