unit IntervalArithmetic32and64;

// Delphi unit (version 5.3) for 32-bit and 64-bit Windows environments
// (C) Copyright 1998-2025 by Andrzej Marciniak
// Poznan University of Technology, Institute of Computing Science

// Note: Do not use this unit for any Delphi's compiler older than XE!

interface

{$IFDEF WIN64}
// Delphi's 64 bit compiler do not support 80-bit Extended floating point
// values on Win64 (Extended = Double on Win64).
// The uTExtendedX87 unit provides for Win64 a replacement FPU-backed 80-bit
// Extended floating point type called TExtendedX87. This unit is available from
// http://cc.embarcadero.com/Item/28488
// Be sure that one of the defines EnableHelperRoutines or
// EnableFWAITsEverywhere is define within this unit by the $DEFINE compiler
// directive (both these defines are given as comments in uTExtendedX87 unit
// - see lines 126 and 128 in uTExtendedX87)
uses uTExtendedX87;
type Extended = TExtendedX87;
{$ENDIF}

// Basic interval type with definitions of overloading operators for proper
// intervals
type interval = record
                  var a, b : Extended;
                  class operator Implicit (x : Extended) : interval;
                  class operator Negative (x : interval) : interval;
                  class operator Positive (x : interval) : interval;
                  class operator Add (x, y : interval) : interval;
                  class operator Subtract (x, y : interval) : interval;
                  class operator Multiply (x, y : interval) : interval;
                  class operator Divide (x, y : interval) : interval;
                end;

// Functions for basic arithmetic operations for proper intervals (one can use
// these functions instead of overloading operators in interval record type
// what significantly reduces computational time)
function int_width (const x : interval) : Extended;
function iadd (const x, y : interval) : interval;
function isub (const x, y : interval) : interval;
function imul (const x, y : interval) : interval;
function idiv (const x, y : interval) : interval;

// Basic interval type with definitions of overloading operators for directed
// (improper) intervals
// For a theory see
// http://www.cs.put.poznan.pl/amarciniak/KONF-referaty/DirectedArithmetic.pdf
type dinterval = record
                   var a, b : Extended;
                   class operator Implicit (x : Extended) : dinterval;
                   class operator Negative (x : dinterval) : dinterval;
                   class operator Positive (x : dinterval) : dinterval;
                   class operator Add (x, y : dinterval) : dinterval;
                   class operator Subtract (x, y : dinterval) : dinterval;
                   class operator Multiply (x, y : dinterval) : dinterval;
                   class operator Divide (x, y : dinterval) : dinterval;
                 end;

// Basic arithmetic operations for directed (improper) intervals (one can use
// these functions instead of overloading operators in dinterval record type
// what significantly reduces computational time)
function dint_width (const x : dinterval) : Extended;
function projection (const x : dinterval) : interval;
function opposite (const x : dinterval) : dinterval;
function inverse (const x : dinterval) : dinterval;
function diadd (const x, y : dinterval) : dinterval;
function disub (const x, y : dinterval) : dinterval;
function dimul (const x, y : dinterval) : dinterval;
function didiv (const x, y : dinterval) : dinterval;

// Data reading functions for proper intervals
function int_read (const sa : string) : interval;
function left_read (const sa : string) : Extended;
function right_read (const sa : string) : Extended;

// Data reading functions for directed (improper) intervals
function dleft_read (const sa : string) : Extended;
function dright_read (const sa : string) : Extended;

// A procedure for transforming ends of proper intervals into strings
procedure iends_to_strings (const x         : interval;
                            out left, right : string);

// Basic functions for proper intervals.
// Note: Before using the values returned by functions, check the value of st.
// interval absolute value
function iabs (const x : interval) : interval;
// interval sine
function isin (const x : interval;
               out st  : Integer) : interval;
// interval cosine
function icos (const x : interval;
               out st  : Integer) : interval;
// interval tangent (|x|<Pi/2)
function itan (const x : interval;
               out st  : Integer) : interval;
// interval cotangent (0<|x|<Pi)
function icot (const x : interval;
               out st  : Integer) : interval;
// interval arcsin (interval inverse of the sine function, |x|<1)
function iarcsin (const x : interval;
                  out st  : Integer) : interval;
// interval arccos (interval inverse of the cosine function, |x|<1)
function iarccos (const x : interval;
                  out st  : Integer) : interval;
// interval arctan (interval inverse of the tangent function)
function iarctan (const x : interval;
                  out st  : Integer) : interval;
// interval arccot (interval inverse of the cotangent function)
function iarccot (const x : interval;
                  out st  : Integer) : interval;
// interval hiperbolic sine
function isinh (const x : interval;
                out st  : Integer) : interval;
// interval hiperbolic cosine
function icosh (const x : interval;
                out st  : Integer) : interval;
// interval hiperbolic tangent (|x|<Pi/2)
function itanh (const x : interval;
                out st  : Integer) : interval;
// interval hiperbolic cotangent (0<|x|<Pi)
function icoth (const x : interval;
                out st  : Integer) : interval;
// interval inverse hiperbolic sine
function iarcsinh (const x : interval;
                   out st  : Integer) : interval;
// interval inverse hiperbolic cosine (x>=1)
function iarccosh (const x : interval;
                   out st  : Integer) : interval;
// interval inverse hiperbolic tangent (|x|<1)
function iarctanh (const x : interval;
                   out st  : Integer) : interval;
// interval inverse hiperbolic cotangent (|x|>1)
function iarccoth (const x : interval;
                   out st  : Integer) : interval;
// interval natural logarithm (|x|>0)
function iln (const x : interval;
              out st  : Integer) : interval;
// interval exponential function
function iexp (const x : interval;
               out st  : Integer) : interval;
// interval square function
function isqr (const x : interval;
               out st  : Integer) : interval;
// interval square root function (x>=0)
function isqrt (const x : interval;
                out st  : Integer) : interval;
// interval power function (a^x, a>0, a<>1)
function iax (const a, x : interval;
              out st     : Integer) : interval;
// interval (1+x)^m (|x|<=1, m>0)
function ionepxtom (const x, m : interval;
                    out st     : Integer) : interval;
// interval (1-x)^m (|x|<=1, m>0)
function ionemxtom (const x, m : interval;
                    out st     : Integer) : interval;
// interval (1+x)^(-m) (|x|<1, m>0)
function ionepxtomm (const x, m : interval;
                     out st     : Integer) : interval;
// interval (1-m)^(-m) (|x|<1, m>0)
function ionemxtomm (const x, m : interval;
                     out st     : Integer) : interval;

// Interval constants (in the form of proper intervals)
function isqrt2 : interval;
function isqrt3 : interval;
function isqrt5 : interval;
function isqrt6 : interval;
function isqrt7 : interval;
function isqrt8 : interval;
function isqrt10 : interval;
function ipi : interval;

