source: git/kernel/fglmvec.cc @ bc669c

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