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
30 class fglmVectorRep
31 {
32 private:
33 int ref_count;
34 int N;
35 number *elems;
36 public:
fglmVectorRep()37 fglmVectorRep ():ref_count (1), N (0), elems (0)
38 {
39 }
fglmVectorRep(int n,number * e)40 fglmVectorRep (int n, number * e):ref_count (1), N (n), elems (e)
41 {
42 }
fglmVectorRep(int n)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 }
~fglmVectorRep()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
clone() const65 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 }
deleteObject()78 BOOLEAN deleteObject ()
79 {
80 return --ref_count == 0;
81 }
copyObject()82 fglmVectorRep *copyObject ()
83 {
84 ref_count++;
85 return this;
86 }
refcount() const87 int refcount () const
88 {
89 return ref_count;
90 }
isUnique() const91 BOOLEAN isUnique () const
92 {
93 return ref_count == 1;
94 }
95
size() const96 int size () const
97 {
98 return N;
99 }
isZero() const100 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 }
numNonZeroElems() const110 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 }
setelem(int i,number n)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 }
ejectelem(int i,number n)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 }
getelem(int i)134 number & getelem (int i)
135 {
136 fglmASSERT (0 < i && i <= N, "getelem: wrong index");
137 return elems[i - 1];
138 }
getconstelem(int i) const139 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
fglmVector(fglmVectorRep * r)152 fglmVector::fglmVector (fglmVectorRep * r):rep (r)
153 {
154 }
155
fglmVector()156 fglmVector::fglmVector ():rep (new fglmVectorRep ())
157 {
158 }
159
fglmVector(int size)160 fglmVector::fglmVector (int size):rep (new fglmVectorRep (size))
161 {
162 }
163
fglmVector(int size,int basis)164 fglmVector::fglmVector (int size, int basis):rep (new fglmVectorRep (size))
165 {
166 rep->setelem (basis, nInit (1));
167 }
168
fglmVector(const fglmVector & v)169 fglmVector::fglmVector (const fglmVector & v)
170 {
171 rep = v.rep->copyObject ();
172 }
173
~fglmVector()174 fglmVector::~fglmVector ()
175 {
176 if(rep->deleteObject ())
177 delete rep;
178 }
179
180 #ifndef HAVE_EXPLICIT_CONSTR
mac_constr(const fglmVector & v)181 void fglmVector::mac_constr (const fglmVector & v)
182 {
183 rep = v.rep->copyObject ();
184 }
185
mac_constr_i(int size)186 void fglmVector::mac_constr_i (int size)
187 {
188 rep = new fglmVectorRep (size);
189 }
190
clearelems()191 void fglmVector::clearelems ()
192 {
193 if(rep->deleteObject ())
194 delete rep;
195 }
196 #endif
197
makeUnique()198 void fglmVector::makeUnique ()
199 {
200 if(rep->refcount () != 1)
201 {
202 rep->deleteObject ();
203 rep = rep->clone ();
204 }
205 }
206
size() const207 int fglmVector::size () const
208 {
209 return rep->size ();
210 }
211
numNonZeroElems() const212 int fglmVector::numNonZeroElems () const
213 {
214 return rep->numNonZeroElems ();
215 }
216
217 void
nihilate(const number fac1,const number fac2,const fglmVector v)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 oder 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
operator =(const fglmVector & v)261 fglmVector & 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
operator ==(const fglmVector & v)272 int 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
operator !=(const fglmVector & v)290 int fglmVector::operator != (const fglmVector & v)
291 {
292 return !(*this == v);
293 }
294
isZero()295 int fglmVector::isZero ()
296 {
297 return rep->isZero ();
298 }
299
elemIsZero(int i)300 int fglmVector::elemIsZero (int i)
301 {
302 return nIsZero (rep->getconstelem (i));
303 }
304
operator +=(const fglmVector & v)305 fglmVector & 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
operator -=(const fglmVector & v)328 fglmVector & 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
operator *=(const number & n)350 fglmVector & 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
operator /=(const number & n)371 fglmVector & 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
operator -(const fglmVector & v)398 fglmVector 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
operator +(const fglmVector & lhs,const fglmVector & rhs)412 fglmVector operator + (const fglmVector & lhs, const fglmVector & rhs)
413 {
414 fglmVector temp = lhs;
415 temp += rhs;
416 return temp;
417 }
418
operator -(const fglmVector & lhs,const fglmVector & rhs)419 fglmVector operator - (const fglmVector & lhs, const fglmVector & rhs)
420 {
421 fglmVector temp = lhs;
422 temp -= rhs;
423 return temp;
424 }
425
operator *(const fglmVector & v,const number n)426 fglmVector operator * (const fglmVector & v, const number n)
427 {
428 fglmVector temp = v;
429 temp *= n;
430 return temp;
431 }
432
operator *(const number n,const fglmVector & v)433 fglmVector operator * (const number n, const fglmVector & v)
434 {
435 fglmVector temp = v;
436 temp *= n;
437 return temp;
438 }
439
getelem(int i)440 number & fglmVector::getelem (int i)
441 {
442 makeUnique ();
443 return rep->getelem (i);
444 }
445
getconstelem(int i) const446 number fglmVector::getconstelem (int i) const
447 {
448 return rep->getconstelem (i);
449 }
450
setelem(int i,number & n)451 void fglmVector::setelem (int i, number & n)
452 {
453 makeUnique ();
454 rep->setelem (i, n);
455 n = n_Init (0, currRing->cf);
456 }
457
gcd() const458 number 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
clearDenom()502 number 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: ***
542