source: git/kernel/fglm/fglmvec.cc @ 4f8fd1d

fieker-DuValspielwiese
Last change on this file since 4f8fd1d was 4f8fd1d, checked in by Frédéric Chapoton <chapoton@…>, 20 months ago
trying to add a codespell linter for kernel/
  • Property mode set to 100644
File size: 10.8 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
15
16
17#include "kernel/mod2.h"
18
19#include "kernel/structs.h"
20#include "coeffs/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  {
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    {
105      if(!nIsZero (getconstelem (k)))
106          return 0;
107    }
108    return 1;
109  }
110  int numNonZeroElems () const
111  {
112    int num = 0;
113    int k;
114    for(k = N; k > 0; k--)
115    {
116      if(!nIsZero (getconstelem (k)))
117          num++;
118    }
119    return num;
120  }
121  void setelem (int i, number n)
122  {
123    fglmASSERT (0 < i && i <= N, "setelem: wrong index");
124    nDelete (elems + i - 1);
125    elems[i - 1] = n;
126  }
127  number ejectelem (int i, number n)
128  {
129    fglmASSERT (isUnique (), "should only be called if unique!");
130    number temp = elems[i - 1];
131    elems[i - 1] = n;
132    return temp;
133  }
134  number & getelem (int i)
135  {
136    fglmASSERT (0 < i && i <= N, "getelem: wrong index");
137    return elems[i - 1];
138  }
139  number getconstelem (int i) const
140  {
141    fglmASSERT (0 < i && i <= N, "getconstelem: wrong index");
142    return elems[i - 1];
143  }
144  friend class fglmVector;
145};
146
147
148///--------------------------------------------------------------------------------
149/// Implementation of class fglmVector
150///--------------------------------------------------------------------------------
151
152fglmVector::fglmVector (fglmVectorRep * r):rep (r)
153{
154}
155
156fglmVector::fglmVector ():rep (new fglmVectorRep ())
157{
158}
159
160fglmVector::fglmVector (int size):rep (new fglmVectorRep (size))
161{
162}
163
164fglmVector::fglmVector (int size, int basis):rep (new fglmVectorRep (size))
165{
166  rep->setelem (basis, nInit (1));
167}
168
169fglmVector::fglmVector (const fglmVector & v)
170{
171  rep = v.rep->copyObject ();
172}
173
174fglmVector::~fglmVector ()
175{
176  if(rep->deleteObject ())
177    delete rep;
178}
179
180#ifndef HAVE_EXPLICIT_CONSTR
181void fglmVector::mac_constr (const fglmVector & v)
182{
183  rep = v.rep->copyObject ();
184}
185
186void fglmVector::mac_constr_i (int size)
187{
188  rep = new fglmVectorRep (size);
189}
190
191void fglmVector::clearelems ()
192{
193  if(rep->deleteObject ())
194    delete rep;
195}
196#endif
197
198void fglmVector::makeUnique ()
199{
200  if(rep->refcount () != 1)
201  {
202    rep->deleteObject ();
203    rep = rep->clone ();
204  }
205}
206
207int fglmVector::size () const
208{
209  return rep->size ();
210}
211
212int fglmVector::numNonZeroElems () const
213{
214  return rep->numNonZeroElems ();
215}
216
217void
218  fglmVector::nihilate (const number fac1, const number fac2,
219                        const fglmVector v)
220{
221  int i;
222  int vsize = v.size ();
223  number term1, term2;
224  fglmASSERT (vsize <= rep->size (), "v has to be smaller or equal");
225  if(rep->isUnique ())
226  {
227    for(i = vsize; i > 0; i--)
228    {
229      term1 = nMult (fac1, rep->getconstelem (i));
230      term2 = nMult (fac2, v.rep->getconstelem (i));
231      rep->setelem (i, nSub (term1, term2));
232      nDelete (&term1);
233      nDelete (&term2);
234    }
235    for(i = rep->size (); i > vsize; i--)
236    {
237      rep->setelem (i, nMult (fac1, rep->getconstelem (i)));
238    }
239  }
240  else
241  {
242    number *newelems;
243    newelems = (number *) omAlloc (rep->size () * sizeof (number));
244    for(i = vsize; i > 0; i--)
245    {
246      term1 = nMult (fac1, rep->getconstelem (i));
247      term2 = nMult (fac2, v.rep->getconstelem (i));
248      newelems[i - 1] = nSub (term1, term2);
249      nDelete (&term1);
250      nDelete (&term2);
251    }
252    for(i = rep->size (); i > vsize; i--)
253    {
254      newelems[i - 1] = nMult (fac1, rep->getconstelem (i));
255    }
256    rep->deleteObject ();
257    rep = new fglmVectorRep (rep->size (), newelems);
258  }
259}
260
261fglmVector & fglmVector::operator = (const fglmVector & v)
262{
263  if(this != &v)
264  {
265    if(rep->deleteObject ())
266      delete rep;
267    rep = v.rep->copyObject ();
268  }
269  return *this;
270}
271
272int fglmVector::operator == (const fglmVector & v)
273{
274  if(rep->size () == v.rep->size ())
275  {
276    if(rep == v.rep)
277      return 1;
278    else
279    {
280      int i;
281      for(i = rep->size (); i > 0; i--)
282        if(!nEqual (rep->getconstelem (i), v.rep->getconstelem (i)))
283          return 0;
284      return 1;
285    }
286  }
287  return 0;
288}
289
290int fglmVector::operator != (const fglmVector & v)
291{
292  return !(*this == v);
293}
294
295int fglmVector::isZero ()
296{
297  return rep->isZero ();
298}
299
300int fglmVector::elemIsZero (int i)
301{
302  return nIsZero (rep->getconstelem (i));
303}
304
305fglmVector & fglmVector::operator += (const fglmVector & v)
306{
307  fglmASSERT (size () == v.size (), "incompatible vectors");
308  // ACHTUNG : Das Verhalten hier mit gcd genau ueberpruefen!
309  int i;
310  if(rep->isUnique ())
311  {
312    for(i = rep->size (); i > 0; i--)
313      rep->setelem (i, nAdd (rep->getconstelem (i), v.rep->getconstelem (i)));
314  }
315  else
316  {
317    int n = rep->size ();
318    number *newelems;
319    newelems = (number *) omAlloc (n * sizeof (number));
320    for(i = n; i > 0; i--)
321      newelems[i - 1] = nAdd (rep->getconstelem (i), v.rep->getconstelem (i));
322    rep->deleteObject ();
323    rep = new fglmVectorRep (n, newelems);
324  }
325  return *this;
326}
327
328fglmVector & fglmVector::operator -= (const fglmVector & v)
329{
330  fglmASSERT (size () == v.size (), "incompatible vectors");
331  int i;
332  if(rep->isUnique ())
333  {
334    for(i = rep->size (); i > 0; i--)
335      rep->setelem (i, nSub (rep->getconstelem (i), v.rep->getconstelem (i)));
336  }
337  else
338  {
339    int n = rep->size ();
340    number *newelems;
341    newelems = (number *) omAlloc (n * sizeof (number));
342    for(i = n; i > 0; i--)
343      newelems[i - 1] = nSub (rep->getconstelem (i), v.rep->getconstelem (i));
344    rep->deleteObject ();
345    rep = new fglmVectorRep (n, newelems);
346  }
347  return *this;
348}
349
350fglmVector & fglmVector::operator *= (const number & n)
351{
352  int s = rep->size ();
353  int i;
354  if(!rep->isUnique ())
355  {
356    number *temp;
357    temp = (number *) omAlloc (s * sizeof (number));
358    for(i = s; i > 0; i--)
359      temp[i - 1] = nMult (rep->getconstelem (i), n);
360    rep->deleteObject ();
361    rep = new fglmVectorRep (s, temp);
362  }
363  else
364  {
365    for(i = s; i > 0; i--)
366      rep->setelem (i, nMult (rep->getconstelem (i), n));
367  }
368  return *this;
369}
370
371fglmVector & fglmVector::operator /= (const number & n)
372{
373  int s = rep->size ();
374  int i;
375  if(!rep->isUnique ())
376  {
377    number *temp;
378    temp = (number *) omAlloc (s * sizeof (number));
379    for(i = s; i > 0; i--)
380    {
381      temp[i - 1] = nDiv (rep->getconstelem (i), n);
382      nNormalize (temp[i - 1]);
383    }
384    rep->deleteObject ();
385    rep = new fglmVectorRep (s, temp);
386  }
387  else
388  {
389    for(i = s; i > 0; i--)
390    {
391      rep->setelem (i, nDiv (rep->getconstelem (i), n));
392      nNormalize (rep->getelem (i));
393    }
394  }
395  return *this;
396}
397
398fglmVector operator - (const fglmVector & v)
399{
400  fglmVector temp (v.size ());
401  int i;
402  number n;
403  for(i = v.size (); i > 0; i--)
404  {
405    n = nCopy (v.getconstelem (i));
406    n = nInpNeg (n);
407    temp.setelem (i, n);
408  }
409  return temp;
410}
411
412fglmVector operator + (const fglmVector & lhs, const fglmVector & rhs)
413{
414  fglmVector temp = lhs;
415  temp += rhs;
416  return temp;
417}
418
419fglmVector operator - (const fglmVector & lhs, const fglmVector & rhs)
420{
421  fglmVector temp = lhs;
422  temp -= rhs;
423  return temp;
424}
425
426fglmVector operator * (const fglmVector & v, const number n)
427{
428  fglmVector temp = v;
429  temp *= n;
430  return temp;
431}
432
433fglmVector operator * (const number n, const fglmVector & v)
434{
435  fglmVector temp = v;
436  temp *= n;
437  return temp;
438}
439
440number & fglmVector::getelem (int i)
441{
442  makeUnique ();
443  return rep->getelem (i);
444}
445
446number fglmVector::getconstelem (int i) const
447{
448  return rep->getconstelem (i);
449}
450
451void fglmVector::setelem (int i, number & n)
452{
453  makeUnique ();
454  rep->setelem (i, n);
455  n = n_Init (0, currRing->cf);
456}
457
458number fglmVector::gcd () const
459{
460  int i = rep->size ();
461  BOOLEAN found = FALSE;
462  BOOLEAN gcdIsOne = FALSE;
463  number theGcd;
464  number current;
465  while(i > 0 && !found)
466  {
467    current = rep->getconstelem (i);
468    if(!nIsZero (current))
469    {
470      theGcd = nCopy (current);
471      found = TRUE;
472      if(!nGreaterZero (theGcd))
473      {
474        theGcd = nInpNeg (theGcd);
475      }
476      if(nIsOne (theGcd))
477        gcdIsOne = TRUE;
478    }
479    i--;
480  }
481  if(found)
482  {
483    while(i > 0 && !gcdIsOne)
484    {
485      current = rep->getconstelem (i);
486      if(!nIsZero (current))
487      {
488        number temp = n_SubringGcd (theGcd, current, currRing->cf);
489        nDelete (&theGcd);
490        theGcd = temp;
491        if(nIsOne (theGcd))
492          gcdIsOne = TRUE;
493      }
494      i--;
495    }
496  }
497  else
498    theGcd = nInit (0);
499  return theGcd;
500}
501
502number fglmVector::clearDenom ()
503{
504  number theLcm = nInit (1);
505  BOOLEAN isZero = TRUE;
506  int i;
507  for(i = size (); i > 0; i--)
508  {
509    if(!nIsZero (rep->getconstelem (i)))
510    {
511      isZero = FALSE;
512      number temp = n_NormalizeHelper (theLcm, rep->getconstelem (i), currRing->cf);
513      nDelete (&theLcm);
514      theLcm = temp;
515    }
516  }
517  if(isZero)
518  {
519    nDelete (&theLcm);
520    theLcm = nInit (0);
521  }
522  else
523  {
524    if(!nIsOne (theLcm))
525    {
526      *this *= theLcm;
527      for(i = size (); i > 0; i--)
528      {
529        nNormalize (rep->getelem (i));
530      }
531    }
532  }
533  return theLcm;
534}
535
536// ----------------------------------------------------------------------------
537// Local Variables: ***
538// compile-command: "make Singular" ***
539// page-delimiter: "^\\(\\|//!\\)" ***
540// fold-internal-margins: nil ***
541// End: ***
Note: See TracBrowser for help on using the repository browser.