My Project
npolygon.cc
Go to the documentation of this file.
1 // ----------------------------------------------------------------------------
2 // npolygon.cc
3 // begin of file
4 // Stephan Endrass, endrass@mathematik.uni-mainz.de
5 // 23.7.99
6 // ----------------------------------------------------------------------------
7 
8 #define NPOLYGON_CC
9 
10 
11 
12 
13 #include "kernel/mod2.h"
14 
15 #ifdef HAVE_SPECTRUM
16 
17 #ifdef NPOLYGON_PRINT
18 #include <iostream.h>
19 #ifndef NPOLYGON_IOSTREAM
20 #include <stdio.h>
21 #endif
22 #endif
23 
25 #include "polys/monomials/ring.h"
26 
27 #include "kernel/spectrum/GMPrat.h"
30 #include "kernel/structs.h"
31 
32 // ----------------------------------------------------------------------------
33 // Allocate memory for a linear form of k terms
34 // ----------------------------------------------------------------------------
35 
37 {
38  if( k > 0 )
39  {
40  c = new Rational[k];
41 
42  #ifndef NBDEBUG
43  if( c == (Rational*)NULL )
44  {
45  #ifdef NPOLYGON_PRINT
46  #ifdef NPOLYGON_IOSTREAM
47  cerr <<
48  "void linearForm::copy_new( int k ): no memory left ...\n" ;
49  #else
50  fprintf( stderr,
51  "void linearForm::copy_new( int k ): no memory left ...\n");
52  #endif
53  #endif
54  HALT();
55  }
56  #endif
57  }
58  else if( k == 0 )
59  {
60  c = (Rational*)NULL;
61  }
62  else if( k < 0 )
63  {
64  #ifdef NPOLYGON_PRINT
65  #ifdef NPOLYGON_IOSTREAM
66  cerr <<
67  "void linearForm::copy_new( int k ): k < 0 ...\n";
68  #else
69  fprintf( stderr,
70  "void linearForm::copy_new( int k ): k < 0 ...\n" );
71  #endif
72  #endif
73 
74  HALT();
75  }
76 }
77 
78 // ----------------------------------------------------------------------------
79 // Delete the memory of a linear form
80 // ----------------------------------------------------------------------------
81 
83 {
84  if( c != (Rational*)NULL && N > 0 )
85  delete [] c;
86  copy_zero( );
87 }
88 
89 // ----------------------------------------------------------------------------
90 // Initialize deep from another linear form
91 // ----------------------------------------------------------------------------
92 
94 {
95  copy_new( l.N );
96  for( int i=l.N-1; i>=0; i-- )
97  {
98  c[i] = l.c[i];
99  }
100  N = l.N;
101 }
102 
103 // ----------------------------------------------------------------------------
104 // Copy constructor
105 // ----------------------------------------------------------------------------
106 
108 {
109  copy_deep( l );
110 }
111 
112 // ----------------------------------------------------------------------------
113 // Destructor
114 // ----------------------------------------------------------------------------
115 
117 {
118  copy_delete( );
119 }
120 
121 // ----------------------------------------------------------------------------
122 // Define `=` to be a deep copy
123 // ----------------------------------------------------------------------------
124 
126 {
127  copy_delete( );
128  copy_deep( l );
129 
130  return *this;
131 }
132 
133 // ----------------------------------------------------------------------------
134 // ostream for linear form
135 // ----------------------------------------------------------------------------
136 
137 #ifdef NPOLYGON_PRINT
138 
139 ostream & operator << ( ostream &s,const linearForm &l )
140 {
141  for( int i=0; i<l.N; i++ )
142  {
143  if( l.c[i] != (Rational)0 )
144  {
145  if( i > 0 && l.c[i] >= (Rational)0 )
146  {
147  #ifdef NPOLYGON_IOSTREAM
148  s << "+";
149  #else
150  fprintf( stdout,"+" );
151  #endif
152  }
153 
154  s << l.c[i];
155 
156  #ifdef NPOLYGON_IOSTREAM
157  s << "*x" << i+1;
158  #else
159  fprintf( stdout,"*x%d",i+1 );
160  #endif
161  }
162  }
163  return s;
164 }
165 
166 #endif
167 
168 // ----------------------------------------------------------------------------
169 // Equality of linear forms
170 // ----------------------------------------------------------------------------
171 
172 int operator == ( const linearForm &l1,const linearForm &l2 )
173 {
174  if( l1.N!=l2.N )
175  return FALSE;
176  for( int i=l1.N-1; i >=0 ; i-- )
177  {
178  if( l1.c[i]!=l2.c[i] )
179  return FALSE;
180  }
181  return TRUE;
182 }
183 
184 
185 // ----------------------------------------------------------------------------
186 // Newton weight of a monomial
187 // ----------------------------------------------------------------------------
188 
189 Rational linearForm::weight( poly m, const ring r ) const
190 {
191  Rational ret=(Rational)0;
192 
193  for( int i=0,j=1; i<N; i++,j++ )
194  {
195  ret += c[i]*(Rational)p_GetExp( m,j,r );
196  }
197 
198  return ret;
199 }
200 
201 // ----------------------------------------------------------------------------
202 // Newton weight of a polynomial
203 // ----------------------------------------------------------------------------
204 
205 Rational linearForm::pweight( poly m, const ring r ) const
206 {
207  if( m==(poly)NULL )
208  return (Rational)0;
209 
210  Rational ret = weight( m, r );
211  Rational tmp;
212 
213  for( m=pNext(m); m!=(poly)NULL; pIter(m) )
214  {
215  tmp = weight( m, r );
216  if( tmp<ret )
217  {
218  ret = tmp;
219  }
220  }
221 
222  return ret;
223 }
224 
225 // ----------------------------------------------------------------------------
226 // Newton weight of a monomial shifted by the product of the variables
227 // ----------------------------------------------------------------------------
228 
229 Rational linearForm::weight_shift( poly m, const ring r ) const
230 {
231  Rational ret=(Rational)0;
232 
233  for( int i=0,j=1; i<N; i++,j++ )
234  {
235  ret += c[i]*(Rational)( p_GetExp( m,j,r ) + 1 );
236  }
237 
238  return ret;
239 }
240 
241 // ----------------------------------------------------------------------------
242 // Newton weight of a monomial (forgetting the first variable)
243 // ----------------------------------------------------------------------------
244 
245 Rational linearForm::weight1( poly m, const ring r ) const
246 {
247  Rational ret=(Rational)0;
248 
249  for( int i=0,j=2; i<N; i++,j++ )
250  {
251  ret += c[i]*(Rational)p_GetExp( m,j,r );
252  }
253 
254  return ret;
255 }
256 
257 // ----------------------------------------------------------------------------
258 // Newton weight of a monomial shifted by the product of the variables
259 // (forgetting the first variable)
260 // ----------------------------------------------------------------------------
261 
262 Rational linearForm::weight_shift1( poly m, const ring r ) const
263 {
264  Rational ret=(Rational)0;
265 
266  for( int i=0,j=2; i<N; i++,j++ )
267  {
268  ret += c[i]*(Rational)( p_GetExp( m,j,r ) + 1 );
269  }
270 
271  return ret;
272 }
273 
274 
275 // ----------------------------------------------------------------------------
276 // Test if all coefficients of a linear form are positive
277 // ----------------------------------------------------------------------------
278 
280 {
281  for( int i=0; i<N; i++ )
282  {
283  if( c[i] <= (Rational)0 )
284  {
285  return FALSE;
286  }
287  }
288  return TRUE;
289 }
290 
291 
292 // ----------------------------------------------------------------------------
293 // Allocate memory for a Newton polygon of k linear forms
294 // ----------------------------------------------------------------------------
295 
297 {
298  if( k > 0 )
299  {
300  l = new linearForm[k];
301 
302  #ifndef SING_NDEBUG
303  if( l == (linearForm*)NULL )
304  {
305  #ifdef NPOLYGON_PRINT
306  #ifdef NPOLYGON_IOSTREAM
307  cerr <<
308  "void newtonPolygon::copy_new( int k ): no memory left ...\n";
309  #else
310  fprintf( stderr,
311  "void newtonPolygon::copy_new( int k ): no memory left ...\n" );
312  #endif
313  #endif
314 
315  HALT();
316  }
317  #endif
318  }
319  else if( k == 0 )
320  {
321  l = (linearForm*)NULL;
322  }
323  else if( k < 0 )
324  {
325  #ifdef NPOLYGON_PRINT
326  #ifdef NPOLYGON_IOSTREAM
327  cerr << "void newtonPolygon::copy_new( int k ): k < 0 ...\n";
328  #else
329  fprintf( stderr,
330  "void newtonPolygon::copy_new( int k ): k < 0 ...\n" );
331  #endif
332  #endif
333 
334  HALT();
335  }
336 }
337 
338 // ----------------------------------------------------------------------------
339 // Delete the memory of a Newton polygon
340 // ----------------------------------------------------------------------------
341 
343 {
344  if( l != (linearForm*)NULL && N > 0 )
345  delete [] l;
346  copy_zero( );
347 }
348 
349 // ----------------------------------------------------------------------------
350 // Initialize deep from another Newton polygon
351 // ----------------------------------------------------------------------------
352 
354 {
355  copy_new( np.N );
356  for( int i=0; i<np.N; i++ )
357  {
358  l[i] = np.l[i];
359  }
360  N = np.N;
361 }
362 
363 // ----------------------------------------------------------------------------
364 // Copy constructor
365 // ----------------------------------------------------------------------------
366 
368 {
369  copy_deep( np );
370 }
371 
372 // ----------------------------------------------------------------------------
373 // Destructor
374 // ----------------------------------------------------------------------------
375 
377 {
378  copy_delete( );
379 }
380 
381 // ----------------------------------------------------------------------------
382 // We define `=` to be a deep copy
383 // ----------------------------------------------------------------------------
384 
386 {
387  copy_delete( );
388  copy_deep( np );
389 
390  return *this;
391 }
392 
393 // ----------------------------------------------------------------------------
394 // Initialize a Newton polygon from a polynomial
395 // ----------------------------------------------------------------------------
396 
397 newtonPolygon::newtonPolygon( poly f, const ring s )
398 {
399  copy_zero( );
400 
401  int *r=new int[s->N];
402  poly *m=new poly[s->N];
403 
404 
405  KMatrix<Rational> mat(s->N,s->N+1 );
406 
407  int i,j,stop=FALSE;
408  linearForm sol;
409 
410  // ---------------
411  // init counters
412  // ---------------
413 
414  for( i=0; i<s->N; i++ )
415  {
416  r[i] = i;
417  }
418 
419  m[0] = f;
420 
421  for( i=1; i<s->N; i++ )
422  {
423  m[i] = pNext(m[i-1]);
424  }
425 
426  // -----------------------------
427  // find faces (= linear forms)
428  // -----------------------------
429 
430  do
431  {
432  // ---------------------------------------------------
433  // test if monomials p.m[r[0]]m,...,p.m[r[p.vars-1]]
434  // are linearely independent
435  // ---------------------------------------------------
436 
437  for( i=0; i<s->N; i++ )
438  {
439  for( j=0; j<s->N; j++ )
440  {
441  // mat.set( i,j,pGetExp( m[r[i]],j+1 ) );
442  mat.set( i,j,p_GetExp( m[i],j+1,s ) );
443  }
444  mat.set( i,j,1 );
445  }
446 
447  if( mat.solve( &(sol.c),&(sol.N) ) == s->N )
448  {
449  // ---------------------------------
450  // check if linearForm is positive
451  // check if linearForm is extremal
452  // ---------------------------------
453 
454  if( sol.positive( ) && sol.pweight( f,s ) >= (Rational)1 )
455  {
456  // ----------------------------------
457  // this is a face or the polyhedron
458  // ----------------------------------
459 
460  add_linearForm( sol );
461  sol.c = (Rational*)NULL;
462  sol.N = 0;
463  }
464  }
465 
466  // --------------------
467  // increment counters
468  // --------------------
469 
470  for( i=1; r[i-1] + 1 == r[i] && i < s->N; i++ );
471 
472  for( j=0; j<i-1; j++ )
473  {
474  r[j]=j;
475  }
476 
477  if( i>1 )
478  {
479  m[0]=f;
480  for( j=1; j<i-1; j++ )
481  {
482  m[j]=pNext(m[j-1]);
483  }
484  }
485  r[i-1]++;
486  m[i-1]=pNext(m[i-1]);
487 
488  if( m[s->N-1] == (poly)NULL )
489  {
490  stop = TRUE;
491  }
492  } while( stop == FALSE );
493 }
494 
495 #ifdef NPOLYGON_PRINT
496 
497 ostream & operator << ( ostream &s,const newtonPolygon &a )
498 {
499  #ifdef NPOLYGON_IOSTREAM
500  s << "Newton polygon:" << endl;
501  #else
502  fprintf( stdout,"Newton polygon:\n" );
503  #endif
504 
505  for( int i=0; i<a.N; i++ )
506  {
507  s << a.l[i];
508 
509  #ifdef NPOLYGON_IOSTREAM
510  s << endl;
511  #else
512  fprintf( stdout,"\n" );
513  #endif
514  }
515 
516  return s;
517 }
518 
519 #endif
520 
521 // ----------------------------------------------------------------------------
522 // Add a face to the newton polygon
523 // ----------------------------------------------------------------------------
524 
526 {
527  int i;
528  newtonPolygon np;
529 
530  // -------------------------------------
531  // test if linear form is already here
532  // -------------------------------------
533 
534  for( i=0; i<N; i++ )
535  {
536  if( l0==l[i] )
537  {
538  return;
539  }
540  }
541 
542  np.copy_new( N+1 );
543  np.N = N+1;
544 
545  for( i=0; i<N; i++ )
546  {
547  np.l[i].copy_shallow( l[i] );
548  l[i].copy_zero( );
549  }
550 
551  np.l[N] = l0;
552 
553  copy_delete( );
554  copy_shallow( np );
555  np.copy_zero( );
556 
557  return;
558 }
559 
560 // ----------------------------------------------------------------------------
561 // Newton weight of a monomial
562 // ----------------------------------------------------------------------------
563 
564 Rational newtonPolygon::weight( poly m, const ring r ) const
565 {
566  Rational ret = l[0].weight( m,r );
567  Rational tmp;
568 
569  for( int i=1; i<N; i++ )
570  {
571  tmp = l[i].weight( m,r );
572 
573  if( tmp < ret )
574  {
575  ret = tmp;
576  }
577  }
578  return ret;
579 }
580 
581 // ----------------------------------------------------------------------------
582 // Newton weight of a monomial shifted by the product of the variables
583 // ----------------------------------------------------------------------------
584 
585 Rational newtonPolygon::weight_shift( poly m, const ring r ) const
586 {
587  Rational ret = l[0].weight_shift( m, r );
588  Rational tmp;
589 
590  for( int i=1; i<N; i++ )
591  {
592  tmp = l[i].weight_shift( m, r );
593 
594  if( tmp < ret )
595  {
596  ret = tmp;
597  }
598  }
599  return ret;
600 }
601 
602 // ----------------------------------------------------------------------------
603 // Newton weight of a monomial (forgetting the first variable)
604 // ----------------------------------------------------------------------------
605 
606 Rational newtonPolygon::weight1( poly m, const ring r ) const
607 {
608  Rational ret = l[0].weight1( m, r );
609  Rational tmp;
610 
611  for( int i=1; i<N; i++ )
612  {
613  tmp = l[i].weight1( m, r );
614 
615  if( tmp < ret )
616  {
617  ret = tmp;
618  }
619  }
620  return ret;
621 }
622 
623 // ----------------------------------------------------------------------------
624 // Newton weight of a monomial shifted by the product of the variables
625 // (forgetting the first variable)
626 // ----------------------------------------------------------------------------
627 
628 Rational newtonPolygon::weight_shift1( poly m, const ring r ) const
629 {
630  Rational ret = l[0].weight_shift1( m, r );
631  Rational tmp;
632 
633  for( int i=1; i<N; i++ )
634  {
635  tmp = l[i].weight_shift1( m, r );
636 
637  if( tmp < ret )
638  {
639  ret = tmp;
640  }
641  }
642  return ret;
643 }
644 
645 /*
646 // ----------------------------------------------------------------------------
647 // Check if the polynomial belonging to the Newton polygon
648 // is semiquasihomogeneous (sqh)
649 // ----------------------------------------------------------------------------
650 
651 int newtonPolygon::is_sqh( void ) const
652 {
653  return ( N==1 ? TRUE : FALSE );
654 }
655 
656 // ----------------------------------------------------------------------------
657 // If the polynomial belonging to the Newton polygon is sqh,
658 // return its weights vector
659 // ----------------------------------------------------------------------------
660 
661 Rational* newtonPolygon::sqh_weights( void ) const
662 {
663  return ( N==1 ? l[0].c : (Rational*)NULL );
664 }
665 
666 int newtonPolygon::sqh_N( void ) const
667 {
668  return ( N==1 ? l[0].N : 0 );
669 }
670 */
671 
672 #endif /* HAVE_SPECTRUM */
673 // ----------------------------------------------------------------------------
674 // npolygon.cc
675 // end of file
676 // ----------------------------------------------------------------------------
677 
#define TRUE
Definition: auxiliary.h:100
#define FALSE
Definition: auxiliary.h:96
int l
Definition: cfEzgcd.cc:100
int m
Definition: cfEzgcd.cc:128
int i
Definition: cfEzgcd.cc:132
int k
Definition: cfEzgcd.cc:99
FILE * f
Definition: checklibs.c:9
void set(int, int, const K &)
Definition: kmatrix.h:354
int solve(K **, int *)
Definition: kmatrix.h:599
Rational weight_shift(poly, const ring r) const
Definition: npolygon.cc:229
Rational weight(poly, const ring r) const
Definition: npolygon.cc:189
void copy_zero(void)
Definition: npolygon.h:109
linearForm()
Definition: npolygon.h:130
Rational weight1(poly, const ring r) const
Definition: npolygon.cc:245
void copy_new(int)
Definition: npolygon.cc:36
Rational * c
Definition: npolygon.h:22
Rational weight_shift1(poly, const ring r) const
Definition: npolygon.cc:262
void copy_delete(void)
Definition: npolygon.cc:82
linearForm & operator=(const linearForm &)
Definition: npolygon.cc:125
Rational pweight(poly, const ring r) const
Definition: npolygon.cc:205
int positive(void)
Definition: npolygon.cc:279
void copy_shallow(linearForm &)
Definition: npolygon.h:119
void copy_deep(const linearForm &)
Definition: npolygon.cc:93
void add_linearForm(const linearForm &)
Definition: npolygon.cc:525
void copy_new(int)
Definition: npolygon.cc:296
void copy_shallow(newtonPolygon &)
Definition: npolygon.h:154
Rational weight(poly, const ring r) const
Definition: npolygon.cc:564
Rational weight_shift1(poly, const ring r) const
Definition: npolygon.cc:628
void copy_delete(void)
Definition: npolygon.cc:342
newtonPolygon & operator=(const newtonPolygon &)
Definition: npolygon.cc:385
linearForm * l
Definition: npolygon.h:66
Rational weight_shift(poly, const ring r) const
Definition: npolygon.cc:585
void copy_zero(void)
Definition: npolygon.h:144
Rational weight1(poly, const ring r) const
Definition: npolygon.cc:606
void copy_deep(const newtonPolygon &)
Definition: npolygon.cc:353
const CanonicalForm int s
Definition: facAbsFact.cc:51
int j
Definition: facHensel.cc:110
static void HALT(void)
Definition: mod2.h:126
#define pIter(p)
Definition: monomials.h:37
#define pNext(p)
Definition: monomials.h:36
int operator==(const linearForm &l1, const linearForm &l2)
Definition: npolygon.cc:172
#define NULL
Definition: omList.c:12
static long p_GetExp(const poly p, const unsigned long iBitmask, const int VarOffset)
get a single variable exponent @Note: the integer VarOffset encodes:
Definition: p_polys.h:467
ostream & operator<<(ostream &s, const spectrum &spec)
Definition: semic.cc:249