implementation
  uses System.SysUtils, System.Math, Vcl.Dialogs;
  type char_tab        = array [1..80] of Char;
{$IFDEF WIN64}
       interval_double = record
                           a, b : Double
                         end;
{$ENDIF}
  const bit : array [0..7] of Byte = ($01, $02, $04, $08, $10, $20, $40, $80);
        ldi : array [0..63] of string =
{2^0}      ('1.000000000000000000000000000000000000000000000000000000000000000',
            '0.500000000000000000000000000000000000000000000000000000000000000',
            '0.250000000000000000000000000000000000000000000000000000000000000',
            '0.125000000000000000000000000000000000000000000000000000000000000',
            '0.062500000000000000000000000000000000000000000000000000000000000',
            '0.031250000000000000000000000000000000000000000000000000000000000',
            '0.015625000000000000000000000000000000000000000000000000000000000',
            '0.007812500000000000000000000000000000000000000000000000000000000',
            '0.003906250000000000000000000000000000000000000000000000000000000',
            '0.001953125000000000000000000000000000000000000000000000000000000',
{2^(-10)}   '0.000976562500000000000000000000000000000000000000000000000000000',
            '0.000488281250000000000000000000000000000000000000000000000000000',
            '0.000244140625000000000000000000000000000000000000000000000000000',
            '0.000122070312500000000000000000000000000000000000000000000000000',
            '0.000061035156250000000000000000000000000000000000000000000000000',
            '0.000030517578125000000000000000000000000000000000000000000000000',
            '0.000015258789062500000000000000000000000000000000000000000000000',
            '0.000007629394531250000000000000000000000000000000000000000000000',
            '0.000003814697265625000000000000000000000000000000000000000000000',
            '0.000001907348632812500000000000000000000000000000000000000000000',
{2^(-20)}   '0.000000953674316406250000000000000000000000000000000000000000000',
            '0.000000476837158203125000000000000000000000000000000000000000000',
            '0.000000238418579101562500000000000000000000000000000000000000000',
            '0.000000119209289550781250000000000000000000000000000000000000000',
            '0.000000059604644775390625000000000000000000000000000000000000000',
            '0.000000029802322387695312500000000000000000000000000000000000000',
            '0.000000014901161193847656250000000000000000000000000000000000000',
            '0.000000007450580596923828125000000000000000000000000000000000000',
            '0.000000003725290298461914062500000000000000000000000000000000000',
            '0.000000001862645149230957031250000000000000000000000000000000000',
{2^(-30)}   '0.000000000931322574615478515625000000000000000000000000000000000',
            '0.000000000465661287307739257812500000000000000000000000000000000',
            '0.000000000232830643653869628906250000000000000000000000000000000',
            '0.000000000116415321826934814453125000000000000000000000000000000',
            '0.000000000058207660913467407226562500000000000000000000000000000',
            '0.000000000029103830456733703613281250000000000000000000000000000',
            '0.000000000014551915228366851806640625000000000000000000000000000',
            '0.000000000007275957614183425903320312500000000000000000000000000',
            '0.000000000003637978807091712951660156250000000000000000000000000',
            '0.000000000001818989403545856475830078125000000000000000000000000',
{2^(-40)}   '0.000000000000909494701772928237915039062500000000000000000000000',
            '0.000000000000454747350886464118957519531250000000000000000000000',
            '0.000000000000227373675443232059478759765625000000000000000000000',
            '0.000000000000113686837721616029739379882812500000000000000000000',
            '0.000000000000056843418860808014869689941406250000000000000000000',
            '0.000000000000028421709430404007434844970703125000000000000000000',
            '0.000000000000014210854715202003717422485351562500000000000000000',
            '0.000000000000007105427357601001858711242675781250000000000000000',
            '0.000000000000003552713678800500929355621337890625000000000000000',
            '0.000000000000001776356839400250464677810668945312500000000000000',
{2^(-50)}   '0.000000000000000888178419700125232338905334472656250000000000000',
            '0.000000000000000444089209850062616169452667236328125000000000000',
            '0.000000000000000222044604925031308084726333618164062500000000000',
            '0.000000000000000111022302462515654042363166809082031250000000000',
            '0.000000000000000055511151231257827021181583404541015625000000000',
            '0.000000000000000027755575615628913510590791702270507812500000000',
            '0.000000000000000013877787807814456755295395851135253906250000000',
            '0.000000000000000006938893903907228377647697925567626953125000000',
            '0.000000000000000003469446951953614188823848962783813476562500000',
            '0.000000000000000001734723475976807094411924481391906738281250000',
{2^(-60)}   '0.000000000000000000867361737988403547205962240695953369140625000',
            '0.000000000000000000433680868994201773602981120347976684570312500',
            '0.000000000000000000216840434497100886801490560173988342285156250',
            '0.000000000000000000108420217248550443400745280086994171142578125');
{$IFDEF WIN64}
        ldi_double : array [0..52] of string =
{2^0}        ('1.0000000000000000000000000000000000000000000000000000',
              '0.5000000000000000000000000000000000000000000000000000',
              '0.2500000000000000000000000000000000000000000000000000',
              '0.1250000000000000000000000000000000000000000000000000',
              '0.0625000000000000000000000000000000000000000000000000',
              '0.0312500000000000000000000000000000000000000000000000',
              '0.0156250000000000000000000000000000000000000000000000',
              '0.0078125000000000000000000000000000000000000000000000',
              '0.0039062500000000000000000000000000000000000000000000',
              '0.0019531250000000000000000000000000000000000000000000',
{2^(-10)}     '0.0009765625000000000000000000000000000000000000000000',
              '0.0004882812500000000000000000000000000000000000000000',
              '0.0002441406250000000000000000000000000000000000000000',
              '0.0001220703125000000000000000000000000000000000000000',
              '0.0000610351562500000000000000000000000000000000000000',
              '0.0000305175781250000000000000000000000000000000000000',
              '0.0000152587890625000000000000000000000000000000000000',
              '0.0000076293945312500000000000000000000000000000000000',
              '0.0000038146972656250000000000000000000000000000000000',
              '0.0000019073486328125000000000000000000000000000000000',
{2^(-20)}     '0.0000009536743164062500000000000000000000000000000000',
              '0.0000004768371582031250000000000000000000000000000000',
              '0.0000002384185791015625000000000000000000000000000000',
              '0.0000001192092895507812500000000000000000000000000000',
              '0.0000000596046447753906250000000000000000000000000000',
              '0.0000000298023223876953125000000000000000000000000000',
              '0.0000000149011611938476562500000000000000000000000000',
              '0.0000000074505805969238281250000000000000000000000000',
              '0.0000000037252902984619140625000000000000000000000000',
              '0.0000000018626451492309570312500000000000000000000000',
{2^(-30)}     '0.0000000009313225746154785156250000000000000000000000',
              '0.0000000004656612873077392578125000000000000000000000',
              '0.0000000002328306436538696289062500000000000000000000',
              '0.0000000001164153218269348144531250000000000000000000',
              '0.0000000000582076609134674072265625000000000000000000',
              '0.0000000000291038304567337036132812500000000000000000',
              '0.0000000000145519152283668518066406250000000000000000',
              '0.0000000000072759576141834259033203125000000000000000',
              '0.0000000000036379788070917129516601562500000000000000',
              '0.0000000000018189894035458564758300781250000000000000',
{2^(-40)}     '0.0000000000009094947017729282379150390625000000000000',
              '0.0000000000004547473508864641189575195312500000000000',
              '0.0000000000002273736754432320594787597656250000000000',
              '0.0000000000001136868377216160297393798828125000000000',
              '0.0000000000000568434188608080148696899414062500000000',
              '0.0000000000000284217094304040074348449707031250000000',
              '0.0000000000000142108547152020037174224853515625000000',
              '0.0000000000000071054273576010018587112426757812500000',
              '0.0000000000000035527136788005009293556213378906250000',
              '0.0000000000000017763568394002504646778106689453125000',
{2^(-50)}     '0.0000000000000008881784197001252323389053344726562500',
              '0.0000000000000004440892098500626161694526672363281250',
              '0.0000000000000002220446049250313080847263336181640625');
{$ENDIF}

  class operator interval.Implicit (x : Extended) : interval;
  var s : string;
  begin
    Str (x:26, s);
    Result.a:=left_read(s);
    Result.b:=right_read(s)
  end {Implicit};

  class operator interval.Negative (x : interval) : interval;
  var z : interval;
  begin
    z:=x;
    Result.a:=-z.b;
    Result.b:=-z.a
  end {Negative};

  class operator interval.Positive (x : interval) : interval;
  begin
    Result.a:=x.a;
    Result.b:=x.b
  end {Positive};

  function int_width (const x : interval) : Extended;
  begin
    SetRoundMode (rmUp);
    Result:=x.b-x.a;
    SetRoundMode (rmNearest)
  end {int_width};

  function iadd (const x, y : interval) : interval;
  begin
    SetRoundMode (rmDown);
    Result.a:=x.a+y.a;
    SetRoundMode (rmUp);
    Result.b:=x.b+y.b;
    SetRoundMode (rmNearest)
  end {iadd};

  class operator interval.Add (x, y : interval) : interval;
  begin
    Result:=iadd(x, y)
  end {Add};

  function isub (const x, y : interval) : interval;
  begin
    SetRoundMode (rmDown);
    Result.a:=x.a-y.b;
    SetRoundMode (rmUp);
    Result.b:=x.b-y.a;
    SetRoundMode (rmNearest)
  end {isub};

  class operator interval.Subtract (x, y : interval) : interval;
  begin
    Result:=isub(x, y)
  end {Subtract};

  function imul (const x, y : interval) : interval;
  var x1y1, x1y2, x2y1 : Extended;
  begin
    SetRoundMode (rmDown);
    x1y1:=x.a*y.a;
    x1y2:=x.a*y.b;
    x2y1:=x.b*y.a;
    with Result do
      begin
        a:=x.b*y.b;
        if x2y1<a
          then a:=x2y1;
        if x1y2<a
          then a:=x1y2;
        if x1y1<a
          then a:=x1y1
      end;
    SetRoundMode (rmUp);
    x1y1:=x.a*y.a;
    x1y2:=x.a*y.b;
    x2y1:=x.b*y.a;
    with Result do
      begin
        b:=x.b*y.b;
        if x2y1>b
          then b:=x2y1;
        if x1y2>b
          then b:=x1y2;
        if x1y1>b
          then b:=x1y1
      end;
    SetRoundMode (rmNearest)
  end {imul};

  class operator interval.Multiply (x, y : interval) : interval;
  begin
    Result:=imul(x, y)
  end {Multiply};

  function idiv (const x, y : interval) : interval;
  var x1y1, x1y2, x2y1 : Extended;
  begin
    if (y.a<=0.0) and (y.b>=0.0)
      then raise EZeroDivide.Create ('Division by an interval containing 0.')
      else begin
             SetRoundMode (rmDown);
             x1y1:=x.a/y.a;
             x1y2:=x.a/y.b;
             x2y1:=x.b/y.a;
             with Result do
               begin
                 a:=x.b/y.b;
                 if x2y1<a
                   then a:=x2y1;
                 if x1y2<a
                   then a:=x1y2;
                 if x1y1<a
                   then a:=x1y1
               end;
             SetRoundMode (rmUp);
             x1y1:=x.a/y.a;
             x1y2:=x.a/y.b;
             x2y1:=x.b/y.a;
             with Result do
               begin
                 b:=x.b/y.b;
                 if x2y1>b
                   then b:=x2y1;
                 if x1y2>b
                   then b:=x1y2;
                 if x1y1>b
                   then b:=x1y1
               end
           end;
    SetRoundMode (rmNearest)
  end {idiv};

  class operator interval.Divide (x, y : interval) : interval;
  begin
    Result:=idiv(x, y)
  end {Divide};

  class operator dinterval.Implicit (x : Extended) : dinterval;
  var s : string;
  begin
    Str (x:26, s);
    Result.a:=left_read(s);
    Result.b:=right_read(s)
  end {dImplicit};

  class operator dinterval.Negative (x : dinterval) : dinterval;
  var z : dinterval;
  begin
    z:=x;
    Result.a:=-z.b;
    Result.b:=-z.a
  end {dNegative};

  class operator dinterval.Positive (x : dinterval) : dinterval;
  begin
    Result.a:=x.a;
    Result.b:=x.b
  end {dPositive};

  function dint_width (const x : dinterval) : Extended;
  var w1, w2 : Extended;
  begin
    SetRoundMode (rmUp);
    w1:=x.b-x.a;
    if w1<0
      then w1:=-w1;
    SetRoundMode (rmDown);
    w2:=x.b-x.a;
    if w2<0
      then w2:=-w2;
    if w1>w2
      then Result:=w1
      else Result:=w2;
    SetRoundMode (rmNearest)
  end {dint_width};

  function projection (const x : dinterval) : interval;
  var z : dinterval;
  begin
    if x.a>x.b
      then begin
             z:=x;
             Result.a:=z.b;
             Result.b:=z.a
           end
      else begin
             Result.a:=x.a;
             Result.b:=x.b
           end;
  end {projection};

  function opposite (const x : dinterval) : dinterval;
  begin
    Result.a:=-x.a;
    Result.b:=-x.b;
  end {opposite};

  function inverse (const x : dinterval) : dinterval;
  var z1, z2 : dinterval;
  begin
    SetRoundMode (rmDown);
    z1.a:=1/x.a;
    z2.b:=1/x.b;
    SetRoundMode (rmUp);
    z1.b:=1/x.b;
    z2.a:=1/x.a;
    if dint_width(z1)>=dint_width(z2)
      then Result:=z1
      else Result:=z2;
    SetRoundMode (rmNearest)
  end {inverse};

  function diadd (const x, y : dinterval) : dinterval;
  var z1, z2 : dinterval;
  begin
    SetRoundMode (rmDown);
    if (x.a<=x.b) and (y.a<=y.b)
      then begin
             Result.a:=x.a+y.a;
             SetRoundMode (rmUp);
             Result.b:=x.b+y.b
           end
      else begin
             z1.a:=x.a+y.a;
             z2.b:=x.b+y.b;
             SetRoundMode (rmUp);
             z1.b:=x.b+y.b;
             z2.a:=x.a+y.a;
             if dint_width(z1)>=dint_width(z2)
               then Result:=z1
               else Result:=z2
           end;
    SetRoundMode (rmNearest)
  end {diadd};

  class operator dinterval.Add (x, y : dinterval) : dinterval;
  begin
    Result:=diadd(x, y)
  end {dAdd};

  function disub (const x, y : dinterval) : dinterval;
  var z1, z2 : dinterval;
  begin
    SetRoundMode (rmDown);
    if (x.a<=x.b) and (y.a<=y.b)
      then begin
             Result.a:=x.a-y.b;
             SetRoundMode (rmUp);
             Result.b:=x.b-y.a
           end
      else begin
             z1.a:=x.a-y.b;
             z2.b:=x.b-y.a;
             SetRoundMode (rmUp);
             z1.b:=x.b-y.a;
             z2.a:=x.a-y.b;
             if dint_width(z1)>=dint_width(z2)
               then Result:=z1
               else Result:=z2
           end;
    SetRoundMode (rmNearest)
  end {disub};

  class operator dinterval.Subtract (x, y : dinterval) : dinterval;
  begin
    Result:=disub(x, y)
  end {dSubtract};

  function dimul (const x, y : dinterval) : dinterval;
  var z1, z2               : dinterval;
  var x1y1, x1y2, x2y1, z  : Extended;
      xn, xp, yn, yp, zero : Boolean;
  begin
    SetRoundMode (rmDown);
    if (x.a<=x.b) and (y.a<=y.b)
      then begin
             x1y1:=x.a*y.a;
             x1y2:=x.a*y.b;
             x2y1:=x.b*y.a;
             with Result do
               begin
                 a:=x.b*y.b;
                 if x2y1<a
                   then a:=x2y1;
                 if x1y2<a
                   then a:=x1y2;
                 if x1y1<a
                   then a:=x1y1
               end;
             SetRoundMode (rmUp);
             x1y1:=x.a*y.a;
             x1y2:=x.a*y.b;
             x2y1:=x.b*y.a;
             with Result do
               begin
                 b:=x.b*y.b;
                 if x2y1>b
                   then b:=x2y1;
                 if x1y2>b
                   then b:=x1y2;
                 if x1y1>b
                   then b:=x1y1
               end
           end
      else begin
             xn:=(x.a<0) and (x.b<0);
             xp:=(x.a>0) and (x.b>0);
             yn:=(y.a<0) and (y.b<0);
             yp:=(y.a>0) and (y.b>0);
             zero:=False;
