source: git/kernel/fglmvec.cc @ fbc7cb

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