Bill Daly on Thu, 27 May 1999 11:59:44 -0400

[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

 Problem with symbolic calculation

```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