// A, B in H-T
             if (xn or xp) and (yn or yp)
               then if xp and yp
                      then begin
                             z1.a:=x.a*y.a;
                             z2.b:=x.b*y.b;
                             SetRoundMode (rmUp);
                             z1.b:=x.b*y.b;
                             z2.a:=x.a*y.a
                           end
                      else if xp and yn
                             then begin
                                    z1.a:=x.b*y.a;
                                    z2.b:=x.a*y.b;
                                    SetRoundMode (rmUp);
                                    z1.b:=x.a*y.b;
                                    z2.a:=x.b*y.a
                                  end
                             else if xn and yp
                                    then begin
                                           z1.a:=x.a*y.b;
                                           z2.b:=x.b*y.a;
                                           SetRoundMode (rmUp);
                                           z1.b:=x.b*y.a;
                                           z2.a:=x.a*y.b
                                         end
                                    else begin
                                           z1.a:=x.b*y.b;
                                           z2.b:=x.a*y.a;
                                           SetRoundMode (rmUp);
                                           z1.b:=x.a*y.a;
                                           z2.a:=x.b*y.b
                                         end
// A in H-T, B in T
               else if (xn or xp)
                        and ((y.a<=0) and (y.b>=0) or (y.a>=0) and (y.b<=0))
                      then if xp and (y.a<=y.b)
                             then begin
                                    z1.a:=x.b*y.a;
                                    z2.b:=x.b*y.b;
                                    SetRoundMode (rmUp);
                                    z1.b:=x.b*y.b;
                                    z2.a:=x.b*y.a
                                  end
                             else if xp and (y.a>y.b)
                                    then begin
                                           z1.a:=x.a*y.a;
                                           z2.b:=x.a*y.b;
                                           SetRoundMode (rmUp);
                                           z1.b:=x.a*y.b;
                                           z2.a:=x.a*y.a
                                         end
                                    else if xn and (y.a<=y.b)
                                           then begin
                                                  z1.a:=x.a*y.b;
                                                  z2.b:=x.a*y.a;
                                                  SetRoundMode (rmUp);
                                                  z1.b:=x.a*y.a;
                                                  z2.a:=x.a*y.b
                                                end
                                           else begin
                                                  z1.a:=x.b*y.b;
                                                  z2.b:=x.b*y.a;
                                                  SetRoundMode (rmUp);
                                                  z1.b:=x.b*y.a;
                                                  z2.a:=x.b*y.b
                                                end
// A in T, B in H-T
                      else if ((x.a<=0) and (x.b>=0) or (x.a>=0) and (x.b<=0))
                               and (yn or yp)
                             then if (x.a<=x.b) and yp
                                    then begin
                                           z1.a:=x.a*y.b;
                                           z2.b:=x.b*y.b;
                                           SetRoundMode (rmUp);
                                           z1.b:=x.b*y.b;
                                           z2.a:=x.a*y.b
                                         end
                                    else if (x.a<=0) and yn
                                           then begin
                                                  z1.a:=x.b*y.a;
                                                  z2.b:=x.a*y.a;
                                                  SetRoundMode (rmUp);
                                                  z1.b:=x.a*y.a;
                                                  z2.a:=x.b*y.a
                                                end
                                           else if (x.a>x.b) and yp
                                                  then begin
                                                         z1.a:=x.a*y.a;
                                                         z2.b:=x.b*y.a;
                                                         SetRoundMode (rmUp);
                                                         z1.b:=x.b*y.a;
                                                         z2.a:=x.a*y.a
                                                       end
                                                  else begin
                                                         z1.a:=x.b*y.b;
                                                         z2.b:=x.a*y.b;
                                                         SetRoundMode (rmUp);
                                                         z1.b:=x.a*y.b;
                                                         z2.a:=x.b*y.b
                                                       end
// A, B in Z-
                             else if (x.a>=0) and (x.b<=0) and (y.a>=0)
                                      and (y.b<=0)
                                   then begin
                                          z1.a:=x.a*y.a;
                                          z:=x.b*y.b;
                                          if z1.a<z
                                            then z1.a:=z;
                                          z2.b:=x.a*y.b;
                                          z:=x.b*y.a;
                                          if z<z2.b
                                            then z2.b:=z;
                                          SetRoundMode (rmUp);
                                          z1.b:=x.a*y.b;
                                          z:=x.b*y.a;
                                          if z<z1.b
                                            then z1.b:=z;
                                          z2.a:=x.a*y.a;
                                          z:=x.b*y.b;
                                          if z2.a<z
                                            then z2.a:=z
                                        end
// A in Z and B in Z- or A in Z- and B in Z
                                   else zero:=True;
             if zero
               then begin
                      Result.a:=0;
                      Result.b:=0
                    end
               else if dint_width(z1)>=dint_width(z2)
                      then Result:=z1
                      else Result:=z2
           end;
    SetRoundMode (rmNearest)
  end {dimul};

  class operator dinterval.Multiply (x, y : dinterval) : dinterval;
  begin
    Result:=dimul(x, y)
  end {dMultiply};

  function didiv (const x, y : dinterval) : dinterval;
  var x1y1, x1y2, x2y1     : Extended;
      z1, z2               : dinterval;
      xn, xp, yn, yp, zero : Boolean;
  begin
    SetRoundMode (rmDown);
    if (x.a<=x.b) and (y.a<=y.b)
      then begin
             if (y.a<=0.0) and (y.b>=0.0)
               then begin
                      SetRoundMode (rmNearest);
                      raise EZeroDivide.Create ('Division by an interval '
                                                +'containing 0.')
                    end
               else begin
                      x1y1:=x.a/y.a;
                      x1y2:=x.a/y.b;
                      x2y1:=x.b/y.a;
                      with Result do
                        begin
                          a:=x.b/y.b;
                          if x2y1<a
                            then a:=x2y1;
                          if x1y2<a
                            then a:=x1y2;
                          if x1y1<a
                            then a:=x1y1
                        end;
                      SetRoundMode (rmUp);
                      x1y1:=x.a/y.a;
                      x1y2:=x.a/y.b;
                      x2y1:=x.b/y.a;
                      with Result do
                        begin
                          b:=x.b/y.b;
                          if x2y1>b
                            then b:=x2y1;
                          if x1y2>b
                            then b:=x1y2;
                          if x1y1>b
                            then b:=x1y1
                        end
                    end
           end
      else begin
             xn:=(x.a<0) and (x.b<0);
             xp:=(x.a>0) and (x.b>0);
             yn:=(y.a<0) and (y.b<0);
             yp:=(y.a>0) and (y.b>0);
             zero:=False;
// A, B in H-T
             if (xn or xp) and (yn or yp)
               then if xp and yp
                      then begin
                             z1.a:=x.a/y.b;
                             z2.b:=x.b/y.a;
                             SetRoundMode (rmUp);
                             z1.b:=x.b/y.a;
                             z2.a:=x.a/y.b
                           end
                      else if xp and yn
                             then begin
                                    z1.a:=x.b/y.b;
                                    z2.b:=x.a/y.a;
                                    SetRoundMode (rmUp);
                                    z1.b:=x.a/y.a;
                                    z2.a:=x.b/y.b
                                  end
                           else if xn and yp
                                  then begin
                                         z1.a:=x.a/y.a;
                                         z2.b:=x.b/y.b;
                                         SetRoundMode (rmUp);
                                         z1.b:=x.b/y.b;
                                         z2.a:=x.a/y.a
                                       end
                                  else begin
                                         z1.a:=x.b/y.a;
                                         z2.b:=x.a/y.b;
                                         SetRoundMode (rmUp);
                                         z1.b:=x.a/y.b;
                                         z2.a:=x.b/y.a
                                       end
