program TimeP(output);
{ Timing of a pascal program composed mainly of shortReal arithmetic }
{ P.S.Williams }

var

 x1,x2,x3,j11,j12,j13,j21,j22,j23,j31,j32,j33 : shortReal;
 f1,f2,f3,i11,i12,i13,i21,i22,i23,i31,i32,i33 : shortReal;
 newX1,newX2,newX3,c1,c2,c3,det : shortReal;
 t1,t2,t : integer;

 count,maxNoIts :integer;


function monoTime  : integer ;
{Read the monotonic centi-clock}
var r0 : integer;
begin  
 *SWI_16_42;
 *STR_0,r0;
 monoTime := r0;
end;


begin


 writeln('Generalised Newton Raphson in Pascal');

 t1 := monoTime;

 maxNoIts := 10000;

 x1 := 0.25; x2 := 0.5; x3 := 0.75;

 c1 := (x1+x2+x3)*(x1+x2+x3);
 c2 := (x1+2*x2+3*x3)*(x1+2*x2+3*x3);
 c3 := (2*x1-x2-x3)*(2*x1-x2-x3);

 x1 := 0.2; x2 := 0.45; x3 := 0.7;

 for count := 0 to maxNoIts do
 begin
 
  f1 := (x1+x2+x3)*(x1+x2+x3)         - c1;
  f2 := (x1+2*x2+3*x3)*(x1+2*x2+3*x3) - c2;
  f3 := (2*x1-x2-x3)*(2*x1-x2-x3)     - c3;

  j11 := 2*(x1+x2+x3);     j12 := 2*(x1+x2+x3);     j13 := 2*(x1+x2+x3);
  j21 := 2*(x1+2*x2+3*x3); j22 := 4*(x1+2*x2+3*x3); j23 := 6*(x1+2*x2+3*x3);
  j31 := 4*(2*x1-x2-x3);   j32 := -(2*x1-x2-x3);    j33 := -(2*x1-x2-x3);



  det :=  j11*j22*j33 - j11*j23*j32 
       - j12*j21*j33 + j12*j23*j31
       + j13*j21*j32 - j13*j22*j31;


  i11 :=  (j22*j33 - j23*j32)/det;
  i12 := -(j12*j33 - j13*j32)/det;   
  i13 :=  (j12*j23 - j13*j22)/det;  
                                                                           
  i21 := -(j21*j33 - j23*j31)/det;
  i22 :=  (j11*j33 - j13*j31)/det;
  i23 := -(j11*j23 - j13*j21)/det;
                                                                           
  i31 :=  (j21*j32 - j22*j31)/det;
  i32 := -(j11*j32 - j12*j31)/det;
  i33 :=  (j11*j22 - j12*j21)/det;

  newX1 := x1 - (i11*f1+i12*f2+i13*f3);
  newX2 := x2 - (i21*f1+i22*f2+i23*f3);  
  newX3 := x3 - (i31*f1+i32*f2+i33*f3);

  x1 := newX1; x2 := newX2; x3 := newX3;
                         
     
 end;

 t2 := monoTime;
 t := t2-t1;
 writeln(maxNoIts,', iterations took ',t,' centi-secs');
 writeln('x1 = ',x1,' x2 = ',x2,' x3 = ',x3);


end.
