source: git/kernel/npolygon.cc @ abe5c8

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