// A in T, B in H-T
               else if (x.a<=0) and (x.b>=0) or (x.a>=0) and (x.b<=0)
                        and (yn or yp)
                      then if (x.a<=x.b) and yp
                             then begin
                                    z1.a:=x.a/y.a;
                                    z2.b:=x.b/y.a;
                                    SetRoundMode (rmUp);
                                    z1.b:=x.b/y.a;
                                    z2.a:=x.a/y.a
                                  end
                             else if (x.a<=x.b) and yn
                                    then begin
                                           z1.a:=x.b/y.b;
                                           z2.b:=x.a/y.b;
                                           SetRoundMode (rmUp);
                                           z1.b:=x.a/y.b;
                                           z2.a:=x.b/y.b
                                         end
                                    else if (x.a>x.b) and yp
                                           then begin
                                                  z1.a:=x.a/y.b;
                                                  z2.b:=x.b/y.b;
                                                  SetRoundMode (rmUp);
                                                  z1.b:=x.b/y.b;
                                                  z2.a:=x.a/y.b
                                                end
                                           else begin
                                                  z1.a:=x.b/y.a;
                                                  z2.b:=x.a/y.a;
                                                  SetRoundMode (rmUp);
                                                  z1.b:=x.a/y.a;
                                                  z2.a:=x.b/y.a
                                                end
                      else zero:=True;
             if zero
               then begin
                      SetRoundMode (rmNearest);
                      raise EZeroDivide.Create ('Division by an interval '
                                                +'containing 0.')
                    end
               else if dint_width(z1)>=dint_width(z2)
                      then Result:=z1
                      else Result:=z2
           end;
    SetRoundMode (rmNearest)
  end {didiv};

  class operator dinterval.Divide (x, y : dinterval) : dinterval;
  begin
    Result:=didiv(x, y)
  end {dDivide};

  procedure to_fixed_point (const awzi      : char_tab;
                            var significand : string);
  var exponent              : SmallInt;
      i, j, k, code         : Integer;
      remember, s1, s2, sum : Byte;
      short_sumz            : ShortString;
      sumz                  : string;
  begin
    exponent:=0;
    j:=1;
    for i:=16 downto 2 do
      begin
        if awzi[i]='1'
          then exponent:=exponent+j;
        j:=2*j
      end;
    exponent:=exponent-16383;
    for i:=80 downto 17 do
      if awzi[i]='1'
        then begin
               remember:=0;
               for j:=65 downto 3 do
                 begin
                   Val (significand[j], s1, code);
                   Val (ldi[i-17,j], s2, code);
                   sum:=s1+s2+remember;
                   Str (sum, short_sumz);
                   sumz:=string(short_sumz);
                   if sum>9
                     then begin
                            significand[j]:=sumz[2];
                            Val (sumz[1], remember, code);
                            if j=3
                              then begin
                                     Val (significand[1], s1, code);
                                     sum:=s1+remember;
                                     Str (sum, short_sumz);
                                     sumz:=string(short_sumz);
                                     significand[1]:=sumz[1]
                                   end
                          end
                     else begin
                            significand[j]:=sumz[1];
                            remember:=0
                          end
                 end;
               Val (significand[1], s1, code);
               Val (ldi[i-17,1], s2, code);
               sum:=s1+s2;
               Str (sum, short_sumz);
               sumz:=string(short_sumz);
               significand[1]:=sumz[1];
             end;
    if exponent>0
      then for i:=1 to exponent do
             begin
               j:=Length(significand);
               remember:=0;
               for k:=j downto j-62 do
                 begin
                   Val (significand[k], s1, code);
                   sum:=2*s1+remember;
                   Str (sum, short_sumz);
                   sumz:=string(short_sumz);
                   if sum>9
                     then begin
                            significand[k]:=sumz[2];
                            Val (sumz[1], remember, code)
                          end
                     else begin
                            significand[k]:=sumz[1];
                            remember:=0
                          end
                 end;
               for k:=j-64 downto 1 do
                 begin
                   Val (significand[k], s1, code);
                   sum:=2*s1+remember;
                   Str (sum, short_sumz);
                   sumz:=string(short_sumz);
                   if sum>9
                     then begin
                            significand[k]:=sumz[2];
                            Val (sumz[1], remember, code);
                            if k=1
                              then significand:=sumz[1]+significand
                          end
                     else begin
                            significand[k]:=sumz[1];
                            remember:=0
                          end
                 end
             end
      else if exponent<0
             then for i:=1 to -exponent do
                    begin
                      j:=Length(significand);
                      if significand[1]='1'
                        then begin
                               significand[1]:='0';
                               remember:=10
                             end
                        else remember:=0;
                      for k:=3 to j do
                        begin
                          Val (significand[k], s1, code);
                          sum:=remember+s1;
                          s1:=sum div 2;
                          Str (s1, short_sumz);
                          sumz:=string(short_sumz);
                          significand[k]:=sumz[1];
                          remember:=10*(sum mod 2);
                          if (k=j) and (remember<>0)
                            then significand:=significand+'5'
                        end
                    end;
    if awzi[1]='1'
      then significand:='-'+significand
      else significand:='+'+significand;
    if FormatSettings.DecimalSeparator=','
      then while (significand[Length(significand)]='0')
               and (significand[Length(significand)-1]<>',') do
             significand:=Copy(significand, 1, Length(significand)-1)
      else while (significand[Length(significand)]='0')
               and (significand[Length(significand)-1]<>'.') do
             significand:=Copy(significand, 1, Length(significand)-1)
  end {to_fixed_point};

  function int_read (const sa : string) : interval;
  var x, px, nx          : Extended;
      sx, sa1            : string;
      i, j               : Integer;
      tab                : array [1..10] of Byte absolute x;
      eps                : array [1..10] of Byte;
      epsx               : Extended absolute eps;
      epsw               : Word absolute eps;
      digits, rev_digits : char_tab;
      ix                 : interval;
      sep                : Char;
  begin
    sa1:=sa;
    if FormatSettings.DecimalSeparator=','
      then sep:=','
      else sep:='.';
    if (Pos('.', sa1)>0) and (FormatSettings.DecimalSeparator=',')
      then sa1[Pos('.', sa1)]:=',';
    x:=StrToFloat(sa1);
    if Pos('e', sa1)>0
      then sa1[Pos('e', sa1)]:='E';
    while sa1[1]=' ' do
      Delete (sa1, 1, 1);
    while sa1[Length(sa1)]=' ' do
      Delete (sa1, Length(sa1), 1);
    if (sa1[1]<>'-') and (sa1[1]<>'+')
      then Insert ('+', sa1, 1);
    while (sa1[2]='0') and (Length(sa1)>2) and (sa1[3]<>'e') and (sa1[3]<>'E')
        and (sa1[3]<>sep) do
      Delete (sa1, 2, 1);
    if (sa1[Length(sa1)]='E') or (sa1[Length(sa1)]='+')
        or (sa1[Length(sa1)]='-')
      then sa1:=sa1+'0'
      else if Pos('E', sa1)=0
             then sa1:=sa1+'E0';
    if Pos(sep, sa1)=0
      then Insert (sep+'0', sa1, Pos('E', sa1));
    sx:=Copy(sa1, Pos('E', sa1)+1, Length(sa1)-Pos('E', sa1));
    sa1:=Copy(sa1, 1, Pos('E', sa1)-1);
    j:=StrToInt(sx);
    if j>0
      then for i:=1 to j do
             begin
               Insert (sep, sa1, Pos(sep, sa1)+2);
               Delete (sa1, Pos(sep, sa1), 1);
               if Pos(sep, sa1)=Length(sa1)
                 then sa1:=sa1+'0'
             end
      else if j<0
             then for i:=j to -1 do
                    begin
                      Insert (sep, sa1, Pos(sep, sa1)-1);
                      Delete (sa1, Pos(sep, sa1)+2, 1);
                      if sa1[2]=sep
                        then Insert ('0', sa1, 2)
                    end;
    while (sa1[Length(sa1)]='0') and (sa1[Length(sa1)-1]<>sep) do
      sa1:=Copy(sa1, 1, Length(sa1)-1);
    for i:=1 to 10 do
      for j:=7 downto 0 do
        if tab[i] and bit[j] = bit[j]
          then digits[8*i-j]:='1'
          else digits[8*i-j]:='0';
    for i:=1 to 10 do
      for j:=1 to 8 do
        rev_digits[8*(i-1)+j]:=digits[80-8*i+j];
    sx:='0'+sep
        +'000000000000000000000000000000000000000000000000000000000000000';
    to_fixed_point (rev_digits, sx);
    if sa1=sx
      then begin
             ix.a:=x;
             ix.b:=x
           end
      else begin
             for i:=18 to 80 do
               rev_digits[i]:='0';
             rev_digits[17]:='1';
             rev_digits[1]:='0';
             for i:=1 to 2 do
               begin
                 eps[i]:=0;
                 for j:=1 to 8 do
                   if rev_digits[8*(i-1)+j]='1'
                     then eps[i]:=eps[i] or bit[8-j]
               end;
             epsw:=Swap(epsw);
             epsw:=epsw-63;
             epsw:=Swap(epsw);
             for i:=1 to 2 do
               for j:=7 downto 0 do
                 if eps[i] and bit[j] = bit[j]
                   then rev_digits[8*i-j]:='1'
                   else rev_digits[8*i-j]:='0';
             for i:=1 to 10 do
               for j:=1 to 8 do
                 digits[8*(i-1)+j]:=rev_digits[80-8*i+j];
             for i:=1 to 10 do
               begin
                 eps[i]:=0;
                 for j:=1 to 8 do
                   if digits[8*(i-1)+j]='1'
                     then eps[i]:=eps[i] or bit[8-j]
               end;
             px:=x-epsx;
             nx:=x+epsx;
             i:=Length(sa1)-Pos(sep, sa1);
             j:=Length(sx)-Pos(sep, sx);
             if j>i
               then i:=j;
             while Length(sa1)-Pos(sep, sa1)<i do
               sa1:=sa1+'0';
             while Length(sx)-Pos(sep, sx)<i do
               sx:=sx+'0';
             i:=Pos(sep, sa1);
             j:=Pos(sep, sx);
             if j>i
               then i:=j;
             while Pos(sep, sa1)<i do
               Insert ('0', sa1, 2);
             while Pos(sep, sx)<i do
               Insert ('0', sx, 2);
             if sx[1]='+'
               then if sa1<sx
                      then begin
                             ix.a:=px;
                             ix.b:=x
                           end
                      else begin
                             ix.a:=x;
                             ix.b:=nx
                           end
               else if sa1<sx
                      then begin
                             ix.a:=x;
                             ix.b:=nx
                           end
                      else begin
                             ix.a:=px;
                             ix.b:=x
                           end
           end;
    Result.a:=ix.a;
    Result.b:=ix.b
  end {int_read};

  function left_read (const sa : string) : Extended;
  var int_number : interval;
  begin
    int_number:=int_read(sa);
    Result:=int_number.a
  end {left_read};

  function right_read (const sa : string) : Extended;
  var int_number : interval;
  begin
    int_number:=int_read(sa);
    Result:=int_number.b
  end {right_read};

  function dleft_read (const sa : string) : Extended;
  var int_number : interval;
  begin
    int_number:=int_read(sa);
    Result:=int_number.b
  end {dleft_read};

  function dright_read (const sa : string) : Extended;
  var int_number : interval;
  begin
    int_number:=int_read(sa);
    Result:=int_number.a
  end {dright_read};

