
{$IFDEF FPC}
  {$MODE DELPHI}{$H+}
{$ENDIF}

{$DEFINE arith87}
{$S-}
{$X+}
Unit arithmet;

interface

{***********}
{type Number}
{***********}

type
   PNumber = ^Number;
   number  = object
              sig: {$IFDEF FPC_HAS_TYPE_EXTENDED}extended{$ELSE}double{$ENDIF};
              exp: smallint;
           procedure init(p:PNumber);
           procedure initzero;
           procedure initone;
           function sign:integer;
      end;


const
  {$IFDEF FPC_HAS_TYPE_EXTENDED}
   places=17;
   IntegralLimit=100000000000000000.;
   MaxExp=4932-Places+1;
   MinExp=-4932-Places;
  {$ELSE}
   places=15;
   IntegralLimit=1000000000000000.;
   MaxExp=308-Places+1;
   MinExp=-324-Places;
  {$ENDIF}
   MaxNumber:number =( sig: IntegralLimit - 1.0 ; exp: maxExp  );


   MaxExpnDecimal=99-places+1;
   MinExpnDecimal=-100-places+1;

  // MaxNum:Number=(sig: IntegralLimit; exp: MaxExpnDecimal);
  // NegMaxNum:number=(sig: -IntegralLimit; exp:MaxExpnDecimal);
  constMAXNUM:number =( sig:IntegralLimit/10.0 ; exp:maxExpnDecimal+1 );
  constNEgMAXNUM:number =( sig:-IntegralLimit/10.0 ; exp:maxExpnDecimal+1 );

type
  unaryoperation =procedure(var x:number);
  binaryoperation = procedure (var a,b:Number; var x:Number);


procedure subst(var p:PNumber; var n:number);
procedure disposenumber(var p:PNumber);

function checkRangeDecimal(var n:number; extyp:integer):boolean;

function IntegerVal(a:number; var c:integer):Integer;
                                   {c=0:nomal ; else c:sign}
function sgn(n:PNumber):integer;
function isZero(n:PNumber):boolean;
function isinteger(var n:number):boolean;
function compare(var a,b:number):integer;

procedure add(var a,b:Number; var x:Number);
procedure sbt(var a,b:Number; var x:Number);
procedure mlt(var a,b:Number; var x:Number);
procedure qtt(var a,b:Number; var x:Number);
procedure divide(var a,b:number; var x,y:number);

procedure qtt2(var n:number);
procedure oposite(var n:Number);
procedure oppose(var n:Number);
procedure absolute(var n:number);
procedure opposite(var n:Number);

procedure intpart(var n:number);
procedure fractpart(var n:number);

procedure intround(var n:number);    {round to integer}
procedure ceil(var n:number);
procedure BasicInt(var n:number);
procedure EpsDecimal(var n:number);
procedure EpsNative(var n:number);
procedure BasicMod(var a,b:Number; var x:number);
procedure Remainder(var  a,b:Number; var x:number);
//procedure eps(var n:number);
procedure min(var a,b:number; var n:number);
procedure max(var a,b:number; var n:number);

//function  min( a,b:number; var n:number):boolean;
//function  max( a,b:number; var n:number):boolean;

procedure round(var x,n:number; var y:number);
procedure Truncate(var x,n:number; var y:number);

//procedure round2digit(var n:number);
procedure round9( var n:number);
procedure sqrlong(var a:number);
procedure square(var n:number);

procedure initdecimal(var n:number; x:extended);
procedure initinteger(var n:number; b:integer);
procedure initlongint(var n:number; i:longint);

function LongintVal(x:extended; var c:integer):longint;
function LongintVal(const a:number; var c:integer):longint;
                                   {c=0:normal ; else c:sign}
type
    roundprocedure =procedure (var n:number);
var       RoundExpression:roundprocedure;
procedure roundprecision(var n:number);
procedure RoundVariable(var n:number);

procedure convert(a:extended; var n:number);
procedure convert16(a:extended; var n:number);
function doubleVal(a:Number):double;
function logN(var a:number):double;
procedure power(var a,b:Number; var x:Number);
//procedure intpower(var a,b:number; var n:number);

procedure  NumericRep(var x:extended;var code:integer;var line:string;var cp:integer);
procedure NVal(s:string; var x:extended);

function DStr(n:Number):string;
function DStrF(x:Double; signiwidth:integer):ansistring;
procedure ConvertToString(const n:number;var digits:ansistring;var exp1:integer; places0:integer);
procedure roundstring(var s:ansistring; n:integer; var exp:integer);
procedure setOpModeNative;
procedure setOpModeDecimal;
procedure setOpModeHigh;
procedure setOpModeRational;

