My Project
Loading...
Searching...
No Matches
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 4559 of file ipshell.cc.

4560{
4561 res->data= (void*)loNewtonPolytope( (ideal)arg1->Data() );
4562 return FALSE;
4563}
#define FALSE
Definition auxiliary.h:97
void * Data()
Definition subexpr.cc:1192
CanonicalForm res
Definition facAbsFact.cc:60
ideal loNewtonPolytope(const ideal id)
Definition mpr_base.cc:3191

◆ loSimplex()

BOOLEAN loSimplex ( leftv res,
leftv args )

Implementation of the Simplex Algorithm.

For args, see class simplex.

Definition at line 4565 of file ipshell.cc.

4566{
4567 if ( !(rField_is_long_R(currRing)) )
4568 {
4569 WerrorS("Ground field not implemented!");
4570 return TRUE;
4571 }
4572
4573 simplex * LP;
4574 matrix m;
4575
4576 leftv v= args;
4577 if ( v->Typ() != MATRIX_CMD ) // 1: matrix
4578 return TRUE;
4579 else
4580 m= (matrix)(v->CopyD());
4581
4582 LP = new simplex(MATROWS(m),MATCOLS(m));
4583 LP->mapFromMatrix(m);
4584
4585 v= v->next;
4586 if ( v->Typ() != INT_CMD ) // 2: m = number of constraints
4587 return TRUE;
4588 else
4589 LP->m= (int)(long)(v->Data());
4590
4591 v= v->next;
4592 if ( v->Typ() != INT_CMD ) // 3: n = number of variables
4593 return TRUE;
4594 else
4595 LP->n= (int)(long)(v->Data());
4596
4597 v= v->next;
4598 if ( v->Typ() != INT_CMD ) // 4: m1 = number of <= constraints
4599 return TRUE;
4600 else
4601 LP->m1= (int)(long)(v->Data());
4602
4603 v= v->next;
4604 if ( v->Typ() != INT_CMD ) // 5: m2 = number of >= constraints
4605 return TRUE;
4606 else
4607 LP->m2= (int)(long)(v->Data());
4608
4609 v= v->next;
4610 if ( v->Typ() != INT_CMD ) // 6: m3 = number of == constraints
4611 return TRUE;
4612 else
4613 LP->m3= (int)(long)(v->Data());
4614
4615#ifdef mprDEBUG_PROT
4616 Print("m (constraints) %d\n",LP->m);
4617 Print("n (columns) %d\n",LP->n);
4618 Print("m1 (<=) %d\n",LP->m1);
4619 Print("m2 (>=) %d\n",LP->m2);
4620 Print("m3 (==) %d\n",LP->m3);
4621#endif
4622
4623 LP->compute();
4624
4625 lists lres= (lists)omAlloc( sizeof(slists) );
4626 lres->Init( 6 );
4627
4628 lres->m[0].rtyp= MATRIX_CMD; // output matrix
4629 lres->m[0].data=(void*)LP->mapToMatrix(m);
4630
4631 lres->m[1].rtyp= INT_CMD; // found a solution?
4632 lres->m[1].data=(void*)(long)LP->icase;
4633
4634 lres->m[2].rtyp= INTVEC_CMD;
4635 lres->m[2].data=(void*)LP->posvToIV();
4636
4637 lres->m[3].rtyp= INTVEC_CMD;
4638 lres->m[3].data=(void*)LP->zrovToIV();
4639
4640 lres->m[4].rtyp= INT_CMD;
4641 lres->m[4].data=(void*)(long)LP->m;
4642
4643 lres->m[5].rtyp= INT_CMD;
4644 lres->m[5].data=(void*)(long)LP->n;
4645
4646 res->data= (void*)lres;
4647
4648 return FALSE;
4649}
#define TRUE
Definition auxiliary.h:101
int m
Definition cfEzgcd.cc:128
Linear Programming / Linear Optimization using Simplex - Algorithm.
intvec * zrovToIV()
BOOLEAN mapFromMatrix(matrix m)
void compute()
matrix mapToMatrix(matrix m)
intvec * posvToIV()
int rtyp
Definition subexpr.h:91
void * data
Definition subexpr.h:88
Definition lists.h:24
sleftv * m
Definition lists.h:46
INLINE_THIS void Init(int l=0)
#define Print
Definition emacs.cc:80
const Variable & v
< [in] a sqrfree bivariate poly
Definition facBivar.h:39
void WerrorS(const char *s)
Definition feFopen.cc:24
@ MATRIX_CMD
Definition grammar.cc:287
ip_smatrix * matrix
Definition matpol.h:43
#define MATROWS(i)
Definition matpol.h:26
#define MATCOLS(i)
Definition matpol.h:27
slists * lists
#define omAlloc(size)
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition polys.cc:13
static BOOLEAN rField_is_long_R(const ring r)
Definition ring.h:548
sleftv * leftv
Definition structs.h:53
@ INTVEC_CMD
Definition tok.h:101
@ INT_CMD
Definition tok.h:96

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