{$IFDEF WIN64}
  procedure to_fixed_point_Win64 (const awzi      : char_tab;
                                  var significand : string);
  var exponent              : SmallInt;
      i, j, k, code         : Integer;
      remember, s1, s2, sum : Byte;
      short_sumz            : ShortString;
      sumz                  : string;
      exponent_zero         : Boolean;
  begin
    exponent:=0;
    j:=1;
    for i:=12 downto 2 do
      begin
        if awzi[i]='1'
          then exponent:=exponent+j;
        j:=2*j
      end;
    if exponent=0
      then begin
             exponent_zero:=True;
             exponent:=-1022
           end
      else begin
             exponent_zero:=False;
             exponent:=exponent-1023
           end;
    for i:=64 downto 13 do
      if awzi[i]='1'
        then begin
               remember:=0;
               for j:=54 downto 3 do
                 begin
                   Val (significand[j], s1, code);
                   Val (ldi_double[i-12,j], s2, code);
                   sum:=s1+s2+remember;
                   Str (sum, short_sumz);
                   sumz:=string(short_sumz);
                   if sum>9
                     then begin
                            significand[j]:=sumz[2];
                            Val (sumz[1], remember, code);
                            if j=3
                              then begin
                                     Val (significand[1], s1, code);
                                     sum:=s1+remember;
                                     Str (sum, short_sumz);
                                     sumz:=string(short_sumz);
                                     significand[1]:=sumz[1]
                                   end
                          end
                     else begin
                            significand[j]:=sumz[1];
                            remember:=0
                          end
                 end;
               Val (significand[1], s1, code);
               Val (ldi_double[i-12,1], s2, code);
               sum:=s1+s2;
               Str (sum, short_sumz);
               sumz:=string(short_sumz);
               significand[1]:=sumz[1];
               if (i=13) and not exponent_zero
                 then begin
                        Val (ldi[0], s2, code);
                        sum:=sum+s2;
                        Str (sum, short_sumz);
                        sumz:=string(short_sumz);
                        significand[1]:=sumz[1]
                      end
             end;
    if (significand[1]='0') and not exponent_zero
      then significand[1]:='1';
    if exponent>0
      then for i:=1 to exponent do
             begin
               j:=Length(significand);
               remember:=0;
               for k:=j downto j-51 do
                 begin
                   Val (significand[k], s1, code);
                   sum:=2*s1+remember;
                   Str (sum, short_sumz);
                   sumz:=string(short_sumz);
                   if sum>9
                     then begin
                            significand[k]:=sumz[2];
                            Val (sumz[1], remember, code)
                          end
                     else begin
                            significand[k]:=sumz[1];
                            remember:=0
                          end
                 end;
               for k:=j-53 downto 1 do
                 begin
                   Val (significand[k], s1, code);
                   sum:=2*s1+remember;
                   Str (sum, short_sumz);
                   sumz:=string(short_sumz);
                   if sum>9
                     then begin
                            significand[k]:=sumz[2];
                            Val (sumz[1], remember, code);
                            if k=1
                              then significand:=sumz[1]+significand
                          end
                     else begin
                            significand[k]:=sumz[1];
                            remember:=0
                          end
                 end
             end
      else if exponent<0
             then for i:=1 to -exponent do
                    begin
                      j:=Length(significand);
                      if significand[1]='1'
                        then begin
                               significand[1]:='0';
                               remember:=10
                             end
                        else remember:=0;
                      for k:=3 to j do
                        begin
                          Val (significand[k], s1, code);
                          sum:=remember+s1;
                          s1:=sum div 2;
                          Str (s1, short_sumz);
                          sumz:=string(short_sumz);
                          significand[k]:=sumz[1];
                          remember:=10*(sum mod 2);
                          if (k=j) and (remember<>0)
                            then significand:=significand+'5'
                        end
                    end;
    if awzi[1]='1'
      then significand:='-'+significand
      else significand:='+'+significand;
    if FormatSettings.DecimalSeparator=','
      then while (significand[Length(significand)]='0')
               and (significand[Length(significand)-1]<>',') do
             significand:=Copy(significand, 1, Length(significand)-1)
      else while (significand[Length(significand)]='0')
               and (significand[Length(significand)-1]<>'.') do
             significand:=Copy(significand, 1, Length(significand)-1)
  end {to_fixed_point_Win64};

  function int_read_Win64 (const sa : string) : interval_double;
  var x, px, nx          : Double;
      sx, sa1            : string;
      i, j               : Integer;
      tab                : array [1..8] of Byte absolute x;
      eps                : array [1..8] of Byte;
      epsx               : Double absolute eps;
      epsw               : Word;
      digits, rev_digits : char_tab;
      ix                 : interval_double;
      sep                : Char;
  begin
    sa1:=sa;
    if FormatSettings.DecimalSeparator=','
      then sep:=','
      else sep:='.';
    if (Pos('.', sa1)>0) and (FormatSettings.DecimalSeparator=',')
      then sa1[Pos('.', sa1)]:=',';
    x:=StrToFloat(sa1);
    if (sa1[1]<>' ')
      then Insert ('+', sa1, 1);
    sx:=Copy(sa1, Pos('E', sa1)+1, Length(sa1)-Pos('E', sa1));
    sa1:=Copy(sa1, 1, Pos('E', sa1)-1);
    j:=StrToInt(sx);
    if j>0
      then for i:=1 to j do
             begin
               Insert (sep, sa1, Pos(sep, sa1)+2);
               Delete (sa1, Pos(sep, sa1), 1);
               if Pos(sep, sa1)=Length(sa1)
                 then sa1:=sa1+'0'
             end
      else if j<0
             then for i:=j to -1 do
                    begin
                      Insert (sep, sa1, Pos(sep, sa1)-1);
                      Delete (sa1, Pos(sep, sa1)+2, 1);
                      if sa1[2]=sep
                        then Insert ('0', sa1, 2)
                    end;
    while (sa1[Length(sa1)]='0') and (sa1[Length(sa1)-1]<>sep) do
      sa1:=Copy(sa1, 1, Length(sa1)-1);
    for i:=1 to 8 do
      for j:=7 downto 0 do
        if tab[i] and bit[j] = bit[j]
          then digits[8*i-j]:='1'
          else digits[8*i-j]:='0';
    for i:=1 to 8 do
      for j:=1 to 8 do
        rev_digits[8*(i-1)+j]:=digits[64-8*i+j];
    sx:='0'+sep
        +'0000000000000000000000000000000000000000000000000000';
    to_fixed_point_Win64 (rev_digits, sx);
    if sa1=sx
      then begin
             ix.a:=x;
             ix.b:=x
           end
      else begin
             for i:=13 to 64 do
               rev_digits[i]:='0';
             rev_digits[1]:='0';
             epsw:=0;
             j:=1;
             for i:=10 downto 2 do
               begin
                 if rev_digits[i]='1'
                   then epsw:=epsw+j;
                 j:=2*j
               end;
             epsw:=epsw-13;
             for i:=2 to 9 do
               begin
                 j:=j div 2;
                 if epsw div j =1
                   then begin
                          rev_digits[i]:='1';
                          epsw:=epsw-j
                        end
                   else rev_digits[i]:='0'
               end;
             if epsw=1
               then rev_digits[10]:='1'
               else rev_digits[10]:='0';
             for i:=1 to 8 do
               for j:=1 to 8 do
                 digits[8*(i-1)+j]:=rev_digits[64-8*i+j];
             for i:=1 to 8 do
               begin
                 eps[i]:=0;
                 for j:=1 to 8 do
                   if digits[8*(i-1)+j]='1'
                     then eps[i]:=eps[i] or bit[8-j]
               end;
             px:=x-epsx;
             nx:=x+epsx;
             i:=Length(sa1)-Pos(sep, sa1);
             j:=Length(sx)-Pos(sep, sx);
             if j>i
               then i:=j;
             while Length(sa1)-Pos(sep, sa1)<i do
               sa1:=sa1+'0';
             while Length(sx)-Pos(sep, sx)<i do
               sx:=sx+'0';
             i:=Pos(sep, sa1);
             j:=Pos(sep, sx);
             if j>i
               then i:=j;
             while Pos(sep, sa1)<i do
               Insert ('0', sa1, 2);
             while Pos(sep, sx)<i do
               Insert ('0', sx, 2);
             if sx[1]='+'
               then if sa1<sx
                      then begin
                             ix.a:=px;
                             ix.b:=x
                           end
                      else begin
                             ix.a:=x;
                             ix.b:=nx
                           end
               else if sa1<sx
                      then begin
                             ix.a:=x;
                             ix.b:=nx
                           end
                      else begin
                             ix.a:=px;
                             ix.b:=x
                           end
           end;
    Result.a:=ix.a;
    Result.b:=ix.b
  end {int_read_Win64};

  function left_read_Win64 (const sa : string) : Double;
  var int_number : interval_double;
  begin
    int_number:=int_read_Win64(sa);
    Result:=int_number.a
  end {left_read_Win64};

  function right_read_Win64 (const sa : string) : Double;
  var int_number : interval_double;
  begin
    int_number:=int_read_Win64(sa);
    Result:=int_number.b
  end {right_read_Win64};
{$ENDIF}

  procedure iends_to_strings (const x         : interval;
                              out left, right : string);
  procedure modify_mantissa (const i      : Integer;
                             var mantissa : string);
  var s, s1    : string;
      short_s1 : ShortString;
  begin
    if i>=0
{$IFDEF WIN64}
      then Insert ('+', mantissa, 18)
      else Insert ('-', mantissa, 18);
{$ELSE}
      then Insert ('+', mantissa, 21)
      else Insert ('-', mantissa, 21);
{$ENDIF}
    Str (System.Abs(i), short_s1);
    s1:=string(short_s1);
    if i<10
      then s:=string('000')+s1
      else if i<100
             then s:=string('00')+s1
             else if i<1000
                    then s:=string('0')+s1
                    else s:=s1;
{$IFDEF WIN64}
    Insert (s, mantissa, 19)
{$ELSE}
    Insert (s, mantissa, 22)
{$ENDIF}
  end;
  function take_up (var fl_str : string) : string;
  var s, s1      : string;
      short_s    : ShortString;
      code, i, k : Integer;
      finished   : Boolean;
  begin
    finished:=False;
{$IFDEF WIN64}
    k:=16;
{$ELSE}
    k:=19;
{$ENDIF}
    repeat
      s:=Copy(fl_str, k, 1);
      Delete (fl_str, k, 1);
      Val (s, i, code);
      i:=i+1;
      if i<10
        then begin
               Str (i, short_s);
               s:=string(short_s);
               Insert (s, fl_str, k);
               finished:=True
             end
        else begin
               Insert ('0', fl_str, k);
               k:=k-1
             end
    until finished or (k<4);
    if not finished
      then begin
             s:=Copy(fl_str, 2, 1);
             Delete (fl_str, 2, 1);
             Val (s, i, code);
             i:=i+1;
             if i<10
               then begin
                      Str (i, short_s);
                      s:=string(short_s);
                      Insert (s, fl_str, 2)
                    end
               else begin
                      Insert ('1', fl_str, 2);
                      s:='0';
{$IFDEF WIN64}
                      for k:=4 to 16 do
{$ELSE}
                      for k:=4 to 19 do
{$ENDIF}
                        begin
                          s1:=Copy(fl_str, k, 1);
                          Delete (fl_str, k, 1);
                          Insert (s, fl_str, k);
                          s:=s1
                        end;
{$IFDEF WIN64}
                      s:=Copy(fl_str, 18, 5);
                      Delete (fl_str, 18, 5);
{$ELSE}
                      s:=Copy(fl_str, 21, 5);
                      Delete (fl_str, 21, 5);
{$ENDIF}
                      Val (s, i, code);
                      i:=i-1;
                      modify_mantissa (i, fl_str)
                    end
           end;
    Result:=fl_str
  end;
  function take_down (var fl_str : string) : string;
  var s          : string;
      short_s    : ShortString;
      code, i, k : Integer;
      finished   : Boolean;
  begin
    finished:=False;
{$IFDEF WIN64}
    k:=16;
{$ELSE}
    k:=19;
{$ENDIF}
    repeat
      s:=Copy(fl_str, k, 1);
      Delete (fl_str, k, 1);
      Val (s, i, code);
      i:=i-1;
      if i>-1
        then begin
               Str (i, short_s);
               s:=string(short_s);
               Insert (s, fl_str, k);
               finished:=True
             end
        else begin
               Insert ('9', fl_str, k);
               k:=k-1
             end
    until finished or (k<4);
    if not finished
      then begin
             s:=Copy(fl_str, 2, 1);
             Delete (fl_str, 2, 1);
             Val (s, i, code);
             i:=i-1;
             if i>0
               then begin
                      Str (i, short_s);
                      s:=string(short_s);
                      Insert (s, fl_str, 2)
                    end
               else begin
                      s:=Copy(fl_str, 4, 1);
                      Insert (s, fl_str, 2);
{$IFDEF WIN64}
                      for k:=4 to 15 do
{$ELSE}
                      for k:=4 to 18 do
{$ENDIF}
                        begin
                          s:=Copy(fl_str, k+1, 1);
                          Delete (fl_str, k+1, 1);
                          Insert (s, fl_str, k)
                        end;
{$IFDEF WIN64}
                      Delete (fl_str, 16, 1);
                      Insert ('9', fl_str, 16);
                      s:=Copy(fl_str, 18, 5);
                      Delete (fl_str, 18, 5);
{$ELSE}
                      Delete (fl_str, 19, 1);
                      Insert ('9', fl_str, 19);
                      s:=Copy(fl_str, 21, 5);
                      Delete (fl_str, 21, 5);
{$ENDIF}
                      Val (s, i, code);
                      i:=i-1;
                      modify_mantissa (i, fl_str)
                    end
           end;
    Result:=fl_str
  end;
  var code                    : Integer;
{$IFDEF WIN64}
      y, z                    : Double;
{$ELSE}
      y, z                    : Extended;
{$ENDIF}
      short_left, short_right : ShortString;
  begin
    if x.a<=x.b
      then if x.a>=0
             then begin
{$IFDEF WIN64}
                    Str (x.a:23, short_left);
{$ELSE}
                    Str (x.a:26, short_left);
{$ENDIF}
                    left:=string(short_left);
{$IFDEF WIN64}
                    Delete (left, 17, 1);
                    y:=right_read(left);
{$ELSE}
                    Delete (left, 20, 1);
{$ENDIF}
                    Val (left, z, code);
{$IFDEF WIN64}
                    if (x.a<z) or (y<>x.a)
{$ELSE}
                    if x.a<z
{$ENDIF}
                      then left:=take_down(left);
{$IFDEF WIN64}
                    Str (x.b:22, short_right);
{$ELSE}
                    Str (x.b:25, short_right);
{$ENDIF}
                    right:=string(short_right);
{$IFDEF WIN64}
                    y:=left_read_Win64(right);
{$ELSE}
                    y:=left_read(right);
{$ENDIF}
                    Val (right, z, code);
                    if (x.b>=z) and (x.a<>x.b) and (y<>x.b)
                      then right:=take_up(right)
                  end
             else if x.b<=0
                    then begin
{$IFDEF WIN64}
                           Str (x.a:22, short_left);
{$ELSE}
                           Str (x.a:25, short_left);
{$ENDIF}
                           left:=string(short_left);
{$IFDEF WIN64}
                           y:=right_read_Win64(left);
{$ELSE}
                           y:=right_read(left);
{$ENDIF}
                           Val (left, z, code);
                           if (x.a<=z) and (x.a<>x.b) and (y<>x.a)
                             then left:=take_up(left);
{$IFDEF WIN64}
                           Str (x.b:23, short_right);
{$ELSE}
                           Str (x.b:26, short_right);
{$ENDIF}
                           right:=string(short_right);
{$IFDEF WIN64}
                           Delete (right, 17, 1);
                           y:=left_read(right);
{$ELSE}
                           Delete (right, 20, 1);
{$ENDIF}
                           Val (right, z, code);
{$IFDEF WIN64}
                           if (x.b>z) or (y<>x.b)
{$ELSE}
                           if x.b>z
{$ENDIF}
                             then right:=take_down(right)
                         end
                    else begin
{$IFDEF WIN64}
                           Str (x.a:22, short_left);
{$ELSE}
                           Str (x.a:25, short_left);
{$ENDIF}
                           left:=string(short_left);
{$IFDEF WIN64}
                           y:=right_read_Win64(left);
{$ELSE}
                           y:=right_read(left);
{$ENDIF}
                           Val (left, z, code);
                           if (x.a<=z) and (y<>x.a)
                             then left:=take_up(left);
{$IFDEF WIN64}
                           Str (x.b:22, short_right);
{$ELSE}
                           Str (x.b:25, short_right);
{$ENDIF}
                           right:=string(short_right);
{$IFDEF WIN64}
                           y:=left_read_Win64(right);
{$ELSE}
                           y:=left_read(right);
{$ENDIF}
                           Val (right, z, code);
                           if (x.b>=z) and (y<>x.b)
                             then right:=take_up(right)
                         end
  end {iends_to_strings};

  function iabs (const x : interval) : interval;
  begin
    if Abs(x.b)>=Abs(x.a)
      then Result.b:=Abs(x.b)
      else Result.b:=Abs(x.a);
    if x.b<=0
      then Result.a:=Abs(x.b)
      else if x.a<=0
             then if Abs(x.a)<x.b
                    then Result.a:=Abs(x.a)
                    else Result.a:=x.b
  end {iabs};

  function isin (const x : interval;
                 out st  : Integer) : interval;
  var finished, is_even : Boolean;
      k                 : Integer;
      d, s, w, w1, x2   : interval;
  begin
    if x.a>x.b
      then st:=1
      else begin
             s:=x;
             w:=x;
             x2:=isqr(x,st);
             k:=1;
             is_even:=True;
             finished:=False;
             st:=0;
             repeat
               d.a:=(k+1)*(k+2);
               d.b:=d.a;
               s:=imul(s,idiv(x2,d));
               if is_even
                 then w1:=isub(w,s)
                 else w1:=iadd(w,s);
               if (w.a<>0) and (w.b<>0)
                 then if (Abs(w.a-w1.a)/Abs(w.a)<1e-18)
                          and (Abs(w.b-w1.b)/Abs(w.b)<1e-18)
                        then finished:=True
                        else
                 else if (w.a=0) and (w.b<>0)
                        then if (Abs(w.a-w1.a)<1e-18)
                                 and (Abs(w.b-w1.b)/Abs(w.b)<1e-18)
                               then finished:=True
                               else
                         else if w.a<>0
                                then if (Abs(w.a-w1.a)/Abs(w.a)<1e-18)
                                         and (Abs(w.b-w1.b)<1e-18)
                                       then finished:=True
                                       else
                                else if (Abs(w.a-w1.a)<1e-18)
                                         and (Abs(w.b-w1.b)<1e-18)
                                       then finished:=True;
               if finished
                 then begin
                        if w1.b>1
                          then begin
                                 w1.b:=1;
                                 if w1.a>1
                                   then w1.a:=1
                               end;
                        if w1.a<-1
                          then begin
                                 w1.a:=-1;
                                 if w1.b<-1
                                   then w1.b:=-1
                               end;
                        Result:=w1
                      end
                 else begin
                        w:=w1;
                        k:=k+2;
                        is_even:=not is_even
                      end
             until finished or (k>MaxInt/2);
             if not finished
               then st:=2
           end
  end {isin};

  function icos (const x : interval;
                 out st  : Integer) : interval;
  var finished, is_even : Boolean;
      k                 : Integer;
      c, d, w, w1, x2   : interval;
  begin
    if x.a>x.b
      then st:=1
      else begin
             c.a:=1;
             c.b:=1;
             w:=c;
             x2:=isqr(x,st);
             k:=1;
             is_even:=True;
             finished:=False;
             st:=0;
             repeat
               d.a:=k*(k+1);
               d.b:=d.a;
               c:=imul(c,idiv(x2,d));
               if is_even
                 then w1:=isub(w,c)
                 else w1:=iadd(w,c);
               if (w.a<>0) and (w.b<>0)
                 then if (Abs(w.a-w1.a)/Abs(w.a)<1e-18)
                          and (Abs(w.b-w1.b)/Abs(w.b)<1e-18)
                        then finished:=True
                        else
                 else if (w.a=0) and (w.b<>0)
                        then if (Abs(w.a-w1.a)<1e-18)
                                 and (Abs(w.b-w1.b)/Abs(w.b)<1e-18)
                               then finished:=True
                               else
                         else if w.a<>0
                                then if (Abs(w.a-w1.a)/Abs(w.a)<1e-18)
                                         and (Abs(w.b-w1.b)<1e-18)
                                       then finished:=True
                                       else
                                else if (Abs(w.a-w1.a)<1e-18)
                                         and (Abs(w.b-w1.b)<1e-18)
                                       then finished:=True;
               if finished
                 then begin
                        if w1.b>1
                          then begin
                                 w1.b:=1;
                                 if w1.a>1
                                   then w1.a:=1
                               end;
                        if w1.a<-1
                          then begin
                                 w1.a:=-1;
                                 if w1.b<-1
                                   then w1.b:=-1
                               end;
                        Result:=w1
                      end
                 else begin
                        w:=w1;
                        k:=k+2;
                        is_even:=not is_even
                      end
             until finished or (k>MaxInt/2);
             if not finished
               then st:=2
           end
  end {icos};

  function itan (const x : interval;
                 out st  : Integer) : interval;