var
   signiwidth:smallint=places;



var
  zero:PNumber;
  one:PNumber;
  ten:PNumber;
  half:PNumber;
  MAXNUM:PNumber;
  decimalPI:PNumber;


{************}
implementation
{************}
uses math, base,sconsts,float ;


//var
   //zero :number;
   //one  :number;
   //decimalPI:number ;

procedure number.init(p:PNumber);
begin
  self:=p^
end;

procedure number.initzero;
begin
  init(zero)
end;

procedure number.initone;
begin
  init(one)
end;

function number.sign:integer;
begin
  result:=sgn(@self)
end;


function roundext(x:double):double;
begin
   result:=system.Round(x)
end;


{***********}
{type Number}
{***********}

const
   UpperLimit:array[places..places+1]of extended=(integralLimit,integralLimit*10);

var
   DecimalMode:boolean=true;





function Normalize(var a:number):boolean;
var
   x:extended;
begin
  result:=true;
  with a do
     begin
        if sig=0. then
           exp:=minexp
        else if abs(sig)<integralLimit/10.0 then
          begin
             x:=sig*10.0;
             while (abs(x)<IntegralLimit) and (exp>minexp) do
                begin
                    sig:=x;
                    x:=sig*10.0 ;
                    dec(exp)
                end;
          end
        else if abs(sig)>=IntegralLimit then
             while (abs(sig)>=IntegralLimit) and (exp<=maxexp) do
               begin
                    sig:=sig/10.0 ;
                    inc(exp);
               end;

        if DecimalMode then
          begin
             sig:=roundext(sig/2)*2;             //最下位桁を偶数に丸める
             if abs(sig)=IntegralLimit then
                  begin
                      sig:=sig/10.0;
                      inc(exp);
                  end;
          end
        else
           //sig:=Int(sig);                    // 切り捨て
           begin
           sig:=roundext(sig);             //最下位桁を丸める
           if abs(sig)=IntegralLimit then
                begin
                    sig:=sig/10.0;
                    inc(exp);
                end;
           end;


         if exp>maxexp then
               begin
                   a:=maxnumber;
                   setexception(1002);
                   result:=false ;
               end;
     end;
end;


procedure  add(var a,b:Number; var x:Number);
var
   aa,bb:number;
begin
   if a.exp<b.exp then
      add(b,a,x)
   else
      begin
         aa.init(@a);
         bb.init(@b);
         if aa.exp<=bb.exp+(places+1) then
              begin
                 if aa.exp=bb.exp+1 then
                      begin
                          aa.sig:=aa.sig*10.;
                          dec(aa.exp)
                      end
                 else
                      while aa.exp>bb.exp do
                         begin
                            bb.sig:=bb.sig/10. ;
                            inc(bb.exp)
                         end;
                 aa.sig:=aa.sig+bb.sig;
              end;
          normalize(aa);
          x.init(@aa)
      end;
end;


procedure sbt(var a,b:Number; var x:Number);
var
   bb:Number;
begin
     bb.init(@b);
     bb.sig:=-bb.sig;
     add(a,bb,x);
end;

procedure mlt(var a,b:Number; var x:Number);
begin
     x.sig:=a.sig*b.sig/IntegralLimit;
     x.exp:=a.exp+b.exp+places;
     normalize(x);
end;


procedure qtt(var a,b:Number; var x:Number);
begin
     if b.sig=0. then
        begin
             a:=maxnumber;
             setexception(3001);
        end;
     x.sig:=a.sig*IntegralLimit/b.sig;
     x.exp:=a.exp-b.exp-places;
     normalize(x);
end;



function IntegerVal(a:number; var c:integer):Integer;
var
   x  :double;
begin
    x:=doubleval(a);
    c:=0;
    if x<MinInt  {-2147483648.0} then
       c:=-1
    else if x>MaxInt {2147483647.0} then
       c:=1
    else
       IntegerVal:=system.round(x) ;
end;






function isZero(n:PNumber):boolean;
begin
   iszero:=n^.sig=0.0
end;

function isinteger(var n:number):boolean;
var
   r:number;
begin
   basicmod(n,one^,r);
   result:=iszero(@r)
end;

{**********}
{arithmetic}
{**********}


function sgn(n:Pnumber):integer;
begin
   with n^ do
     if sig>0 then sgn:=1
     else if sig<0 then sgn:=-1
     else sgn:=0
end;


procedure oposite(var n:Number);
begin
    n.sig:=-n.sig
end;

