Re: Problem with symbolic calculation

[Karim BELABAS <Karim.Belabas@math.u-psud.fr>]

[Bill Daly:]
>    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.

Polynomial arithmetic is highly inefficient in PARI when the number of
variables gets large (25 is very large there). If the columns had been
oriented the other way round (easy first step, the way you do it by hand
later), it would have succeeded, after a major computational effort (some
kind of sparse representation for multinomials is sorely needed there...)

>    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

Internal bug in the generic kernel: one subfunction assumed it would only
treat simplified objects (as in GP), and it received (0 / some denominator).
First step was to compute numerator content (0), and cancel it (0 / 0). Ouch.

[of course the bug would not occur when pasting the two polynomials you sent:
the internal structure (including those stupid 0s...) has been simplified out
of the way]

The patch below also corrects two weird generic bugs:

1)  [1]~ * [[1]] --> SEGV

2)  complicated multiplication between polmods with incompatible intmod
coefficients could produce 0 instead of the expected result.

  Karim.

*** src/basemath/gen1.c.orig	Wed Jun  2 16:43:17 1999
--- src/basemath/gen1.c	Thu May 27 18:37:20 1999
***************
*** 737,745 ****
  
    if (gcmp0(x)) return gcopy(x);
  
!   z = cgetg(3, t_RFRAC);
!   y1=(GEN)y[1];
    y2=(GEN)y[2]; tx = typ(x);
    if (is_const_t(tx) || varn(x) > mingvar(y1,y2)) { cx = x; x = gun; }
    else
    {
--- 737,745 ----
  
    if (gcmp0(x)) return gcopy(x);
  
!   y1=(GEN)y[1]; if (gcmp0(y1)) return gcopy(y1);
    y2=(GEN)y[2]; tx = typ(x);
+   z = cgetg(3, t_RFRAC);
    if (is_const_t(tx) || varn(x) > mingvar(y1,y2)) { cx = x; x = gun; }
    else
    {
***************
*** 828,836 ****
      {
        case t_INTMOD: break;
        case t_INT:
!         if (*p) { y[i] = lmodii(p1, *p); continue; }
        default:
!         return Q? 1: 0;
      }
      p2 = (GEN)p1[1];
      if (pr==NULL) pr = p2;
--- 828,837 ----
      {
        case t_INTMOD: break;
        case t_INT:
!         if (*p) p1 = modii(p1, *p);
!         y[i] = (long)p1; continue;
        default:
!         return (Q && !pr)? 1: 0;
      }
      p2 = (GEN)p1[1];
      if (pr==NULL) pr = p2;
***************
*** 839,845 ****
        if (!egalii(p2, pr))
        {
          if (DEBUGMEM) err(warner,"different modulus in ff_poltype");
!         return Q? 1: 0;
        }
        if (DEBUGMEM > 2) err(warner,"different pointers in ff_poltype");
      }
--- 840,846 ----
        if (!egalii(p2, pr))
        {
          if (DEBUGMEM) err(warner,"different modulus in ff_poltype");
!         return 0;
        }
        if (DEBUGMEM > 2) err(warner,"different pointers in ff_poltype");
      }
***************
*** 1165,1171 ****
          {
            case t_VEC:
              z=cgetg(ly,t_MAT);
!             for (i=1; i<ly; i++) z[i]=lmul((GEN)y[i],x);
              return z;
  
            case t_MAT:
--- 1166,1177 ----
          {
            case t_VEC:
              z=cgetg(ly,t_MAT);
!             for (i=1; i<ly; i++)
!             {
!               p1 = gmul((GEN)y[i],x);
!               if (typ(p1) != t_COL) err(gmuleri,tx,ty);
!               z[i]=(long)p1;
!             }
              return z;
  
            case t_MAT:
__
Karim Belabas                    email: Karim.Belabas@math.u-psud.fr
Dep. de Mathematiques, Bat. 425
Universite Paris-Sud             Tel: (00 33) 1 69 15 57 48
F-91405 Orsay (France)           Fax: (00 33) 1 69 15 60 19
--
PARI/GP Home Page: http://hasse.mathematik.tu-muenchen.de/ntsw/pari/