source: git/kernel/npolygon.cc @ fbc7cb

jengelh-datetimespielwiese
Last change on this file since fbc7cb was ba5e9e, checked in by Oleksandr Motsak <motsak@…>, 9 years ago
Changed configure-scripts to generate individual public config files for each package: resources, libpolys, singular (main) fix: sources should include correct corresponding config headers.
  • 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#ifdef HAVE_CONFIG_H
11#include "singularconfig.h"
12#endif /* HAVE_CONFIG_H */
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>
25#include <polys/monomials/p_polys.h>
26#include <polys/monomials/ring.h>
27
28#include <kernel/GMPrat.h>
29#include <kernel/kmatrix.h>
30#include <kernel/npolygon.h>
31#include <kernel/structs.h>
32
33// ----------------------------------------------------------------------------
34//  Allocate memory for a linear form of  k  terms
35// ----------------------------------------------------------------------------
36
37void    linearForm::copy_new( int k )
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
83void    linearForm::copy_delete( void )
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
94void    linearForm::copy_deep( const linearForm &l )
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
108linearForm::linearForm( const linearForm &l )
109{
110  copy_deep( l );
111}
112
113// ----------------------------------------------------------------------------
114//  Destructor
115// ----------------------------------------------------------------------------
116
117linearForm::~linearForm( )
118{
119  copy_delete( );
120}
121
122// ----------------------------------------------------------------------------
123//  Define `=` to be a deep copy
124// ----------------------------------------------------------------------------
125
126linearForm & linearForm::operator = ( const linearForm &l )
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
140ostream & 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
173int     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
190Rational 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
206Rational 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
230Rational linearForm::weight_shift( poly m, const ring r ) const
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
246Rational 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
263Rational linearForm::weight_shift1( poly m, const ring r ) const
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
280int linearForm::positive( void )
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
297void    newtonPolygon::copy_new( int k )
298{
299  if( k > 0 )
300  {
301    l = new linearForm[k];
302
303    #ifndef 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
343void    newtonPolygon::copy_delete( void )
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
354void    newtonPolygon::copy_deep( const newtonPolygon &np )
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
368newtonPolygon::newtonPolygon( const newtonPolygon &np )
369{
370  copy_deep( np );
371}
372
373// ----------------------------------------------------------------------------
374//  Destructor
375// ----------------------------------------------------------------------------
376
377newtonPolygon::~newtonPolygon( )
378{
379  copy_delete( );
380}
381
382// ----------------------------------------------------------------------------
383//  We define `=` to be a deep copy
384// ----------------------------------------------------------------------------
385
386newtonPolygon & newtonPolygon::operator = ( const newtonPolygon &np )
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
398template class KMatrix<Rational>;
399
400newtonPolygon::newtonPolygon( poly f, const ring s )
401{
402  copy_zero( );
403
404  int *r=new int[s->N];
405  poly *m=new poly[s->N];
406
407
408  KMatrix<Rational> mat(s->N,s->N+1 );
409
410  int i,j,stop=FALSE;
411  linearForm sol;
412
413  // ---------------
414  //  init counters
415  // ---------------
416
417  for( i=0; i<s->N; i++ )
418  {
419    r[i] = i;
420  }
421
422  m[0] = f;
423
424  for( i=1; i<s->N; i++ )
425  {
426    m[i] = pNext(m[i-1]);
427  }
428
429  // -----------------------------
430  //  find faces (= linear forms)
431  // -----------------------------
432
433  do
434  {
435    // ---------------------------------------------------
436    //  test if monomials p.m[r[0]]m,...,p.m[r[p.vars-1]]
437    //  are linearely independent
438    // ---------------------------------------------------
439
440    for( i=0; i<s->N; i++ )
441    {
442      for( j=0; j<s->N; j++ )
443      {
444        //    mat.set( i,j,pGetExp( m[r[i]],j+1 ) );
445        mat.set( i,j,p_GetExp( m[i],j+1,s ) );
446      }
447      mat.set( i,j,1 );
448    }
449
450    if( mat.solve( &(sol.c),&(sol.N) ) == s->N )
451    {
452      // ---------------------------------
453      //  check if linearForm is positive
454      //  check if linearForm is extremal
455      // ---------------------------------
456
457      if( sol.positive( ) && sol.pweight( f,s ) >= (Rational)1 )
458      {
459        // ----------------------------------
460        //  this is a face or the polyhedron
461        // ----------------------------------
462
463        add_linearForm( sol );
464        sol.c = (Rational*)NULL;
465        sol.N = 0;
466      }
467    }
468
469    // --------------------
470    //  increment counters
471    // --------------------
472
473    for( i=1; r[i-1] + 1 == r[i] && i < s->N; i++ );
474
475    for( j=0; j<i-1; j++ )
476    {
477      r[j]=j;
478    }
479
480    if( i>1 )
481    {
482      m[0]=f;
483      for( j=1; j<i-1; j++ )
484      {
485        m[j]=pNext(m[j-1]);
486      }
487    }
488    r[i-1]++;
489    m[i-1]=pNext(m[i-1]);
490
491    if( m[s->N-1] == (poly)NULL )
492    {
493      stop = TRUE;
494    }
495  } while( stop == FALSE );
496}
497
498#ifdef  NPOLYGON_PRINT
499
500ostream & operator << ( ostream &s,const newtonPolygon &a )
501{
502  #ifdef NPOLYGON_IOSTREAM
503    s << "Newton polygon:" << endl;
504  #else
505    fprintf( stdout,"Newton polygon:\n" );
506  #endif
507
508  for( int i=0; i<a.N; i++ )
509  {
510    s << a.l[i];
511
512    #ifdef NPOLYGON_IOSTREAM
513      s << endl;
514    #else
515      fprintf( stdout,"\n" );
516    #endif
517  }
518
519  return s;
520}
521
522#endif
523
524// ----------------------------------------------------------------------------
525//  Add a face to the newton polygon
526// ----------------------------------------------------------------------------
527
528void    newtonPolygon::add_linearForm( const linearForm &l0 )
529{
530  int           i;
531  newtonPolygon np;
532
533  // -------------------------------------
534  //  test if linear form is already here
535  // -------------------------------------
536
537  for( i=0; i<N; i++ )
538  {
539    if( l0==l[i] )
540    {
541      return;
542    }
543  }
544
545  np.copy_new( N+1 );
546  np.N = N+1;
547
548  for( i=0; i<N; i++ )
549  {
550    np.l[i].copy_shallow( l[i] );
551    l[i].copy_zero( );
552  }
553
554  np.l[N] = l0;
555
556  copy_delete( );
557  copy_shallow( np );
558  np.copy_zero( );
559
560  return;
561}
562
563// ----------------------------------------------------------------------------
564//  Newton weight of a monomial
565// ----------------------------------------------------------------------------
566
567Rational  newtonPolygon::weight( poly m, const ring r ) const
568{
569  Rational ret = l[0].weight( m,r );
570  Rational tmp;
571
572  for( int i=1; i<N; i++ )
573  {
574    tmp = l[i].weight( m,r );
575
576    if( tmp < ret )
577    {
578      ret = tmp;
579    }
580  }
581  return ret;
582}
583
584// ----------------------------------------------------------------------------
585//  Newton weight of a monomial shifted by the product of the variables
586// ----------------------------------------------------------------------------
587
588Rational  newtonPolygon::weight_shift( poly m, const ring r ) const
589{
590  Rational ret = l[0].weight_shift( m, r );
591  Rational tmp;
592
593  for( int i=1; i<N; i++ )
594  {
595    tmp = l[i].weight_shift( m, r );
596
597    if( tmp < ret )
598    {
599      ret = tmp;
600    }
601  }
602  return ret;
603}
604
605// ----------------------------------------------------------------------------
606//  Newton weight of a monomial (forgetting the first variable)
607// ----------------------------------------------------------------------------
608
609Rational  newtonPolygon::weight1( poly m, const ring r ) const
610{
611  Rational ret = l[0].weight1( m, r );
612  Rational tmp;
613
614  for( int i=1; i<N; i++ )
615  {
616    tmp = l[i].weight1( m, r );
617
618    if( tmp < ret )
619    {
620      ret = tmp;
621    }
622  }
623  return ret;
624}
625
626// ----------------------------------------------------------------------------
627//  Newton weight of a monomial shifted by the product of the variables
628//  (forgetting the first variable)
629// ----------------------------------------------------------------------------
630
631Rational  newtonPolygon::weight_shift1( poly m, const ring r ) const
632{
633  Rational ret = l[0].weight_shift1( m, r );
634  Rational tmp;
635
636  for( int i=1; i<N; i++ )
637  {
638    tmp = l[i].weight_shift1( m, r );
639
640    if( tmp < ret )
641    {
642      ret = tmp;
643    }
644  }
645  return ret;
646}
647
648/*
649// ----------------------------------------------------------------------------
650//  Chcek if the polynomial belonging to the Newton polygon
651//  is semiquasihomogeneous (sqh)
652// ----------------------------------------------------------------------------
653
654int     newtonPolygon::is_sqh( void ) const
655{
656  return ( N==1 ? TRUE : FALSE );
657}
658
659// ----------------------------------------------------------------------------
660//  If the polynomial belonging to the Newton polygon is sqh,
661//  return its weights vector
662// ----------------------------------------------------------------------------
663
664Rational*   newtonPolygon::sqh_weights( void ) const
665{
666  return ( N==1 ? l[0].c : (Rational*)NULL );
667}
668
669int    newtonPolygon::sqh_N( void ) const
670{
671  return ( N==1 ? l[0].N : 0 );
672}
673*/
674
675#endif /* HAVE_SPECTRUM */
676// ----------------------------------------------------------------------------
677//  npolygon.cc
678//  end of file
679// ----------------------------------------------------------------------------
680
Note: See TracBrowser for help on using the repository browser.