procedure oppose(var n:Number);
begin
    n.sig:=-n.sig
end;

procedure qtt2(var n:number);
begin
  with n do sig:=sig/2.0;
end;

{**********}
{arithmetic}
{**********}



function compare(var a,b:number):integer; {do not change a,b}
var
   nonNegA,nonNegB:boolean;
   sign:integer;
begin
   nonNegA:=(a.sig>=0);
   nonNegB:=(b.sig>=0);
   if nonNegA and not nonNegB then
      compare:=1
   else if not NonNegA and NonNegB then
      compare:=-1
   else
      begin
         if nonNegA then sign:=1 else sign:=-1;
         if a.exp>b.exp then
               compare:=sign
         else if a.exp<b.exp then
               compare:=-sign
         else if a.sig>b.sig then
                 compare:=1
         else if a.sig<b.sig then
                 compare:=-1
         else
                 compare:=0
      end;
end;




procedure denormalize(var n:number);
begin
    n.sig:=doubleval(n);
    n.sig:=0;
end;



procedure intpart(var n:number);
begin
   if n.exp<0 then
       convert(int(doubleval(n)),n)
end;

procedure fractpart(var n:number);
var
   m:number;
begin
   m:=n;
   intpart(m);
   sbt(n,m,n);
end;

procedure  BasicInt(var n:number);
var
    m:number;
begin
    if n.sig>=0 then
       intpart(n)
    else
       begin
            m:=n;
            intpart(m);
            if compare(m,n)=0 then
               n:=m
            else
               sbt(m,one^,n)
       end
end;



procedure ceil(var n:number);
begin
    oposite(n);
    BasicInt(n);
    oposite(n);
end;

procedure BasicMod(var a,b:Number; var x:Number);
var
    aa,bb:number;
begin
     if isZero(@b) then
        begin
            setexception(3006);
        end
     else
        begin
           aa.init(@a);
           bb.init(@b);
           qtt(aa,bb,x);
           basicInt(x);
           mlt(x,bb,x);
           sbt(aa,x,x);
        end
end;

procedure Remainder(var a,b:Number; var x:Number);
var
    aa,bb:number;
begin
     if isZero(@b) then
        setexception(3006)
     else
        begin
           aa.init(@a);
           bb.init(@b);
           qtt(aa,bb,x);
           Intpart(x);
           mlt(x,bb,x);
           sbt(aa,x,x);
        end
end;



procedure divide(var a,b:number; var x,y:number);
var
  q,r:number;
begin
    qtt(a,b,q);
    basicint(q);
    basicmod(a,b,r);
    x.init(@q);
    y.init(@r);
    if not isinteger(q) then setexception(SystemErr)
end;

{*******}
{convert}
{*******}

var
   arrayPower10: array[0..15] of double;

procedure  power10init;
var
  i:integer;
begin
  arraypower10[0]:=1;
  for i:=1 to 15 do
      arraypower10[i]:=arraypower10[i-1]*10;
end;

function power10(i:integer):double;
var
   x,y:double;
   h:integer;
begin
    h:=i div 16;
    i:=i mod 16;
    if h=0 then
       power10:=arrayPower10[i]
    else
       begin
          x:=1E16;
          y:=1.;
          while h>0 do
              begin
                 if h mod 2 =1 then
                    y:=y*x;
                 h:=h div 2;
                 if h>0 then x:=x*x;
              end;
          power10:=y*arrayPower10[i]
       end;
end;



const IntegralLimit2=IntegralLimit*IntegralLimit;
const IntegralLimit20=IntegralLimit*IntegralLimit*10;


function normalize1(var a:number):boolean;
begin
    while (abs(a.sig)>IntegralLimit20) and (a.exp<=maxexp-places) do
         begin
              a.sig:=a.sig/IntegralLimit ;
              a.exp:=a.exp+places;
         end;
    if a.sig=0. then a.exp:=minexp;
    while (abs(a.sig)<=1.0) and (a.exp>=minexp+places) do
         begin
              a.sig:=a.sig*IntegralLimit ;
              a.exp:=a.exp-places;
         end;
    normalize1:=normalize(a)
end;

procedure NoRound(var n:number);
begin

end;

var
    RoundConv:roundprocedure=NoRound;
    RoundVari:roundprocedure=NoRound;

procedure RoundVariable(var n:number);
begin
    RoundVari(n);
    checkrangedecimal(n,1002);
end;

procedure RoundConvert(var n:number; extyp:integer);
begin
    Roundconv(n);
    if DecimalMode then
       checkrangedecimal(n,extyp);
end;

