mpr_inout.h File Reference

Go to the source code of this file.

Macros
#define	DEFAULT_DIGITS   30
#define	MPR_DENSE   1
#define	MPR_SPARSE   2

Functions
BOOLEAN	nuUResSolve (leftv res, leftv args)
	solve a multipolynomial system using the u-resultant Input ideal must be 0-dimensional and (currRing->N) == IDELEMS(ideal). More...
BOOLEAN	nuMPResMat (leftv res, leftv arg1, leftv arg2)
	returns module representing the multipolynomial resultant matrix Arguments 2: ideal i, int k k=0: use sparse resultant matrix of Gelfand, Kapranov and Zelevinsky k=1: use resultant matrix of Macaulay (k=0 is default) More...
BOOLEAN	nuLagSolve (leftv res, leftv arg1, leftv arg2, leftv arg3)
	find the (complex) roots an univariate polynomial Determines the roots of an univariate polynomial using Laguerres' root-solver. More...
BOOLEAN	nuVanderSys (leftv res, leftv arg1, leftv arg2, leftv arg3)
	COMPUTE: polynomial p with values given by v at points p1,..,pN derived from p; more precisely: consider p as point in K^n and v as N elements in K, let p1,..,pN be the points in K^n obtained by evaluating all monomials of degree 0,1,...,N at p in lexicographical order, then the procedure computes the polynomial f satisfying f(pi) = v[i] RETURN: polynomial f of degree d. More...
BOOLEAN	loNewtonP (leftv res, leftv arg1)
	compute Newton Polytopes of input polynomials More...
BOOLEAN	loSimplex (leftv res, leftv args)
	Implementation of the Simplex Algorithm. More...

Macro Definition Documentation

#define DEFAULT_DIGITS   30

Definition at line 13 of file mpr_inout.h.

#define MPR_DENSE   1

Definition at line 15 of file mpr_inout.h.

#define MPR_SPARSE   2

Definition at line 16 of file mpr_inout.h.

Function Documentation

BOOLEAN loNewtonP	(	leftv	res,
		leftv	arg1
	)

compute Newton Polytopes of input polynomials

Definition at line 4269 of file ipshell.cc.

{
  res->data= (void*)loNewtonPolytope( (ideal)arg1->Data() );
  return FALSE;
}

BOOLEAN loSimplex	(	leftv	res,
		leftv	args
	)

Implementation of the Simplex Algorithm.

For args, see class simplex.

Definition at line 4275 of file ipshell.cc.

{
  if ( !(rField_is_long_R(currRing)) )
  {
    WerrorS("Ground field not implemented!");
    return TRUE;
  }

  simplex * LP;
  matrix m;

  leftv v= args;
  if ( v->Typ() != MATRIX_CMD ) // 1: matrix
    return TRUE;
  else
    m= (matrix)(v->CopyD());

  LP = new simplex(MATROWS(m),MATCOLS(m));
  LP->mapFromMatrix(m);

  v= v->next;
  if ( v->Typ() != INT_CMD )    // 2: m = number of constraints
    return TRUE;
  else
    LP->m= (int)(long)(v->Data());

  v= v->next;
  if ( v->Typ() != INT_CMD )    // 3: n = number of variables
    return TRUE;
  else
    LP->n= (int)(long)(v->Data());

  v= v->next;
  if ( v->Typ() != INT_CMD )    // 4: m1 = number of <= constraints
    return TRUE;
  else
    LP->m1= (int)(long)(v->Data());

  v= v->next;
  if ( v->Typ() != INT_CMD )    // 5: m2 = number of >= constraints
    return TRUE;
  else
    LP->m2= (int)(long)(v->Data());

  v= v->next;
  if ( v->Typ() != INT_CMD )    // 6: m3 = number of == constraints
    return TRUE;
  else
    LP->m3= (int)(long)(v->Data());

#ifdef mprDEBUG_PROT
  Print("m (constraints) %d\n",LP->m);
  Print("n (columns) %d\n",LP->n);
  Print("m1 (<=) %d\n",LP->m1);
  Print("m2 (>=) %d\n",LP->m2);
  Print("m3 (==) %d\n",LP->m3);
#endif

  LP->compute();

  lists lres= (lists)omAlloc( sizeof(slists) );
  lres->Init( 6 );

  lres->m[0].rtyp= MATRIX_CMD; // output matrix
  lres->m[0].data=(void*)LP->mapToMatrix(m);

  lres->m[1].rtyp= INT_CMD;   // found a solution?
  lres->m[1].data=(void*)(long)LP->icase;

  lres->m[2].rtyp= INTVEC_CMD;
  lres->m[2].data=(void*)LP->posvToIV();

  lres->m[3].rtyp= INTVEC_CMD;
  lres->m[3].data=(void*)LP->zrovToIV();

  lres->m[4].rtyp= INT_CMD;
  lres->m[4].data=(void*)(long)LP->m;

  lres->m[5].rtyp= INT_CMD;
  lres->m[5].data=(void*)(long)LP->n;

  res->data= (void*)lres;

  return FALSE;
}

