unit float;
{$IFDEF FPC}
  {$MODE DELPHI}{$H+}
{$ENDIF}
(***************************************)
(* Copyright (C) 2003, SHIRAISHI Kazuo *)
(***************************************)
interface
uses SysUtils,math;
type
    floatFunction1=procedure (var x:double);
    floatFunction2=procedure (var x,y:double);
function LongIntRound(x:double):longint;
procedure  add(var x,y:double);
procedure  sbt(var x,y:double);
procedure  mlt(var x,y:double);
procedure  qtt(var x,y:double);
procedure  power(var x,y:double);
procedure  opposite(var x:double);
procedure  basicmod(var x,y:double);
procedure  square(var x:double);
procedure  FMAX(var x,y:double);
procedure  FMIN(var x,y:double);
procedure  FABS  (var x:double);
procedure  FCEIL (var x:double);
procedure  FFLOOR(var x:double);
procedure  FSQRT (var x:double);
procedure  FROUND(var x:double);
procedure  FEPS  (var x:double);
procedure  FSIN  (var x:double);
procedure  FCOS  (var x:double);
procedure  FTAN  (var x:double);
procedure  FCOT  (var x:double);
function fcompare(var x,y:double):integer;
function fsign(var x:double):integer;
function NPXpower(x,y:double):double;
function NPXpower1plus(x,y:double):double;
//procedure invalidoperation;
implementation
 uses
      base;

function LongIntRound(x:double):longint;
begin
   if (x>maxint) or (x<minint) then
           raise EInValidOp.create('');
   result:=System.Round(x);
end;
procedure  opposite(var x:double);
begin
  x:=-x
end;
procedure  add(var x,y:double);
begin
   x:=x+y
end;
procedure  sbt(var x,y:double);
begin
   x:=x-y
end;
procedure  mlt(var x,y:double);
begin
   x:=x*y
end;
procedure  qtt(var x,y:double);
begin
   x:=x/y
end;
procedure square(var x:double);
begin
 x:=sqr(x)
end;
function fsign(var x:double):integer;
begin
 if x>0 then
    result:=1
 else if x<0 then
    result:=-1
 else
    result:=0
end;
function fcompare(var x,y:double):integer;
var
   z:double;
begin
   z:=x-y;
   result:=fsign(z)
end;
function round(x:double):double;
begin
 if (x>-$7fffffff) and (x<$7fffffff) then
    result:=system.round(x)
 else
    result:=x
end;
function floor(x:double):double;
begin
  result:=system.int(x);
  if (x<0) and (result<>x) then
     result:=result-1
end;
function ceil(x:double):double;
begin
    result:=-floor(-x)
end;
procedure  BasicMod(var x,y:double);
var
   q:double;
begin
   q:=floor(x/y);
   x:=x-y*q;
end;
procedure  FMAX(var x,y:double);
begin
  if x<y then x:=y
end;
procedure  FMIN(var x,y:double);
begin
 if x>y then x:=y
end;
procedure  FABS  (var x:double);
begin
  x:=abs(x)
end;
procedure  FCEIL (var x:double);
begin
  x:=ceil(x)
end;
procedure FFLOOR  (var x:double);
begin
  x:=Floor(x)
end;
procedure  FROUND(var x:double);
begin
  x:=round(x)
end;
var
   number2:array[0..3]of word=($ffff,$ffff,$ffff,$7fef);
var
   number0:double absolute number2;
procedure  FEPS(var x:double);
//const
//   number2:array[0..3]of word=($ffff,$ffff,$ffff,$7fef);
//var
//   number0:double absolute number2;
var
   e:word;
begin
    number0:=x;
    e:=(number2[3] and $7ff0) div $10;
    if e>0 then
       begin
         number2[3]:=e*$10 ;
         number2[2]:=0;
         number2[1]:=0;
         number2[0]:=0;
         x:=number0/4503599627370496.
       end
    else
       begin
         number2[3]:=0;
         number2[2]:=0;
         number2[1]:=0;
         number2[0]:=1;
         x:=number0
       end;
end;
procedure  FSQRT (var x:double);
begin
 x:=sqrt(x)
end;
procedure  FSIN (var x:double);
begin
  x:=sin(x)
end;
procedure  FCOS (var x:double);
begin
 x:=cos(x)
end;
procedure  FTAN (var x:double);
begin
 x:=math.tan(x)
end;
procedure  FCOT (var x:double);
begin
 x:=math.cot(x)
end;
procedure FASIN (var x:double);
begin
  x:=math.arcsin(x);
end;
procedure FACOS (var x:double);
begin
  x:=math.arccos(x);
end;

function log1plus(x:double):double;
var
  x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20:double;
