My Project
Loading...
Searching...
No Matches
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
26
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
139ostream & 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
172int 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
189Rational 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
205Rational 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
229Rational 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
245Rational 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
262Rational 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
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
497ostream & 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
564Rational 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
585Rational 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
606Rational 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
628Rational 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
651int 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
661Rational* newtonPolygon::sqh_weights( void ) const
662{
663 return ( N==1 ? l[0].c : (Rational*)NULL );
664}
665
666int 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