BOOLEAN nuLagSolve	(	leftv	res,
		leftv	arg1,
		leftv	arg2,
		leftv	arg3
	)

find the (complex) roots an univariate polynomial Determines the roots of an univariate polynomial using Laguerres' root-solver.

Good for polynomials with low and middle degree (<40). Arguments 3: poly arg1 , int arg2 , int arg3 arg2>0: defines precision of fractional part if ground field is Q arg3: number of iterations for approximation of roots (default=2) Returns a list of all (complex) roots of the polynomial arg1

Definition at line 4384 of file ipshell.cc.

{

  poly gls;
  gls= (poly)(arg1->Data());
  int howclean= (int)(long)arg3->Data();

  if ( !(rField_is_R(currRing) ||
         rField_is_Q(currRing) ||
         rField_is_long_R(currRing) ||
         rField_is_long_C(currRing)) )
  {
    WerrorS("Ground field not implemented!");
    return TRUE;
  }

  if ( !(rField_is_R(currRing) || rField_is_long_R(currRing) || \
                          rField_is_long_C(currRing)) )
  {
    unsigned long int ii = (unsigned long int)arg2->Data();
    setGMPFloatDigits( ii, ii );
  }

  if ( gls == NULL || pIsConstant( gls ) )
  {
    WerrorS("Input polynomial is constant!");
    return TRUE;
  }

  int ldummy;
  int deg= currRing->pLDeg( gls, &ldummy, currRing );
  //  int deg= pDeg( gls );
  //  int len= pLength( gls );
  int i,vpos=0;
  poly piter;
  lists elist;
  lists rlist;

  elist= (lists)omAlloc( sizeof(slists) );
  elist->Init( 0 );

  if ( rVar(currRing) > 1 )
  {
    piter= gls;
    for ( i= 1; i <= rVar(currRing); i++ )
      if ( pGetExp( piter, i ) )
      {
        vpos= i;
        break;
      }
    while ( piter )
    {
      for ( i= 1; i <= rVar(currRing); i++ )
        if ( (vpos != i) && (pGetExp( piter, i ) != 0) )
        {
          WerrorS("The input polynomial must be univariate!");
          return TRUE;
        }
      pIter( piter );
    }
  }

  rootContainer * roots= new rootContainer();
  number * pcoeffs= (number *)omAlloc( (deg+1) * sizeof( number ) );
  piter= gls;
  for ( i= deg; i >= 0; i-- )
  {
    //if ( piter ) Print("deg %d, pDeg(piter) %d\n",i,pTotaldegree(piter));
    if ( piter && pTotaldegree(piter) == i )
    {
      pcoeffs[i]= nCopy( pGetCoeff( piter ) );
      //nPrint( pcoeffs[i] );PrintS("  ");
      pIter( piter );
    }
    else
    {
      pcoeffs[i]= nInit(0);
    }
  }

#ifdef mprDEBUG_PROT
  for (i=deg; i >= 0; i--)
  {
    nPrint( pcoeffs[i] );PrintS("  ");
  }
  PrintLn();
#endif

  roots->fillContainer( pcoeffs, NULL, 1, deg, rootContainer::onepoly, 1 );
  roots->solver( howclean );

  int elem= roots->getAnzRoots();
  char *dummy;
  int j;

  rlist= (lists)omAlloc( sizeof(slists) );
  rlist->Init( elem );

  if (rField_is_long_C(currRing))
  {
    for ( j= 0; j < elem; j++ )
    {
      rlist->m[j].rtyp=NUMBER_CMD;
      rlist->m[j].data=(void *)nCopy((number)(roots->getRoot(j)));
      //rlist->m[j].data=(void *)(number)(roots->getRoot(j));
    }
  }
  else
  {
    for ( j= 0; j < elem; j++ )
    {
      dummy = complexToStr( (*roots)[j], gmp_output_digits, currRing->cf );
      rlist->m[j].rtyp=STRING_CMD;
      rlist->m[j].data=(void *)dummy;
    }
  }

  elist->Clean();
  //omFreeSize( (ADDRESS) elist, sizeof(slists) );

  // this is (via fillContainer) the same data as in root
  //for ( i= deg; i >= 0; i-- ) nDelete( &pcoeffs[i] );
  //omFreeSize( (ADDRESS) pcoeffs, (deg+1) * sizeof( number ) );

  delete roots;

  res->rtyp= LIST_CMD;
  res->data= (void*)rlist;

  return FALSE;
}