4675{
4676 poly gls;
4677 gls= (poly)(arg1->Data());
4678 int howclean= (int)(long)arg3->Data();
4679
4680 if ( gls == NULL || pIsConstant( gls ) )
4681 {
4682 WerrorS("Input polynomial is constant!");
4683 return TRUE;
4684 }
4685
4687 {
4688 int* r=Zp_roots(gls, currRing);
4689 lists rlist;
4690 rlist= (lists)omAlloc( sizeof(slists) );
4691 rlist->Init( r[0] );
4692 for(int i=r[0];i>0;i--)
4693 {
4694 rlist->m[i-1].data=n_Init(r[i],currRing->cf);
4695 rlist->m[i-1].rtyp=NUMBER_CMD;
4696 }
4697 omFree(r);
4698 res->data=rlist;
4699 res->rtyp= LIST_CMD;
4700 return FALSE;
4701 }
4702 if ( !(rField_is_R(currRing) ||
4706 {
4707 WerrorS("Ground field not implemented!");
4708 return TRUE;
4709 }
4710
4713 {
4714 unsigned long int ii = (unsigned long int)arg2->Data();
4715 setGMPFloatDigits( ii, ii );
4716 }
4717
4718 int ldummy;
4719 int deg= currRing->pLDeg( gls, &ldummy, currRing );
4720 int i,vpos=0;
4721 poly piter;
4722 lists elist;
4723
4724 elist= (lists)omAlloc( sizeof(slists) );
4725 elist->Init( 0 );
4726
4727 if ( rVar(currRing) > 1 )
4728 {
4729 piter= gls;
4730 for ( i= 1; i <= rVar(currRing); i++ )
4731 if ( pGetExp( piter, i ) )
4732 {
4733 vpos= i;
4734 break;
4735 }
4736 while ( piter )
4737 {
4738 for ( i= 1; i <= rVar(currRing); i++ )
4739 if ( (vpos != i) && (pGetExp( piter, i ) != 0) )
4740 {
4741 WerrorS("The input polynomial must be univariate!");
4742 return TRUE;
4743 }
4744 pIter( piter );
4745 }
4746 }
4747
4748 rootContainer * roots= new rootContainer();
4749 number * pcoeffs= (number *)omAlloc( (deg+1) * sizeof( number ) );
4750 piter= gls;
4751 for ( i= deg; i >= 0; i-- )
4752 {
4753 if ( piter && pTotaldegree(piter) == i )
4754 {
4755 pcoeffs[i]= nCopy( pGetCoeff( piter ) );
4756 //nPrint( pcoeffs[i] );PrintS(" ");
4757 pIter( piter );
4758 }
4759 else
4760 {
4761 pcoeffs[i]= nInit(0);
4762 }
4763 }
4764
4765#ifdef mprDEBUG_PROT
4766 for (i=deg; i >= 0; i--)
4767 {
4768 nPrint( pcoeffs[i] );PrintS(" ");
4769 }
4770 PrintLn();
4771#endif
4772
4773 roots->fillContainer( pcoeffs, NULL, 1, deg, rootContainer::onepoly, 1 );
4774 roots->solver( howclean );
4775
4776 int elem= roots->getAnzRoots();
4777 char *dummy;
4778 int j;
4779
4780 lists rlist;
4781 rlist= (lists)omAlloc( sizeof(slists) );
4782 rlist->Init( elem );
4783
4785 {
4786 for ( j= 0; j < elem; j++ )
4787 {
4788 rlist->m[j].rtyp=NUMBER_CMD;
4789 rlist->m[j].data=(void *)nCopy((number)(roots->getRoot(j)));
4790 //rlist->m[j].data=(void *)(number)(roots->getRoot(j));
4791 }
4792 }
4793 else
4794 {
4795 for ( j= 0; j < elem; j++ )
4796 {
4797 dummy = complexToStr( (*roots)[j], gmp_output_digits, currRing->cf );
4798 rlist->m[j].rtyp=STRING_CMD;
4799 rlist->m[j].data=(void *)dummy;
4800 }
4801 }
4802
4803 elist->Clean();
4804 //omFreeSize( (ADDRESS) elist, sizeof(slists) );
4805
4806 // this is (via fillContainer) the same data as in root
4807 //for ( i= deg; i >= 0; i-- ) nDelete( &pcoeffs[i] );
4808 //omFreeSize( (ADDRESS) pcoeffs, (deg+1) * sizeof( number ) );
4809
4810 delete roots;
4811
4812 res->data= (void*)rlist;
4813
4814 return FALSE;
4815}
int i
Definition cfEzgcd.cc:132
int * Zp_roots(poly p, const ring r)
Definition clapsing.cc:2191
complex root finder for univariate polynomials based on laguers algorithm
Definition mpr_numeric.h:66
gmp_complex * getRoot(const int i)
Definition mpr_numeric.h:88
void fillContainer(number *_coeffs, number *_ievpoint, const int _var, const int _tdg, const rootType _rt, const int _anz)
int getAnzRoots()
Definition mpr_numeric.h:97
bool solver(const int polishmode=PM_NONE)
void Clean(ring r=currRing)
Definition lists.h:26
static FORCE_INLINE number n_Init(long i, const coeffs r)
a number representing i in the given coeff field/ring r
Definition coeffs.h:539
int j
Definition facHensel.cc:110
@ NUMBER_CMD
Definition grammar.cc:289
#define pIter(p)
Definition monomials.h:37
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition monomials.h:44
EXTERN_VAR size_t gmp_output_digits
Definition mpr_base.h:115
char * complexToStr(gmp_complex &c, const unsigned int oprec, const coeffs src)
void setGMPFloatDigits(size_t digits, size_t rest)
Set size of mantissa digits - the number of output digits (basis 10) the size of mantissa consists of...
#define nCopy(n)
Definition numbers.h:15
#define nPrint(a)
only for debug, over any initialized currRing
Definition numbers.h:46
#define nInit(i)
Definition numbers.h:24
#define omFree(addr)
#define NULL
Definition omList.c:12
static long pTotaldegree(poly p)
Definition polys.h:283
#define pIsConstant(p)
like above, except that Comp must be 0
Definition polys.h:239
#define pGetExp(p, i)
Exponent.
Definition polys.h:42
void PrintS(const char *s)
Definition reporter.cc:284
void PrintLn()
Definition reporter.cc:310
static BOOLEAN rField_is_R(const ring r)
Definition ring.h:524
static BOOLEAN rField_is_Zp(const ring r)
Definition ring.h:506
static BOOLEAN rField_is_long_C(const ring r)
Definition ring.h:551
static BOOLEAN rField_is_Q(const ring r)
Definition ring.h:512
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition ring.h:598
@ LIST_CMD
Definition tok.h:118
@ STRING_CMD
Definition tok.h:187

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

4652{
4653 ideal gls = (ideal)(arg1->Data());
4654 int imtype= (int)(long)arg2->Data();
4655
4656 uResultant::resMatType mtype= determineMType( imtype );
4657
4658 // check input ideal ( = polynomial system )
4659 if ( mprIdealCheck( gls, arg1->Name(), mtype, true ) != mprOk )
4660 {
4661 return TRUE;
4662 }
4663
4664 uResultant *resMat= new uResultant( gls, mtype, false );
4665 if (resMat!=NULL)
4666 {
4667 res->rtyp = MODUL_CMD;
4668 res->data= (void*)resMat->accessResMat()->getMatrix();
4669 if (!errorreported) delete resMat;
4670 }
4671 return errorreported;
4672}
virtual ideal getMatrix()
Definition mpr_base.h:31
const char * Name()
Definition subexpr.h:120
Base class for solving 0-dim poly systems using u-resultant.
Definition mpr_base.h:63
resMatrixBase * accessResMat()
Definition mpr_base.h:78
VAR short errorreported
Definition feFopen.cc:23
@ MODUL_CMD
Definition grammar.cc:288
@ mprOk
Definition mpr_base.h:98
uResultant::resMatType determineMType(int imtype)
mprState mprIdealCheck(const ideal theIdeal, const char *name, uResultant::resMatType mtype, BOOLEAN rmatrix=false)

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

4919{
4920 leftv v= args;
4921
4922 ideal gls;
4923 int imtype;
4924 int howclean;
4925
4926 // get ideal
4927 if ( v->Typ() != IDEAL_CMD )
4928 return TRUE;
4929 else gls= (ideal)(v->Data());
4930 v= v->next;
4931
4932 // get resultant matrix type to use (0,1)
4933 if ( v->Typ() != INT_CMD )
4934 return TRUE;
4935 else imtype= (int)(long)v->Data();
4936 v= v->next;
4937
4938 if (imtype==0)
4939 {
4940 ideal test_id=idInit(1,1);
4941 int j;
4942 for(j=IDELEMS(gls)-1;j>=0;j--)
4943 {
4944 if (gls->m[j]!=NULL)
4945 {
4946 test_id->m[0]=gls->m[j];
4947 intvec *dummy_w=id_QHomWeight(test_id, currRing);
4948 if (dummy_w!=NULL)
4949 {
4950 WerrorS("Newton polytope not of expected dimension");
4951 delete dummy_w;
4952 return TRUE;
4953 }
4954 }
4955 }
4956 }
4957
4958 // get and set precision in digits ( > 0 )
4959 if ( v->Typ() != INT_CMD )
4960 return TRUE;
4961 else if ( !(rField_is_R(currRing) || rField_is_long_R(currRing) || \
4963 {
4964 unsigned long int ii=(unsigned long int)v->Data();
4965 setGMPFloatDigits( ii, ii );
4966 }
4967 v= v->next;
4968
4969 // get interpolation steps (0,1,2)
4970 if ( v->Typ() != INT_CMD )
4971 return TRUE;
4972 else howclean= (int)(long)v->Data();
4973
4974 uResultant::resMatType mtype= determineMType( imtype );
4975 int i,count;
4976 lists listofroots= NULL;
4977 number smv= NULL;
4978 BOOLEAN interpolate_det= (mtype==uResultant::denseResMat)?TRUE:FALSE;
4979
4980 //emptylist= (lists)omAlloc( sizeof(slists) );
4981 //emptylist->Init( 0 );
4982
4983 //res->rtyp = LIST_CMD;
4984 //res->data= (void *)emptylist;
4985
4986 // check input ideal ( = polynomial system )
4987 if ( mprIdealCheck( gls, args->Name(), mtype ) != mprOk )
4988 {
4989 return TRUE;
4990 }
4991
4992 uResultant * ures;
4993 rootContainer ** iproots;
4994 rootContainer ** muiproots;
4995 rootArranger * arranger;
4996
4997 // main task 1: setup of resultant matrix
4998 ures= new uResultant( gls, mtype );
4999 if ( ures->accessResMat()->initState() != resMatrixBase::ready )
5000 {
5001 WerrorS("Error occurred during matrix setup!");
5002 return TRUE;
5003 }
5004
5005 // if dense resultant, check if minor nonsingular
5006 if ( mtype == uResultant::denseResMat )
5007 {
5008 smv= ures->accessResMat()->getSubDet();
5009#ifdef mprDEBUG_PROT
5010 PrintS("// Determinant of submatrix: ");nPrint(smv);PrintLn();
5011#endif
5012 if ( nIsZero(smv) )
5013 {
5014 WerrorS("Unsuitable input ideal: Minor of resultant matrix is singular!");
5015 return TRUE;
5016 }
5017 }
5018
5019 // main task 2: Interpolate specialized resultant polynomials
5020 if ( interpolate_det )
5021 iproots= ures->interpolateDenseSP( false, smv );
5022 else
5023 iproots= ures->specializeInU( false, smv );
5024
5025 // main task 3: Interpolate specialized resultant polynomials
5026 if ( interpolate_det )
5027 muiproots= ures->interpolateDenseSP( true, smv );
5028 else
5029 muiproots= ures->specializeInU( true, smv );
5030
5031#ifdef mprDEBUG_PROT
5032 int c= iproots[0]->getAnzElems();
5033 for (i=0; i < c; i++) pWrite(iproots[i]->getPoly());
5034 c= muiproots[0]->getAnzElems();
5035 for (i=0; i < c; i++) pWrite(muiproots[i]->getPoly());
5036#endif
5037
5038 // main task 4: Compute roots of specialized polys and match them up
5039 arranger= new rootArranger( iproots, muiproots, howclean );
5040 arranger->solve_all();
5041
5042 // get list of roots
5043 if ( arranger->success() )
5044 {
5045 arranger->arrange();
5046 listofroots= listOfRoots(arranger, gmp_output_digits );
5047 }
5048 else
5049 {
5050 WerrorS("Solver was unable to find any roots!");
5051 return TRUE;
5052 }
5053
5054 // free everything
5055 count= iproots[0]->getAnzElems();
5056 for (i=0; i < count; i++) delete iproots[i];
5057 omFreeSize( (ADDRESS) iproots, count * sizeof(rootContainer*) );
5058 count= muiproots[0]->getAnzElems();
5059 for (i=0; i < count; i++) delete muiproots[i];
5060 omFreeSize( (ADDRESS) muiproots, count * sizeof(rootContainer*) );
5061
5062 delete ures;
5063 delete arranger;
5064 if (smv!=NULL) nDelete( &smv );
5065
5066 res->data= (void *)listofroots;
5067
5068 //emptylist->Clean();
5069 // omFreeSize( (ADDRESS) emptylist, sizeof(slists) );
5070
5071 return FALSE;
5072}
int BOOLEAN
Definition auxiliary.h:88
void * ADDRESS
Definition auxiliary.h:120
virtual number getSubDet()
Definition mpr_base.h:37
virtual IStateType initState() const
Definition mpr_base.h:41
int getAnzElems()
Definition mpr_numeric.h:95
rootContainer ** specializeInU(BOOLEAN matchUp=false, const number subDetVal=NULL)
Definition mpr_base.cc:3060
rootContainer ** interpolateDenseSP(BOOLEAN matchUp=false, const number subDetVal=NULL)
Definition mpr_base.cc:2922
@ denseResMat
Definition mpr_base.h:65
@ IDEAL_CMD
Definition grammar.cc:285
lists listOfRoots(rootArranger *self, const unsigned int oprec)
Definition ipshell.cc:5075
#define nDelete(n)
Definition numbers.h:16
#define nIsZero(n)
Definition numbers.h:19
#define omFreeSize(addr, size)
void pWrite(poly p)
Definition polys.h:309
int status int void size_t count
Definition si_signals.h:69
ideal idInit(int idsize, int rank)
initialise an ideal / module
intvec * id_QHomWeight(ideal id, const ring r)
#define IDELEMS(i)

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

4818{
4819 int i;
4820 ideal p,w;
4821 p= (ideal)arg1->Data();
4822 w= (ideal)arg2->Data();
4823
4824 // w[0] = f(p^0)
4825 // w[1] = f(p^1)
4826 // ...
4827 // p can be a vector of numbers (multivariate polynom)
4828 // or one number (univariate polynom)
4829 // tdg = deg(f)
4830
4831 int n= IDELEMS( p );
4832 int m= IDELEMS( w );
4833 int tdg= (int)(long)arg3->Data();
4834
4835 res->data= (void*)NULL;
4836
4837 // check the input
4838 if ( tdg < 1 )
4839 {
4840 WerrorS("Last input parameter must be > 0!");
4841 return TRUE;
4842 }
4843 if ( n != rVar(currRing) )
4844 {
4845 Werror("Size of first input ideal must be equal to %d!",rVar(currRing));
4846 return TRUE;
4847 }
4848 if ( m != (int)pow((double)tdg+1,(double)n) )
4849 {
4850 Werror("Size of second input ideal must be equal to %d!",
4851 (int)pow((double)tdg+1,(double)n));
4852 return TRUE;
4853 }
4854 if ( !(rField_is_Q(currRing) /* ||
4855 rField_is_R() || rField_is_long_R() ||
4856 rField_is_long_C()*/ ) )
4857 {
4858 WerrorS("Ground field not implemented!");
4859 return TRUE;
4860 }
4861
4862 number tmp;
4863 number *pevpoint= (number *)omAlloc( n * sizeof( number ) );
4864 for ( i= 0; i < n; i++ )
4865 {
4866 pevpoint[i]=nInit(0);
4867 if ( (p->m)[i] )
4868 {
4869 tmp = pGetCoeff( (p->m)[i] );
4870 if ( nIsZero(tmp) || nIsOne(tmp) || nIsMOne(tmp) )
4871 {
4872 omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4873 WerrorS("Elements of first input ideal must not be equal to -1, 0, 1!");
4874 return TRUE;
4875 }
4876 } else tmp= NULL;
4877 if ( !nIsZero(tmp) )
4878 {
4879 if ( !pIsConstant((p->m)[i]))
4880 {
4881 omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4882 WerrorS("Elements of first input ideal must be numbers!");
4883 return TRUE;
4884 }
4885 pevpoint[i]= nCopy( tmp );
4886 }
4887 }
4888
4889 number *wresults= (number *)omAlloc( m * sizeof( number ) );
4890 for ( i= 0; i < m; i++ )
4891 {
4892 wresults[i]= nInit(0);
4893 if ( (w->m)[i] && !nIsZero(pGetCoeff((w->m)[i])) )
4894 {
4895 if ( !pIsConstant((w->m)[i]))
4896 {
4897 omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4898 omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4899 WerrorS("Elements of second input ideal must be numbers!");
4900 return TRUE;
4901 }
4902 wresults[i]= nCopy(pGetCoeff((w->m)[i]));
4903 }
4904 }
4905
4906 vandermonde vm( m, n, tdg, pevpoint, FALSE );
4907 number *ncpoly= vm.interpolateDense( wresults );
4908 // do not free ncpoly[]!!
4909 poly rpoly= vm.numvec2poly( ncpoly );
4910
4911 omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4912 omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4913
4914 res->data= (void*)rpoly;
4915 return FALSE;
4916}
Rational pow(const Rational &a, int e)
Definition GMPrat.cc:411
int p
Definition cfModGcd.cc:4086
vandermonde system solver for interpolating polynomials from their values
Definition mpr_numeric.h:29
const CanonicalForm & w
Definition facAbsFact.cc:51
#define nIsMOne(n)
Definition numbers.h:26
#define nIsOne(n)
Definition numbers.h:25
void Werror(const char *fmt,...)
Definition reporter.cc:189