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