BOOLEAN nuMPResMat	(	leftv	res,
		leftv	arg1,
		leftv	arg2
	)

returns module representing the multipolynomial resultant matrix Arguments 2: ideal i, int k k=0: use sparse resultant matrix of Gelfand, Kapranov and Zelevinsky k=1: use resultant matrix of Macaulay (k=0 is default)

Definition at line 4361 of file ipshell.cc.

{
  ideal gls = (ideal)(arg1->Data());
  int imtype= (int)(long)arg2->Data();

  uResultant::resMatType mtype= determineMType( imtype );

  // check input ideal ( = polynomial system )
  if ( mprIdealCheck( gls, arg1->Name(), mtype, true ) != mprOk )
  {
    return TRUE;
  }

  uResultant *resMat= new uResultant( gls, mtype, false );
  if (resMat!=NULL)
  {
    res->rtyp = MODUL_CMD;
    res->data= (void*)resMat->accessResMat()->getMatrix();
    if (!errorreported) delete resMat;
  }
  return errorreported;
}

BOOLEAN nuUResSolve	(	leftv	res,
		leftv	args
	)

solve a multipolynomial system using the u-resultant Input ideal must be 0-dimensional and (currRing->N) == IDELEMS(ideal).

Resultant method can be MPR_DENSE, which uses Macaulay Resultant (good for dense homogeneous polynoms) or MPR_SPARSE, which uses Sparse Resultant (Gelfand, Kapranov, Zelevinsky). Arguments 4: ideal i, int k, int l, int m k=0: use sparse resultant matrix of Gelfand, Kapranov and Zelevinsky k=1: use resultant matrix of Macaulay (k=0 is default) l>0: defines precision of fractional part if ground field is Q m=0,1,2: number of iterations for approximation of roots (default=2) Returns a list containing the roots of the system.

Definition at line 4617 of file ipshell.cc.