// |x|<Pi/2
  var intc, ints : interval;
  begin
    if x.a>x.b
      then st:=1
      else if (x.b>=Pi/2) or (x.a<=-Pi/2)
             then st:=3
             else begin
                    ints:=isin(x,st);
                    if st=0
                      then begin
                             intc:=icos(x,st);
                             if st=0
                               then Result:=idiv(ints,intc)
                           end
                  end
  end {itan};

  function icot (const x : interval;
                 out st  : Integer) : interval;
// 0<|x|<Pi
  var intc, ints : interval;
  begin
    if x.a>x.b
      then st:=1
      else if (x.a<=0) and (x.b>=0) or (x.a<=-Pi) or (x.b>=Pi)
             then st:=3
             else begin
                    intc:=icos(x,st);
                    if st=0
                      then begin
                             ints:=isin(x,st);
                             if st=0
                               then Result:=idiv(intc,ints)
                           end
                  end
  end {icot};

  function iarcsin (const x : interval;
                    out st  : Integer) : interval;
// |x|<1
// Important note: For x.b>0.998 it is assumed Result.b=Pi/2,
//                 and for x.a<-0.998 it is assumed Result.a=-Pi/2.
//                 In both cases st=4 and the resulting interval
//                 may not be satisfactory.
//                 For x close to 1 the calculations are time expensive.
  var finished                 : Boolean;
      k                        : Integer;
      d1, d2, d3, s, w, w1, x2 : interval;
  begin
    if x.a>x.b
      then st:=1
      else if (x.a<=-1) or (x.b>=1)
             then st:=3
             else begin
                    s:=x;
                    w:=s;
                    if x.a<-0.998
                      then w.a:=-Pi/2
                      else if x.b>0.998
                             then w.b:=Pi/2;
                      x2:=isqr(x,st);
                      k:=1;
                      finished:=False;
                      st:=0;
                      repeat
                        d1.a:=k;
                        d1.b:=d1.a;
                        d2.a:=k+1;
                        d2.b:=d2.a;
                        d3.a:=k+2;
                        d3.b:=d3.a;
                        s:=imul(s,idiv(imul(d1,x2),d2));
                        w1:=iadd(w,idiv(s,d3));
                        if x.a<-0.998
                          then w1.a:=-Pi/2
                          else if x.b>0.998
                                 then w1.b:=Pi/2;
                        if (w.a<>0) and (w.b<>0)
                          then if (Abs(w.a-w1.a)/Abs(w.a)<1e-18)
                                   and (Abs(w.b-w1.b)/Abs(w.b)<1e-18)
                                 then finished:=True
                                 else
                          else if (w.a=0) and (w.b<>0)
                                 then if (Abs(w.a-w1.a)<1e-18)
                                          and (Abs(w.b-w1.b)/Abs(w.b)<1e-18)
                                        then finished:=True
                                        else
                                 else if w.a<>0
                                        then if (Abs(w.a-w1.a)/Abs(w.a)<1e-18)
                                                 and (Abs(w.b-w1.b)<1e-18)
                                               then finished:=True
                                               else
                                        else if (Abs(w.a-w1.a)<1e-18)
                                                 and (Abs(w.b-w1.b)<1e-18)
                                               then finished:=True;
                        if finished
                          then begin
                                 Result:=w1;
                                 if (w1.b=Pi/2) or (w1.a=-Pi/2)
                                   then st:=4
                               end
                          else begin
                                 w:=w1;
                                 k:=k+2
                               end
                      until finished or (k>MaxInt/2);
                    if not finished
                      then st:=2
                  end
  end {iarcsin};

  function iarccos (const x : interval;
                    out st  : Integer) : interval;
