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