{
  leftv v= args;

  ideal gls;
  int imtype;
  int howclean;

  // get ideal
  if ( v->Typ() != IDEAL_CMD )
    return TRUE;
  else gls= (ideal)(v->Data());
  v= v->next;

  // get resultant matrix type to use (0,1)
  if ( v->Typ() != INT_CMD )
    return TRUE;
  else imtype= (int)(long)v->Data();
  v= v->next;

  if (imtype==0)
  {
    ideal test_id=idInit(1,1);
    int j;
    for(j=IDELEMS(gls)-1;j>=0;j--)
    {
      if (gls->m[j]!=NULL)
      {
        test_id->m[0]=gls->m[j];
        intvec *dummy_w=id_QHomWeight(test_id, currRing);
        if (dummy_w!=NULL)
        {
          WerrorS("Newton polytope not of expected dimension");
          delete dummy_w;
          return TRUE;
        }
      }
    }
  }

  // get and set precision in digits ( > 0 )
  if ( v->Typ() != INT_CMD )
    return TRUE;
  else if ( !(rField_is_R(currRing) || rField_is_long_R(currRing) || \
                          rField_is_long_C(currRing)) )
  {
    unsigned long int ii=(unsigned long int)v->Data();
    setGMPFloatDigits( ii, ii );
  }
  v= v->next;

  // get interpolation steps (0,1,2)
  if ( v->Typ() != INT_CMD )
    return TRUE;
  else howclean= (int)(long)v->Data();

  uResultant::resMatType mtype= determineMType( imtype );
  int i,count;
  lists listofroots= NULL;
  number smv= NULL;
  BOOLEAN interpolate_det= (mtype==uResultant::denseResMat)?TRUE:FALSE;

  //emptylist= (lists)omAlloc( sizeof(slists) );
  //emptylist->Init( 0 );

  //res->rtyp = LIST_CMD;
  //res->data= (void *)emptylist;

  // check input ideal ( = polynomial system )
  if ( mprIdealCheck( gls, args->Name(), mtype ) != mprOk )
  {
    return TRUE;
  }

  uResultant * ures;
  rootContainer ** iproots;
  rootContainer ** muiproots;
  rootArranger * arranger;

  // main task 1: setup of resultant matrix
  ures= new uResultant( gls, mtype );
  if ( ures->accessResMat()->initState() != resMatrixBase::ready )
  {
    WerrorS("Error occurred during matrix setup!");
    return TRUE;
  }

  // if dense resultant, check if minor nonsingular
  if ( mtype == uResultant::denseResMat )
  {
    smv= ures->accessResMat()->getSubDet();
#ifdef mprDEBUG_PROT
    PrintS("// Determinant of submatrix: ");nPrint(smv);PrintLn();
#endif
    if ( nIsZero(smv) )
    {
      WerrorS("Unsuitable input ideal: Minor of resultant matrix is singular!");
      return TRUE;
    }
  }

  // main task 2: Interpolate specialized resultant polynomials
  if ( interpolate_det )
    iproots= ures->interpolateDenseSP( false, smv );
  else
    iproots= ures->specializeInU( false, smv );

  // main task 3: Interpolate specialized resultant polynomials
  if ( interpolate_det )
    muiproots= ures->interpolateDenseSP( true, smv );
  else
    muiproots= ures->specializeInU( true, smv );

#ifdef mprDEBUG_PROT
  int c= iproots[0]->getAnzElems();
  for (i=0; i < c; i++) pWrite(iproots[i]->getPoly());
  c= muiproots[0]->getAnzElems();
  for (i=0; i < c; i++) pWrite(muiproots[i]->getPoly());
#endif

  // main task 4: Compute roots of specialized polys and match them up
  arranger= new rootArranger( iproots, muiproots, howclean );
  arranger->solve_all();

  // get list of roots
  if ( arranger->success() )
  {
    arranger->arrange();
    listofroots= listOfRoots(arranger, gmp_output_digits );
  }
  else
  {
    WerrorS("Solver was unable to find any roots!");
    return TRUE;
  }

  // free everything
  count= iproots[0]->getAnzElems();
  for (i=0; i < count; i++) delete iproots[i];
  omFreeSize( (ADDRESS) iproots, count * sizeof(rootContainer*) );
  count= muiproots[0]->getAnzElems();
  for (i=0; i < count; i++) delete muiproots[i];
  omFreeSize( (ADDRESS) muiproots, count * sizeof(rootContainer*) );

  delete ures;
  delete arranger;
  nDelete( &smv );

  res->data= (void *)listofroots;

  //emptylist->Clean();
  //  omFreeSize( (ADDRESS) emptylist, sizeof(slists) );

  return FALSE;
}