procedure convert(a:extended; var n:number);
         {convert double to Decimal}
begin
    n.sig:=a;
    n.exp:=0;
    normalize1(n);
    RoundConvert(n,1003);
end;

procedure convert1002(a:extended; var n:number);
begin
  n.sig:=a;
  n.exp:=0;
  normalize1(n);
  RoundConvert(n,1002);
end;

function normalize2(var a:number):boolean;forward;

procedure convert16(a:extended; var n:number);
         {convert double to 16 digit or more Decimal}
begin
    n.sig:=a;
    n.exp:=0;
    normalize2(n);
end;
procedure initdecimal(var n:number; x:extended);
begin
    convert(x,n);
end;

procedure initinteger(var n:number ;b:integer);
begin
    initDecimal(n,b)
end;

procedure initlongint(var n:number; i:longint);
begin
    initDecimal(n,i)
end;


function doubleVal(a:Number):double;
begin
    if a.sig=0.0 then
        doubleval:=0.0
    else if a.exp>=0 then
        doubleval:=a.sig*power10(a.exp)
    else
        begin
            while a.exp<-8 do
                  begin
                      a.exp:=a.exp+8;
                      a.sig:=a.sig/power10(8)
                  end;
            doubleval:=a.sig/power10(-a.exp)
        end;
end;

(*
function LongintVal(var a:number; var c:integer):longint;
                                   {c=0:normal ; else c:sign}
begin
  try
    result:=math.floor(ExtendedVal(a)+0.5);
    c:=0;
  except
    c:=sgn(@a);
  end;
end;
*)

function LongintVal(const a:number; var c:integer):longint;
                                   {c=0:normal ; else c:sign}
var
  x:extended;
begin
  x:=doubleVal(a);
  result:= LongintVal(x,c)
end;

function LongintVal(x:extended; var c:integer):longint;
                                   {c=0:normal ; else c:sign}
begin
  if x>maxint then begin c:=1; result:=maxint end
  else if x<minint then begin c:=-1; result:=minint end
  else begin c:=0; result:=system.round(x)  end
end;

{**********}
{arithmetic}
{**********}

procedure sqrlong(var a:number);
begin
  convert(sqrt(doubleval(a)),a)
end;

procedure square(var n:number);
begin
  convert(sqr(doubleval(n)),n)
end;

const ln2:double=0.693147180559945309417;


function log1plus(x:double):double;
var
  x2,x3,x4:double;
begin
  //result:=lnxp1(x);
  x2:=x*x;
  x3:=x2*x;
  x4:=x2*x2;
  result:=-x4/4+x3/3-x2/2+x;
end;

function nearly1(var n:number):boolean;
begin
 nearly1:=abs(doubleval(n)-1)<1E-6
end;


function logN(var a:number):double;
var
   x:number;
begin
  if sgn(@a)<=0 then
       begin
          setexceptionwith(s_InvalidArgInLOG,3004);
          logN:=0
       end
  else if nearly1(a) then
       begin
          sbt(a,one^,x);
          logN:=log1plus(doubleval(x))
       end
  else
           LogN:=ln(doubleval(a));
end;



procedure LongintPower(var a:Number; b:longint; var x:Number);     //2010.3.28
var
   y,xx:number;
label
   L1;
begin
    xx.initone;
    if b=0 then goto L1;
    y.init(@a);
    if b>0 then
       begin
          while b<>0 do
              begin
                 if b mod 2<>0 then mlt(xx,y,xx);
                 b:=b div 2;
                 if b<>0 then mlt(y,y,y);
              end;
       end
    else
       try
          while b<>0 do
              begin
                 if b mod 2<>0 then mlt(xx,y,xx);
                 b:=b div 2;
                 if b<>0 then mlt(y,y,y);
              end;
          qtt(one^,xx,xx);     {y:=1/y}
       except
         on E:EExtype do
           if extype=1002 then
              begin
                extype:=0;
                xx.initzero;
              end
           else if extype=3001 then
              setexception(1002)
           else
              raise E;
       end;
  L1:
    x.init(@xx);
end;


procedure RegularPower(var a,b:Number; var x:number);
var
    i:longint;
    c:integer;
    a1,m,n:number;
    d1,d2:double;