begin
  //result:=lnxp1(x);
  x2:=x*x;
  x3:=x2*x;
  x4:=x2*x2;
  x5:=x4*x;
  x6:=x4*x2;
  x7:=x4*x3;
  x8:=x4*x4;
  x9:=x8*x;
  x10:=x8*x2;
  x11:=x8*x3;
  x12:=x8*x4;
  x13:=x8*x5;
  x14:=x8*x6;
  x15:=x8*x7;
  x16:=x8*x8;
  x17:=x16*x;
  x18:=x16*x2;
  x19:=x16*x3;
  x20:=x16*x4;
  {
  x21:=x16*x5;
  x22:=x16*x6;
  x23:=x16*x7;
  x24:=x16*x8;
  x25:=x16*x9;
  x26:=x16*x10;
  x27:=x16*x11;
  x28:=x16*x12;
  x29:=x16*x13;
  x30:=x16*x14;
  }
  result:= {-x30/30 +x29/29 -x28/28 +x27/27 -x26/26 +x25/25 -x24/24 +x23/23 -x22/22 +x21/21}
           -x20/20 +x19/19 -x18/18 +x17/17 -x16/16 +x15/15 -x14/14 +x13/13 -x12/12 +x11/11
           -x10/10 +x9/9 -x8/8 +x7/7 -x6/6 +x5/5 -x4/4 +x3/3 -x2/2 +x;
end;

function NPXpower1plus(x,y:double):double;
begin
  // result:=exp(y*lnxp1(x))
  result:=exp(y* log1plus(x))
end;

{$IFDEF FPC_HAS_TYPE_EXTENDED}
function NPXpower(x,y:double):double;
begin
   result:=exp(y*ln(x))
end;
{$ELSE}
 function NPXpower(x,y:double):double;  inline;
 begin
     result:=math.power(x,y)
 end;
{$ENDIF}

procedure  power(var x,y:double);
var
   t:double;
begin
   if x>0 then
      begin
         t:=x-1;
         if abs(t)<0.125 then
            x:=NPXpower1plus(t,y)
         else
            x:=NPXPower(x,y)
      end
   else if x=0 then
      if y>0 then
         x:=0
      else if y=0 then
         x:=1
      else
         setexception(3003)     //2022.2.11 修正
   else
      begin
         if int(y)=y then
            begin
               x:=-x;
               power(x,y);
               t:=y/2;
               if int(t)<>t then
                   x:=-x;
            end
         else
            setexception(3002)  //2022.2.11 修正
      end;
end;
function FStr(x:double):ansistring;
var
     s          :string[21];
     sign,sign1 :string[1];
     exrad      :string[6];
     i,e        :integer;
const places=18;
    function pureInt(x :double):double;
    var
       i:double;
    begin
       i:=int(x);
       if x>=0 then
              pureInt:=i
       else
              if i=x then
                 pureint:=i
              else
                 pureint:=i-1
    end;
begin
 if x<>0 then
   begin
    e:=LongintRound(pureint(system.ln(abs(x)) / system.ln(10)))  ;
    if (-5<=e) and (e<places) then
        begin
            if e>=-2 then str(x:1:17,s)
                     else str(x:1:16,s);
            i:=length(s);
            while s[i]='0' do dec(i);
            if s[i]='.' then dec(i);
            s:=copy(s,1,i);
            if s[1]='-' then
                begin
                   s:=copy(s,2,19);
                   sign:='-'
                end
            else
                sign:='';
            if s[1]='0' then  s:=copy(s,2,19);
            s:=sign+s
        end
    else
        begin
           if (e>=-999) and (e<=999) then
              str(x:20,s)
           else
              str(x:19,s);
           i:=pos('E',s);
           sign1:=copy(s,i+1,1);
           exrad:=copy(s,i+2,4);
           s:=copy(s,1,i-1);
            i:=length(s);
            while s[i]='0' do dec(i);
            s:=copy(s,1,i);
            i:=1;
            while s[i]=' ' do inc(i);
            s:=copy(s,i,19);
           if sign1='+' then sign1:='';
           i:=1;
           while exrad[i]='0' do inc(i);
           exrad:=copy(exrad,i,4);
           s:=s+'E'+sign1+exrad
        end;
   end
 else if (x=0) then
        s:='0';
 if s[1]<>'-' then s:=' '+s;
   FStr:=s
end;
(*
procedure invalidoperation;
var
   x:double;
begin
   x:=sqrt(-1)
end;
*)

initialization
  SetExceptionMask([exDenormalized, exZeroDivide, exOverflow,   exUnderflow,  exPrecision]);
  SetRoundMode(rmNearest);
end.