BOOLEAN nuVanderSys	(	leftv	res,
		leftv	arg1,
		leftv	arg2,
		leftv	arg3
	)

COMPUTE: polynomial p with values given by v at points p1,..,pN derived from p; more precisely: consider p as point in K^n and v as N elements in K, let p1,..,pN be the points in K^n obtained by evaluating all monomials of degree 0,1,...,N at p in lexicographical order, then the procedure computes the polynomial f satisfying f(pi) = v[i] RETURN: polynomial f of degree d.

Definition at line 4516 of file ipshell.cc.

{
  int i;
  ideal p,w;
  p= (ideal)arg1->Data();
  w= (ideal)arg2->Data();

  // w[0] = f(p^0)
  // w[1] = f(p^1)
  // ...
  // p can be a vector of numbers (multivariate polynom)
  //   or one number (univariate polynom)
  // tdg = deg(f)

  int n= IDELEMS( p );
  int m= IDELEMS( w );
  int tdg= (int)(long)arg3->Data();

  res->data= (void*)NULL;

  // check the input
  if ( tdg < 1 )
  {
    WerrorS("Last input parameter must be > 0!");
    return TRUE;
  }
  if ( n != rVar(currRing) )
  {
    Werror("Size of first input ideal must be equal to %d!",rVar(currRing));
    return TRUE;
  }
  if ( m != (int)pow((double)tdg+1,(double)n) )
  {
    Werror("Size of second input ideal must be equal to %d!",
      (int)pow((double)tdg+1,(double)n));
    return TRUE;
  }
  if ( !(rField_is_Q(currRing) /* ||
         rField_is_R() || rField_is_long_R() ||
         rField_is_long_C()*/ ) )
         {
    WerrorS("Ground field not implemented!");
    return TRUE;
  }

  number tmp;
  number *pevpoint= (number *)omAlloc( n * sizeof( number ) );
  for ( i= 0; i < n; i++ )
  {
    pevpoint[i]=nInit(0);
    if (  (p->m)[i] )
    {
      tmp = pGetCoeff( (p->m)[i] );
      if ( nIsZero(tmp) || nIsOne(tmp) || nIsMOne(tmp) )
      {
        omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
        WerrorS("Elements of first input ideal must not be equal to -1, 0, 1!");
        return TRUE;
      }
    } else tmp= NULL;
    if ( !nIsZero(tmp) )
    {
      if ( !pIsConstant((p->m)[i]))
      {
        omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
        WerrorS("Elements of first input ideal must be numbers!");
        return TRUE;
      }
      pevpoint[i]= nCopy( tmp );
    }
  }

  number *wresults= (number *)omAlloc( m * sizeof( number ) );
  for ( i= 0; i < m; i++ )
  {
    wresults[i]= nInit(0);
    if ( (w->m)[i] && !nIsZero(pGetCoeff((w->m)[i])) )
    {
      if ( !pIsConstant((w->m)[i]))
      {
        omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
        omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
        WerrorS("Elements of second input ideal must be numbers!");
        return TRUE;
      }
      wresults[i]= nCopy(pGetCoeff((w->m)[i]));
    }
  }

  vandermonde vm( m, n, tdg, pevpoint, FALSE );
  number *ncpoly= vm.interpolateDense( wresults );
  // do not free ncpoly[]!!
  poly rpoly= vm.numvec2poly( ncpoly );

  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
  omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );

  res->data= (void*)rpoly;
  return FALSE;
}