program interval(input,output);{dataflow simulator for real interval arithmetic}const	IMem = 500; {number of instructions}	DMem = 200;  {number of interval variables allowed}	Par = 3;   {max number of paramters for an operator}	Maxexp = 10; Minexp = -9;	Maxinf =10000; Mininf =-10000; {Tied to values of Minman and Maxman}	Maxman = 9999; Minman = -9999; {cannot exceed sqrt(maxint)}	Splitman = 1000; {Smallest normalized mantisa}	Digits = 4; {number of digits in mantissa}type	Positive= 0..maxint;	State  = -1..99;  {Used for holding state of operator -1:done}	OpType = (print,pr,tr,soln,readr,halve,halves,linh,mult,add,intgr,less,leq,noteq,sqrr,minr,maxr,modu,absr,trig,expr,lb,ub,copy,stop); {!!}	Ptr    = 1..DMem;	Loc    = 1..IMem;	Loc0   = 0..IMem;	EdgeT  = (hout,lin,hin,lout); {Warning this order is important in}				      {predicates such as gtS,geS}	CardT  = (finite,infinite);	ExpT   = Minexp..Maxexp;	ManT   = Mininf..Maxinf; 	Pflag  = (PNull,PSoln,PTrace,PPrint);	Sreal  = record		    edge:EdgeT;		    cardinality:CardT;		    exp:ExpT; {exponent}		    mantissa:ManT;		 end;	Int    = record		    hi:Sreal;		    lo:Sreal;	 end;	Instr  = record		    Code:OpType;		    Pars: array[0..Par] of 0..DMem;		 end;	DataMem= record		    D        :array [Ptr] of Int;		    S        :array [Loc] of State;		    LastHalve:Loc;		    RHalve   :array [Loc] of real;		 end;	DataFlags=record		    PF	     :array [Ptr] of Pflag;		 end;var	Debug  : (none,activity,post,trace,dump);	Cut    : (once,all);	GlobalEnd,Verifiable:boolean;	HalveThreshold:real;	I      : array [Loc] of Instr; {Memory holding instructions}	End    : Loc; {last instruction in I}	ParN   : array [OpType] of -1..Par; {number of parameters for each 			opcode. -1 means no result}        ParIntersect : array [OpType] of boolean ;	DInit  : DataMem; {initial memory which is cleared and 				used in first call}	DF     : DataFlags; {hold flags for variables, e.g. print/trace}	MaxDMem:0..DMem;	Shift  : array[0..Digits] of 1..maxint;{array of constant multipliers}						{used for alignment etc.}	Dummy  :Positive;	{constant intervals and Sreals}	PlusInfS,MinusInfS,PlusSmallS,MinusSmallS,ZeroS,	PlusFiniteS,MinusFiniteS:Sreal;	Zero,All,AllFinite:Int;procedure deblank;var Ch:char;begin   while (not eof) and (input^ in [' ','	']) do read(Ch);end;procedure InitialOptions;#include '/user/profs/cleary/bin/options.i';   procedure Option;   begin      case Opt of      'a','A':Debug:=activity;      'd','D':Debug:=dump;      'h','H':HalveThreshold:=StringNum/100;      'n','N':Debug:=none;      'p','P':Debug:=post;      't','T':Debug:=trace;      'v','V':Verifiable:=true;      end;   end;begin   Debug:=trace;   Verifiable:=false;   HalveThreshold:=67/100;   Options;   writeln(Debug);   writeln('Verifiable:',Verifiable);   writeln('Halve threshold',HalveThreshold);end;{InitialOptions}procedure NormalizeUp(E,M:integer;var S:Sreal;var Closed:boolean);beginwith S dobegin   if M=0 then S:=ZeroS else   if M>0 then   begin      while M>=Maxinf do      begin 	 if M mod 10 > 0 then begin Closed:=false;M:=(M div 10)+1 end	 else M:=M div 10;	 E:=E+1;      end;	       while M < Maxinf div 10 do      begin M:=M*10; E:=E-1;       end;      if E > Maxexp then {overflow-set to infinity}      begin 	 S:=PlusInfS;	 Closed:=false;      end else      if E < Minexp then {underflow-set to smallest positive value}      begin          S:=PlusSmallS;         Closed:=false;      end else      begin cardinality:=finite;exp:=E;mantissa:=M;      end;   end else	    if M < 0 then   begin      while M <= Mininf do      begin 	 if M mod 10 < 0 then Closed:=false else	 if M mod 10 > 0 then halt;	 M:=M div 10;	 E:=E+1;      end;	       while M > (Mininf div 10) do      begin M:=M*10; E:=E-1;       end;      if E > Maxexp then {overflow-set to most negative value}      begin          S:=MinusFiniteS;         Closed:=false;      end       else      if E < Minexp then {underflow-set to zero}      begin         S:=ZeroS;         Closed:=false;      end else      begin         cardinality:=finite;exp:=E;mantissa:=M;      end;   end;end;end;{NormalizeUp}procedure NormalizeDn(E,M:integer;var S:Sreal;var Closed:boolean);beginwith S dobegin   if M=0 then S:=ZeroS else   if M>0 then   begin      while M >= Maxinf do      begin 	 if M mod 10 > 0 then Closed:=false else	 if M mod 10 < 0 then halt;	 M:=M div 10;	 E:=E+1;      end;	       while (M < Maxinf div 10) do      begin M:=M*10; E:=E-1;       end;      if E > Maxexp then {overflow-set to largest positive value}      begin 	 S:=PlusFiniteS;	 Closed:=false;      end else      if E < Minexp then {underflow-set to zero}      begin S:=ZeroS; Closed:=false;      end else      begin cardinality:=finite;exp:=E;mantissa:=M;      end;   end else	    if M < 0 then   begin      while M <= Mininf do      begin 	 if M mod 10 < 0 then 	 begin Closed:=false; M:=M div 10 -1;end	 else 	 if M mod 10 = 0 then M:=M div 10 	 else halt;	 E:=E+1;      end;	       while (M>Mininf div 10) do      begin M:=M*10; E:=E-1;       end;      if E > Maxexp then {overflow}      begin          S:=MinusInfS;         Closed:=false;      end       else      if E < Minexp then {underflow}      begin S:=MinusSmallS; Closed:=false;      end else      begin         cardinality:=finite;exp:=E;mantissa:=M;      end;   end;end;end;{NormalizeDn}procedure WriteS(X:Sreal);var E,M:integer;beginwith X dobegin   case edge of   lin: write('[');   lout: write('(');   hin,hout:   end;      case cardinality of   infinite: write('inf':Digits+4);    finite:       if mantissa = 0 then write(0:Digits+1,' ':3)      else begin         M:=mantissa;	 E:=exp;          while (M mod 10 = 0) do	 begin M:=M div 10; E:=E+1;	 end;         write(M:Digits+1,'e',E-Digits:2);      end;   end;      case edge of    hin: write(']');   hout:write(')');   lin,lout:   end;end;end;{WriteS}procedure WriteInt(I:Int);begin   with I do begin WriteS(lo); write(','); WriteS(hi); end;end;{WriteInt}   procedure DumpS(X:Sreal);beginwith X do   write(edge:4,cardinality:9,mantissa:7,exp:3);end;{DumpS}procedure DumpInt(I:Int);begin   with I do begin DumpS(lo); write(' || '); DumpS(hi); end;end;{DumpInt}   procedure ReadInt(var I:Int);var   Ch:char;      Cll,Clu:boolean;	   procedure ReadSUp(var X:Sreal; var Closed:boolean);   var E,M:integer;   begin      with X do      begin         deblank;         case input^ of         '~':begin X:=PlusInfS;read(Ch);	     end;         '-','+','0','1','2','3','4','5','6','7','8','9':	 begin	    cardinality:=finite;   	    read(M);	    read(E); E:=E+Digits;	    NormalizeUp(E,M,X,Closed);	 end;	 end;{case}      end;   end;{ReadSUp}   procedure ReadSDn(var X:Sreal; var Closed:boolean);   var E,M:integer;       Ch:char;   begin      with X do      begin         deblank;         case input^ of         '~':begin X:=MinusInfS;read(Ch);	     end;         '-','+','0','1','2','3','4','5','6','7','8','9':	 begin	    cardinality:=finite;   	    read(M);	    read(E); E:=E+Digits;	    NormalizeDn(E,M,X,Closed);	 end;	 end;{case}      end;   end;{ReadSDn}begin{ReadInt}   with I do    begin       deblank; read(Ch);       case Ch of      '[':Cll:=true;      '(':Cll:=false;      end;      ReadSDn(lo,Cll);if Cll then lo.edge:=lin else lo.edge:=lout;      deblank;      read(Ch); assert(Ch=',');      Clu:=true;      ReadSUp(hi,Clu);      deblank;      read(Ch);      case Ch of      ']':if Clu then hi.edge:=hin else hi.edge:=hout;      ')':hi.edge:=hout;      end;   end;end;{ReadInt}   procedure DumpTables;var tL:Loc; tPar:0..Par; tOp:OpType;begin	for tOp := print to stop do	   writeln(tOp:6,ParN[tOp]:2);	writeln;	for tL := 1 to End do	with I[tL] do	begin	   write(Code:5);	   for tPar := 0 to Par do	      if Pars[tPar] <> 0 then write(Pars[tPar]:4);	   writeln;	end;	writeln('number of memory locations used:',MaxDMem:0);	writeln;end;{DumpTables}	procedure AlignUp   (E0:ExpT;M0:ManT;E1:ExpT;M1:ManT;var E,N0,N1:integer;var Closed:boolean);{Align mantissas M0,M1 preserving accuracy and rounding up wherever possible}{common resulting exponents in E, and mantissas in N0,N1}var D:Positive;begin   if M0=0 then begin E:=E1;N0:=0;N1:=M1;end else   if M1=0 then begin E:=E0;N0:=M0;N1:=0;end else   if E0=E1 then   begin E:=E0; N0:=M0; N1:=M1;   end else   if (E0>E1) then AlignUp(E1,M1,E0,M0,E,N1,N0,Closed) else   begin      D:=E1-E0;      if D>= 2*Digits then      begin          N1:=M1*Maxinf; E:=E1-Digits;	 if M0<0 then N0:=0 else N0:=1;	 Closed:=false;      end else      if D > Digits then      begin          N1:=M1*Maxinf; E:=E1-Digits; 	 if (M0 mod Shift[D-Digits]) = 0 	 then N0:=(M0 div Shift[D-Digits])	 else	    if M0 > 0 then N0:=(M0 div Shift[D-Digits])+1	              else N0:=(M0 div Shift[D-Digits]);      end else      {Digits>=D>=0}      begin N1:=M1*Shift[D]; E:=E1-D; N0:=M0;      end;   end;end;{AlignUp}function gtS(X,Y:Sreal):boolean;{X>Y  careful need to be able to compare x] and (x etc.}var gt:boolean;begin   if (X.exp=Y.exp)and(X.mantissa=Y.mantissa) then gt:=X.edge>Y.edge else   if X.exp = Y.exp then gt:= (X.mantissa > Y.mantissa) else   if X.mantissa = 0 then gt:= 0 > Y.mantissa else   if Y.mantissa = 0 then gt:= X.mantissa > 0 else   if (X.mantissa>0) and (Y.mantissa>0) then gt:= (X.exp > Y.exp) else   if (X.mantissa>0) and (Y.mantissa<0) then gt:= true else   if (X.mantissa<0) and (Y.mantissa>0) then gt:= false else   if (X.mantissa<0) and (Y.mantissa<0) then gt:= (X.exp < Y.exp)    else  writeln('error in gtS');      gtS:=gt;end;{gtS}   function geS(X,Y:Sreal):boolean;{X>=Y  careful need to be able to compare x] and (x etc.}begin   if (X.exp=Y.exp)and(X.mantissa=Y.mantissa) then geS:=X.edge>=Y.edge else   if X.exp = Y.exp then geS:= (X.mantissa >= Y.mantissa) else   if X.mantissa = 0 then geS:= 0 >= Y.mantissa else   if Y.mantissa = 0 then geS:= X.mantissa >= 0 else   if (X.mantissa>0) and (Y.mantissa>0) then geS:= (X.exp > Y.exp) else   if (X.mantissa>0) and (Y.mantissa<0) then geS:= true else   if (X.mantissa<0) and (Y.mantissa>0) then geS:= false else   if (X.mantissa<0) and (Y.mantissa<0) then geS:= (X.exp < Y.exp)    else  writeln('error in geS');end;{geS}   function Point(X:Int):boolean;{X=[x,x]}beginwith X do   Point:=(lo.edge=lin)and (hi.edge=hin) and 	  (lo.mantissa=hi.mantissa) and	  (lo.exp=hi.exp);end;{Point}procedure maxS(X,Y:Sreal;var max:Sreal);begin	if gtS(X,Y) then max:=X else max:=Y;end;procedure minS(X,Y:Sreal;var min:Sreal);begin	if gtS(X,Y) then min:=Y else min:=X;end;procedure Inter(P,Q:Int;var R:Int);begin   minS(P.hi,Q.hi,R.hi);   maxS(P.lo,Q.lo,R.lo);end;function CheckHi(X:Sreal):boolean;var OK:boolean;begin   OK:=true;   with X do   begin      case cardinality of      infinite:         if (exp=Maxexp)and(mantissa=Maxinf) then	 else writeln('**Invalid hi infinity');      finite:      begin         if (mantissa=Maxinf) or (mantissa=Mininf) then	 begin OK:=false; writeln('**Invalid finite value - hi');	 end;	          if mantissa = 0 then	    if (exp=0) then 	    else 	    begin OK:=false; writeln('**Invalid zero - hi')	    end	 else	 begin	    if (mantissa > 0) then	       if mantissa >= (Maxinf div 10) then {OK}	       else 	       begin OK:=false; writeln('**Incorrect normalization - hi') 	       end	    else{mantissa<0}	       if mantissa > (Mininf div 10) then	       begin OK:=false; writeln('**Incorrect normalization - hi') 	       end;	 end;      end;      end;{case}