begin
  a1.init(@a);
  if sgn(@a)<=0 then
      begin
        setexception(3004);
        x.initzero
      end
  else if nearly1(a) then
      begin
        sbt(a,one^,a1);
        d1:=doubleval(a1);  d2:=doubleval(b);
        convert1002(NPXpower1plus(d1,d2),x)
        //convert1002(math.power(doubleval(a),doubleval(b)),x)
      end
  else //if compareabs(b,PNumber(@const1024)^)<0 then
    begin
      m.init(@b);
      intpart(m);
      sbt(b,m,n);
      i:=longintval(m,c);
      if c=0 then
         begin
           LongintPower(a,i,x);
           if iszero(@n) then
             begin
               RoundConv(x);
               checkrangedecimal(x,1002)
             end
           else
             convert1002(doubleVal(x)*NPXPower(doubleval(a),doubleval(n)),x);
         end
      else
         convert1002(NPXpower(doubleVal(a),doubleVal(b)),x);
    end;
end;


procedure power(var a,b:Number; var x:number);
var
   n,y:number;
begin
   if iszero(@b) then
       x.initone
   else if a.sign>0 then
      begin
         regularPower(a,b,y);
         x.init(@y)
      end
   else if a.sign=0 then
       begin
               if b.sign>0 then
                   x.initzero
               else if b.sign=0 then
                   x.initone
               else
                   setexception(3003)
       end
   else {if a<0 then}
       begin
            if isinteger(b) then
                begin
                      n.init(@a);
                      oppose(n);
                      power(n,b,y);
                      n.init(@b);
                      qtt2(n);
                      if not isinteger(n) then
                           oppose(y);
                      x.init(@y);
                end
            else
                      setexception(3002);
       end;
end;




(*
function mini(a,b:longint):longint;
begin
   if a<=b then  result:=a
           else  result:=b
end;
*)
(*
procedure IncremPowerLongint(var a:Number; b:longint; var x:Number);    //2010.3.28
var
   svlimit:LongInt;
   bb:longint;
   y,xx:number;
label
   L1;
begin
    svlimit:=limit;
    limit:=mini(limit + 2 ,maxplace-1);
    try
      xx.initone;
      if b=0 then goto L1;
      y.init(@a);
      if b>0 then
         begin
            while b<>0 do
                begin
                   if b mod 2<>0 then mlt(xx,y,xx);
                   b:=b div 2;
                   if b<>0 then mlt(y,y,y);
                end;
         end
      else
        if iszero(@a) then
           setexception(3003)
        else
           try
              while b<>0 do
                  begin
                     if b mod 2<>0 then mlt(xx,y,xx);
                     b:=b div 2;
                     if b<>0 then mlt(y,y,y);
                  end;
              qtt(one^,xx,xx);     {y:=1/y}
           except
           on E:EExtype do
             if extype=1002 then
                begin
                  extype:=0;
                  xx.initzero;
                end
             else if extype=3001 then
                setexception(1002)
             else
                raise E;
           end;
    L1:
      x.init(@xx);
    finally
      limit:=svlimit;
    end;
end;



procedure IncremPowerComp(var a:Number; b:int64; var x:Number);     //2010.3.28
var
   svlimit:LongInt;
   y,xx:number;
   c:int64;
begin
    svlimit:=limit;
    limit:=mini(limit + 2 ,maxplace-1);
    try
        xx.initone;
        y.init(@a);
        if b>0 then
           begin
              while b<>0 do
                  begin
                     c:=b div 2;
                     if b-2*c<>0 then mlt(xx,y,xx);
                     b:=c;
                     if b<>0 then mlt(y,y,y);
                  end;
           end
        else if b<0 then
           try
              while b<>0 do
                  begin
                     c:=b div 2;
                     if b-2*c<>0 then mlt(xx,y,xx);
                     b:=c;
                     if b<>0 then mlt(y,y,y);
                  end;
               qtt(one^,xx,xx);     {y:=1/y}
           except
             on E:EExtype do
               if extype=1002 then
                  begin
                    extype:=0;
                    xx.initzero;
                  end
               else if extype=3001 then
                  setexception(1002)
               else
                  raise E;
           end;
        x.init(@xx);
    finally
        limit:=svlimit;
    end;
end;




procedure intpower(var a,b:number; var n:number);
var
   i:longint;
   c:integer;
   ii:int64;
begin


   if isinteger(b) then
      begin
        i:=LongintVal(b,c);
        if c=0 then
           begin
             IncremPowerLongint(a,i,n);
             exit
           end
        else
           begin
               c:=0;
               try
                  ii:=system.round(doubleVal(b));
               except
                  c:=1
               end;
               if c=0 then
                  begin
                    IncremPowerComp(a,ii,n);
                    exit
                  end;
           end;
      end;

   if UseTranscendentalFunction then
         power(a,b,n)
   else
         setexceptionwith(s_PowerIndex,1000);       //2010.3.28
end;

procedure incrempower(var a:Number; b:LongInt; var x:Number);
                                  {assume b>=0}
var
   n:number;
   z:number;
begin
    z.initone;
    n.init(@a);
    while b<>0 do
         begin
            if b mod 2 <>0 then mlt(z,n,z);
            b:=b div 2;
            if b<>0 then mlt(n,n,n);
         end;
   x.init(@z);
end;
*)


