source: git/Singular/fglmvec.cc @ 634dab0

spielwiese
Last change on this file since 634dab0 was 958e16, checked in by Hans Schönemann <hannes@…>, 23 years ago
*hannes: merged number-ops git-svn-id: file:///usr/local/Singular/svn/trunk@5035 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 11.8 KB
RevLine 
[0e1846]1// emacs edit mode for this file is -*- C++ -*-
[958e16]2// $Id: fglmvec.cc,v 1.15 2001-01-09 15:40:05 Singular Exp $
[8bafbf0]3
4/****************************************
5*  Computer Algebra System SINGULAR     *
6****************************************/
[5812c69]7/*
[8bafbf0]8* ABSTRACT - The FGLM-Algorithm
9*   Implementation of number-vectors for the fglm algorithm.
10*   (See fglm.cc). Based on a letter-envelope implementation, mainly
[5812c69]11*   written to be used by the fglm algorithm. Hence they are
[8bafbf0]12*   specialized for this purpose.
13*/
14
15#include "mod2.h"
[34ab5de]16
17#ifdef HAVE_FGLM
[512a2b]18#include "omalloc.h"
[5812c69]19#include "tok.h"
[8bafbf0]20#include "structs.h"
[0e1846]21#include "numbers.h"
[8bafbf0]22#include "fglm.h"
23#include "fglmvec.h"
[0e1846]24
[d37f27]25#define PROT(msg)
26#define STICKYPROT(msg) if (BTEST1(OPT_PROT)) Print(msg)
27#define PROT2(msg,arg)
28#define STICKYPROT2(msg,arg) if (BTEST1(OPT_PROT)) Print(msg,arg)
29#define fglmASSERT(ignore1,ignore2)
30
[8bafbf0]31class fglmVectorRep
[0e1846]32{
33private:
[8bafbf0]34    int ref_count;
[0e1846]35    int N;
36    number * elems;
37public:
[8bafbf0]38    fglmVectorRep() : ref_count(1), N(0), elems(0) {}
39    fglmVectorRep( int n, number * e ) : ref_count(1), N(n), elems(e) {}
[5812c69]40    fglmVectorRep( int n ) : ref_count(1), N(n)
[0e1846]41    {
[5812c69]42        fglmASSERT( N >= 0, "illegal Vector representation" );
43        if ( N == 0 )
44            elems= 0;
45        else {
[c232af]46            elems= (number *)omAlloc( N*sizeof( number ) );
[5812c69]47            for ( int i= N-1; i >= 0; i-- )
48                elems[i]= nInit( 0 );
49        }
[0e1846]50    }
51    ~fglmVectorRep()
52    {
[5812c69]53        if ( N > 0 ) {
54            for ( int i= N-1; i >= 0; i-- )
55                nDelete( elems + i );
[c232af]56            omFreeSize( (ADDRESS)elems, N*sizeof( number ) );
[5812c69]57        }
[0e1846]58    }
59
60    fglmVectorRep* clone() const
61    {
[5812c69]62        if ( N > 0 ) {
63            number * elems_clone;
[c232af]64            elems_clone= (number *)omAlloc( N*sizeof( number ) );
[5812c69]65            for ( int i= N-1; i >= 0; i-- )
66                elems_clone[i] = nCopy( elems[i] );
67            return new fglmVectorRep( N, elems_clone );
68        } else
69            return new fglmVectorRep( N, 0 );
[0e1846]70    }
[8bafbf0]71    BOOLEAN deleteObject() { return --ref_count == 0; }
72    fglmVectorRep * copyObject() { ref_count++; return this; }
73    int refcount() const { return ref_count; }
74    BOOLEAN isUnique() const { return ref_count == 1; }
75
[0e1846]76    int size() const { return N; }
[5812c69]77    int isZero() const
[8bafbf0]78    {
[5812c69]79        int k;
80        for ( k= N; k > 0; k-- )
81            if ( ! nIsZero( getconstelem( k ) ) )
82                return 0;
83        return 1;
[8bafbf0]84    }
[5812c69]85    int numNonZeroElems() const
[0e1846]86    {
[5812c69]87        int num = 0;
88        int k;
89        for ( k= N; k > 0; k-- )
90            if ( ! nIsZero( getconstelem( k ) ) ) num++;
91        return num;
[0e1846]92    }
93    void setelem( int i, number n )
94    {
[5812c69]95        fglmASSERT( 0 < i && i <= N, "setelem: wrong index" );
96        nDelete( elems + i-1 );
97        elems[i-1]= n;
[0e1846]98    }
[5812c69]99    number ejectelem( int i, number n )
[8bafbf0]100    {
[5812c69]101        fglmASSERT( isUnique(), "should only be called if unique!" );
102        number temp= elems[i-1];
103        elems[i-1]= n;
104        return temp;
[8bafbf0]105    }
[0e1846]106    number & getelem( int i )
107    {
[5812c69]108        fglmASSERT( 0 < i && i <= N, "getelem: wrong index" );
109        return elems[i-1];
[0e1846]110    }
111    const number getconstelem( int i) const
112    {
[5812c69]113        fglmASSERT( 0 < i && i <= N, "getconstelem: wrong index" );
114        return elems[i-1];
[0e1846]115    }
116    friend class fglmVector;
117};
118
119
120///--------------------------------------------------------------------------------
121/// Implementation of class fglmVector
122///--------------------------------------------------------------------------------
123
124fglmVector::fglmVector( fglmVectorRep * r ) : rep(r) {}
125
126fglmVector::fglmVector() : rep( new fglmVectorRep() ) {}
127
128fglmVector::fglmVector( int size ) : rep( new fglmVectorRep( size ) ) {}
129
130fglmVector::fglmVector( int size, int basis ) : rep( new fglmVectorRep( size ) )
131{
132    rep->setelem( basis, nInit( 1 ) );
133}
134
135fglmVector::fglmVector( const fglmVector & v )
136{
[8bafbf0]137    rep= v.rep->copyObject();
[0e1846]138}
139
140fglmVector::~fglmVector()
141{
142    if ( rep->deleteObject() )
[5812c69]143        delete rep;
[0e1846]144}
[b9624d6]145
[c7db8e]146#ifndef HAVE_EXPLICIT_CONSTR
[0e584f]147void
148fglmVector::mac_constr( const fglmVector & v)
149{
150    rep= v.rep->copyObject();
151}
[b9624d6]152
[0e584f]153void
154fglmVector::mac_constr_i( int size )
155{
156    rep= new fglmVectorRep( size );
157}
[b9624d6]158
159void
[5812c69]160fglmVector::clearelems()
[b9624d6]161{
162    if ( rep->deleteObject() )
163      delete rep;
164}
[0e584f]165#endif
[b9624d6]166
[5812c69]167void
[0e1846]168fglmVector::makeUnique()
169{
170    if ( rep->refcount() != 1 ) {
[5812c69]171        rep->deleteObject();
172        rep= rep->clone();
[0e1846]173    }
174}
175
[5812c69]176int
[0e1846]177fglmVector::size() const
178{
179    return rep->size();
180}
181
182int
183fglmVector::numNonZeroElems() const
184{
185    return rep->numNonZeroElems();
186}
187
[8bafbf0]188void
[5812c69]189fglmVector::nihilate( const number fac1, const number fac2, const fglmVector v )
[8bafbf0]190{
191    int i;
192    int vsize= v.size();
193    number term1, term2;
194    fglmASSERT( vsize <= rep->size(), "v has to be smaller oder equal" );
195    if ( rep->isUnique() ) {
[5812c69]196        for ( i= vsize; i > 0; i-- ) {
197            term1= nMult( fac1, rep->getconstelem( i ) );
198            term2= nMult( fac2, v.rep->getconstelem( i ) );
199            rep->setelem( i, nSub( term1, term2 ) );
200            nDelete( &term1 );
201            nDelete( &term2 );
202        }
203        for ( i= rep->size(); i > vsize; i-- ) {
204            rep->setelem( i, nMult( fac1, rep->getconstelem( i ) ) );
205        }
[8bafbf0]206    }
[5812c69]207    else
[8bafbf0]208    {
[5812c69]209        number* newelems;
[c232af]210        newelems= (number *)omAlloc( rep->size()*sizeof( number ) );
[5812c69]211        for ( i= vsize; i > 0; i-- ) {
212            term1= nMult( fac1, rep->getconstelem( i ) );
213            term2= nMult( fac2, v.rep->getconstelem( i ) );
214            newelems[i-1]= nSub( term1, term2 );
215            nDelete( &term1 );
216            nDelete( &term2 );
217        }
218        for ( i= rep->size(); i > vsize; i-- ) {
219            newelems[i-1]= nMult( fac1, rep->getconstelem( i ) );
220        }
221        rep->deleteObject();
222        rep= new fglmVectorRep( rep->size(), newelems );
[8bafbf0]223    }
224}
225
[5812c69]226fglmVector &
[0e1846]227fglmVector::operator = ( const fglmVector & v )
228{
229    if ( this != &v ) {
[5812c69]230        if ( rep->deleteObject() )
231            delete rep;
232        rep = v.rep->copyObject();
[0e1846]233    }
234    return *this;
235}
236
237int
238fglmVector::operator == ( const fglmVector & v )
239{
240    if ( rep->size() == v.rep->size() ) {
[5812c69]241        if ( rep == v.rep ) return 1;
242        else {
243            int i;
244            for ( i= rep->size(); i > 0; i-- )
245                if ( ! nEqual( rep->getconstelem( i ), v.rep->getconstelem( i ) ) )
246                    return 0;
247            return 1;
248        }
[0e1846]249    }
250    return 0;
251}
252
253int
254fglmVector::operator != ( const fglmVector & v )
255{
256    return !( *this == v );
257}
258
[5812c69]259int
[0e1846]260fglmVector::isZero()
261{
[8bafbf0]262    return rep->isZero();
263}
264
[5812c69]265int
[8bafbf0]266fglmVector::elemIsZero( int i )
267{
268    return nIsZero( rep->getconstelem( i ) );
[0e1846]269}
270
271fglmVector &
272fglmVector::operator += ( const fglmVector & v )
273{
[8bafbf0]274    fglmASSERT( size() == v.size(), "incompatible vectors" );
[0e1846]275    // ACHTUNG : Das Verhalten hier mit gcd genau ueberpruefen!
276    int i;
277    if ( rep->isUnique() ) {
[5812c69]278        for ( i= rep->size(); i > 0; i-- )
279            rep->setelem( i, nAdd( rep->getconstelem( i ), v.rep->getconstelem( i ) ) );
[0e1846]280    }
[5812c69]281    else
[0e1846]282    {
[5812c69]283        int n = rep->size();
284        number* newelems;
[c232af]285        newelems= (number *)omAlloc( n*sizeof( number ) );
[5812c69]286        for ( i= n; i > 0; i-- )
287            newelems[i-1]= nAdd( rep->getconstelem( i ), v.rep->getconstelem( i ) );
288        rep->deleteObject();
289        rep= new fglmVectorRep( n, newelems );
[0e1846]290    }
291    return *this;
292}
293
294fglmVector &
295fglmVector::operator -= ( const fglmVector & v )
296{
[8bafbf0]297    fglmASSERT( size() == v.size(), "incompatible vectors" );
[0e1846]298    int i;
299    if ( rep->isUnique() ) {
[5812c69]300        for ( i= rep->size(); i > 0; i-- )
301            rep->setelem( i, nSub( rep->getconstelem( i ), v.rep->getconstelem( i ) ) );
[0e1846]302    }
[5812c69]303    else
[0e1846]304    {
[5812c69]305        int n = rep->size();
306        number* newelems;
[c232af]307        newelems= (number *)omAlloc( n*sizeof( number ) );
[5812c69]308        for ( i= n; i > 0; i-- )
309            newelems[i-1]= nSub( rep->getconstelem( i ), v.rep->getconstelem( i ) );
310        rep->deleteObject();
311        rep= new fglmVectorRep( n, newelems );
[0e1846]312    }
313    return *this;
314}
315
316fglmVector &
317fglmVector::operator *= ( const number & n )
318{
319    int s= rep->size();
320    int i;
321    if ( ! rep->isUnique() ) {
[5812c69]322        number * temp;
[c232af]323        temp= (number *)omAlloc( s*sizeof( number ) );
[5812c69]324        for ( i= s; i > 0; i-- )
325            temp[i-1]= nMult( rep->getconstelem( i ), n );
326        rep->deleteObject();
327        rep= new fglmVectorRep( s, temp );
[0e1846]328    }
[5812c69]329    else
[0e1846]330    {
[5812c69]331        for (i= s; i > 0; i-- )
332            rep->setelem( i, nMult( rep->getconstelem( i ), n ) );
[0e1846]333    }
334    return *this;
335}
336
337fglmVector &
338fglmVector::operator /= ( const number & n )
339{
340    int s= rep->size();
341    int i;
342    if ( ! rep->isUnique() ) {
[5812c69]343        number * temp;
[c232af]344        temp= (number *)omAlloc( s*sizeof( number ) );
[5812c69]345        for ( i= s; i > 0; i-- ) {
346            temp[i-1]= nDiv( rep->getconstelem( i ), n );
347            nNormalize( temp[i-1] );
348        }
349        rep->deleteObject();
350        rep= new fglmVectorRep( s, temp );
[0e1846]351    }
[5812c69]352    else
[0e1846]353    {
[5812c69]354        for (i= s; i > 0; i-- ) {
355            rep->setelem( i, nDiv( rep->getconstelem( i ), n ) );
356            nNormalize( rep->getelem( i ) );
357        }
[0e1846]358    }
359    return *this;
360}
361
[5812c69]362fglmVector
363operator - ( const fglmVector & v )
[0e1846]364{
365    fglmVector temp( v.size() );
366    int i;
367    number n ;
368    for ( i= v.size(); i > 0; i-- ) {
[5812c69]369        n= nCopy( v.getconstelem( i ) );
370        n= nNeg( n );
371        temp.setelem( i, n );
[0e1846]372    }
373    return temp;
374}
375
[5812c69]376fglmVector
377operator + ( const fglmVector & lhs, const fglmVector & rhs )
[0e1846]378{
379    fglmVector temp= lhs;
380    temp+= rhs;
381    return temp;
382}
383
[5812c69]384fglmVector
[0e1846]385operator - ( const fglmVector & lhs, const fglmVector & rhs )
386{
387    fglmVector temp= lhs;
388    temp-= rhs;
389    return temp;
390}
391
[5812c69]392fglmVector
[0e1846]393operator * ( const fglmVector & v, const number n )
394{
395    fglmVector temp= v;
396    temp*= n;
397    return temp;
398}
399
400fglmVector
401operator * ( const number n, const fglmVector & v )
402{
403    fglmVector temp= v;
404    temp*= n;
405    return temp;
406}
407
408number &
409fglmVector::getelem( int i )
410{
411    makeUnique();
412    return rep->getelem( i );
413}
414
[5812c69]415const number
[0e1846]416fglmVector::getconstelem( int i ) const
417{
418    return rep->getconstelem( i );
419}
420
421void
422fglmVector::setelem( int i, number & n )
423{
424    makeUnique();
425    rep->setelem( i, n );
426    n= nNULL;
427}
428
[5812c69]429number
[0e1846]430fglmVector::gcd() const
[5812c69]431{
[0e1846]432    int i= rep->size();
[8bafbf0]433    BOOLEAN found = FALSE;
434    BOOLEAN gcdIsOne = FALSE;
[0e1846]435    number theGcd;
436    number current;
437    while( i > 0 && ! found ) {
[5812c69]438        current= rep->getconstelem( i );
439        if ( ! nIsZero( current ) ) {
440            theGcd= nCopy( current );
441            found= TRUE;
442            if ( ! nGreaterZero( theGcd ) ) {
443                theGcd= nNeg( theGcd );
444            }
445            if ( nIsOne( theGcd ) ) gcdIsOne= TRUE;
446        }
447        i--;
[0e1846]448    }
449    if ( found ) {
[5812c69]450        while ( i > 0 && ! gcdIsOne ) {
451            current= rep->getconstelem( i );
452            if ( ! nIsZero( current ) ) {
[958e16]453                number temp= nGcd( theGcd, current, currRing );
[5812c69]454                nDelete( &theGcd );
455                theGcd= temp;
456                if ( nIsOne( theGcd ) ) gcdIsOne= TRUE;
457            }
458            i--;
459        }
[0e1846]460    }
[5812c69]461    else
462        theGcd= nInit( 0 );
[0e1846]463    return theGcd;
464}
465
466number
[5812c69]467fglmVector::clearDenom()
[0e1846]468{
469    number theLcm = nInit( 1 );
470    number current;
[8bafbf0]471    BOOLEAN isZero = TRUE;
[0e1846]472    int i;
473    for ( i= size(); i > 0; i-- ) {
[5812c69]474        if ( ! nIsZero( rep->getconstelem(i) ) ) {
475            isZero= FALSE;
[958e16]476            number temp= nLcm( theLcm, rep->getconstelem( i ), currRing );
[5812c69]477            nDelete( &theLcm );
478            theLcm= temp;
479        }
[0e1846]480    }
481    if ( isZero ) {
[5812c69]482        nDelete( &theLcm );
483        theLcm= nInit( 0 );
[0e1846]484    }
485    else {
[5812c69]486        if ( ! nIsOne( theLcm ) ) {
487            *this *= theLcm;
488            for ( i= size(); i > 0; i-- ) {
489                nNormalize( rep->getelem( i ) );
490            }
491        }
[0e1846]492    }
493    return theLcm;
494}
495
[34ab5de]496#endif
[8bafbf0]497// ----------------------------------------------------------------------------
498// Local Variables: ***
499// compile-command: "make Singular" ***
500// page-delimiter: "^\\(\\|//!\\)" ***
501// fold-internal-margins: nil ***
502// End: ***
Note: See TracBrowser for help on using the repository browser.