mpr_inout.cc
Go to the documentation of this file.
1 /****************************************
2 * Computer Algebra System SINGULAR *
3 ****************************************/
4 
5 
6 /*
7 * ABSTRACT - multipolynomial resultant
8 */
9 
10 
11 #include <kernel/mod2.h>
12 #include <misc/auxiliary.h>
13 
14 //#ifdef HAVE_MPR
15 
16 //-> includes
17 #include <omalloc/omalloc.h>
18 
19 #include <misc/mylimits.h>
20 #include <misc/options.h>
21 #include <misc/intvec.h>
22 
23 #include <coeffs/numbers.h>
24 #include <coeffs/mpr_global.h>
25 
26 #include <polys/matpol.h>
27 
28 
29 #include <kernel/structs.h>
30 #include <kernel/polys.h>
31 #include <kernel/ideals.h>
32 
33 
34 #include <math.h>
35 
36 #include "mpr_base.h"
37 #include "mpr_numeric.h"
38 #include "mpr_inout.h"
39 
40 // to get detailed timigs, define MPR_TIMING
41 #ifdef MPR_TIMING
42 #define TIMING
43 #endif
44 #include <factory/timing.h>
45 TIMING_DEFINE_PRINT(mpr_overall)
46 TIMING_DEFINE_PRINT(mpr_check)
47 TIMING_DEFINE_PRINT(mpr_constr)
48 TIMING_DEFINE_PRINT(mpr_ures)
49 TIMING_DEFINE_PRINT(mpr_mures)
50 TIMING_DEFINE_PRINT(mpr_arrange)
51 TIMING_DEFINE_PRINT(mpr_solver)
52 
53 #define TIMING_EPR(t,msg) TIMING_END_AND_PRINT(t,msg);TIMING_RESET(t);
54 
55 //<-
56 
57 //------------------------------------------------------------------------------
58 
59 //-> void mprPrintError( mprState state )
60 void mprPrintError( mprState state, const char * name )
61 {
62  switch (state)
63  {
64  case mprWrongRType:
65  WerrorS("Unknown chosen resultant matrix type!");
66  break;
67  case mprHasOne:
68  Werror("One element of the ideal %s is constant!",name);
69  break;
70  case mprInfNumOfVars:
71  Werror("Wrong number of elements in given ideal %s, should be %d resp. %d!",
72  name,(currRing->N)+1,(currRing->N));
73  break;
74  case mprNotZeroDim:
75  Werror("The given ideal %s must be 0-dimensional!",name);
76  break;
77  case mprNotHomog:
78  Werror("The given ideal %s has to be homogeneous in the first ring variable!",
79  name);
80  break;
81  case mprNotReduced:
82  Werror("The given ideal %s has to reduced!",name);
83  break;
84  case mprUnSupField:
85  WerrorS("Ground field not implemented!");
86  break;
87  default:
88  break;
89  }
90 }
91 //<-
92 
93 //-> mprState mprIdealCheck()
94 mprState mprIdealCheck( const ideal theIdeal,
95  const char * /*name*/,
97  BOOLEAN rmatrix )
98 {
99  mprState state = mprOk;
100  // int power;
101  int k;
102 
103  int numOfVars= mtype == uResultant::denseResMat?(currRing->N)-1:(currRing->N);
104  if ( rmatrix ) numOfVars++;
105 
106  if ( mtype == uResultant::none )
107  state= mprWrongRType;
108 
109  if ( IDELEMS(theIdeal) != numOfVars )
110  state= mprInfNumOfVars;
111 
112  for ( k= IDELEMS(theIdeal) - 1; (state == mprOk) && (k >= 0); k-- )
113  {
114  poly p = (theIdeal->m)[k];
115  if ( pIsConstant(p) ) state= mprHasOne;
116  else
117  if ( (mtype == uResultant::denseResMat) && !p_IsHomogeneous(p, currRing) )
118  state=mprNotHomog;
119  }
120 
121  if ( !(rField_is_R(currRing)||
125  (rmatrix && rField_is_Q_a(currRing))) )
126  state= mprUnSupField;
127 
128  if ( state != mprOk ) mprPrintError( state, "" /* name */ );
129 
130  return state;
131 }
132 //<-
133 
134 //-> uResultant::resMatType determineMType( int imtype )
136 {
137  switch ( imtype )
138  {
139  case MPR_DENSE:
141  case 0:
142  case MPR_SPARSE:
144  default:
145  return uResultant::none;
146  }
147 }
148 //<-
149 
150 //-> function u_resultant_det
151 poly u_resultant_det( ideal gls, int imtype )
152 {
153  uResultant::resMatType mtype= determineMType( imtype );
154  poly resdet;
155  poly emptypoly= pInit();
156  number smv= NULL;
157 
158  TIMING_START(mpr_overall);
159 
160  // check input ideal ( = polynomial system )
161  if ( mprIdealCheck( gls, "", mtype ) != mprOk )
162  {
163  return emptypoly;
164  }
165 
166  uResultant *ures;
167 
168  // main task 1: setup of resultant matrix
169  TIMING_START(mpr_constr);
170  ures= new uResultant( gls, mtype );
171  TIMING_EPR(mpr_constr,"construction");
172 
173  // if dense resultant, check if minor nonsingular
174  if ( mtype == uResultant::denseResMat )
175  {
176  smv= ures->accessResMat()->getSubDet();
177 #ifdef mprDEBUG_PROT
178  PrintS("// Determinant of submatrix: ");nPrint(smv); PrintLn();
179 #endif
180  if ( nIsZero(smv) )
181  {
182  WerrorS("Unsuitable input ideal: Minor of resultant matrix is singular!");
183  return emptypoly;
184  }
185  }
186 
187  // main task 2: Interpolate resultant polynomial
188  TIMING_START(mpr_ures);
189  resdet= ures->interpolateDense( smv );
190  TIMING_EPR(mpr_ures,"ures");
191 
192  // free mem
193  delete ures;
194  nDelete( &smv );
195  pDelete( &emptypoly );
196 
197  TIMING_EPR(mpr_overall,"overall");
198 
199  return ( resdet );
200 }
201 //<-
202 
203 //-----------------------------------------------------------------------------
204 
205 //#endif // HAVE_MPR
206 
207 // local Variables: ***
208 // folded-file: t ***
209 // compile-command-1: "make installg" ***
210 // compile-command-2: "make install" ***
211 // End: ***
212 
213 // in folding: C-c x
214 // leave fold: C-c y
215 // foldmode: F10
void PrintLn()
Definition: reporter.cc:327
Base class for solving 0-dim poly systems using u-resultant.
Definition: mpr_base.h:62
#define MPR_DENSE
Definition: mpr_inout.h:15
Compatiblity layer for legacy polynomial operations (over currRing)
return P p
Definition: myNF.cc:203
BOOLEAN p_IsHomogeneous(poly p, const ring r)
Definition: p_polys.cc:3187
static BOOLEAN rField_is_R(const ring r)
Definition: ring.h:467
resMatrixBase * accessResMat()
Definition: mpr_base.h:78
static BOOLEAN rField_is_Q_a(const ring r)
Definition: ring.h:488
uResultant::resMatType determineMType(int imtype)
void WerrorS(const char *s)
Definition: feFopen.cc:24
int k
Definition: cfEzgcd.cc:93
Definition: mpr_base.h:98
#define nPrint(a)
only for debug, over any initalized currRing
Definition: numbers.h:46
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:12
#define MPR_SPARSE
Definition: mpr_inout.h:16
All the auxiliary stuff.
#define pIsConstant(p)
like above, except that Comp might be != 0
Definition: polys.h:209
#define TIMING_EPR(t, msg)
void PrintS(const char *s)
Definition: reporter.cc:294
static BOOLEAN rField_is_Q(const ring r)
Definition: ring.h:461
#define IDELEMS(i)
Definition: simpleideals.h:24
mprState mprIdealCheck(const ideal theIdeal, const char *name, uResultant::resMatType mtype, BOOLEAN rmatrix=false)
#define nDelete(n)
Definition: numbers.h:16
static BOOLEAN rField_is_long_C(const ring r)
Definition: ring.h:494
TIMING_DEFINE_PRINT(mpr_overall) TIMING_DEFINE_PRINT(mpr_check) TIMING_DEFINE_PRINT(mpr_constr) TIMING_DEFINE_PRINT(mpr_ures) TIMING_DEFINE_PRINT(mpr_mures) TIMING_DEFINE_PRINT(mpr_arrange) TIMING_DEFINE_PRINT(mpr_solver) void mprPrintError(mprState state
char name(const Variable &v)
Definition: variable.h:95
#define nIsZero(n)
Definition: numbers.h:19
#define NULL
Definition: omList.c:10
#define TIMING_START(t)
Definition: timing.h:87
static BOOLEAN rField_is_long_R(const ring r)
Definition: ring.h:491
#define pInit()
allocates a new monomial and initializes everything to 0
Definition: polys.h:61
poly interpolateDense(const number subDetVal=NULL)
Definition: mpr_base.cc:2770
#define pDelete(p_ptr)
Definition: polys.h:157
mprState
Definition: mpr_base.h:96
polyrec * poly
Definition: hilb.h:10
int BOOLEAN
Definition: auxiliary.h:131
void Werror(const char *fmt,...)
Definition: reporter.cc:199
virtual number getSubDet()
Definition: mpr_base.h:37