{*************}
{ EPS function}
{*************}
const
   lowerbound=integralLimit/10.0;

procedure eps(var n:number);
begin
 with n do
   if iszero(@n) then
      begin
          sig:=lowerbound;
          exp:=minexp;
      end
   else
       begin
          sig:=lowerbound;
          exp:=exp-places+1;
          if exp<minexp then exp:=minexp;
      end;
end;

procedure EpsDecimal(var n:number);
begin
 with n do
   if iszero(@n) then
      begin
          sig:=lowerbound;
          exp:=MinExpnDecimal;
      end
    else
      begin
          sig:=lowerbound;
          exp:=exp-places+3;
          if exp<MinExpnDecimal then exp:=MinExpnDecimal;
      end;
end;

procedure EpsNative(var n:number);
begin
    eps(n)
end;



procedure  min(var a,b:number; var n:number);
begin
     n:=zero^; {error value}
     if compare(a,b)<=0 then n:=a else n:=b
end;

procedure max(var a,b:number; var n:number);
begin
     n:=zero^;
     if compare(a,b)>=0 then n:=a else n:=b
end;

{
procedure  tenfold(var x:number ; n:integer);
var
   i:integer;
begin
   if n>0 then
      for i:=1 to n do mlt(x,ten,x)
   else if n<0 then
      for i:=-1 downto n do qtt(x,ten,x)
end;
}


procedure intround(var n:number);
begin
   if n.exp>=0 then exit;
   if n.sig>=0 then                             //2021.06.18
      convert(int(doubleval(n)+0.5),n)
   else
      convert(int(doubleval(n)-0.5),n);
end;


procedure round(var x,n:number; var y:number);
var
   e:integer;
   c:integer;
begin
   e:=IntegerVal(n,c);
   x.exp:=x.exp+e;
   intround(x);
   x.exp:=x.exp-e;
   normalize(x);
   y:=x;
end;


procedure truncate(var x,n:number; var y:number);
var
   e:integer;
   c:integer;
begin
   e:=IntegerVal(n,c);
   x.exp:=x.exp+e;
   intpart(x);
   x.exp:=x.exp-e;
   y:=x;
end;



{**********************}
{numeric representation}
{**********************}

procedure roundprecision(var n:number);
begin
  //要修正
end;

{**********************}
{numeric representation}
{**********************}


procedure  NumericRep(var x:extended; var code:integer;var line:string;var cp:integer);

var
      cpintpart,cpfractpart,cpexrad:integer;
      intpartlen,fractpartlen,exradlen:integer;
      scaledrep:boolean;
      s:string{string31};

   function isDigit:boolean;
   begin
       case line[cp] of
           '0'..'9':
              isDigit:=true
            else
              isDigit:=false
       end
   end;

begin
      while (line[cp]=' ') do inc(cp); {spacecut}

      {intPart}
      cpintpart:=cp;
      while isDigit do inc(cp);
      intpartlen:=cp-cpintpart;

      {fractpart}
      if line[cp]='.' then inc(cp);
      cpfractpart:=cp;
      while isDigit do inc(cp);
      fractpartlen:=cp-cpfractpart;

      {exrad}
      if (line[cp]='E') or (line[cp]='e') then
         begin
            cpExrad:=cp;
            inc(cp);
            if (line[cp]='+')  or (line[cp]='-') then inc(cp);
            while isDigit do inc(cp);
            exradlen:=cp-cpexrad;
            scaledrep:=true;
          end
      else
          begin
             cpExrad:=cp;
             exradlen:=0;
             scaledrep:=false;
          end;



          {give a value}
          s:='0'+ copy(line,cpintpart,  intpartlen)
                +'.'
                +copy(line,cpfractpart,fractpartlen)+'0'
                +copy(line,cpexrad,    exradlen);
          val(s,x,code);

         if (intpartlen=0) and (fractpartlen=0)
             or (code<>0) then
               begin
                   code:=1001;
                   x:=base.maxnumber;
               end
         else
                begin
                    extype:=0;
                end;

end;

procedure NVal(s:string; var x:extended);
var
   c,cp:integer;
   m:boolean;