// |x|<1
// Important note: For x.b>0.998 it is assumed Result.a=0,
//                 and for x.a<-0.998 it is assumed Result.b=Pi.
//                 In both cases st=4 and the resulting interval
//                 may not be satisfactory.
//                 For x close to 1 the calculations are time expensive.
  const d : interval = (a:2; b:2);
  var s : interval;
  begin
    if x.a>x.b
      then st:=1
      else if (x.a<=-1) or (x.b>=1)
             then st:=3
             else begin
                    s:=iarcsin(x,st);
                    if (st=0) or (st=4)
                      then Result:=isub(idiv(ipi,d),s);
                    if st=4
                      then if x.b>0.998
                             then Result.a:=0
                             else Result.b:=Pi
                  end
  end {iarccos};

  function iarctan (const x : interval;
                    out st  : Integer) : interval;
  var s : interval;
  begin
    if x.a>x.b
      then st:=1
      else begin
             s.a:=1;
             s.b:=s.a;
             s:=iadd(s,isqr(x,st));
             s:=isqrt(s,st);
             s:=idiv(x,s);
             Result:=iarcsin(s,st)
           end
  end {iarctan};

  function iarccot (const x : interval;
                    out st  : Integer) : interval;
  var c : interval;
  begin
    if x.a>x.b
      then st:=1
      else begin
             c.a:=1;
             c.b:=c.a;
             c:=iadd(c,isqr(x,st));
             c:=idiv(x,isqrt(c,st));
             Result:=iarccos(c,st)
           end
  end {iarccot};

  function isinh (const x : interval;
                  out st  : Integer) : interval;
  var finished        : Boolean;
      k               : Integer;
      d, s, w, w1, x2 : interval;
  begin
    if x.a>x.b
      then st:=1
      else begin
             s:=x;
             w:=x;
             x2:=isqr(x,st);
             k:=1;
             finished:=False;
             st:=0;
             repeat
               d.a:=(k+1)*(k+2);
               d.b:=d.a;
               s:=imul(s,idiv(x2,d));
               w1:=iadd(w,s);
               if (w.a<>0) and (w.b<>0)
                 then if (Abs(w.a-w1.a)/Abs(w.a)<1e-18)
                          and (Abs(w.b-w1.b)/Abs(w.b)<1e-18)
                        then finished:=True
                        else
                 else if (w.a=0) and (w.b<>0)
                        then if (Abs(w.a-w1.a)<1e-18)
                                 and (Abs(w.b-w1.b)/Abs(w.b)<1e-18)
                               then finished:=True
                               else
                         else if w.a<>0
                                then if (Abs(w.a-w1.a)/Abs(w.a)<1e-18)
                                         and (Abs(w.b-w1.b)<1e-18)
                                       then finished:=True
                                       else
                                else if (Abs(w.a-w1.a)<1e-18)
                                         and (Abs(w.b-w1.b)<1e-18)
                                       then finished:=True;
               if finished
                 then Result:=w1
                 else begin
                        w:=w1;
                        k:=k+2
                      end
             until finished or (k>MaxInt/2);
             if not finished
               then st:=2
           end
  end {isinh};

  function icosh (const x : interval;
                  out st  : Integer) : interval;
  var finished        : Boolean;
      k               : Integer;
      c, d, w, w1, x2 : interval;
  begin
    if x.a>x.b
      then st:=1
      else begin
             c.a:=1;
             c.b:=1;
             w:=c;
             x2:=isqr(x,st);
             k:=1;
             finished:=False;
             st:=0;
             repeat
               d.a:=k*(k+1);
               d.b:=d.a;
               c:=imul(c,idiv(x2,d));
               w1:=iadd(w,c);
               if (w.a<>0) and (w.b<>0)
                 then if (Abs(w.a-w1.a)/Abs(w.a)<1e-18)
                          and (Abs(w.b-w1.b)/Abs(w.b)<1e-18)
                        then finished:=True
                        else
                 else if (w.a=0) and (w.b<>0)
                        then if (Abs(w.a-w1.a)<1e-18)
                                 and (Abs(w.b-w1.b)/Abs(w.b)<1e-18)
                               then finished:=True
                               else
                         else if w.a<>0
                                then if (Abs(w.a-w1.a)/Abs(w.a)<1e-18)
                                         and (Abs(w.b-w1.b)<1e-18)
                                       then finished:=True
                                       else
                                else if (Abs(w.a-w1.a)<1e-18)
                                         and (Abs(w.b-w1.b)<1e-18)
                                       then finished:=True;
               if finished
                 then begin
                        if w1.a<1
                          then w1.a:=1;
                        Result:=w1
                      end
                 else begin
                        w:=w1;
                        k:=k+2
                      end
             until finished or (k>MaxInt/2);
             if not finished
               then st:=2
           end
  end {icosh};

  function itanh (const x : interval;
                  out st  : Integer) : interval;
  var c, s : interval;
  begin
    if x.a>x.b
      then st:=1
      else begin
             s:=isinh(x,st);
             if st=0
               then begin
                      c:=icosh(x,st);
                      if st=0
                        then Result:=idiv(s,c)
                    end
           end
  end {itanh};

  function icoth (const x : interval;
                  out st  : Integer) : interval;
  var c, s : interval;
  begin
    if x.a>x.b
      then st:=1
      else begin
             c:=icosh(x,st);
             if st=0
               then begin
                      s:=isinh(x,st);
                      if st=0
                        then Result:=idiv(c,s)
                    end
           end
  end {icoth};

  function iarcsinh (const x : interval;
                     out st  : Integer) : interval;
  const d : interval = (a:1; b:1);
  var s : interval;
  begin
    if x.a>x.b
      then st:=1
      else begin
             s:=isqr(x,st);
             if st=0
               then begin
                      s:=isqrt(iadd(s,d),st);
                      if st=0
                        then begin
                               s:=iln(iadd(x,s),st);
                               if st=0
                                 then Result:=s
                             end
                    end
           end
  end {iarcsinh};

  function iarccosh (const x : interval;
                     out st  : Integer) : interval;
// x>=1
  const d : interval = (a:1; b:1);
  var c : interval;
  begin
    if x.a>x.b
      then st:=1
      else if x.a<1
             then st:=3
             else begin
                    c:=isqr(x,st);
                    if st=0
                      then begin
                             c:=isqrt(isub(c,d),st);
                             if st=0
                               then begin
                                      c:=iln(iadd(x,c),st);
                                      if st=0
                                        then Result:=c
                                    end
                           end
                  end
  end {iarccosh};

  function iarctanh (const x : interval;
                     out st  : Integer) : interval;
// |x|<1
  const d1 : interval = (a:1; b:1);
        d2 : interval = (a:2; b:2);
  var t : interval;
  begin
    if x.a>x.b
      then st:=1
      else if (x.a<=-1) or (x.b>=1)
             then st:=3
             else begin
                    t:=idiv(iadd(d1,x),isub(d1,x));
                    t:=iln(t,st);
                    if st=0
                      then Result:=idiv(t,d2)
                  end
  end {iarctanh};

  function iarccoth (const x : interval;
                     out st  : Integer) : interval;
// |x|>1
  const d1 : interval = (a:1; b:1);
        d2 : interval = (a:2; b:2);
  var c : interval;
  begin
    if x.a>x.b
      then st:=1
      else if (x.a<=1) or (x.b>=-1)
             then st:=3
             else begin
                    c:=idiv(iadd(x,d1),isub(x,d1));
                    c:=iln(c,st);
                    if st=0
                      then Result:=idiv(c,d2)
                  end
  end {iarccoth};

  function iln (const x : interval;
                out st  : Integer) : interval;
// x>0
  var finished         : Boolean;
      k                : Integer;
      d, w, w1, w2, x1 : interval;
  begin
    if x.a>x.b
      then st:=1
      else if x.a<=0
             then st:=3
             else begin
                    d.a:=1;
                    d.b:=d.a;
                    x1:=idiv(isub(x,d),iadd(x,d));
                    w:=x1;
                    w2:=w;
                    k:=1;
                    finished:=False;
                    st:=0;
                    repeat
                      d.a:=k+2;
                      d.b:=d.a;
                      w2:=imul(w2,imul(x1,x1));
                      w1:=iadd(w,idiv(w2,d));
                      if (w.a<>0) and (w.b<>0)
                        then if (Abs(w.a-w1.a)/Abs(w.a)<1e-18)
                                 and (Abs(w.b-w1.b)/Abs(w.b)<1e-18)
                               then finished:=True
                               else
                        else if (w.a=0) and (w.b<>0)
                               then if (Abs(w.a-w1.a)<1e-18)
                                        and (Abs(w.b-w1.b)/Abs(w.b)<1e-18)
                                      then finished:=True
                                      else
                               else if w.a<>0
                                      then if (Abs(w.a-w1.a)/Abs(w.a)<1e-18)
                                               and (Abs(w.b-w1.b)<1e-18)
                                             then finished:=True
                                             else
                                      else if (Abs(w.a-w1.a)<1e-18)
                                               and (Abs(w.b-w1.b)<1e-18)
                                             then finished:=True;
                      if finished
                        then begin
                               d.a:=2;
                               d.b:=d.a;
                               Result:=imul(d,w1)
                             end
                        else begin
                               w:=w1;
                               k:=k+2
                             end
                    until finished or (k>MaxInt/2);
                    if not finished
                      then st:=2
                  end
  end {iln};

  function iexp (const x : interval;
                 out st  : Integer) : interval;
  var finished    : Boolean;
      k           : Integer;
      d, e, w, w1 : interval;
  begin
    if x.a>x.b
      then st:=1
      else begin
             e.a:=1;
             e.b:=1;
             w:=e;
             k:=1;
             finished:=False;
             st:=0;
             repeat
               d.a:=k;
               d.b:=k;
               e:=imul(e,idiv(x,d));
               w1:=iadd(w,e);
               if (Abs(w.a-w1.a)/Abs(w.a)<1e-18)
                   and (Abs(w.b-w1.b)/Abs(w.b)<1e-18)
                 then begin
                        finished:=True;
                        Result:=w1
                      end
                 else begin
                        w:=w1;
                        k:=k+1
                      end
             until finished or (k>MaxInt/2);
             if not finished
               then st:=2
           end
  end {iexp};

  function isqr (const x : interval;
                 out st  : Integer) : interval;
  var maxx, minx : Extended;
  begin
    if x.a>x.b
      then st:=1
      else begin
             st:=0;
             if (x.a<=0) and (x.b>=0)
               then minx:=0
               else if x.a>0
                      then minx:=x.a
                      else minx:=x.b;
             if Abs(x.a)>Abs(x.b)
               then maxx:=Abs(x.a)
               else maxx:=Abs(x.b);
             SetRoundMode (rmDown);
             Result.a:=minx*minx;
             SetRoundMode (rmUp);
             Result.b:=maxx*maxx;
             SetRoundMode (rmNearest)
           end
  end {isqr};

  function isqrt (const x : interval;
                  out st  : Integer) : interval;
// x>=0
  begin
    if x.a>x.b
      then st:=1
      else if x.a<0
             then st:=3
             else begin
                    st:=0;
                    SetRoundMode (rmDown);
                    Result.a:=Sqrt(x.a);
                    SetRoundMode (rmUp);
                    Result.b:=Sqrt(x.b);
                    SetRoundMode (rmNearest)
                  end
  end {isqrt};

  function iax (const a, x : interval;
                out st     : Integer) : interval;
