source: git/kernel/fglmvec.cc @ abe5c8

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