begin
  cp:=1;
  while (cp<=length(s)) and (s[cp]=' ') do inc(cp);
  m:=false;
  if (cp<=length(s)) then
    begin
       if s[cp]='+' then
          inc(cp)
       else if (s[cp]='-') then
       begin
          m:=true;
          inc(cp)
       end;
    end;
  NumericRep(x,c,s,cp) ;
  if c<>0 then setexception(c);
  if m then x:=-x;
  while (cp<=length(s)) and (s[cp]=' ') do inc(cp);
  if cp<=length(s) then setexception(4001);
end;



function ConvertToStringSub(x:double; places0:integer):ansistring;
var
   i:integer;
   s:ansistring;
   p,q,r:double;
begin
    s:='';
    p:=x;
    for i:=places0 downto 1 do
        begin
            q:=INT(p/10);
            r:=p-q*10;
            s:=char(byte('0')+system.round(r))+s;
            p:=q;
        end;
    result:=s;
end;

procedure ConvertToString(const n:number; var digits:ansistring; var exp1:integer; places0:integer);
var
   i:integer;
begin
  //{$IFDEF FPC_HAS_TYPE_EXTENDED}
  // str(abs(n.sig):places:0,digits);
  //{$ELSE}
   digits:=ConvertToStringSub(abs(n.sig),places0);
  //{$ENDIF}
   exp1:=n.exp +places0 ;

   i:=1;
   while (i<=length(digits)) and (digits[i]='0') do
                                            begin inc(i); dec(exp1) end;
   delete(digits,1,i-1) {digits:=copy(digits,i,255)} ;
   while (length(digits)>0) and (digits[length(digits)]='0') do
                                       setlength(digits,length(digits)-1);
end;





procedure roundstring(var s:ansistring; n:integer; var exp:integer);
var
    carry:boolean;
    t:char;
begin

    if n<0 then begin s:='';  exit end;

    if length(s)>n then
       begin
           t:=s[n+1];
           setlength(s,n);
           if t>='5' then
              begin
                 carry:=true;
                 while carry and (length(s)>0) do
                 begin
                     s[length(s)]:=succ(s[length(s)]);
                     if s[length(s)]<='9' then
                        carry:=false
                     else
                        setlength(s,length(s)-1);
                 end;
                 if length(s)=0 then
                    begin
                        s:='1';
                        inc(exp)
                    end;
              end;
       end;

    while (length(s)>0) and (s[length(s)]='0') do
                                       setlength(s,length(s)-1);

end;


function DStrSub(n:Number; signiwidth:integer):ansistring;
var
     sign    :string[1];
     digits  :ansistring;
     exp     :integer;
     exrad   :ansistring;
     e       :integer;
begin
   {roundvariable(n);}
   if sgn(@n)=0 then begin DstrSub:=' 0' ; exit end;
   if sgn(@n)>0 then sign:=' ' else sign:='-';
   if signiwidth>places then
       ConvertToString(n,digits,exp,places+1)
   else
       ConvertToString(n,digits,exp,places);
   roundstring(digits,signiwidth,exp);

   if (exp>0) and (exp<=signiwidth) then
      begin
         if exp>=length(digits) then        {wothout fraction part}
            begin
               while length(digits)<exp do digits:=digits + '0';
               DstrSub:=sign + digits;
            end
         else
             begin
                DstrSub:=sign + copy(digits,1,exp) + '.' +copy(digits,exp+1,maxint{255});
             end
      end
   else if (exp<=0) and ((length(digits)-exp<=signiwidth)) then
      begin
          e:=exp;
          while e<0 do
                    begin
                        digits:='0' + digits;
                        inc(e);
                    end;
          DstrSub:=sign + '.' + digits;
      end
   else
      begin
          str(exp-1,exrad);
          if exp-1>0 then exrad:='+'+exrad;    //2021.12.30
          DstrSub:=sign + copy(digits,1,1) + '.' + copy(digits,2,maxint{255})
                                        + 'E' +exrad;
      end;
end;

function DStr(n:Number):ansistring;
begin
 DStr:=DStrSub(n, signiwidth)
end;

{$IFDEF FPC_HAS_EXTENDED}
function normalize2(var a:number):boolean; inline;
begin
   result:=Normalize(a);
end;
{$ELSE}
//const IntegralLimit200=IntegralLimit*IntegralLimit*100;


const  IntegralLimit10=IntegralLimit*10;

function NormalizeLarge(var a:number ):boolean;
var
     x:extended;