// a>0, a<>1
  var finished          : Boolean;
      k                 : Integer;
      d, w, w1, w2, xln : interval;
  begin
    if (a.a>a.b) or (x.a>x.b)
      then st:=1
      else if (a.b<=0) or (a.a<=1) and (a.b>=1)
             then st:=3
             else begin
                    w:=iln(a,st);
                    if st=0
                      then begin
                             xln:=imul(x,w);
                             w2:=xln;
                             d.a:=1;
                             d.b:=d.a;
                             w:=iadd(d,xln);
                             k:=1;
                             finished:=False;
                             st:=0;
                             repeat
                               d.a:=k+1;
                               d.b:=d.a;
                               w2:=idiv(imul(w2,xln),d);
                               w1:=iadd(w,w2);
                               if (w.a<>0) and (w.b<>0)
                                 then if (Abs(w.a-w1.a)/Abs(w.a)<1e-18)
                                          and (Abs(w.b-w1.b)/Abs(w.b)<1e-18)
                                        then finished:=True
                                        else
                                 else if (w.a=0) and (w.b<>0)
                                        then if (Abs(w.a-w1.a)<1e-18)
                                                 and (Abs(w.b-w1.b)/Abs(w.b)
                                                      <1e-18)
                                               then finished:=True
                                               else
                                        else if w.a<>0
                                               then if (Abs(w.a-w1.a)/Abs(w.a)
                                                        <1e-18)
                                                        and (Abs(w.b-w1.b)
                                                             <1e-18)
                                                      then finished:=True
                                                      else
                                               else if (Abs(w.a-w1.a)<1e-18)
                                                        and (Abs(w.b-w1.b)
                                                             <1e-18)
                                                      then finished:=True;
                               if finished
                                 then Result:=w1
                                 else begin
                                        w:=w1;
                                        k:=k+1
                                      end
                             until finished or (k>MaxInt/2);
                             if not finished
                               then st:=2
                           end
                  end
  end {iax};

  function ionepxtom (const x, m : interval;
                      out st     : Integer) : interval;
// |x|<=1, m>0
  var finished         : Boolean;
      k                : Integer;
      d, km, w, w1, w2 : interval;
  begin
    if (x.a>x.b) or (m.a>m.b)
      then st:=1
      else if (x.a<-1) or (x.b>1) or (m.b<=0)
             then st:=3
             else begin
                    d.a:=1;
                    d.b:=d.a;
                    w:=d;
                    w2:=d;
                    km:=m;
                    k:=1;
                    finished:=False;
                    st:=0;
                    repeat
                      d.a:=k;
                      d.b:=d.a;
                      w2:=imul(idiv(imul(km,w2),d),x);
                      w1:=iadd(w,w2);
                      if (w.a<>0) and (w.b<>0)
                        then if (Abs(w.a-w1.a)/Abs(w.a)<1e-18)
                                 and (Abs(w.b-w1.b)/Abs(w.b)<1e-18)
                               then finished:=True
                               else
                        else if (w.a=0) and (w.b<>0)
                               then if (Abs(w.a-w1.a)<1e-18)
                                        and (Abs(w.b-w1.b)/Abs(w.b)<1e-18)
                                      then finished:=True
                                      else
                               else if w.a<>0
                                      then if (Abs(w.a-w1.a)/Abs(w.a)<1e-18)
                                               and (Abs(w.b-w1.b)<1e-18)
                                             then finished:=True
                                             else
                                      else if (Abs(w.a-w1.a)<1e-18)
                                               and (Abs(w.b-w1.b)<1e-18)
                                             then finished:=True;
                      if finished
                        then Result:=w1
                        else begin
                               w:=w1;
                               d.a:=1;
                               d.b:=d.a;
                               km:=isub(km,d);
                               k:=k+1
                             end
                    until finished or (k>MaxInt/2);
                    if not finished
                      then st:=2
                  end
  end {ionepxtom};

  function ionemxtom (const x, m : interval;
                      out st     : Integer) : interval;
// |x|<=1, m>0
  var finished, is_even : Boolean;
      k                 : Integer;
      d, km, w, w1, w2  : interval;
  begin
    if (x.a>x.b) or (m.a>m.b)
      then st:=1
      else if (x.a<-1) or (x.b>1) or (m.a<=0)
             then st:=3
             else begin
                    d.a:=1;
                    d.b:=d.a;
                    w:=d;
                    w2:=d;
                    km:=m;
                    k:=1;
                    finished:=False;
                    is_even:=True;
                    st:=0;
                    repeat
                      d.a:=k;
                      d.b:=d.a;
                      w2:=imul(idiv(imul(km,w2),d),x);
                      if is_even
                        then w1:=isub(w,w2)
                        else w1:=iadd(w,w2);
                      if (w.a<>0) and (w.b<>0)
                        then if (Abs(w.a-w1.a)/Abs(w.a)<1e-18)
                                 and (Abs(w.b-w1.b)/Abs(w.b)<1e-18)
                               then finished:=True
                               else
                        else if (w.a=0) and (w.b<>0)
                               then if (Abs(w.a-w1.a)<1e-18)
                                        and (Abs(w.b-w1.b)/Abs(w.b)<1e-18)
                                      then finished:=True
                                      else
                               else if w.a<>0
                                      then if (Abs(w.a-w1.a)/Abs(w.a)<1e-18)
                                               and (Abs(w.b-w1.b)<1e-18)
                                             then finished:=True
                                             else
                                      else if (Abs(w.a-w1.a)<1e-18)
                                               and (Abs(w.b-w1.b)<1e-18)
                                             then finished:=True;
                      if finished
                        then Result:=w1
                        else begin
                               w:=w1;
                               d.a:=1;
                               d.b:=d.a;
                               km:=isub(km,d);
                               k:=k+1;
                               is_even:=not is_even
                             end
                    until finished or (k>MaxInt/2);
                    if not finished
                      then st:=2
                  end
  end {ionemxtom};

  function ionepxtomm (const x, m : interval;
                       out st     : Integer) : interval;
// |x|<1, m>0
  var finished, is_even : Boolean;
      k                 : Integer;
      d, km, w, w1, w2  : interval;
  begin
    if (x.a>x.b) or (m.a>m.b)
      then st:=1
      else if (x.a<=-1) or (x.b>=1) or (m.a<=0)
             then st:=3
             else begin
                    d.a:=1;
                    d.b:=d.a;
                    w:=d;
                    w2:=d;
                    km:=m;
                    k:=1;
                    finished:=False;
                    is_even:=True;
                    st:=0;
                    repeat
                      d.a:=k;
                      d.b:=d.a;
                      w2:=imul(idiv(imul(km,w2),d),x);
                      if is_even
                        then w1:=isub(w,w2)
                        else w1:=iadd(w,w2);
                      if (w.a<>0) and (w.b<>0)
                        then if (Abs(w.a-w1.a)/Abs(w.a)<1e-18)
                                 and (Abs(w.b-w1.b)/Abs(w.b)<1e-18)
                               then finished:=True
                               else
                        else if (w.a=0) and (w.b<>0)
                               then if (Abs(w.a-w1.a)<1e-18)
                                        and (Abs(w.b-w1.b)/Abs(w.b)<1e-18)
                                      then finished:=True
                                      else
                               else if w.a<>0
                                      then if (Abs(w.a-w1.a)/Abs(w.a)<1e-18)
                                               and (Abs(w.b-w1.b)<1e-18)
                                             then finished:=True
                                             else
                                      else if (Abs(w.a-w1.a)<1e-18)
                                               and (Abs(w.b-w1.b)<1e-18)
                                             then finished:=True;
                      if finished
                        then Result:=w1
                        else begin
                               w:=w1;
                               d.a:=1;
                               d.b:=d.a;
                               km:=iadd(km,d);
                               k:=k+1;
                               is_even:=not is_even
                             end
                    until finished or (k>MaxInt/2);
                    if not finished
                      then st:=2
                  end
  end {ionepxtomm};

  function ionemxtomm (const x, m : interval;
                       out st     : Integer) : interval;
// |x|<1, m>0
  var finished         : Boolean;
      k                : Integer;
      d, km, w, w1, w2 : interval;
  begin
    if (x.a>x.b) or (m.a>m.b)
      then st:=1
      else if (x.a<=-1) or (x.b>=1) or (m.a<=0)
            then st:=3
            else begin
                    d.a:=1;
                    d.b:=d.a;
                    w:=d;
                    w2:=d;
                    km:=m;
                    k:=1;
                    finished:=False;
                    st:=0;
                    repeat
                      d.a:=k;
                      d.b:=d.a;
                      w2:=imul(idiv(imul(km,w2),d),x);
                      w1:=iadd(w,w2);
                      if (w.a<>0) and (w.b<>0)
                        then if (Abs(w.a-w1.a)/Abs(w.a)<1e-18)
                                 and (Abs(w.b-w1.b)/Abs(w.b)<1e-18)
                               then finished:=True
                               else
                        else if (w.a=0) and (w.b<>0)
                               then if (Abs(w.a-w1.a)<1e-18)
                                        and (Abs(w.b-w1.b)/Abs(w.b)<1e-18)
                                      then finished:=True
                                      else
                               else if w.a<>0
                                      then if (Abs(w.a-w1.a)/Abs(w.a)<1e-18)
                                               and (Abs(w.b-w1.b)<1e-18)
                                             then finished:=True
                                             else
                                      else if (Abs(w.a-w1.a)<1e-18)
                                               and (Abs(w.b-w1.b)<1e-18)
                                             then finished:=True;
                      if finished
                        then Result:=w1
                        else begin
                               w:=w1;
                               d.a:=1;
                               d.b:=d.a;
                               km:=iadd(km,d);
                               k:=k+1
                             end
                    until finished or (k>MaxInt/2);
                    if not finished
                      then st:=2
                 end
  end {ionemxtomm};

  function isqrt2 : interval;
  var i2 : string;
  begin
    i2:='1.414213562373095048';
    Result.a:=left_read(i2);
    i2:='1.414213562373095049';
    Result.b:=right_read(i2)
  end {isqrt2};

  function isqrt3 : interval;
  var i3 : string;
  begin
    i3:='1.732050807568877293';
    Result.a:=left_read(i3);
    i3:='1.732050807568877294';
    Result.b:=right_read(i3)
  end {isqrt3};

  function isqrt5 : interval;
  var i5 : string;
  begin
    i5:='2.236067977499789696';
    Result.a:=left_read(i5);
    i5:='2.236067977499789697';
    Result.b:=right_read(i5)
  end {isqrt5};

  function isqrt6 : interval;
  var i6 : string;
  begin
    i6:='2.449489742783178098';
    Result.a:=left_read(i6);
    i6:='2.449489742783178099';
    Result.b:=right_read(i6)
  end {isqrt6};

  function isqrt7 : interval;
  var i7 : string;
  begin
    i7:='2.645751311064590590';
    Result.a:=left_read(i7);
    i7:='2.645751311064590591';
    Result.b:=right_read(i7)
  end {isqrt7};

  function isqrt8 : interval;
  var i8 : string;
  begin
    i8:='2.828427124746190097';
    Result.a:=left_read(i8);
    i8:='2.828427124746190098';
    Result.b:=right_read(i8)
  end {isqrt8};

  function isqrt10 : interval;
  var i10 : string;
  begin
    i10:='3.162277660168379331';
    Result.a:=left_read(i10);
    i10:='3.162277660168379332';
    Result.b:=right_read(i10)
  end {isqrt10};

  function ipi : interval;
  var ipistr : string;
  begin
    ipistr:='3.141592653589793238';
    Result.a:=left_read(ipistr);
    ipistr:='3.141592653589793239';
    Result.b:=right_read(ipistr)
  end {ipi};

{$IFDEF WIN64}
initialization
  ShowMessage ('Although on Win64 environment all internal calculations will '
               +'be executed using the TExtendedX87 type (a replacement type '
               +'for Win32''s Extended on Win64), do not enter any data that '
               +'exceed the Double type range.');
{$ENDIF}
end.