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).
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)
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.
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.
BOOLEAN	loNewtonP (leftv res, leftv arg1)
	compute Newton Polytopes of input polynomials
BOOLEAN	loSimplex (leftv res, leftv args)
	Implementation of the Simplex Algorithm.

Macro Definition Documentation

DEFAULT_DIGITS

#define DEFAULT_DIGITS   30

Definition at line 13 of file mpr_inout.h.

MPR_DENSE

#define MPR_DENSE   1

Definition at line 15 of file mpr_inout.h.

MPR_SPARSE

#define MPR_SPARSE   2

Definition at line 16 of file mpr_inout.h.

Function Documentation

loNewtonP()

BOOLEAN loNewtonP	(	leftv	res,
		leftv	arg1
	)

compute Newton Polytopes of input polynomials

Definition at line 4647 of file ipshell.cc.

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

loSimplex()

BOOLEAN loSimplex	(	leftv	res,
		leftv	args
	)

Implementation of the Simplex Algorithm.

For args, see class simplex.

Definition at line 4653 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;
}

nuLagSolve()

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 4762 of file ipshell.cc.

{
  poly gls;
  gls= (poly)(arg1->Data());
  int howclean= (int)(long)arg3->Data();
 
  if ( gls == NULL || pIsConstant( gls ) )
  {
    WerrorS("Input polynomial is constant!");
    return TRUE;
  }
 
  if (rField_is_Zp(currRing))
  {
    int* r=Zp_roots(gls, currRing);
    lists rlist;
    rlist= (lists)omAlloc( sizeof(slists) );
    rlist->Init( r[0] );
    for(int i=r[0];i>0;i--)
    {
      rlist->m[i-1].data=n_Init(r[i],currRing);
      rlist->m[i-1].rtyp=NUMBER_CMD;
    }
    omFree(r);
    res->data=rlist;
    res->rtyp= LIST_CMD;
    return FALSE;
  }
  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 );
  }
 
  int ldummy;
  int deg= currRing->pLDeg( gls, &ldummy, currRing );
  int i,vpos=0;
  poly piter;
  lists elist;
 
  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 && 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;
 
  lists rlist;
  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->data= (void*)rlist;
 
  return FALSE;
}

nuMPResMat()

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 4739 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;
}

nuUResSolve()

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 5006 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;
}

nuVanderSys()

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 4905 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;
}