Problem with symbolic calculation

[bill.daly@tradition.co.uk (Bill Daly)]

I tried the following GP program. The comments should explain the
problem.

/*
   I want to solve (symbolically) the problem of defining a conic

      y^2 + r*x*y + s*x^2 + t*y + u*x + v = 0

   passing through a given set of 5 points [x1,y1] through [x5,y5]. This
is
   a straightforward linear algebra problem which PARI should be able to
handle.

   I start by defining a matrix m whose i-th row is

      [xi*yi,xi^2,yi,xi,1]

   and a column vector v = [-y1^2,-y2^2,-y3^2,-y4^2,-y5^2]~.
*/

m = matid(5);
v = vectorv(5,j,0);
{
   for (j = 1, 5,
      sx = eval(Str("x"j));
      sy = eval(Str("y"j));
      m[j,1] = sx*sy;
      m[j,2] = sx^2;
      m[j,3] = sy;
      m[j,4] = sx;
      m[j,5] = 1;
      v[j] = -sy^2;
   );
}
/*
   At this point, I should be able to solve for [r,s,t,u,v]~ by
computing
   m^-1*v, but when I tried to invert m, PARI failed to give an answer
   even with a 64 Meg stack. This doesn't look like such a difficult
problem,
   so I tried to solve the system of equations longhand.
*/
{
   forstep (j = 5, 2, -1,
      for (k = 1, j-1,
         c = m[k,j]/m[j,j];
         v[k] -= v[j]*c;
         m[k,] -= m[j,]*c;
      );
   );
}
/*
   The above loop should work, but PARI gives a divide-by-zero error in
   calculating c when k = 1 and j = 2. At this point, we have

      m[1,2] = x1^2 + ((((y4 - y5)/(x4 - x5))*x3^2 + (((-x4^2 +
x5^2)/(x4 - x5))
         *y3 + ((y5*x4^2 - x5^2*y4)/(x4 - x5))))/(((-y4 + y5)/(x4 -
x5))*x3 + 
         (y3 + ((-y5*x4 + x5*y4)/(x4 - x5)))))*x1 + (((-x3^2 + ((x4^2 -
x5^2)/
         (x4 - x5))*x3 + ((-x5*x4^2 + x5^2*x4)/(x4 - x5)))/(((-y4 + y5)/
         (x4 - x5))*x3 + (y3 + ((-y5*x4 + x5*y4)/(x4 - x5)))))*y1 +
((((y5*x4 - 
         x5*y4)/(x4 - x5))*x3^2 + ((-y5*x4^2 + x5^2*y4)/(x4 - x5))*x3 + 
         ((x5*x4^2 - x5^2*x4)/(x4 - x5))*y3)/(((-y4 + y5)/(x4 - x5))*x3
+ 
         (y3 + ((-y5*x4 + x5*y4)/(x4 - x5))))))

      m[2,2] = ((((y4 - y5)/(x4 - x5))*x2 + (-y2 + ((y5*x4 - x5*y4)/(x4
- x5))))
         *x3^2 + (((-y4 + y5)/(x4 - x5))*x2^2 + (((x4^2 - x5^2)/(x4 -
x5))*y2 + 
         ((-y5*x4^2 + x5^2*y4)/(x4 - x5))))*x3 + ((x2^2 + ((-x4^2 +
x5^2)/
         (x4 - x5))*x2 + ((x5*x4^2 - x5^2*x4)/(x4 - x5)))*y3 + 
         (((-y5*x4 + x5*y4)/(x4 - x5))*x2^2 + ((y5*x4^2 - x5^2*y4)/(x4 -
x5))
         *x2 + ((-x5*x4^2 + x5^2*x4)/(x4 - x5))*y2)))/(((-y4 + y5)/(x4 -
x5))
         *x3 + (y3 + ((-y5*x4 + x5*y4)/(x4 - x5))))

   and m[1,2]/m[2,2] gives a divide-by-zero error. This can't be right.
*/


Regards,

Bill
**********************************************************************
This email and any files transmitted with it are confidential and
intended solely for the use of the individual or entity to whom
they are addressed.
This footnote also confirms that this email message has been swept by
MIMEsweeper for the presence of computer viruses.
**********************************************************************