program DirectedIntervalPoisson4thOrder_Example_1;
{$APPTYPE CONSOLE}
// Tis program differs from the IntervalPoisson4thOrder_Example_1 only by the
// directedfullintervalPoisson procedure, where the directed interval arithmetic
// is applied to solve the Poisson equation, and by using dinterval type instead
// of interval type
uses SysUtils, DateUtils, IntervalArithmetic32and64;
type intervalmatrix        = array of array of dinterval;
     intervalrightfunction = function (const X, Y : dinterval;
                                       out st     : Integer) : dinterval;
     intervalboundary      = function (const X : dinterval;
                                       out st  : Integer) : dinterval;

procedure directedfullintervalPoisson (const n, m        : Integer;
                                       const ALPHA, BETA : dinterval;
                                       const eps         : Extended;
                                       const F           : intervalrightfunction;
                                       const PHI1, PHI2,
                                             PHI3, PHI4  : intervalboundary;
                                       const THETA, KSI  : intervalrightfunction;
                                       const Pconst,
                                             Qconst      : Extended;
                                       out U             : intervalmatrix;
                                       out st            : Integer);
{------------------------------------------------------------------------------}
{                                                                              }
{  The procedure fullintervalPoisson solves the second order partial           }
{  differential equation u|xx + u|yy = f(x,y), where u=u(x,y), u|xx denotes    }
{  the second order partial derivative of u with respect to x,                 }
{  0 <= x <= alpha, 0 <= y <= beta, with the boundary conditions u = phi1(y)   }
{  for x = 0, u = phi2(x) for y = 0, u = phi3(y) for x = alpha, u = phi4(x)    }
{  for y = beta, where phi1(0) = phi2(0), phi2(alpha) = phi3(0),               }
{  phi3(beta) = phi4(alpha), phi4(0) = phi1(beta).                             }
{                                                                              }
{  The problem is solved by an interval difference method of fourth order      }
{  using directed interval arithmetic.                                         }
{                                                                              }
{  Data:                                                                       }
{    n           - the number of subintervals in the interval [0,alpha],       }
{    m           - the number of subintervals in the interval [0,beta],        }
{    ALPHA       - an interval representation of maximum value of x,           }
{    BETA        - an interval representation of maximum value of y,           }
{    eps         - the maximum difference for which the values of PHI1([0,0])  }
{                  and PHI2(0), PHI2(ALPHA) and PHI3([0,0]), PHI3(BETA) and    }
{                  PHI4(ALPHA), PHI4([0,0]) and PHI1(BETA), where PHI1, PHI2,  }
{                  PHI3 and PHI4 are interval extensions of phi1, phi2, phi3   }
{                  and phi4, respectively, can be considered as equal, i.e.    }
{                  the following inequalities are fulfilled:                   }
{                    Abs(PHI([0,0]).z - PHI2([0,0].z))/Abs(PHI1([0,0]).z) < eps}
{                               for PHI1([0,0]).z <> 0 and PHI2([0,0]).z <> 0, }
{                    Abs(PHI2([0,0].z)) < eps for PHI1([0,0]) = 0,             }
{                    Abs(PHI1([0,0].z)) < eps for PHI2([0,0]) = 0,             }
{                    Abs(PHI2(ALPHA).z-PHI3([0,0].z))/Abs(PHI2(ALPHA).z) < eps }
{                               for PHI2(ALPHA).z <> 0 and PHI3([0,0]).z <> 0, }
{                    Abs(PHI3([0,0])).z < eps for PHI2(ALPHA).z = 0,           }
{                    Abs(PHI2(ALPHA)).z < eps for PHI3([0,0]).z = 0,           }
{                    Abs(PHI3(BETA).z-PHI4(ALPHA).z)/Abs(PHI3(BETA).z)         }
{                               for PHI3(BETA).z <> 0 and PHI4(ALPHA).z <> 0,  }
{                    Abs(PHI4(ALPHA).z) < eps for PHI3(BETA).z = 0,            }
{                    Abs(PHI3(BETA).z) < eps for PHI4(ALPHA).z = 0,            }
{                    Abs(PHI4([0,0].z)-PHI1(BETA).z)/Abs(PHI([0,0].z)) < eps   }
{                               for PHI4([0,0]).z <> 0 and PHI1(BETA).z <> 0,  }
{                    Abs(PHI1(BETA).z) < eps for PHI4([0,0]).z = 0,            }
{                    Abs(PHI4([0,0]).z) < eps for PHI1(BETA).z = 0,            }
{                  where z denotes either left or right end of interval),      }
{    F           - a Delphi Pascal function which calculates the value of      }
{                  interval extension of function f = f(x,y) occurring in the  }
{                  differential equation,                                      }
{    PHI1, PHI2,                                                               }
{    PHI3, PHI4  - Delphi Pascal functions which calculate the values of       }
{                  interval extensions of functions phi1 = phi1(y),            }
{                  phi2 = phi2(x), phi3 = phi3(y) and phi4 = phi4(x) occurring }
{                  in the boundary conditions,                                 }
{    THETA       - a Delphi Pascal function which calculates the value of      }
{                  interval extension of function f|xxxx (the fourth order     }
{                  partial derivative of f with respect to x),                 }
{    KSI         - a Delphi Pascal function which calculates the value of      }
{                  interval extension of function f|yyyy (the fourth order     }
{                  partial derivative of f with respect to y),                 }
{    Pconst      - the value such that Abs(u|xxxxyy) <= Pconst for each        }
{                  0 <= x <= alpha and 0 <= y <= beta (u|xxxxyy denotes the    }
{                  sixth order mixed partial derivative of u with respect to   }
{                  x and y),                                                   }
{    Qconst      - the value such that Abs(u|xxyyyy) <= Qconst for each        }
{                  0 <= x <= alpha and 0 <= y <= beta (u|xxyyyy denotes the    }
{                  sixth order mixed partial derivative of u with respect to   }
{                  x and y).                                                   }
{  Result:                                                                     }
{    U - an interval array containing the solution (the element                }
{        [u[i,j].left, u[i,j].right] contains the solution at the point        }
{        (ih, jk), where i = 0, 1, ... , n, j = 0, 1, ... , m, h = alpha/n,    }
{        k = beta/m.                                                           }
{  Other parameters:                                                           }
{    st   - a variable which within the procedure diffmethod is assigned the   }
{           value of:                                                          }
{             1, if n < 2 or m < 2,                                            }
{             2, if ALPHA.left <= 0 or BETA.left <= 0,                         }
{             3, if one of standard interval functions cannot be evaluated     }
{             4, if at least one of the conditions given above in the          }
{                description of eps is not fulfilled (that means the boundary  }
{                conditions are wrong),                                        }
{             5, if the matrix of adequate linear system is singular,          }
{             0, otherwise.                                                    }
{           Note: If st<>0, then the elements of array u are not calculated.   }
{  Unlocal identifiers:                                                        }
{    intervalmatrix        - a type identifier of interval dynamical array     }
{                            defined as follows:                               }
{                              type intervalmatrix                             }
{                                   = array of array of interval;              }
{                            which corresponds to the array [0..n, 0..m],      }
{    intervalrightfunction - a procedural-type identifier defined as follows   }
{                              type intervalrightfunction                      }
{                                   = function (const X, Y : interval)         }
{                                               : interval;                    }
{    intervalboundary    - a procedural-type identifier defined as follows     }
{                            type intervalboundary                             }
{                                 = function (const X : interval) : interval;  }
{  Note:                                                                       }
{    Within the procedure it is assigned a memory for the U parameter          }
{    (dynamical array). After using the procedure you have to release this     }
{    memory by calling the procedure Finalize(X), where X denotes your         }
{    argument corresponding to the U parameter.                                }
{                                                                              }
{------------------------------------------------------------------------------}
const izero   : dinterval = (a:0;  b:0);
      ione    : dinterval = (a:1;  b:1);
      itwo    : dinterval = (a:2;  b:2);
      ifive   : dinterval = (a:5;  b:5);
      ieight  : dinterval = (a:8;  b:8);
      itwelve : dinterval = (a:12; b:12);
      itwenty : dinterval = (a:20; b:20);
var i, j, jh, j1, k, kh, l, lh, l1, l2, n1, n2, p, q, rh : Integer;
    z                                                    : Extended;
    AF, BF, CF, H1, H2, H4, HH, iHH, II, iminus1HH,
    iplus1HH, JJ, jKK, jminus1KK, jplus1KK, K1, K2, K4,
    KK, MAX, MM, NN, PPconst, PQ, QQconst, S, S1, S2     : dinterval;
    A1, B1, X                                            : array of dinterval;
    r                                                    : array of Integer;
function boundconds (const B1, B2 : dinterval;
                     const eps    : Extended) : Integer;
begin
  Result:=0;
  if (B1.a<>0) and (B2.a<>0)
    then if Abs(B1.a-B2.a)/Abs(B1.a)>=eps
           then Result:=4
           else
    else if B1.a=0
            then if Abs(B2.a)>=eps
                   then Result:=4
                   else
            else if B2.a=0
                   then if Abs(B1.a)>=eps
                          then Result:=4;
  if (B1.b<>0) and (B2.b<>0)
    then if Abs(B1.b-B2.b)/Abs(B1.b)>=eps
           then Result:=4
           else
    else if B1.b=0
            then if Abs(B2.b)>=eps
                   then Result:=4
                   else
            else if B2.b=0
                   then if Abs(B1.b)>=eps
                          then Result:=4;
end {boundconds};
begin
  st:=0;
  NN.a:=n;
  NN.b:=n;
  MM.a:=m;
  MM.b:=m;
  PPconst.a:=Pconst;
  PPconst.b:=-Pconst;
  QQconst.a:=Qconst;
  QQconst.b:=-Qconst;
  if (n<2) or (m<2)
    then st:=1
    else if (ALPHA.a<=0) or (BETA.a<=0)
           then st:=2;
  if st=0
    then
      begin
        S:=PHI1(izero,st);
        if st=0
          then
            begin
              S1:=PHI2(izero,st);
              if st=0
                then
                  begin
                    st:=boundconds(S,S1,eps);
                    if st=0
                      then
                        begin
                          S:=PHI2(ALPHA,st);
                          if st=0
                            then
                              begin
                                S1:=PHI3(izero,st);
                                if st=0
                                  then
                                    begin
                                      st:=boundconds(S,S1,eps);
                                      if st=0
                                        then
                                          begin
                                            S:=PHI3(BETA,st);
                                            if st=0
                                              then
                                                begin
                                                  S1:=PHI4(ALPHA,st);
                                                  if st=0
                                                    then
{<--- 20}
                                  begin
                                    st:=boundconds(S,S1,eps);
                                    if st=0
                                      then
                                        begin
                                          S:=PHI4(izero,st);
                                          if st=0
                                            then
                                              begin
                                                S1:=PHI1(BETA,st);
                                                if st=0
                                                  then st:=boundconds(S,S1,eps)
                                              end
                                        end
                                  end
{---> 20}
                                                end
                                          end
                                    end
                              end
                        end
                  end
            end
      end;
  if st=0
    then begin
           H1:=didiv(ALPHA,NN);  // interval h
           K1:=didiv(BETA,MM);   // interval k
           H2:=dimul(H1,H1);     // interval h^2
           K2:=dimul(K1,K1);     // interval k^2
           H4:=dimul(H2,H2);     // interval h^4
           K4:=dimul(K2,K2);     // interval k^4
           HH.a:=-H1.b;
           HH.b:=H1.a;          // HH = [-h,h]
           KK.a:=-K1.b;
           KK.b:=K1.a;          // KK = [-k,k]
           AF:=diadd(H2,K2);     // interval h^2+k^2
           BF:=disub(dimul(ifive,H2),K2);
           BF:=dimul(itwo,BF);   // interval 2*(5*h^2-k^2)
           CF:=disub(dimul(ifive,K2),H2);
           CF:=dimul(itwo,CF);   // interval 2*(5*k^2-h^2);
           PQ:=dimul(H4,PPconst);
           PQ:=diadd(PQ,dimul(K4,QQconst));
           PQ:=didiv(PQ,itwenty);
           PPconst:=diadd(PPconst,QQconst);
           PPconst:=dimul(dimul(H2,K2),PPconst);
           PPconst:=didiv(PPconst,itwelve);
           PQ:=diadd(PQ,PPconst);
                                // PQ = (h^4*[-P,P]+k^4*[-Q,Q])/20
                                //       +h^2*k^2*([-P,P]+{-Q,Q])/12
           n1:=(n-1)*(m-1);
           n2:=n1+1;
           p:=n2;
           SetLength (r,n1+1);
           for i:=1 to n2 do
             r[i-1]:=0;
           k:=0;
           j:=0;
           SetLength (A1,n1+1);
           SetLength (B1,n1+1);
           SetLength (X,(n*m-n-m+4)*(n*m-n-m+4) div 4);
           repeat
             k:=k+1;
             for i:=1 to n1 do
               A1[i-1]:=izero;
             j:=j+1;
             i:=Trunc((k-1)/(m-1))+1;
             l1:=(i-2)*(m-1)+j;
             l2:=l1+m-1;
             II.a:=i;
             II.b:=i;
             JJ.a:=j;
             JJ.b:=j;
             iHH:=dimul(II,H1);  // interval ih
             jKK:=dimul(JJ,K1);  // interval jk
             iminus1HH:=dimul(disub(II,ione),H1);
             iplus1HH:=dimul(diadd(II,ione),H1);
             jminus1KK:=dimul(disub(JJ,ione),K1);
             jplus1KK:=dimul(diadd(JJ,ione),K1);
             if i>1
               then begin
                      A1[l1-1]:=CF;
                      if j>1
                        then A1[l1-2]:=AF;
                      if j<m-1
                        then A1[l1]:=AF
                    end;
             S:=dimul(itwenty,AF);
             A1[l2-1].a:=-S.b;
             A1[l2-1].b:=-S.a;
             if j>1
               then A1[l2-2]:=BF;
             if j<m-1
               then A1[l2]:=BF;
             l1:=l2+m-1;
             if i<n-1
               then begin
                      A1[l1-1]:=CF;
                      if j>1
                        then A1[l1-2]:=AF;
                      if j<m-1
                        then A1[l1]:=AF
                    end;
             S:=F(iminus1HH,jKK,st);
             if st=0
               then begin
                      S1:=F(iplus1HH,jKK,st);
                      if st=0
                        then begin
                               S:=diadd(S,S1);
                               S1:=F(iHH,jKK,st);
                               if st=0
                                 then begin
                                        S:=diadd(S,dimul(ieight,S1));
                                        S1:=F(iHH,jminus1KK,st);
                                        if st=0
                                          then begin
                                                 S:=diadd(S,S1);
                                                 S1:=F(iHH,jplus1KK,st);
                                                 if st=0
                                                   then S:=diadd(S,S1);
                                               end
                                      end
                             end
                    end;  // S = F[i-1,j]+F[i+1,j]+8*F[i,j]+F[i,j-1]+F[i,j+1]
             if st=0
               then begin
                      S1:=THETA(diadd(iHH,HH),jKK,st);
                      if st=0
                        then begin
                               S1:=dimul(H4,S1);
                               S2:=KSI(iHH,diadd(jKK,KK),st);
                               if st=0
                                 then begin
                                        S1:=diadd(S1,dimul(K4,S2));
                                        S1:=didiv(S1,itwenty)
                                      end
                             end
                    end;  // S1 = (h^4*THETA(X[i]+[-h,h],Y[j])
                          //       +k^4*KSI(X[i],Y[j]+[-k,k]))/20
             if st=0
               then begin
                      S:=diadd(disub(S,S1),PQ);
                      S:=dimul(dimul(H2,K2),S); // S = LAMBDA[i,j]
                      if i=1
                        then begin
                               S1:=PHI1(jminus1KK,st);
                               if st=0
                                 then
                                   begin
                                     S2:=PHI1(jplus1KK,st);
                                     if st=0
                                       then begin
                                              S1:=dimul(AF,diadd(S1,S2));
                                              S2:=PHI1(jKK,st);
                                              if st=0
                                                then S1:=diadd(S1,dimul(CF,S2))
                                            end;
                                     S:=disub(S,S1);
                                     if j=1
                                       then
                                         begin
                                           S1:=PHI2(iplus1HH,st);
                                           if st=0
                                             then
                                               begin
                                                 S:=disub(S,dimul(AF,S1));
                                                 S1:=PHI2(iHH,st);
                                                 if st=0
                                                   then S:=disub(S,dimul(BF,S1))
                                               end
                                         end;
                                     if j=m-1
                                       then
                                         begin
                                           S1:=PHI4(iplus1HH,st);
                                           if st=0
                                             then
                                               begin
                                                 S:=disub(S,dimul(AF,S1));
                                                 S1:=PHI4(iHH,st);
                                                 if st=0
                                                   then S:=disub(S,dimul(BF,S1))
                                               end
                                         end
                                   end
                             end
                        else if i=n-1
                               then
                                 begin
                                   S1:=PHI3(jminus1KK,st);
                                   if st=0
                                     then
                                       begin
                                         S2:=PHI3(jplus1KK,st);
                                         if st=0
                                           then
                                             begin
                                               S1:=dimul(AF,diadd(S1,S2));
                                               S2:=PHI3(jKK,st);
                                               if st=0
                                                 then S1:=diadd(S1,dimul(CF,S2))
                                             end
                                       end;
                                    S:=disub(S,S1);
                                    if j=1
                                      then
                                        begin
                                          S1:=PHI2(iminus1HH,st);
                                          if st=0
                                            then
                                              begin
                                                S:=disub(S,dimul(AF,S1));
                                                S1:=PHI2(iHH,st);
                                                if st=0
                                                  then S:=disub(S,dimul(BF,S1))
                                              end;
                                        end;
                                    if j=m-1
                                      then
                                        begin
                                          S1:=PHI4(iminus1HH,st);
                                          if st=0
                                            then
                                              begin
                                                S:=disub(S,dimul(AF,S1));
                                                S1:=PHI4(iHH,st);
                                                if st=0
                                                  then S:=disub(S,dimul(BF,S1))
                                              end
                                        end
                                 end
                               else
                                 begin
                                   if j=1
                                     then
                                       begin
                                         S1:=PHI2(iminus1HH,st);
                                         if st=0
                                           then
                                             begin
                                               S2:=PHI2(iplus1HH,st);
                                               if st=0
                                                 then
                                                   begin
                                                     S1:=diadd(S1,S2);
                                                     S:=disub(S,dimul(AF,S1));
                                                     S1:=PHI2(iHH,st);
                                                     if st=0
                                                       then S:=disub(S,
                                                                   dimul(BF,S1))
                                                   end
                                             end
                                       end;
                                   if j=m-1
                                     then
                                       begin
                                         S1:=PHI4(iminus1HH,st);
                                         if st=0
                                           then
                                             begin
                                               S2:=PHI4(iplus1HH,st);
                                               if st=0
                                                 then
                                                   begin
                                                     S1:=diadd(S1,S2);
                                                     S:=disub(S,dimul(AF,S1));
                                                     S1:=PHI4(iHH,st);
                                                     if st=0
                                                       then S:=disub(S,
                                                                   dimul(BF,S1))
                                                   end
                                             end
                                       end
                                 end
                    end;
             if st=0
               then begin
                      A1[n2-1]:=S;
                      for i:=1 to n1 do
                        begin
                          rh:=r[i-1];
                          if rh<>0
                            then B1[rh-1]:=A1[i-1]
                        end;
                      kh:=k-1;
                      l:=0;
                      MAX.a:=0;
                      MAX.b:=0;
                      for j1:=1 to n2 do
                        if r[j1-1]=0
                          then begin
                                 S:=A1[j1-1];
                                 l:=l+1;
                                 q:=l;
                                 for i:=1 to kh do
                                   begin
                                     S:=disub(S,dimul(B1[i-1],X[q-1]));
                                     q:=q+p
                                   end;
                                 A1[l-1]:=S;
                                 if S.a<0
                                   then S.a:=Abs(S.a);
                                 if S.b<0
                                   then S.b:=Abs(S.b);
                                 if S.b<S.a
                                   then begin
                                          z:=S.a;
                                          S.a:=S.b;
                                          S.b:=z
                                        end;
                                 if (j1<n2) and (S.b>MAX.a)
                                   then begin
                                          MAX:=S;
                                          jh:=j1;
                                          lh:=l
                                        end
                               end;
                      if (MAX.a=0) and (MAX.b=0)
                        then st:=5
                        else begin
                               MAX:=didiv(ione,A1[lh-1]);
                               r[jh-1]:=k;
                               for i:=1 to p do
                                 A1[i-1]:=dimul(MAX,A1[i-1]);
                               jh:=0;
                               q:=0;
                               for j1:=1 to kh do
                                 begin
                                   S:=X[q+lh-1];
                                   for i:=1 to p do
                                   if i<>lh
                                     then begin
                                            jh:=jh+1;
                                            X[jh-1]:=disub(X[q+i-1],
                                                           dimul(S,A1[i-1]))
                                          end;
                                   q:=q+p
                                 end;
                               for i:=1 to p do
                                 if i<>lh
                                   then begin
                                          jh:=jh+1;
                                          X[jh-1]:=A1[i-1]
                                        end;
                               p:=p-1
                             end;
                      if j=m-1
                        then j:=0
                    end
           until (k=n1) or (st<>0);
           Finalize (A1);
           Finalize (B1);
           if st=0
             then begin
                    SetLength (U,n+1);
                    for i:=0 to n do
                      SetLength (U[i],m+1);
                    for k:=1 to n1 do
                      begin
                        rh:=r[k-1];
                        if rh<>k
                          then begin
                                 S:=X[k-1];
                                 X[k-1]:=X[rh-1];
                                 i:=r[rh-1];
                                 while i<>k do
                                   begin
                                     X[rh-1]:=X[i-1];
                                     r[rh-1]:=rh;
                                     rh:=i;
                                     i:=r[rh-1]
                                   end;
                                 X[rh-1]:=S;
                                 r[rh-1]:=rh
                               end
                      end;
                    for i:=1 to n-1 do
                      for j:=1 to m-1 do
                        U[i,j]:=X[(i-1)*(m-1)+j-1];
                    i:=0;
                    repeat
                      i:=i+1;
                      II.a:=i;
                      II.b:=i;
                      iHH:=dimul(II,H1);
                      S:=PHI2(iHH,st);
                      if st=0
                        then begin
                               U[i,0]:=S;
                               S:=PHI4(iHH,st);
                               if st=0
                                 then U[i,m]:=S
                             end
                    until (i=n-1) or (st<>0);
                    if st=0
                      then begin
                             j:=-1;
                             repeat
                               j:=j+1;
                               JJ.a:=j;
                               JJ.b:=j;
                               jKK:=dimul(JJ,K1);
                               S:=PHI1(jKK,st);
                               if st=0
                                 then begin
                                        U[0,j]:=S;
                                        S:=PHI3(jKK,st);
                                        if st=0
                                          then U[n,j]:=S
                                      end
                             until (j=m) or (st<>0)
                           end
                  end;
           Finalize (X);
           Finalize (r)
         end
end {directedfullintervalPoisson};

//  The Delphi Pascal interval functions for the equation u|xx + u|yy = 0 with
//  the boundary conditions phi1(y) = cos(3y) for x = 0, phi2(x) = exp(3x) for
//  y = 0, phi3(y) = exp(3)cos(3y) for x = 1, phi4(x) = exp(3x)cos(3) for y = 1
//  (alpha = beta = 1).

function intf (const ix, iy : dinterval;
               out st       : Integer) : dinterval;
begin
  Result.a:=0;
  Result.b:=0;
  st:=0
end {intf};

function intphi1 (const iy : dinterval;
                  out st   : Integer) : dinterval;
const ithree : interval = (a:3; b:3);
var status     : Integer;
    cos3y, diy : interval;
begin
  diy.a:=iy.a;
  diy.b:=iy.b;
  cos3y:=icos(imul(ithree,diy),status);
  if status=0
    then begin
           Result.a:=cos3y.a;
           Result.b:=cos3y.b;
           st:=0
         end
    else begin
           Result.a:=0;
           Result.b:=0;
           st:=3
         end
end {intphi1};

function intphi2 (const ix : dinterval;
                  out st   : Integer) : dinterval;
const ithree : interval = (a:3; b:3);
var status     : Integer;
    exp3x, pix : interval;
begin
  pix.a:=ix.a;
  pix.b:=ix.b;
  exp3x:=iexp(imul(ithree,pix),status);
  if status=0
    then begin
           Result.a:=exp3x.a;
           Result.b:=exp3x.b;
           st:=0
         end
    else begin
           Result.a:=0;
           Result.b:=0;
           st:=3
         end
end {intphi2};

function intphi3 (const iy : dinterval;
                  out st   : Integer) : dinterval;
const ithree : interval = (a:3; b:3);
var status           : Integer;
    cos3y, exp3, piy : interval;
begin
  exp3:=iexp(ithree,status);
  if status=0
    then begin
           piy.a:=iy.a;
           piy.b:=iy.b;
           cos3y:=icos(imul(ithree,piy),status);
           if status=0
             then begin
                    piy:=imul(exp3,cos3y);
                    Result.a:=piy.a;
                    Result.b:=piy.b;
                    st:=0
                  end
             else begin
                    Result.a:=0;
                    Result.b:=0;
                    st:=3
                  end
         end
    else begin
           Result.a:=0;
           Result.b:=0;
           st:=3
         end
end {intphi3};

function intphi4 (const ix : dinterval;
                  out st   : Integer) : dinterval;
const ithree : interval = (a:3; b:3);
var status           : Integer;
    cos3, exp3x, pix : interval;
begin
  cos3:=icos(ithree,status);
  if status=0
    then begin
           pix.a:=ix.a;
           pix.b:=ix.b;
           exp3x:=iexp(imul(ithree,pix),status);
           if status=0
             then begin
                    pix:=imul(cos3,exp3x);
                    Result.a:=pix.a;
                    Result.b:=pix.b;
                    st:=0
                  end
             else begin
                    Result.a:=0;
                    Result.b:=0;
                    st:=3
                  end
         end
    else begin
           Result.a:=0;
           Result.b:=0;
           st:=3
         end
end {intphi4};

function THETA (const ix, iy : dinterval;
                out st       : Integer) : dinterval;
begin
  Result.a:=0;
  Result.b:=0;
  st:=0
end {THETA};

function KSI (const ix, iy : dinterval;
              out st       : Integer) : dinterval;
begin
  Result.a:=0;
  Result.b:=0;
  st:=0
end {KSI};

function modify_time (const time1, time2 : TDateTime) : string;
var code, days, hour_numbers : Integer;
    hours, time              : string;
begin
  time:=TimeToStr(time1-time2);
  days:=DaysBetween(time1,time2);
  if days>0
    then begin
           hours:=Copy(time,1,2);
           Val (hours, hour_numbers, code);
           hour_numbers:=24*days+hour_numbers;
           Str (hour_numbers, hours);
           Result:=hours+Copy(time,3,6)
         end
    else Result:=time
end {modify_time};

var i, imod, j, jmod, l, n, m, st, ui, uj            : Integer;
    alpha, beta, eps, exact, h, k, Pconst, Qconst, w : Extended;
    file_name, left, right, time                     : string;
    z, z1                                            : Char;
    OK, OK1                                          : Boolean;
    plik                                             : TextFile;
    time1, time2                                     : TDateTime;
    puu                                              : interval;
    intalpha, intbeta                                : dinterval;
    uu                                               : intervalmatrix;

begin
  Writeln ('SOLVING THE POISSON EQUATION WITH GIVEN BOUNDARY CONDITIONS');
  Writeln ('BY AN INTERVAL DIFFERENCE METHOD OF FOURTH ORDER');
  Writeln ('WITH DIRECTED INTERVAL ARITHMETIC');
  alpha:=1;
  intalpha.a:=1;
  intalpha.b:=1;
  beta:=1;
  intbeta.a:=1;
  intbeta.b:=1;
  Writeln;
  Writeln ('Insert data');
  repeat
    OK:=false;
    try
      Write ('n = m = ');
      Readln (n);
      if (n>=20) and (n mod 10 = 0)
        then OK:=True
        else Writeln ('Error: the value of n = m must be an integer number '
                      +'equal to 20, 30, 40, ... .');
      except
        Writeln ('Error: the value of n = m must be an integer number equal to '
                 +'20, 30, 40, ... .')
    end
  until OK;
  m:=n;

// The appropriate value of Mconst is very important.
  Pconst:=14643;
  Qconst:=14643;

// Don't change the value of eps!
  eps:=1.0E-16;

  Writeln;
  repeat
    Write ('Output of results (f - file, s - screen): ');
    Readln (z)
  until (z='f') or (z='F') or (z='s') or (z='S');
  if (z='f') or (z='F')
    then begin
           Writeln;
           repeat
             OK:=false;
             try
               Writeln ('Insert the name of file');
               Write ('(the extension .txt will be added automatically): ');
               Readln (file_name);
               file_name:=file_name+'.txt';
               OK1:=False;
               if FileExists(file_name)
                 then begin
                        Writeln ('The file '+file_name+' already exists. ');
                        repeat
                          Writeln ('Do you want to delete this file and create '
                                   +'a new one in its place?');
                          Write ('(y - yes, n - no): ');
                          Readln (z1)
                        until (z1='y') or (z1='Y') or (z1='n') or (z1='N');
                        if (z1='y') or (z1='Y')
                          then OK1:=True
                      end
                 else OK1:=True;
               if OK1
                 then begin
                        AssignFile (plik, file_name);
                        Rewrite (plik);
                        OK:=True;
                        Writeln (plik, 'SOLVING THE POISSON EQUATION WITH GIVEN'
                                 +' BOUNDARY CONDITIONS');
                        Writeln (plik, 'BY AN INTERVAL DIFFERENCE METHOD OF'
                                 +' FOURTH ORDER');
                        Writeln (plik, 'WITH DIRECTED INTERVAL ARITHMETIC')
                      end
             except
               Writeln ('Error: incorrect name of the .txt file.')
             end
           until OK
         end;
  Writeln;
  Writeln ('Calculating ... Please, wait!');
  Writeln;
  time1:=Now;
  directedfullintervalPoisson (n, m, intalpha, intbeta, eps, intf, intphi1,
                               intphi2, intphi3, intphi4, THETA, KSI, Pconst,
                               Qconst, uu, st);
  time2:=Now;
  time:=modify_time(time1,time2);
  Writeln ('status = ', st, ', time = ', time);
  if (z='f') or (z='F')
    then Writeln (plik, 'status = ', st, ', time = ', time);
  if (z='s') or (z='S')
    then begin
           Writeln ('  u - a solution obtained by interval method,');
           Writeln (' eu - the exact (approximate) solution');
           Readln
         end
    else begin
           Writeln (plik, '  u - a solution obtained by interval method,');
           Writeln (plik, ' eu - the exact (approximate) solution')
         end;
  if st=0
    then begin
           h:=alpha/n;
           k:=beta/m;
           l:=0;
           imod:=n div 10;
           jmod:=m div 10;
           ui:=n div 2;
           uj:=m div 2;
           Writeln;
           for j:=0 to m do
             if j mod jmod = 0
               then begin
                      puu:=projection(uu[ui,j]);
                      w:=dint_width(uu[ui,j]);
                      iends_to_strings (puu, left, right);
// Exact solution
                      exact:=Exp(3*ui*h)*Cos(3*j*k);
                      if (z='s') or (z='S')
                        then begin
                               Writeln (' u(', ui*h:3:1, ',', j*k:3:1, ') = [',
                                        left, ',', right, ']');
                               Writeln ('      width =  ', w:25);
                               Writeln ('eu(', ui*h:3:1, ',', j*k:3:1, ') =  ',
                                        exact:25);
                               l:=l+1;
                               if l=4
                                 then begin
                                        l:=0;
                                        Readln
                                      end
                                 else Writeln
                             end
                        else begin
                               Writeln (plik, ' ');
                               Writeln (plik, ' u(', ui*h:3:1, ',', j*k:3:1,
                                        ') = [', left, ',', right, ']');
                               Writeln (plik, '      width =  ', w:25);
                               Writeln (plik,'eu(', ui*h:3:1, ',', j*k:3:1,
                                        ') =  ', exact:25)
                             end
                    end;
           for i:=0 to n do
             if i mod imod = 0
               then begin
                      puu:=projection(uu[i,uj]);
                      w:=dint_width(uu[i,uj]);
                      iends_to_strings (puu, left, right);
// Exact solution
                      exact:=Exp(3*i*h)*Cos(3*uj*k);
                      if (z='s') or (z='S')
                        then begin
                               Writeln (' u(', i*h:3:1, ',', uj*k:3:1, ') = [',
                                        left, ',', right, ']');
                               Writeln ('      width =  ', w:25);
                               Writeln ('eu(', i*h:3:1, ',', uj*k:3:1, ') =  ',
                                        exact:25);
                               l:=l+1;
                               if l=4
                                 then begin
                                        l:=0;
                                        Readln
                                      end
                                 else Writeln
                             end
                        else begin
                               Writeln (plik, ' ');
                               Writeln (plik, ' u(', i*h:3:1, ',', uj*k:3:1,
                                        ') = [', left, ',', right, ']');
                               Writeln (plik, '      width =  ', w:25);
                               Writeln (plik, 'eu(', i*h:3:1, ',', uj*k:3:1,
                                        ') =  ', exact:25)
                             end
                    end;
           Finalize (uu)
         end;
  if (z='f') or (z='F')
    then CloseFile(plik);
  Writeln ('THE END');
  Readln
end.