// $Id: solver.lib,v 1.3 1999-07-06 15:32:58 Singular Exp $
////////////////////////////////////////////////////////////////////////////////
// solver.lib                                                                //
// algorithms for solving algebraic system of dimension zero                              //
// written by Dietmar Hillebrand and ...                                                 //
//                                                                            //
////////////////////////////////////////////////////////////////////////////////

version="$Id: solver.lib,v 1.3 1999-07-06 15:32:58 Singular Exp $";
info="
LIBRARY: solver.lib: PROCEDURES FOR SOLVING ZERO-DIMENSIONAL ALGEBRAIC SYSTEMS
AUTHOR: Dietmar Hillebrand

PROCEDURES:
  triMNewton(G,a,err,itb); computes an improved approximation for
                           a common zero of the triangular system G
                           via multivariate Newton.
";


///////////////////////////////////////////////////////////////////////////////
//
//   Multivariate Newton for triangular systems
//
//
///////////////////////////////////////////////////////////////////////////////

proc triMNewton (ideal G, list a, number err, int itb)
"USAGE:  triMNewtion(G,a,err,itb);
             ideal G=g1,..,gn, a triangular system
                 in n vars, i.e gi=gi(var(n-i+1),..,var(n));
             list of numbers a, an approximation of a common zero of G,
                 (a[i] to be substituted in var(i));
             number err, an error bound;
             int itb, the maximal number of iterations performed.
RETURN:  an improved approximation for a common zero of G;
EXAMPLE: example triMNewton; shows an example
"
{
    if (itb == 0)
    {
        dbprint("iteration bound performed!");
        return(a);
    }

    int i,j,k;
    ideal p=G;
    matrix J=jacob(G);
    list h;
    poly hh;
    int fertig=1;
    int n=nvars(basering);

    for (i = 1; i <= n; i++)
    {
        for (j = n; j >= n-i+1; j--)
        {
            p[i] = subst(p[i],var(j),a[j]);
            for (k = n; k >= n-i+1; k--)
            {
                J[i,k] = subst(J[i,k],var(j),a[j]);
            }
        }
        if (J[i,n-i+1] == 0)
        {
            print("Error: ideal not radical!");
            return();
        }

        // solve linear equations
        hh = -p[i];
        for (j = n; j >= n-i+2; j--)
        {
            hh = hh - J[i,j]*h[j];
        }
        h[n-i+1] = number(hh/J[i,n-i+1]);
    }

    for (i = 1; i <= n; i++)
    {
        if ( abs(h[i]) > err)
        {
            fertig = 0;
            break;
        }
    }
    if ( not fertig )
    {
        for (i = 1; i <= n; i++)
        {
            a[i] = a[i] + h[i];
        }
        return(triMNewton(G,a,err,itb-1));
    }
    else
    {
        return(a);
    }
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=real,(z,y,x),(lp);
   ideal i=x^2-1,y^2+x4-3,z2-y4+x-1;
   list a=2,3,4;
   number e=1.0e-10;

   list l = triMNewton(i,a,e,20);
   l;
}


////////////////////////////////////////////////////////////////////////////////

static proc abs( number r)
{
   if (r >= 0)
   {
       return(r);
   }
   else
   {
       return(-r);
   }
}