begin
    result:=true;
    with a do
       begin
          if sig=0. then
             exp:=minexp
          else if abs(sig)<integralLimit10/10.0 then
            begin
               x:=sig*10.0;
               while (abs(x)<IntegralLimit10) and (exp>minexp) do
                  begin
                      sig:=x;
                      x:=sig*10.0 ;
                      dec(exp)
                  end;
            end
          else if abs(sig)>=IntegralLimit10 then
               while (abs(sig)>=IntegralLimit10) and (exp<=maxexp) do
                 begin
                      sig:=sig/10.0 ;
                      inc(exp);
                 end;

          if DecimalMode then
            begin
               sig:=roundext(sig/2)*2;             //最下位桁を偶数に丸める
               if abs(sig)=IntegralLimit10 then
                    begin
                        sig:=sig/10.0;
                        inc(exp);
                    end;
            end
          else
             //sig:=Int(sig);                    // 切り捨て
             begin
             sig:=roundext(sig);             //最下位桁を丸める
             if abs(sig)=IntegralLimit10 then
                  begin
                      sig:=sig/10.0;
                      inc(exp);
                  end;
             end;

           if exp>maxexp then
                 begin
                     a:=maxnumber;
                     setexception(1002);
                     result:=false ;
                 end;
       end;
end;

function normalize2(var a:number):boolean;
const IntegralLimit10_2=IntegralLimit10*IntegralLimit10*10;
begin
    while (abs(a.sig)>IntegralLimit10_2) and (a.exp<=maxexp-places-1) do
         begin
              a.sig:=a.sig/(IntegralLimit10) ;
              a.exp:=a.exp+(places+1);
         end;
    if a.sig=0. then a.exp:=minexp;
    while (abs(a.sig)<=1.0) and (a.exp>=minexp+(places+1)) do
         begin
              a.sig:=a.sig*(IntegralLimit10) ;
              a.exp:=a.exp-(places+1);
         end;
    normalize2:=normalizeLarge(a)
end;
{$ENDIF}

function DStrF(x:Double; signiwidth:integer):ansistring;
var
     n:Number;
begin
     n.sig:=x;
     n.exp:=0;
    if signiwidth>places then
       Normalize2(n)
    else
       Normalize1(n);
    DStrF:=DStrSub(n,signiwidth);
end;


{round decimal}

procedure round9( var n:number);
begin
  with n do
      sig:=RoundExt(sig/100000000)*100000000;
  Normalize(n)
end;


procedure round2digit( var n:number);
begin
  with n do
      sig:=RoundExt(sig/100)*100;
  Normalize(n)
end;



procedure disposenumber(var p:PNumber);
begin
    if p<>nil then dispose(p);
    p:=nil
end;


procedure subst(var p:PNumber; var n:number);

begin
   if p=nil then
        New(p);
   p^.init(@n)

end;


function checkRangeDecimal(var n:number; extyp:integer):boolean;
begin

    //round2digit(n);
    with n do
         if (exp>maxExpnDecimal) and
            ((compare(n,constMaxNum)>0) or (compare(n,constNegMaxNum)<0)) then
                 setexception(extyp)
          else if (exp<minExpnDecimal) then
                 n:=zero^;

      checkRangeDecimal:=(extype=0);

end;



procedure absolute(var n:Number);
begin
  with n do sig:=abs(sig)
end;

procedure opposite(var n:Number);
begin
  oppose(n)
end;

procedure setOpModeNative;
begin
  DecimalMode:=false;
  RoundExpression:=NoRound;
  RoundVari:=NoRound;
    if signiwidthMore then
        signiwidth:=17
  else
       signiwidth:=15;
 end;

procedure setOpModeDecimal;
begin
  DecimalMode:=true;
  RoundExpression:=NoRound;
  RoundVari:=round2digit;

  if signiwidthMore then
        signiwidth:=places
  else
       signiwidth:=places-2;
end;

procedure setOpModeHigh;
begin
  DecimalMode:=true;
end;

procedure setOpModeRational;
begin
  DecimalMode:=false;
end;

var
   constzero:Number=(sig:0;exp:minexp);
   constone: Number;
   constten: Number;
   consthalf:Number;
   constdecimalPI:number;


procedure initconsts;
begin
   convert(1.0, constone);
   convert(10.0,constten);
   convert(0.5, consthalf);
   convert(pi,constdecimalPI);

   pointer(zero):=@constzero;
   pointer(one):=@constone;
   pointer(ten):=@constten;
   pointer(half):=@consthalf;
   pointer(MAXNUM):=@constMAXNUM;
   pointer(decimalPI):=@constdecimalPI;
end;


begin
     initconsts;
     power10init;
     SetOpModeDecimal;
end.





