My Project
Loading...
Searching...
No Matches
sparsmat.cc
Go to the documentation of this file.
1/****************************************
2* Computer Algebra System SINGULAR *
3****************************************/
4
5/*
6* ABSTRACT: operations with sparse matrices (bareiss, ...)
7*/
8
9#include "misc/auxiliary.h"
10
11#include "misc/mylimits.h"
12
13#include "misc/options.h"
14
15#include "reporter/reporter.h"
16#include "misc/intvec.h"
17
18#include "coeffs/numbers.h"
19
20#include "monomials/ring.h"
21#include "monomials/p_polys.h"
22
23#include "simpleideals.h"
24
25#include "sparsmat.h"
26#include "prCopy.h"
27
28#include "templates/p_Procs.h"
29
30#include "kbuckets.h"
31#include "operations/p_Mult_q.h"
32
33// define SM_NO_BUCKETS, if sparsemat stuff should not use buckets
34// #define SM_NO_BUCKETS
35
36// this is also influenced by TEST_OPT_NOTBUCKETS
37#ifndef SM_NO_BUCKETS
38// buckets do carry a small additional overhead: only use them if
39// min-length of polys is >= SM_MIN_LENGTH_BUCKET
40#define SM_MIN_LENGTH_BUCKET MIN_LENGTH_BUCKET - 5
41#else
42#define SM_MIN_LENGTH_BUCKET MAX_INT_VAL
43#endif
44
45typedef struct smprec sm_prec;
46typedef sm_prec * smpoly;
47struct smprec
48{
49 smpoly n; // the next element
50 int pos; // position
51 int e; // level
52 poly m; // the element
53 float f; // complexity of the element
54};
55
56
57/* declare internal 'C' stuff */
58static void sm_ExactPolyDiv(poly, poly, const ring);
59static BOOLEAN sm_IsNegQuot(poly, const poly, const poly, const ring);
60static void sm_ExpMultDiv(poly, const poly, const poly, const ring);
61static void sm_PolyDivN(poly, const number, const ring);
62static BOOLEAN smSmaller(poly, poly);
63static void sm_CombineChain(poly *, poly, const ring);
64static void sm_FindRef(poly *, poly *, poly, const ring);
65
66static void sm_ElemDelete(smpoly *, const ring);
68static float sm_PolyWeight(smpoly, const ring);
69static smpoly sm_Poly2Smpoly(poly, const ring);
70static poly sm_Smpoly2Poly(smpoly, const ring);
71static BOOLEAN sm_HaveDenom(poly, const ring);
72static number sm_Cleardenom(ideal, const ring);
73
75
76static poly pp_Mult_Coeff_mm_DivSelect_MultDiv(poly p, int &lp, poly m,
77 poly a, poly b, const ring currRing)
78{
79 if (rOrd_is_Comp_dp(currRing) && currRing->ExpL_Size > 2)
80 {
81 // pp_Mult_Coeff_mm_DivSelectMult only works for (c/C,dp) and
82 // ExpL_Size > 2
83 // should be generalized, at least to dp with ExpL_Size == 2
84 // (is the case for 1 variable)
85 int shorter;
86 p = currRing->p_Procs->pp_Mult_Coeff_mm_DivSelectMult(p, m, a, b,
87 shorter, currRing);
88 lp -= shorter;
89 }
90 else
91 {
94 }
95 return p;
96}
97
98static poly sm_SelectCopy_ExpMultDiv(poly p, poly m, poly a, poly b, const ring currRing)
99{
100 int lp = 0;
102}
103
104
105/* class for sparse matrix:
106* 3 parts of matrix during the algorithm
107* m_act[cols][pos(rows)] => m_row[rows][pos(cols)] => m_res[cols][pos(rows)]
108* input pivotcols as rows result
109* pivot like a stack from pivot and pivotcols
110* elimination rows reordered according
111* to pivot choise
112* stored in perm
113* a step is as follows
114* - search of pivot (smPivot)
115* - swap pivot column and last column (smSwapC)
116* select pivotrow to piv and red (smSelectPR)
117* consider sign
118* - elimination (smInitElim, sm0Elim, sm1Elim)
119* clear zero column as result of elimination (smZeroElim)
120* - tranfer from
121* piv and m_row to m_res (smRowToCol)
122* m_act to m_row (smColToRow)
123*/
125private:
126 int nrows, ncols; // dimension of the problem
127 int sign; // for determinant (start: 1)
128 int act; // number of unreduced columns (start: ncols)
129 int crd; // number of reduced columns (start: 0)
130 int tored; // border for rows to reduce
131 int inred; // unreducable part
132 int rpiv, cpiv; // position of the pivot
133 int normalize; // Normalization flag
134 int *perm; // permutation of rows
135 float wpoints; // weight of all points
136 float *wrw, *wcl; // weights of rows and columns
137 smpoly * m_act; // unreduced columns
138 smpoly * m_res; // reduced columns (result)
139 smpoly * m_row; // reduced part of rows
140 smpoly red; // row to reduce
141 smpoly piv, oldpiv; // pivot and previous pivot
142 smpoly dumm; // allocated dummy
143 ring _R;
144
145 void smColToRow();
146 void smRowToCol();
147 void smFinalMult();
149 void smWeights();
150 void smPivot();
151 void smNewWeights();
152 void smNewPivot();
153 void smZeroElim();
154 void smToredElim();
155 void smCopToRes();
156 void smSelectPR();
157 void sm1Elim();
158 void smHElim();
159 void smMultCol();
160 poly smMultPoly(smpoly);
161 void smActDel();
162 void smColDel();
163 void smPivDel();
164 void smSign();
165 void smInitPerm();
166 int smCheckNormalize();
167 void smNormalize();
168public:
169 sparse_mat(ideal, const ring);
170 ~sparse_mat();
171 int smGetSign() { return sign; }
172 smpoly * smGetAct() { return m_act; }
173 int smGetRed() { return tored; }
174 ideal smRes2Mod();
175 poly smDet();
176 void smNewBareiss(int, int);
177 void smToIntvec(intvec *);
178};
179
180/*
181* estimate maximal exponent for det or minors,
182* the module m has di vectors and maximal rank ra,
183* estimate yields for the tXt minors,
184* we have di,ra >= t
185*/
186static void smMinSelect(long *, int, int);
187
188long sm_ExpBound( ideal m, int di, int ra, int t, const ring currRing)
189{
190 poly p;
191 long kr, kc;
192 long *r, *c;
193 int al, bl, i, j, k;
194
195 if (ra==0) ra=1;
196 al = di*sizeof(long);
197 c = (long *)omAlloc(al);
198 bl = ra*sizeof(long);
199 r = (long *)omAlloc0(bl);
200 for (i=di-1;i>=0;i--)
201 {
202 kc = 0;
203 p = m->m[i];
204 while(p!=NULL)
205 {
206 k = p_GetComp(p, currRing)-1;
207 kr = r[k];
208 for (j=currRing->N;j>0;j--)
209 {
210 long t=p_GetExp(p,j, currRing);
211 if(t /*p_GetExp(p,j, currRing)*/ >kc)
212 kc=t; /*p_GetExp(p,j, currRing);*/
213 if(t /*p_GetExp(p,j, currRing)s*/ >kr)
214 kr=t; /*p_GetExp(p,j, currRing);*/
215 }
216 r[k] = kr;
217 pIter(p);
218 }
219 c[i] = kc;
220 }
221 if (t<di) smMinSelect(c, t, di);
222 if (t<ra) smMinSelect(r, t, ra);
223 kr = kc = 0;
224 for (j=t-1;j>=0;j--)
225 {
226 kr += r[j];
227 kc += c[j];
228 }
229 omFreeSize((ADDRESS)c, al);
230 omFreeSize((ADDRESS)r, bl);
231 if (kr<kc) kc = kr;
232 if (kr<1) kr = 1;
233 return kr;
234}
235
236static void smMinSelect(long *c, int t, int d)
237{
238 long m;
239 int pos, i;
240 do
241 {
242 d--;
243 pos = d;
244 m = c[pos];
245 for (i=d-1;i>=0;i--)
246 {
247 if(c[i]<m)
248 {
249 pos = i;
250 m = c[i];
251 }
252 }
253 for (i=pos;i<d;i++) c[i] = c[i+1];
254 } while (d>t);
255}
256
257/* ----------------- ops with rings ------------------ */
258ring sm_RingChange(const ring origR, long bound)
259{
260// *origR =currRing;
261 ring tmpR=rCopy0(origR,FALSE,FALSE);
263 int *block0=(int*)omAlloc0(3*sizeof(int));
264 int *block1=(int*)omAlloc0(3*sizeof(int));
265 ord[0]=ringorder_c;
266 ord[1]=ringorder_dp;
267 tmpR->order=ord;
268 tmpR->OrdSgn=1;
269 block0[1]=1;
270 tmpR->block0=block0;
271 block1[1]=tmpR->N;
272 tmpR->block1=block1;
273 tmpR->bitmask = 2*bound;
274 tmpR->wvhdl = (int **)omAlloc0((3) * sizeof(int*));
275
276// ???
277// if (tmpR->bitmask > currRing->bitmask) tmpR->bitmask = currRing->bitmask;
278
279 rComplete(tmpR,1);
280 if (origR->qideal!=NULL)
281 {
282 tmpR->qideal= idrCopyR_NoSort(origR->qideal, origR, tmpR);
283 }
284 if (TEST_OPT_PROT)
285 Print("[%ld:%d]", (long) tmpR->bitmask, tmpR->ExpL_Size);
286 return tmpR;
287}
288
290{
291 if (r->qideal!=NULL) id_Delete(&(r->qideal),r);
292 for(int i=r->N-1;i>=0;i--) omFree(r->names[i]);
293 omFreeSize(r->names,r->N*sizeof(char*));
295}
296
297/* ----------------- basics (used from 'C') ------------------ */
298/*2
299*returns the determinant of the module I;
300*uses Bareiss algorithm
301*/
302poly sm_CallDet(ideal I,const ring R)
303{
304 if (I->ncols != I->rank)
305 {
306 Werror("det of %ld x %d module (matrix)",I->rank,I->ncols);
307 return NULL;
308 }
309 int r=id_RankFreeModule(I,R);
310 if (I->ncols != r) // some 0-lines at the end
311 {
312 return NULL;
313 }
314 long bound=sm_ExpBound(I,r,r,r,R);
315 number diag,h=n_Init(1,R->cf);
316 poly res;
317 ring tmpR;
318 sparse_mat *det;
319 ideal II;
320
321 tmpR=sm_RingChange(R,bound);
322 II = idrCopyR(I, R, tmpR);
323 diag = sm_Cleardenom(II,tmpR);
324 det = new sparse_mat(II,tmpR);
325 id_Delete(&II,tmpR);
326 if (det->smGetAct() == NULL)
327 {
328 delete det;
330 return NULL;
331 }
332 res=det->smDet();
333 if(det->smGetSign()<0) res=p_Neg(res,tmpR);
334 delete det;
335 res = prMoveR(res, tmpR, R);
337 if (!n_Equal(diag,h,R->cf))
338 {
339 p_Mult_nn(res,diag,R);
341 }
342 n_Delete(&diag,R->cf);
343 n_Delete(&h,R->cf);
344 return res;
345}
346
347void sm_CallBareiss(ideal I, int x, int y, ideal & M, intvec **iv, const ring R)
348{
349 int r=id_RankFreeModule(I,R),t=r;
350 int c=IDELEMS(I),s=c;
351 long bound;
352 ring tmpR;
353 sparse_mat *bareiss;
354
355 if ((x>0) && (x<t))
356 t-=x;
357 if ((y>1) && (y<s))
358 s-=y;
359 if (t>s) t=s;
360 bound=sm_ExpBound(I,c,r,t,R);
361 tmpR=sm_RingChange(R,bound);
362 ideal II = idrCopyR(I, R, tmpR);
363 bareiss = new sparse_mat(II,tmpR);
364 if (bareiss->smGetAct() == NULL)
365 {
366 delete bareiss;
367 *iv=new intvec(1,rVar(tmpR));
368 }
369 else
370 {
371 id_Delete(&II,tmpR);
372 bareiss->smNewBareiss(x, y);
373 II = bareiss->smRes2Mod();
374 *iv = new intvec(bareiss->smGetRed());
375 bareiss->smToIntvec(*iv);
376 delete bareiss;
377 II = idrMoveR(II,tmpR,R);
378 }
380 M=II;
381}
382
383/*
384* constructor
385*/
386sparse_mat::sparse_mat(ideal smat, const ring RR)
387{
388 int i;
389 poly* pmat;
390 _R=RR;
391
392 ncols = smat->ncols;
393 nrows = id_RankFreeModule(smat,RR);
394 if (nrows <= 0)
395 {
396 m_act = NULL;
397 return;
398 }
399 sign = 1;
400 inred = act = ncols;
401 crd = 0;
402 tored = nrows; // without border
403 i = tored+1;
404 perm = (int *)omAlloc(sizeof(int)*(i+1));
405 perm[i] = 0;
406 m_row = (smpoly *)omAlloc0(sizeof(smpoly)*i);
407 wrw = (float *)omAlloc(sizeof(float)*i);
408 i = ncols+1;
409 wcl = (float *)omAlloc(sizeof(float)*i);
410 m_act = (smpoly *)omAlloc(sizeof(smpoly)*i);
411 m_res = (smpoly *)omAlloc0(sizeof(smpoly)*i);
414 m_res[0]->m = NULL;
415 pmat = smat->m;
416 for(i=ncols; i; i--)
417 {
418 m_act[i] = sm_Poly2Smpoly(pmat[i-1], RR);
419 pmat[i-1] = NULL;
420 }
421 this->smZeroElim();
422 oldpiv = NULL;
423}
424
425/*
426* destructor
427*/
429{
430 int i;
431 if (m_act == NULL) return;
434 i = ncols+1;
435 omFreeSize((ADDRESS)m_res, sizeof(smpoly)*i);
436 omFreeSize((ADDRESS)m_act, sizeof(smpoly)*i);
437 omFreeSize((ADDRESS)wcl, sizeof(float)*i);
438 i = nrows+1;
439 omFreeSize((ADDRESS)wrw, sizeof(float)*i);
440 omFreeSize((ADDRESS)m_row, sizeof(smpoly)*i);
441 omFreeSize((ADDRESS)perm, sizeof(int)*(i+1));
442}
443
444/*
445* transform the result to a module
446*/
447
449{
450 ideal res = idInit(crd, crd);
451 int i;
452
453 for (i=crd; i; i--)
454 {
455 res->m[i-1] = sm_Smpoly2Poly(m_res[i],_R);
456 res->rank=si_max(res->rank, p_MaxComp(res->m[i-1],_R));
457 }
458 return res;
459}
460
461/*
462* permutation of rows
463*/
465{
466 int i;
467
468 for (i=v->rows()-1; i>=0; i--)
469 (*v)[i] = perm[i+1];
470}
471/* ---------------- the algorithm's ------------------ */
472/*
473* the determinant (up to sign), uses new Bareiss elimination
474*/
476{
477 poly res = NULL;
478
479 if (sign == 0)
480 {
481 this->smActDel();
482 return NULL;
483 }
484 if (act < 2)
485 {
486 if (act != 0) res = m_act[1]->m;
487 omFreeBin((void *)m_act[1], smprec_bin);
488 return res;
489 }
490 normalize = 0;
491 this->smInitPerm();
492 this->smPivot();
493 this->smSign();
494 this->smSelectPR();
495 this->sm1Elim();
496 crd++;
497 m_res[crd] = piv;
498 this->smColDel();
499 act--;
500 this->smZeroElim();
501 if (sign == 0)
502 {
503 this->smActDel();
504 return NULL;
505 }
506 if (act < 2)
507 {
508 this->smFinalMult();
509 this->smPivDel();
510 if (act != 0) res = m_act[1]->m;
511 omFreeBin((void *)m_act[1], smprec_bin);
512 return res;
513 }
514 loop
515 {
516 this->smNewPivot();
517 this->smSign();
518 this->smSelectPR();
519 this->smMultCol();
520 this->smHElim();
521 crd++;
522 m_res[crd] = piv;
523 this->smColDel();
524 act--;
525 this->smZeroElim();
526 if (sign == 0)
527 {
528 this->smPivDel();
529 this->smActDel();
530 return NULL;
531 }
532 if (act < 2)
533 {
534 if (TEST_OPT_PROT) PrintS(".\n");
535 this->smFinalMult();
536 this->smPivDel();
537 if (act != 0) res = m_act[1]->m;
538 omFreeBin((void *)m_act[1], smprec_bin);
539 return res;
540 }
541 }
542}
543
544/*
545* the new Bareiss elimination:
546* - with x unreduced last rows, pivots from here are not allowed
547* - the method will finish for number of unreduced columns < y
548*/
550{
551 if ((x > 0) && (x < nrows))
552 {
553 tored -= x;
554 this->smToredElim();
555 }
556 if (y < 1) y = 1;
557 if (act <= y)
558 {
559 this->smCopToRes();
560 return;
561 }
562 normalize = this->smCheckNormalize();
563 if (normalize) this->smNormalize();
564 this->smPivot();
565 this->smSelectPR();
566 this->sm1Elim();
567 crd++;
568 this->smColToRow();
569 act--;
570 this->smRowToCol();
571 this->smZeroElim();
572 if (tored != nrows)
573 this->smToredElim();
574 if (act <= y)
575 {
576 this->smFinalMult();
577 this->smCopToRes();
578 return;
579 }
580 loop
581 {
582 if (normalize) this->smNormalize();
583 this->smNewPivot();
584 this->smSelectPR();
585 this->smMultCol();
586 this->smHElim();
587 crd++;
588 this->smColToRow();
589 act--;
590 this->smRowToCol();
591 this->smZeroElim();
592 if (tored != nrows)
593 this->smToredElim();
594 if (act <= y)
595 {
596 if (TEST_OPT_PROT) PrintS(".\n");
597 this->smFinalMult();
598 this->smCopToRes();
599 return;
600 }
601 }
602}
603
604/* ----------------- pivot method ------------------ */
605
606/*
607* prepare smPivot, compute weights for rows and columns
608* and the weight for all points
609*/
611{
612 float wc, wp, w;
613 smpoly a;
614 int i;
615
616 wp = 0.0;
617 for (i=tored; i; i--) wrw[i] = 0.0; // ???
618 for (i=act; i; i--)
619 {
620 wc = 0.0;
621 a = m_act[i];
622 loop
623 {
624 if (a->pos > tored)
625 break;
626 w = a->f = sm_PolyWeight(a,_R);
627 wc += w;
628 wrw[a->pos] += w;
629 a = a->n;
630 if (a == NULL)
631 break;
632 }
633 wp += wc;
634 wcl[i] = wc;
635 }
636 wpoints = wp;
637}
638
639/*
640* compute pivot
641*/
643{
644 float wopt = 1.0e30;
645 float wc, wr, wp, w;
646 smpoly a;
647 int i, copt, ropt;
648
649 this->smWeights();
650 for (i=act; i; i--)
651 {
652 a = m_act[i];
653 loop
654 {
655 if (a->pos > tored)
656 break;
657 w = a->f;
658 wc = wcl[i]-w;
659 wr = wrw[a->pos]-w;
660 if ((wr<0.25) || (wc<0.25)) // row or column with only one point
661 {
662 if (w<wopt)
663 {
664 wopt = w;
665 copt = i;
666 ropt = a->pos;
667 }
668 }
669 else // elimination
670 {
671 wp = w*(wpoints-wcl[i]-wr);
672 wp += wr*wc;
673 if (wp < wopt)
674 {
675 wopt = wp;
676 copt = i;
677 ropt = a->pos;
678 }
679 }
680 a = a->n;
681 if (a == NULL)
682 break;
683 }
684 }
685 rpiv = ropt;
686 cpiv = copt;
687 if (cpiv != act)
688 {
689 a = m_act[act];
690 m_act[act] = m_act[cpiv];
691 m_act[cpiv] = a;
692 }
693}
694
695/*
696* prepare smPivot, compute weights for rows and columns
697* and the weight for all points
698*/
700{
701 float wc, wp, w, hp = piv->f;
702 smpoly a;
703 int i, f, e = crd;
704
705 wp = 0.0;
706 for (i=tored; i; i--) wrw[i] = 0.0; // ???
707 for (i=act; i; i--)
708 {
709 wc = 0.0;
710 a = m_act[i];
711 loop
712 {
713 if (a->pos > tored)
714 break;
715 w = a->f;
716 f = a->e;
717 if (f < e)
718 {
719 w *= hp;
720 if (f) w /= m_res[f]->f;
721 }
722 wc += w;
723 wrw[a->pos] += w;
724 a = a->n;
725 if (a == NULL)
726 break;
727 }
728 wp += wc;
729 wcl[i] = wc;
730 }
731 wpoints = wp;
732}
733
734/*
735* compute pivot
736*/
738{
739 float wopt = 1.0e30, hp = piv->f;
740 float wc, wr, wp, w;
741 smpoly a;
742 int i, copt, ropt, f, e = crd;
743
744 this->smNewWeights();
745 for (i=act; i; i--)
746 {
747 a = m_act[i];
748 loop
749 {
750 if (a->pos > tored)
751 break;
752 w = a->f;
753 f = a->e;
754 if (f < e)
755 {
756 w *= hp;
757 if (f) w /= m_res[f]->f;
758 }
759 wc = wcl[i]-w;
760 wr = wrw[a->pos]-w;
761 if ((wr<0.25) || (wc<0.25)) // row or column with only one point
762 {
763 if (w<wopt)
764 {
765 wopt = w;
766 copt = i;
767 ropt = a->pos;
768 }
769 }
770 else // elimination
771 {
772 wp = w*(wpoints-wcl[i]-wr);
773 wp += wr*wc;
774 if (wp < wopt)
775 {
776 wopt = wp;
777 copt = i;
778 ropt = a->pos;
779 }
780 }
781 a = a->n;
782 if (a == NULL)
783 break;
784 }
785 }
786 rpiv = ropt;
787 cpiv = copt;
788 if (cpiv != act)
789 {
790 a = m_act[act];
791 m_act[act] = m_act[cpiv];
792 m_act[cpiv] = a;
793 }
794}
795
796/* ----------------- elimination ------------------ */
797
798/* first step of elimination */
800{
801 poly p = piv->m; // pivotelement
802 smpoly c = m_act[act]; // pivotcolumn
803 smpoly r = red; // row to reduce
804 smpoly res, a, b;
805 poly w, ha, hb;
806
807 if ((c == NULL) || (r == NULL))
808 {
809 while (r!=NULL) sm_ElemDelete(&r,_R);
810 return;
811 }
812 do
813 {
814 a = m_act[r->pos];
815 res = dumm;
816 res->n = NULL;
817 b = c;
818 w = r->m;
819 loop // combine the chains a and b: p*a + w*b
820 {
821 if (a == NULL)
822 {
823 do
824 {
825 res = res->n = smElemCopy(b);
826 res->m = pp_Mult_qq(b->m, w,_R);
827 res->e = 1;
828 res->f = sm_PolyWeight(res,_R);
829 b = b->n;
830 } while (b != NULL);
831 break;
832 }
833 if (a->pos < b->pos)
834 {
835 res = res->n = a;
836 a = a->n;
837 }
838 else if (a->pos > b->pos)
839 {
840 res = res->n = smElemCopy(b);
841 res->m = pp_Mult_qq(b->m, w,_R);
842 res->e = 1;
843 res->f = sm_PolyWeight(res,_R);
844 b = b->n;
845 }
846 else
847 {
848 ha = pp_Mult_qq(a->m, p,_R);
849 p_Delete(&a->m,_R);
850 hb = pp_Mult_qq(b->m, w,_R);
851 ha = p_Add_q(ha, hb,_R);
852 if (ha != NULL)
853 {
854 a->m = ha;
855 a->e = 1;
856 a->f = sm_PolyWeight(a,_R);
857 res = res->n = a;
858 a = a->n;
859 }
860 else
861 {
862 sm_ElemDelete(&a,_R);
863 }
864 b = b->n;
865 }
866 if (b == NULL)
867 {
868 res->n = a;
869 break;
870 }
871 }
872 m_act[r->pos] = dumm->n;
873 sm_ElemDelete(&r,_R);
874 } while (r != NULL);
875}
876
877/* higher steps of elimination */
879{
880 poly hp = this->smMultPoly(piv);
881 poly gp = piv->m; // pivotelement
882 smpoly c = m_act[act]; // pivotcolumn
883 smpoly r = red; // row to reduce
884 smpoly res, a, b;
885 poly ha, hr, x, y;
886 int e, ip, ir, ia;
887
888 if ((c == NULL) || (r == NULL))
889 {
890 while(r!=NULL) sm_ElemDelete(&r,_R);
891 p_Delete(&hp,_R);
892 return;
893 }
894 e = crd+1;
895 ip = piv->e;
896 do
897 {
898 a = m_act[r->pos];
899 res = dumm;
900 res->n = NULL;
901 b = c;
902 hr = r->m;
903 ir = r->e;
904 loop // combine the chains a and b: (hp,gp)*a(l) + hr*b(h)
905 {
906 if (a == NULL)
907 {
908 do
909 {
910 res = res->n = smElemCopy(b);
911 x = SM_MULT(b->m, hr, m_res[ir]->m,_R);
912 b = b->n;
913 if(ir) SM_DIV(x, m_res[ir]->m,_R);
914 res->m = x;
915 res->e = e;
916 res->f = sm_PolyWeight(res,_R);
917 } while (b != NULL);
918 break;
919 }
920 if (a->pos < b->pos)
921 {
922 res = res->n = a;
923 a = a->n;
924 }
925 else if (a->pos > b->pos)
926 {
927 res = res->n = smElemCopy(b);
928 x = SM_MULT(b->m, hr, m_res[ir]->m,_R);
929 b = b->n;
930 if(ir) SM_DIV(x, m_res[ir]->m,_R);
931 res->m = x;
932 res->e = e;
933 res->f = sm_PolyWeight(res,_R);
934 }
935 else
936 {
937 ha = a->m;
938 ia = a->e;
939 if (ir >= ia)
940 {
941 if (ir > ia)
942 {
943 x = SM_MULT(ha, m_res[ir]->m, m_res[ia]->m,_R);
944 p_Delete(&ha,_R);
945 ha = x;
946 if (ia) SM_DIV(ha, m_res[ia]->m,_R);
947 ia = ir;
948 }
949 x = SM_MULT(ha, gp, m_res[ia]->m,_R);
950 p_Delete(&ha,_R);
951 y = SM_MULT(b->m, hr, m_res[ia]->m,_R);
952 }
953 else if (ir >= ip)
954 {
955 if (ia < crd)
956 {
957 x = SM_MULT(ha, m_res[crd]->m, m_res[ia]->m,_R);
958 p_Delete(&ha,_R);
959 ha = x;
960 SM_DIV(ha, m_res[ia]->m,_R);
961 }
962 y = hp;
963 if(ir > ip)
964 {
965 y = SM_MULT(y, m_res[ir]->m, m_res[ip]->m,_R);
966 if (ip) SM_DIV(y, m_res[ip]->m,_R);
967 }
968 ia = ir;
969 x = SM_MULT(ha, y, m_res[ia]->m,_R);
970 if (y != hp) p_Delete(&y,_R);
971 p_Delete(&ha,_R);
972 y = SM_MULT(b->m, hr, m_res[ia]->m,_R);
973 }
974 else
975 {
976 x = SM_MULT(hr, m_res[ia]->m, m_res[ir]->m,_R);
977 if (ir) SM_DIV(x, m_res[ir]->m,_R);
978 y = SM_MULT(b->m, x, m_res[ia]->m,_R);
979 p_Delete(&x,_R);
980 x = SM_MULT(ha, gp, m_res[ia]->m,_R);
981 p_Delete(&ha,_R);
982 }
983 ha = p_Add_q(x, y,_R);
984 if (ha != NULL)
985 {
986 if (ia) SM_DIV(ha, m_res[ia]->m,_R);
987 a->m = ha;
988 a->e = e;
989 a->f = sm_PolyWeight(a,_R);
990 res = res->n = a;
991 a = a->n;
992 }
993 else
994 {
995 a->m = NULL;
996 sm_ElemDelete(&a,_R);
997 }
998 b = b->n;
999 }
1000 if (b == NULL)
1001 {
1002 res->n = a;
1003 break;
1004 }
1005 }
1006 m_act[r->pos] = dumm->n;
1007 sm_ElemDelete(&r,_R);
1008 } while (r != NULL);
1009 p_Delete(&hp,_R);
1010}
1011
1012/* ----------------- transfer ------------------ */
1013
1014/*
1015* select the pivotrow and store it to red and piv
1016*/
1018{
1019 smpoly b = dumm;
1020 smpoly a, ap;
1021 int i;
1022
1023 if (TEST_OPT_PROT)
1024 {
1025 if ((crd+1)%10)
1026 PrintS(".");
1027 else
1028 PrintS(".\n");
1029 }
1030 a = m_act[act];
1031 if (a->pos < rpiv)
1032 {
1033 do
1034 {
1035 ap = a;
1036 a = a->n;
1037 } while (a->pos < rpiv);
1038 ap->n = a->n;
1039 }
1040 else
1041 m_act[act] = a->n;
1042 piv = a;
1043 a->n = NULL;
1044 for (i=1; i<act; i++)
1045 {
1046 a = m_act[i];
1047 if (a->pos < rpiv)
1048 {
1049 loop
1050 {
1051 ap = a;
1052 a = a->n;
1053 if ((a == NULL) || (a->pos > rpiv))
1054 break;
1055 if (a->pos == rpiv)
1056 {
1057 ap->n = a->n;
1058 a->m = p_Neg(a->m,_R);
1059 b = b->n = a;
1060 b->pos = i;
1061 break;
1062 }
1063 }
1064 }
1065 else if (a->pos == rpiv)
1066 {
1067 m_act[i] = a->n;
1068 a->m = p_Neg(a->m,_R);
1069 b = b->n = a;
1070 b->pos = i;
1071 }
1072 }
1073 b->n = NULL;
1074 red = dumm->n;
1075}
1076
1077/*
1078* store the pivotcol in m_row
1079* m_act[cols][pos(rows)] => m_row[rows][pos(cols)]
1080*/
1082{
1083 smpoly c = m_act[act];
1084 smpoly h;
1085
1086 while (c != NULL)
1087 {
1088 h = c;
1089 c = c->n;
1090 h->n = m_row[h->pos];
1091 m_row[h->pos] = h;
1092 h->pos = crd;
1093 }
1094}
1095
1096/*
1097* store the pivot and the assosiated row in m_row
1098* to m_res (result):
1099* piv + m_row[rows][pos(cols)] => m_res[cols][pos(rows)]
1100*/
1102{
1103 smpoly r = m_row[rpiv];
1104 smpoly a, ap, h;
1105
1106 m_row[rpiv] = NULL;
1107 perm[crd] = rpiv;
1108 piv->pos = crd;
1109 m_res[crd] = piv;
1110 while (r != NULL)
1111 {
1112 ap = m_res[r->pos];
1113 loop
1114 {
1115 a = ap->n;
1116 if (a == NULL)
1117 {
1118 ap->n = h = r;
1119 r = r->n;
1120 h->n = a;
1121 h->pos = crd;
1122 break;
1123 }
1124 ap = a;
1125 }
1126 }
1127}
1128
1129/* ----------------- C++ stuff ------------------ */
1130
1131/*
1132* clean m_act from zeros (after elim)
1133*/
1135{
1136 int i = 0;
1137 int j;
1138
1139 loop
1140 {
1141 i++;
1142 if (i > act) return;
1143 if (m_act[i] == NULL) break;
1144 }
1145 j = i;
1146 loop
1147 {
1148 j++;
1149 if (j > act) break;
1150 if (m_act[j] != NULL)
1151 {
1152 m_act[i] = m_act[j];
1153 i++;
1154 }
1155 }
1156 act -= (j-i);
1157 sign = 0;
1158}
1159
1160/*
1161* clean m_act from cols not to reduced (after elim)
1162* put them to m_res
1163*/
1165{
1166 int i = 0;
1167 int j;
1168
1169 loop
1170 {
1171 i++;
1172 if (i > act) return;
1173 if (m_act[i]->pos > tored)
1174 {
1175 m_res[inred] = m_act[i];
1176 inred--;
1177 break;
1178 }
1179 }
1180 j = i;
1181 loop
1182 {
1183 j++;
1184 if (j > act) break;
1185 if (m_act[j]->pos > tored)
1186 {
1187 m_res[inred] = m_act[j];
1188 inred--;
1189 }
1190 else
1191 {
1192 m_act[i] = m_act[j];
1193 i++;
1194 }
1195 }
1196 act -= (j-i);
1197 sign = 0;
1198}
1199
1200/*
1201* copy m_act to m_res
1202*/
1204{
1205 smpoly a,ap,r,h;
1206 int i,j,k,l;
1207
1208 i = 0;
1209 if (act)
1210 {
1211 a = m_act[act]; // init perm
1212 do
1213 {
1214 i++;
1215 perm[crd+i] = a->pos;
1216 a = a->n;
1217 } while ((a != NULL) && (a->pos <= tored));
1218 for (j=act-1;j;j--) // load all positions of perm
1219 {
1220 a = m_act[j];
1221 k = 1;
1222 loop
1223 {
1224 if (perm[crd+k] >= a->pos)
1225 {
1226 if (perm[crd+k] > a->pos)
1227 {
1228 for (l=i;l>=k;l--) perm[crd+l+1] = perm[crd+l];
1229 perm[crd+k] = a->pos;
1230 i++;
1231 }
1232 a = a->n;
1233 if ((a == NULL) || (a->pos > tored)) break;
1234 }
1235 k++;
1236 if ((k > i) && (a->pos <= tored))
1237 {
1238 do
1239 {
1240 i++;
1241 perm[crd+i] = a->pos;
1242 a = a->n;
1243 } while ((a != NULL) && (a->pos <= tored));
1244 break;
1245 }
1246 }
1247 }
1248 }
1249 for (j=act;j;j--) // renumber m_act
1250 {
1251 k = 1;
1252 a = m_act[j];
1253 while ((a != NULL) && (a->pos <= tored))
1254 {
1255 if (perm[crd+k] == a->pos)
1256 {
1257 a->pos = crd+k;
1258 a = a->n;
1259 }
1260 k++;
1261 }
1262 }
1263 tored = crd+i;
1264 for(k=1;k<=i;k++) // clean this from m_row
1265 {
1266 j = perm[crd+k];
1267 if (m_row[j] != NULL)
1268 {
1269 r = m_row[j];
1270 m_row[j] = NULL;
1271 do
1272 {
1273 ap = m_res[r->pos];
1274 loop
1275 {
1276 a = ap->n;
1277 if (a == NULL)
1278 {
1279 h = ap->n = r;
1280 r = r->n;
1281 h->n = NULL;
1282 h->pos = crd+k;
1283 break;
1284 }
1285 ap = a;
1286 }
1287 } while (r!=NULL);
1288 }
1289 }
1290 while(act) // clean m_act
1291 {
1292 crd++;
1293 m_res[crd] = m_act[act];
1294 act--;
1295 }
1296 for (i=1;i<=tored;i++) // take the rest of m_row
1297 {
1298 if(m_row[i] != NULL)
1299 {
1300 tored++;
1301 r = m_row[i];
1302 m_row[i] = NULL;
1303 perm[tored] = i;
1304 do
1305 {
1306 ap = m_res[r->pos];
1307 loop
1308 {
1309 a = ap->n;
1310 if (a == NULL)
1311 {
1312 h = ap->n = r;
1313 r = r->n;
1314 h->n = NULL;
1315 h->pos = tored;
1316 break;
1317 }
1318 ap = a;
1319 }
1320 } while (r!=NULL);
1321 }
1322 }
1323 for (i=tored+1;i<=nrows;i++) // take the rest of m_row
1324 {
1325 if(m_row[i] != NULL)
1326 {
1327 r = m_row[i];
1328 m_row[i] = NULL;
1329 do
1330 {
1331 ap = m_res[r->pos];
1332 loop
1333 {
1334 a = ap->n;
1335 if (a == NULL)
1336 {
1337 h = ap->n = r;
1338 r = r->n;
1339 h->n = NULL;
1340 h->pos = i;
1341 break;
1342 }
1343 ap = a;
1344 }
1345 } while (r!=NULL);
1346 }
1347 }
1348 while (inred < ncols) // take unreducable
1349 {
1350 crd++;
1351 inred++;
1352 m_res[crd] = m_res[inred];
1353 }
1354}
1355
1356/*
1357* multiply and divide the column, that goes to result
1358*/
1360{
1361 smpoly a = m_act[act];
1362 int e = crd;
1363 poly ha;
1364 int f;
1365
1366 while (a != NULL)
1367 {
1368 f = a->e;
1369 if (f < e)
1370 {
1371 ha = SM_MULT(a->m, m_res[e]->m, m_res[f]->m,_R);
1372 p_Delete(&a->m,_R);
1373 if (f) SM_DIV(ha, m_res[f]->m,_R);
1374 a->m = ha;
1375 if (normalize) p_Normalize(a->m,_R);
1376 }
1377 a = a->n;
1378 }
1379}
1380
1381/*
1382* multiply and divide the m_act finaly
1383*/
1385{
1386 smpoly a;
1387 poly ha;
1388 int i, f;
1389 int e = crd;
1390
1391 for (i=act; i; i--)
1392 {
1393 a = m_act[i];
1394 do
1395 {
1396 f = a->e;
1397 if (f < e)
1398 {
1399 ha = SM_MULT(a->m, m_res[e]->m, m_res[f]->m, _R);
1400 p_Delete(&a->m,_R);
1401 if (f) SM_DIV(ha, m_res[f]->m, _R);
1402 a->m = ha;
1403 }
1404 if (normalize) p_Normalize(a->m, _R);
1405 a = a->n;
1406 } while (a != NULL);
1407 }
1408}
1409
1410/*
1411* check for denominators
1412*/
1414{
1415 int i;
1416 smpoly a;
1417
1418 for (i=act; i; i--)
1419 {
1420 a = m_act[i];
1421 do
1422 {
1423 if(sm_HaveDenom(a->m,_R)) return 1;
1424 a = a->n;
1425 } while (a != NULL);
1426 }
1427 return 0;
1428}
1429
1430/*
1431* normalize
1432*/
1434{
1435 smpoly a;
1436 int i;
1437 int e = crd;
1438
1439 for (i=act; i; i--)
1440 {
1441 a = m_act[i];
1442 do
1443 {
1444 if (e == a->e) p_Normalize(a->m,_R);
1445 a = a->n;
1446 } while (a != NULL);
1447 }
1448}
1449
1450/*
1451* multiply and divide the element, save poly
1452*/
1454{
1455 int f = a->e;
1456 poly r, h;
1457
1458 if (f < crd)
1459 {
1460 h = r = a->m;
1461 h = SM_MULT(h, m_res[crd]->m, m_res[f]->m, _R);
1462 if (f) SM_DIV(h, m_res[f]->m, _R);
1463 a->m = h;
1464 if (normalize) p_Normalize(a->m,_R);
1465 a->f = sm_PolyWeight(a,_R);
1466 return r;
1467 }
1468 else
1469 return NULL;
1470}
1471
1472/*
1473* delete the m_act finaly
1474*/
1476{
1477 smpoly a;
1478 int i;
1479
1480 for (i=act; i; i--)
1481 {
1482 a = m_act[i];
1483 do
1484 {
1485 sm_ElemDelete(&a,_R);
1486 } while (a != NULL);
1487 }
1488}
1489
1490/*
1491* delete the pivotcol
1492*/
1494{
1495 smpoly a = m_act[act];
1496
1497 while (a != NULL)
1498 {
1499 sm_ElemDelete(&a,_R);
1500 }
1501}
1502
1503/*
1504* delete pivot elements
1505*/
1507{
1508 int i=crd;
1509
1510 while (i != 0)
1511 {
1513 i--;
1514 }
1515}
1516
1517/*
1518* the sign of the determinant
1519*/
1521{
1522 int j,i;
1523 if (act > 2)
1524 {
1525 if (cpiv!=act) sign=-sign;
1526 if ((act%2)==0) sign=-sign;
1527 i=1;
1528 j=perm[1];
1529 while(j<rpiv)
1530 {
1531 sign=-sign;
1532 i++;
1533 j=perm[i];
1534 }
1535 while(perm[i]!=0)
1536 {
1537 perm[i]=perm[i+1];
1538 i++;
1539 }
1540 }
1541 else
1542 {
1543 if (cpiv!=1) sign=-sign;
1544 if (rpiv!=perm[1]) sign=-sign;
1545 }
1546}
1547
1549{
1550 int i;
1551 for (i=act;i;i--) perm[i]=i;
1552}
1553
1554/* ----------------- arithmetic ------------------ */
1555#ifdef OLD_DIV
1556/*2
1557* exact division a/b
1558* a destroyed, b NOT destroyed
1559*/
1560void sm_PolyDiv(poly a, poly b, const ring R)
1561{
1562 const number x = pGetCoeff(b);
1563 number y, yn;
1564 poly t, h, dummy;
1565 int i;
1566
1567 if (pNext(b) == NULL)
1568 {
1569 do
1570 {
1571 if (!p_LmIsConstantComp(b,R))
1572 {
1573 for (i=rVar(R); i; i--)
1574 p_SubExp(a,i,p_GetExp(b,i,R),R);
1575 p_Setm(a,R);
1576 }
1577 y = n_Div(pGetCoeff(a),x,R->cf);
1578 n_Normalize(y,R->cf);
1579 p_SetCoeff(a,y,R);
1580 pIter(a);
1581 } while (a != NULL);
1582 return;
1583 }
1584 dummy = p_Init(R);
1585 do
1586 {
1587 for (i=rVar(R); i; i--)
1588 p_SubExp(a,i,p_GetExp(b,i,R),R);
1589 p_Setm(a,R);
1590 y = n_Div(pGetCoeff(a),x,R->cf);
1591 n_Normalize(y,R->cf);
1592 p_SetCoeff(a,y,R);
1593 yn = n_InpNeg(n_Copy(y,R->cf),R->cf);
1594 t = pNext(b);
1595 h = dummy;
1596 do
1597 {
1598 h = pNext(h) = p_Init(R);
1599 //pSetComp(h,0);
1600 for (i=rVar(R); i; i--)
1601 p_SetExp(h,i,p_GetExp(a,i,R)+p_GetExp(t,i,R),R);
1602 p_Setm(h,R);
1603 pSetCoeff0(h,n_Mult(yn, pGetCoeff(t),R->cf));
1604 pIter(t);
1605 } while (t != NULL);
1606 n_Delete(&yn,R->cf);
1607 pNext(h) = NULL;
1608 a = pNext(a) = p_Add_q(pNext(a), pNext(dummy),R);
1609 } while (a!=NULL);
1610 p_LmFree(dummy, R);
1611}
1612#endif
1613
1614//disable that, as it fails with coef buckets
1615//#define X_MAS
1616#ifdef X_MAS
1617// Note: the following in not addapted to SW :(
1618/*
1619/// returns the part of (a*b)/exp(lead(c)) with nonegative exponents
1620poly smMultDiv(poly a, poly b, const poly c)
1621{
1622 poly pa, e, res, r;
1623 BOOLEAN lead;
1624 int la, lb, lr;
1625
1626 if ((c == NULL) || pLmIsConstantComp(c))
1627 {
1628 return pp_Mult_qq(a, b);
1629 }
1630
1631 pqLength(a, b, la, lb, SM_MIN_LENGTH_BUCKET);
1632
1633 // we iter over b, so, make sure b is the shortest
1634 // such that we minimize our iterations
1635 if (lb > la)
1636 {
1637 r = a;
1638 a = b;
1639 b = r;
1640 lr = la;
1641 la = lb;
1642 lb = lr;
1643 }
1644 res = NULL;
1645 e = p_Init();
1646 lead = FALSE;
1647 while (!lead)
1648 {
1649 pSetCoeff0(e,pGetCoeff(b));
1650 if (smIsNegQuot(e, b, c))
1651 {
1652 lead = pLmDivisibleByNoComp(e, a);
1653 r = smSelectCopy_ExpMultDiv(a, e, b, c);
1654 }
1655 else
1656 {
1657 lead = TRUE;
1658 r = pp_Mult__mm(a, e);
1659 }
1660 if (lead)
1661 {
1662 if (res != NULL)
1663 {
1664 smFindRef(&pa, &res, r);
1665 if (pa == NULL)
1666 lead = FALSE;
1667 }
1668 else
1669 {
1670 pa = res = r;
1671 }
1672 }
1673 else
1674 res = p_Add_q(res, r);
1675 pIter(b);
1676 if (b == NULL)
1677 {
1678 pLmFree(e);
1679 return res;
1680 }
1681 }
1682
1683 if (!TEST_OPT_NOT_BUCKETS && lb >= SM_MIN_LENGTH_BUCKET)
1684 {
1685 // use buckets if minimum length is smaller than threshold
1686 spolyrec rp;
1687 poly append;
1688 // find the last monomial before pa
1689 if (res == pa)
1690 {
1691 append = &rp;
1692 pNext(append) = res;
1693 }
1694 else
1695 {
1696 append = res;
1697 while (pNext(append) != pa)
1698 {
1699 assume(pNext(append) != NULL);
1700 pIter(append);
1701 }
1702 }
1703 kBucket_pt bucket = kBucketCreate(currRing);
1704 kBucketInit(bucket, pNext(append), 0);
1705 do
1706 {
1707 pSetCoeff0(e,pGetCoeff(b));
1708 if (smIsNegQuot(e, b, c))
1709 {
1710 lr = la;
1711 r = pp_Mult_Coeff_mm_DivSelect_MultDiv(a, lr, e, b, c);
1712 if (pLmDivisibleByNoComp(e, a))
1713 append = kBucket_ExtractLarger_Add_q(bucket, append, r, &lr);
1714 else
1715 kBucket_Add_q(bucket, r, &lr);
1716 }
1717 else
1718 {
1719 r = pp_Mult__mm(a, e);
1720 append = kBucket_ExtractLarger_Add_q(bucket, append, r, &la);
1721 }
1722 pIter(b);
1723 } while (b != NULL);
1724 pNext(append) = kBucketClear(bucket);
1725 kBucketDestroy(&bucket);
1726 pTest(append);
1727 }
1728 else
1729 {
1730 // use old sm stuff
1731 do
1732 {
1733 pSetCoeff0(e,pGetCoeff(b));
1734 if (smIsNegQuot(e, b, c))
1735 {
1736 r = smSelectCopy_ExpMultDiv(a, e, b, c);
1737 if (pLmDivisibleByNoComp(e, a))
1738 smCombineChain(&pa, r);
1739 else
1740 pa = p_Add_q(pa,r);
1741 }
1742 else
1743 {
1744 r = pp_Mult__mm(a, e);
1745 smCombineChain(&pa, r);
1746 }
1747 pIter(b);
1748 } while (b != NULL);
1749 }
1750 pLmFree(e);
1751 return res;
1752}
1753*/
1754#else
1755
1756/*
1757* returns the part of (a*b)/exp(lead(c)) with nonegative exponents
1758*/
1759poly sm_MultDiv(poly a, poly b, const poly c, const ring R)
1760{
1761 poly pa, e, res, r;
1762 BOOLEAN lead;
1763
1764 if ((c == NULL) || p_LmIsConstantComp(c,R))
1765 {
1766 return pp_Mult_qq(a, b, R);
1767 }
1768 if (smSmaller(a, b))
1769 {
1770 r = a;
1771 a = b;
1772 b = r;
1773 }
1774 res = NULL;
1775 e = p_Init(R);
1776 lead = FALSE;
1777 while (!lead)
1778 {
1780 if (sm_IsNegQuot(e, b, c, R))
1781 {
1782 lead = p_LmDivisibleByNoComp(e, a, R);
1783 r = sm_SelectCopy_ExpMultDiv(a, e, b, c, R);
1784 }
1785 else
1786 {
1787 lead = TRUE;
1788 r = pp_Mult_mm(a, e,R);
1789 }
1790 if (lead)
1791 {
1792 if (res != NULL)
1793 {
1794 sm_FindRef(&pa, &res, r, R);
1795 if (pa == NULL)
1796 lead = FALSE;
1797 }
1798 else
1799 {
1800 pa = res = r;
1801 }
1802 }
1803 else
1804 res = p_Add_q(res, r, R);
1805 pIter(b);
1806 if (b == NULL)
1807 {
1808 p_LmFree(e, R);
1809 return res;
1810 }
1811 }
1812 do
1813 {
1815 if (sm_IsNegQuot(e, b, c, R))
1816 {
1817 r = sm_SelectCopy_ExpMultDiv(a, e, b, c, R);
1818 if (p_LmDivisibleByNoComp(e, a,R))
1819 sm_CombineChain(&pa, r, R);
1820 else
1821 pa = p_Add_q(pa,r,R);
1822 }
1823 else
1824 {
1825 r = pp_Mult_mm(a, e, R);
1826 sm_CombineChain(&pa, r, R);
1827 }
1828 pIter(b);
1829 } while (b != NULL);
1830 p_LmFree(e, R);
1831 return res;
1832}
1833#endif /*else X_MAS*/
1834
1835/*n
1836* exact division a/b
1837* a is a result of smMultDiv
1838* a destroyed, b NOT destroyed
1839*/
1840void sm_SpecialPolyDiv(poly a, poly b, const ring R)
1841{
1842 if (pNext(b) == NULL)
1843 {
1844 sm_PolyDivN(a, pGetCoeff(b),R);
1845 return;
1846 }
1847 sm_ExactPolyDiv(a, b, R);
1848}
1849
1850/* ------------ internals arithmetic ------------- */
1851static void sm_ExactPolyDiv(poly a, poly b, const ring R)
1852{
1853 const number x = pGetCoeff(b);
1854 poly tail = pNext(b), e = p_Init(R);
1855 poly h;
1856 number y, yn;
1857 int lt = pLength(tail);
1858
1860 {
1861 kBucket_pt bucket = kBucketCreate(R);
1862 kBucketInit(bucket, pNext(a), 0);
1863 int lh = 0;
1864 do
1865 {
1866 y = n_Div(pGetCoeff(a), x, R->cf);
1867 n_Normalize(y, R->cf);
1868 p_SetCoeff(a,y, R);
1869 yn = n_InpNeg(n_Copy(y, R->cf), R->cf);
1870 pSetCoeff0(e,yn);
1871 lh = lt;
1872 if (sm_IsNegQuot(e, a, b, R))
1873 {
1874 h = pp_Mult_Coeff_mm_DivSelect_MultDiv(tail, lh, e, a, b, R);
1875 }
1876 else
1877 h = pp_Mult_mm(tail, e, R);
1878 n_Delete(&yn, R->cf);
1879 kBucket_Add_q(bucket, h, &lh);
1880
1881 a = pNext(a) = kBucketExtractLm(bucket);
1882 } while (a!=NULL);
1883 kBucketDestroy(&bucket);
1884 }
1885 else
1886 {
1887 do
1888 {
1889 y = n_Div(pGetCoeff(a), x, R->cf);
1890 n_Normalize(y, R->cf);
1891 p_SetCoeff(a,y, R);
1892 yn = n_InpNeg(n_Copy(y, R->cf), R->cf);
1893 pSetCoeff0(e,yn);
1894 if (sm_IsNegQuot(e, a, b, R))
1895 h = sm_SelectCopy_ExpMultDiv(tail, e, a, b, R);
1896 else
1897 h = pp_Mult_mm(tail, e, R);
1898 n_Delete(&yn, R->cf);
1899 a = pNext(a) = p_Add_q(pNext(a), h, R);
1900 } while (a!=NULL);
1901 }
1902 p_LmFree(e, R);
1903}
1904
1905// obachman --> Wilfried: check the following
1906static BOOLEAN sm_IsNegQuot(poly a, const poly b, const poly c, const ring R)
1907{
1908 if (p_LmDivisibleByNoComp(c, b,R))
1909 {
1910 p_ExpVectorDiff(a, b, c,R);
1911 // Hmm: here used to be a p_Setm(a): but it is unnecessary,
1912 // if b and c are correct
1913 return FALSE;
1914 }
1915 else
1916 {
1917 int i;
1918 for (i=rVar(R); i>0; i--)
1919 {
1920 if(p_GetExp(c,i,R) > p_GetExp(b,i,R))
1921 p_SetExp(a,i,p_GetExp(c,i,R)-p_GetExp(b,i,R),R);
1922 else
1923 p_SetExp(a,i,0,R);
1924 }
1925 // here we actually might need a p_Setm, if a is to be used in
1926 // comparisons
1927 return TRUE;
1928 }
1929}
1930
1931static void sm_ExpMultDiv(poly t, const poly b, const poly c, const ring R)
1932{
1933 p_Test(t,R);
1934 p_LmTest(b,R);
1935 p_LmTest(c,R);
1936 poly bc = p_New(R);
1937
1938 p_ExpVectorDiff(bc, b, c, R);
1939
1940 while(t!=NULL)
1941 {
1942 p_ExpVectorAdd(t, bc, R);
1943 pIter(t);
1944 }
1945 p_LmFree(bc, R);
1946}
1947
1948
1949static void sm_PolyDivN(poly a, const number x, const ring R)
1950{
1951 number y;
1952
1953 do
1954 {
1955 y = n_Div(pGetCoeff(a),x, R->cf);
1956 n_Normalize(y, R->cf);
1957 p_SetCoeff(a,y, R);
1958 pIter(a);
1959 } while (a != NULL);
1960}
1961
1962static BOOLEAN smSmaller(poly a, poly b)
1963{
1964 loop
1965 {
1966 pIter(b);
1967 if (b == NULL) return TRUE;
1968 pIter(a);
1969 if (a == NULL) return FALSE;
1970 }
1971}
1972
1973static void sm_CombineChain(poly *px, poly r, const ring R)
1974{
1975 poly pa = *px, pb;
1976 number x;
1977 int i;
1978
1979 loop
1980 {
1981 pb = pNext(pa);
1982 if (pb == NULL)
1983 {
1984 pa = pNext(pa) = r;
1985 break;
1986 }
1987 i = p_LmCmp(pb, r,R);
1988 if (i > 0)
1989 pa = pb;
1990 else
1991 {
1992 if (i == 0)
1993 {
1994 x = n_Add(pGetCoeff(pb), pGetCoeff(r),R->cf);
1995 p_LmDelete(&r,R);
1996 if (n_IsZero(x,R->cf))
1997 {
1998 p_LmDelete(&pb,R);
1999 pNext(pa) = p_Add_q(pb,r,R);
2000 }
2001 else
2002 {
2003 pa = pb;
2004 p_SetCoeff(pa,x,R);
2005 pNext(pa) = p_Add_q(pNext(pa), r, R);
2006 }
2007 }
2008 else
2009 {
2010 pa = pNext(pa) = r;
2011 pNext(pa) = p_Add_q(pb, pNext(pa),R);
2012 }
2013 break;
2014 }
2015 }
2016 *px = pa;
2017}
2018
2019
2020static void sm_FindRef(poly *ref, poly *px, poly r, const ring R)
2021{
2022 number x;
2023 int i;
2024 poly pa = *px, pp = NULL;
2025
2026 loop
2027 {
2028 i = p_LmCmp(pa, r,R);
2029 if (i > 0)
2030 {
2031 pp = pa;
2032 pIter(pa);
2033 if (pa==NULL)
2034 {
2035 pNext(pp) = r;
2036 break;
2037 }
2038 }
2039 else
2040 {
2041 if (i == 0)
2042 {
2043 x = n_Add(pGetCoeff(pa), pGetCoeff(r),R->cf);
2044 p_LmDelete(&r,R);
2045 if (n_IsZero(x,R->cf))
2046 {
2047 p_LmDelete(&pa,R);
2048 if (pp!=NULL)
2049 pNext(pp) = p_Add_q(pa,r,R);
2050 else
2051 *px = p_Add_q(pa,r,R);
2052 }
2053 else
2054 {
2055 pp = pa;
2056 p_SetCoeff(pp,x,R);
2057 pNext(pp) = p_Add_q(pNext(pp), r, R);
2058 }
2059 }
2060 else
2061 {
2062 if (pp!=NULL)
2063 pp = pNext(pp) = r;
2064 else
2065 *px = pp = r;
2066 pNext(pp) = p_Add_q(pa, pNext(r),R);
2067 }
2068 break;
2069 }
2070 }
2071 *ref = pp;
2072}
2073
2074/* ----------------- internal 'C' stuff ------------------ */
2075
2076static void sm_ElemDelete(smpoly *r, const ring R)
2077{
2078 smpoly a = *r, b = a->n;
2079
2080 p_Delete(&a->m, R);
2081 omFreeBin((void *)a, smprec_bin);
2082 *r = b;
2083}
2084
2086{
2088 memcpy(r, a, sizeof(smprec));
2089/* r->m = pCopy(r->m); */
2090 return r;
2091}
2092
2093/*
2094* from poly to smpoly
2095* do not destroy p
2096*/
2097static smpoly sm_Poly2Smpoly(poly q, const ring R)
2098{
2099 poly pp;
2100 smpoly res, a;
2101 long x;
2102
2103 if (q == NULL)
2104 return NULL;
2106 a->pos = x = p_GetComp(q,R);
2107 a->m = q;
2108 a->e = 0;
2109 loop
2110 {
2111 p_SetComp(q,0,R);
2112 pp = q;
2113 pIter(q);
2114 if (q == NULL)
2115 {
2116 a->n = NULL;
2117 return res;
2118 }
2119 if (p_GetComp(q,R) != x)
2120 {
2121 a = a->n = (smpoly)omAllocBin(smprec_bin);
2122 pNext(pp) = NULL;
2123 a->pos = x = p_GetComp(q,R);
2124 a->m = q;
2125 a->e = 0;
2126 }
2127 }
2128}
2129
2130/*
2131* from smpoly to poly
2132* destroy a
2133*/
2134static poly sm_Smpoly2Poly(smpoly a, const ring R)
2135{
2136 smpoly b;
2137 poly res, pp, q;
2138 long x;
2139
2140 if (a == NULL)
2141 return NULL;
2142 x = a->pos;
2143 q = res = a->m;
2144 loop
2145 {
2146 p_SetComp(q,x,R);
2147 pp = q;
2148 pIter(q);
2149 if (q == NULL)
2150 break;
2151 }
2152 loop
2153 {
2154 b = a;
2155 a = a->n;
2156 omFreeBin((void *)b, smprec_bin);
2157 if (a == NULL)
2158 return res;
2159 x = a->pos;
2160 q = pNext(pp) = a->m;
2161 loop
2162 {
2163 p_SetComp(q,x,R);
2164 pp = q;
2165 pIter(q);
2166 if (q == NULL)
2167 break;
2168 }
2169 }
2170}
2171
2172/*
2173* weigth of a polynomial, for pivot strategy
2174*/
2175static float sm_PolyWeight(smpoly a, const ring R)
2176{
2177 poly p = a->m;
2178 int i;
2179 float res = (float)n_Size(pGetCoeff(p),R->cf);
2180
2181 if (pNext(p) == NULL)
2182 {
2183 for(i=rVar(R); i>0; i--)
2184 {
2185 if (p_GetExp(p,i,R) != 0) return res+1.0;
2186 }
2187 return res;
2188 }
2189 else
2190 {
2191 i = 0;
2192 res = 0.0;
2193 do
2194 {
2195 i++;
2196 res += (float)n_Size(pGetCoeff(p),R->cf);
2197 pIter(p);
2198 }
2199 while (p);
2200 return res+(float)i;
2201 }
2202}
2203
2204static BOOLEAN sm_HaveDenom(poly a, const ring R)
2205{
2206 BOOLEAN sw;
2207 number x;
2208
2209 while (a != NULL)
2210 {
2211 x = n_GetDenom(pGetCoeff(a),R->cf);
2212 sw = n_IsOne(x,R->cf);
2213 n_Delete(&x,R->cf);
2214 if (!sw)
2215 {
2216 return TRUE;
2217 }
2218 pIter(a);
2219 }
2220 return FALSE;
2221}
2222
2223static number sm_Cleardenom(ideal id, const ring R)
2224{
2225 poly a;
2226 number x,y,res=n_Init(1,R->cf);
2227 BOOLEAN sw=FALSE;
2228
2229 for (int i=0; i<IDELEMS(id); i++)
2230 {
2231 a = id->m[i];
2232 sw = sm_HaveDenom(a,R);
2233 if (sw) break;
2234 }
2235 if (!sw) return res;
2236 for (int i=0; i<IDELEMS(id); i++)
2237 {
2238 a = id->m[i];
2239 if (a!=NULL)
2240 {
2241 x = n_Copy(pGetCoeff(a),R->cf);
2242 p_Cleardenom(a, R);
2243 y = n_Div(x,pGetCoeff(a),R->cf);
2244 n_Delete(&x,R->cf);
2245 x = n_Mult(res,y,R->cf);
2246 n_Normalize(x,R->cf);
2247 n_Delete(&res,R->cf);
2248 res = x;
2249 }
2250 }
2251 return res;
2252}
2253
2254/* ----------------- gauss elimination ------------------ */
2255/* in structs.h */
2256typedef struct smnrec sm_nrec;
2257typedef sm_nrec * smnumber;
2258struct smnrec{
2259 smnumber n; // the next element
2260 int pos; // position
2261 number m; // the element
2262};
2263
2265
2266/* declare internal 'C' stuff */
2267static void sm_NumberDelete(smnumber *, const ring R);
2269static smnumber sm_Poly2Smnumber(poly, const ring);
2270static poly sm_Smnumber2Poly(number,const ring);
2271static BOOLEAN smCheckSolv(ideal);
2272
2273/* class for sparse number matrix:
2274*/
2276private:
2277 int nrows, ncols; // dimension of the problem
2278 int act; // number of unreduced columns (start: ncols)
2279 int crd; // number of reduced columns (start: 0)
2280 int tored; // border for rows to reduce
2281 int sing; // indicator for singular problem
2282 int rpiv; // row-position of the pivot
2283 int *perm; // permutation of rows
2284 number *sol; // field for solution
2285 int *wrw, *wcl; // weights of rows and columns
2286 smnumber * m_act; // unreduced columns
2287 smnumber * m_res; // reduced columns (result)
2288 smnumber * m_row; // reduced part of rows
2289 smnumber red; // row to reduce
2290 smnumber piv; // pivot
2291 smnumber dumm; // allocated dummy
2292 ring _R;
2293 void smColToRow();
2294 void smRowToCol();
2295 void smSelectPR();
2296 void smRealPivot();
2297 void smZeroToredElim();
2298 void smGElim();
2299 void smAllDel();
2300public:
2301 sparse_number_mat(ideal, const ring);
2303 int smIsSing() { return sing; }
2304 void smTriangular();
2305 void smSolv();
2306 ideal smRes2Ideal();
2307};
2308
2309/* ----------------- basics (used from 'C') ------------------ */
2310/*2
2311* returns the solution of a linear equation
2312* solution of M*x = r (M has dimension nXn) =>
2313* I = module(transprose(M)) + r*gen(n+1)
2314* uses Gauss-elimination
2315*/
2316ideal sm_CallSolv(ideal I, const ring R)
2317{
2318 sparse_number_mat *linsolv;
2319 ring tmpR;
2320 ideal rr;
2321
2322 if (id_IsConstant(I,R)==FALSE)
2323 {
2324 WerrorS("symbol in equation");
2325 return NULL;
2326 }
2327 I->rank = id_RankFreeModule(I,R);
2328 if (smCheckSolv(I)) return NULL;
2329 tmpR=sm_RingChange(R,1);
2330 rr=idrCopyR(I,R, tmpR);
2331 linsolv = new sparse_number_mat(rr,tmpR);
2332 rr=NULL;
2333 linsolv->smTriangular();
2334 if (linsolv->smIsSing() == 0)
2335 {
2336 linsolv->smSolv();
2337 rr = linsolv->smRes2Ideal();
2338 }
2339 else
2340 WerrorS("singular problem for linsolv");
2341 delete linsolv;
2342 if (rr!=NULL)
2343 rr = idrMoveR(rr,tmpR,R);
2344 sm_KillModifiedRing(tmpR);
2345 return rr;
2346}
2347
2348/*
2349* constructor, destroy smat
2350*/
2352{
2353 int i;
2354 poly* pmat;
2355 _R=R;
2356
2357 crd = sing = 0;
2358 act = ncols = smat->ncols;
2359 tored = nrows = smat->rank;
2360 i = tored+1;
2361 perm = (int *)omAlloc(sizeof(int)*i);
2362 m_row = (smnumber *)omAlloc0(sizeof(smnumber)*i);
2363 wrw = (int *)omAlloc(sizeof(int)*i);
2364 i = ncols+1;
2365 wcl = (int *)omAlloc(sizeof(int)*i);
2366 m_act = (smnumber *)omAlloc(sizeof(smnumber)*i);
2367 m_res = (smnumber *)omAlloc0(sizeof(smnumber)*i);
2369 pmat = smat->m;
2370 for(i=ncols; i; i--)
2371 {
2372 m_act[i] = sm_Poly2Smnumber(pmat[i-1],_R);
2373 }
2374 omFreeSize((ADDRESS)pmat,smat->ncols*sizeof(poly));
2376}
2377
2378/*
2379* destructor
2380*/
2382{
2383 int i;
2385 i = ncols+1;
2386 omFreeSize((ADDRESS)m_res, sizeof(smnumber)*i);
2387 omFreeSize((ADDRESS)m_act, sizeof(smnumber)*i);
2388 omFreeSize((ADDRESS)wcl, sizeof(int)*i);
2389 i = nrows+1;
2390 omFreeSize((ADDRESS)wrw, sizeof(int)*i);
2391 omFreeSize((ADDRESS)m_row, sizeof(smnumber)*i);
2392 omFreeSize((ADDRESS)perm, sizeof(int)*i);
2393}
2394
2395/*
2396* triangularization by Gauss-elimination
2397*/
2399{
2400 tored--;
2401 this->smZeroToredElim();
2402 if (sing != 0) return;
2403 while (act > 1)
2404 {
2405 this->smRealPivot();
2406 this->smSelectPR();
2407 this->smGElim();
2408 crd++;
2409 this->smColToRow();
2410 act--;
2411 this->smRowToCol();
2412 this->smZeroToredElim();
2413 if (sing != 0) return;
2414 }
2415 if (TEST_OPT_PROT) PrintS(".\n");
2416 piv = m_act[1];
2417 rpiv = piv->pos;
2418 m_act[1] = piv->n;
2419 piv->n = NULL;
2420 crd++;
2421 this->smColToRow();
2422 act--;
2423 this->smRowToCol();
2424}
2425
2426/*
2427* solve the triangular system
2428*/
2430{
2431 int i, j;
2432 number x, y, z;
2433 smnumber s, d, r = m_row[nrows];
2434
2435 m_row[nrows] = NULL;
2436 sol = (number *)omAlloc0(sizeof(number)*(crd+1));
2437 while (r != NULL) // expand the rigth hand side
2438 {
2439 sol[r->pos] = r->m;
2440 s = r;
2441 r = r->n;
2443 }
2444 i = crd; // solve triangular system
2445 if (sol[i] != NULL)
2446 {
2447 x = sol[i];
2448 sol[i] = n_Div(x, m_res[i]->m,_R->cf);
2449 n_Delete(&x,_R->cf);
2450 }
2451 i--;
2452 while (i > 0)
2453 {
2454 x = NULL;
2455 d = m_res[i];
2456 s = d->n;
2457 while (s != NULL)
2458 {
2459 j = s->pos;
2460 if (sol[j] != NULL)
2461 {
2462 z = n_Mult(sol[j], s->m,_R->cf);
2463 if (x != NULL)
2464 {
2465 y = x;
2466 x = n_Sub(y, z,_R->cf);
2467 n_Delete(&y,_R->cf);
2468 n_Delete(&z,_R->cf);
2469 }
2470 else
2471 x = n_InpNeg(z,_R->cf);
2472 }
2473 s = s->n;
2474 }
2475 if (sol[i] != NULL)
2476 {
2477 if (x != NULL)
2478 {
2479 y = n_Add(x, sol[i],_R->cf);
2480 n_Delete(&x,_R->cf);
2481 if (n_IsZero(y,_R->cf))
2482 {
2483 n_Delete(&sol[i],_R->cf);
2484 sol[i] = NULL;
2485 }
2486 else
2487 sol[i] = y;
2488 }
2489 }
2490 else
2491 sol[i] = x;
2492 if (sol[i] != NULL)
2493 {
2494 x = sol[i];
2495 sol[i] = n_Div(x, d->m,_R->cf);
2496 n_Delete(&x,_R->cf);
2497 }
2498 i--;
2499 }
2500 this->smAllDel();
2501}
2502
2503/*
2504* transform the result to an ideal
2505*/
2507{
2508 int i, j;
2509 ideal res = idInit(crd, 1);
2510
2511 for (i=crd; i; i--)
2512 {
2513 j = perm[i]-1;
2514 res->m[j] = sm_Smnumber2Poly(sol[i],_R);
2515 }
2516 omFreeSize((ADDRESS)sol, sizeof(number)*(crd+1));
2517 return res;
2518}
2519
2520/* ----------------- pivot method ------------------ */
2521
2522/*
2523* compute pivot
2524*/
2526{
2527 smnumber a;
2528 number x, xo;
2529 int i, copt, ropt;
2530
2531 xo=n_Init(0,_R->cf);
2532 for (i=act; i; i--)
2533 {
2534 a = m_act[i];
2535 while ((a!=NULL) && (a->pos<=tored))
2536 {
2537 x = a->m;
2538 if (n_GreaterZero(x,_R->cf))
2539 {
2540 if (n_Greater(x,xo,_R->cf))
2541 {
2542 n_Delete(&xo,_R->cf);
2543 xo = n_Copy(x,_R->cf);
2544 copt = i;
2545 ropt = a->pos;
2546 }
2547 }
2548 else
2549 {
2550 xo = n_InpNeg(xo,_R->cf);
2551 if (n_Greater(xo,x,_R->cf))
2552 {
2553 n_Delete(&xo,_R->cf);
2554 xo = n_Copy(x,_R->cf);
2555 copt = i;
2556 ropt = a->pos;
2557 }
2558 xo = n_InpNeg(xo,_R->cf);
2559 }
2560 a = a->n;
2561 }
2562 }
2563 rpiv = ropt;
2564 if (copt != act)
2565 {
2566 a = m_act[act];
2567 m_act[act] = m_act[copt];
2568 m_act[copt] = a;
2569 }
2570 n_Delete(&xo,_R->cf);
2571}
2572
2573/* ----------------- elimination ------------------ */
2574
2575/* one step of Gauss-elimination */
2577{
2578 number p = n_Invers(piv->m,_R->cf); // pivotelement
2579 smnumber c = m_act[act]; // pivotcolumn
2580 smnumber r = red; // row to reduce
2581 smnumber res, a, b;
2582 number w, ha, hb;
2583
2584 if ((c == NULL) || (r == NULL))
2585 {
2586 while (r!=NULL) sm_NumberDelete(&r,_R);
2587 return;
2588 }
2589 do
2590 {
2591 a = m_act[r->pos];
2592 res = dumm;
2593 res->n = NULL;
2594 b = c;
2595 w = n_Mult(r->m, p,_R->cf);
2596 n_Delete(&r->m,_R->cf);
2597 r->m = w;
2598 loop // combine the chains a and b: a + w*b
2599 {
2600 if (a == NULL)
2601 {
2602 do
2603 {
2604 res = res->n = smNumberCopy(b);
2605 res->m = n_Mult(b->m, w,_R->cf);
2606 b = b->n;
2607 } while (b != NULL);
2608 break;
2609 }
2610 if (a->pos < b->pos)
2611 {
2612 res = res->n = a;
2613 a = a->n;
2614 }
2615 else if (a->pos > b->pos)
2616 {
2617 res = res->n = smNumberCopy(b);
2618 res->m = n_Mult(b->m, w,_R->cf);
2619 b = b->n;
2620 }
2621 else
2622 {
2623 hb = n_Mult(b->m, w,_R->cf);
2624 ha = n_Add(a->m, hb,_R->cf);
2625 n_Delete(&hb,_R->cf);
2626 n_Delete(&a->m,_R->cf);
2627 if (n_IsZero(ha,_R->cf))
2628 {
2629 sm_NumberDelete(&a,_R);
2630 }
2631 else
2632 {
2633 a->m = ha;
2634 res = res->n = a;
2635 a = a->n;
2636 }
2637 b = b->n;
2638 }
2639 if (b == NULL)
2640 {
2641 res->n = a;
2642 break;
2643 }
2644 }
2645 m_act[r->pos] = dumm->n;
2646 sm_NumberDelete(&r,_R);
2647 } while (r != NULL);
2648 n_Delete(&p,_R->cf);
2649}
2650
2651/* ----------------- transfer ------------------ */
2652
2653/*
2654* select the pivotrow and store it to red and piv
2655*/
2657{
2658 smnumber b = dumm;
2659 smnumber a, ap;
2660 int i;
2661
2662 if (TEST_OPT_PROT)
2663 {
2664 if ((crd+1)%10)
2665 PrintS(".");
2666 else
2667 PrintS(".\n");
2668 }
2669 a = m_act[act];
2670 if (a->pos < rpiv)
2671 {
2672 do
2673 {
2674 ap = a;
2675 a = a->n;
2676 } while (a->pos < rpiv);
2677 ap->n = a->n;
2678 }
2679 else
2680 m_act[act] = a->n;
2681 piv = a;
2682 a->n = NULL;
2683 for (i=1; i<act; i++)
2684 {
2685 a = m_act[i];
2686 if (a->pos < rpiv)
2687 {
2688 loop
2689 {
2690 ap = a;
2691 a = a->n;
2692 if ((a == NULL) || (a->pos > rpiv))
2693 break;
2694 if (a->pos == rpiv)
2695 {
2696 ap->n = a->n;
2697 a->m = n_InpNeg(a->m,_R->cf);
2698 b = b->n = a;
2699 b->pos = i;
2700 break;
2701 }
2702 }
2703 }
2704 else if (a->pos == rpiv)
2705 {
2706 m_act[i] = a->n;
2707 a->m = n_InpNeg(a->m,_R->cf);
2708 b = b->n = a;
2709 b->pos = i;
2710 }
2711 }
2712 b->n = NULL;
2713 red = dumm->n;
2714}
2715
2716/*
2717* store the pivotcol in m_row
2718* m_act[cols][pos(rows)] => m_row[rows][pos(cols)]
2719*/
2721{
2722 smnumber c = m_act[act];
2723 smnumber h;
2724
2725 while (c != NULL)
2726 {
2727 h = c;
2728 c = c->n;
2729 h->n = m_row[h->pos];
2730 m_row[h->pos] = h;
2731 h->pos = crd;
2732 }
2733}
2734
2735/*
2736* store the pivot and the assosiated row in m_row
2737* to m_res (result):
2738* piv + m_row[rows][pos(cols)] => m_res[cols][pos(rows)]
2739*/
2741{
2742 smnumber r = m_row[rpiv];
2743 smnumber a, ap, h;
2744
2745 m_row[rpiv] = NULL;
2746 perm[crd] = rpiv;
2747 piv->pos = crd;
2748 m_res[crd] = piv;
2749 while (r != NULL)
2750 {
2751 ap = m_res[r->pos];
2752 loop
2753 {
2754 a = ap->n;
2755 if (a == NULL)
2756 {
2757 ap->n = h = r;
2758 r = r->n;
2759 h->n = a;
2760 h->pos = crd;
2761 break;
2762 }
2763 ap = a;
2764 }
2765 }
2766}
2767
2768/* ----------------- C++ stuff ------------------ */
2769
2770/*
2771* check singularity
2772*/
2774{
2775 smnumber a;
2776 int i = act;
2777
2778 loop
2779 {
2780 if (i == 0) return;
2781 a = m_act[i];
2782 if ((a==NULL) || (a->pos>tored))
2783 {
2784 sing = 1;
2785 this->smAllDel();
2786 return;
2787 }
2788 i--;
2789 }
2790}
2791
2792/*
2793* delete all smnumber's
2794*/
2796{
2797 smnumber a;
2798 int i;
2799
2800 for (i=act; i; i--)
2801 {
2802 a = m_act[i];
2803 while (a != NULL)
2804 sm_NumberDelete(&a,_R);
2805 }
2806 for (i=crd; i; i--)
2807 {
2808 a = m_res[i];
2809 while (a != NULL)
2810 sm_NumberDelete(&a,_R);
2811 }
2812 if (act)
2813 {
2814 for (i=nrows; i; i--)
2815 {
2816 a = m_row[i];
2817 while (a != NULL)
2818 sm_NumberDelete(&a,_R);
2819 }
2820 }
2821}
2822
2823/* ----------------- internal 'C' stuff ------------------ */
2824
2825static void sm_NumberDelete(smnumber *r, const ring R)
2826{
2827 smnumber a = *r, b = a->n;
2828
2829 n_Delete(&a->m,R->cf);
2831 *r = b;
2832}
2833
2835{
2837 memcpy(r, a, sizeof(smnrec));
2838 return r;
2839}
2840
2841/*
2842* from poly to smnumber
2843* do not destroy p
2844*/
2845static smnumber sm_Poly2Smnumber(poly q, const ring R)
2846{
2847 smnumber a, res;
2848 poly p = q;
2849
2850 if (p == NULL)
2851 return NULL;
2853 a->pos = p_GetComp(p,R);
2854 a->m = pGetCoeff(p);
2855 pGetCoeff(p)=NULL;
2856 loop
2857 {
2858 pIter(p);
2859 if (p == NULL)
2860 {
2861 p_Delete(&q,R);
2862 a->n = NULL;
2863 return res;
2864 }
2865 a = a->n = (smnumber)omAllocBin(smnrec_bin);
2866 a->pos = p_GetComp(p,R);
2867 a->m = pGetCoeff(p);
2868 pGetCoeff(p)=NULL;
2869 }
2870}
2871
2872/*
2873* from smnumber to poly
2874* destroy a
2875*/
2876static poly sm_Smnumber2Poly(number a, const ring R)
2877{
2878 poly res;
2879
2880 if (a == NULL) return NULL;
2881 res = p_Init(R);
2882 pSetCoeff0(res, a);
2883 return res;
2884}
2885
2886/*2
2887* check the input
2888*/
2889static BOOLEAN smCheckSolv(ideal I)
2890{ int i = I->ncols;
2891 if ((i == 0) || (i != I->rank-1))
2892 {
2893 WerrorS("wrong dimensions for linsolv");
2894 return TRUE;
2895 }
2896 for(;i;i--)
2897 {
2898 if(I->m[i-1] == NULL)
2899 {
2900 WerrorS("singular input for linsolv");
2901 return TRUE;
2902 }
2903 }
2904 return FALSE;
2905}
All the auxiliary stuff.
static int si_max(const int a, const int b)
Definition: auxiliary.h:124
int BOOLEAN
Definition: auxiliary.h:87
#define TRUE
Definition: auxiliary.h:100
#define FALSE
Definition: auxiliary.h:96
void * ADDRESS
Definition: auxiliary.h:119
CanonicalForm FACTORY_PUBLIC pp(const CanonicalForm &)
CanonicalForm pp ( const CanonicalForm & f )
Definition: cf_gcd.cc:676
int l
Definition: cfEzgcd.cc:100
int m
Definition: cfEzgcd.cc:128
int i
Definition: cfEzgcd.cc:132
int k
Definition: cfEzgcd.cc:99
Variable x
Definition: cfModGcd.cc:4082
int p
Definition: cfModGcd.cc:4078
CanonicalForm gp
Definition: cfModGcd.cc:4102
CanonicalForm b
Definition: cfModGcd.cc:4103
static CanonicalForm bound(const CFMatrix &M)
Definition: cf_linsys.cc:460
FILE * f
Definition: checklibs.c:9
Definition: intvec.h:23
sparse_mat(ideal, const ring)
Definition: sparsmat.cc:386
int smGetSign()
Definition: sparsmat.cc:171
void smToIntvec(intvec *)
Definition: sparsmat.cc:464
void smFinalMult()
Definition: sparsmat.cc:1384
int smGetRed()
Definition: sparsmat.cc:173
void smInitPerm()
Definition: sparsmat.cc:1548
poly smMultPoly(smpoly)
Definition: sparsmat.cc:1453
void smWeights()
Definition: sparsmat.cc:610
void smCopToRes()
Definition: sparsmat.cc:1203
void smToredElim()
Definition: sparsmat.cc:1164
void smColToRow()
Definition: sparsmat.cc:1081
void smActDel()
Definition: sparsmat.cc:1475
smpoly * m_row
Definition: sparsmat.cc:139
float * wcl
Definition: sparsmat.cc:136
smpoly * m_res
Definition: sparsmat.cc:138
smpoly dumm
Definition: sparsmat.cc:142
void sm1Elim()
Definition: sparsmat.cc:799
smpoly red
Definition: sparsmat.cc:140
smpoly * smGetAct()
Definition: sparsmat.cc:172
void smPivDel()
Definition: sparsmat.cc:1506
float * wrw
Definition: sparsmat.cc:136
smpoly piv
Definition: sparsmat.cc:141
smpoly * m_act
Definition: sparsmat.cc:137
void smColDel()
Definition: sparsmat.cc:1493
poly smDet()
Definition: sparsmat.cc:475
void smNormalize()
Definition: sparsmat.cc:1433
smpoly oldpiv
Definition: sparsmat.cc:141
int smCheckNormalize()
Definition: sparsmat.cc:1413
void smPivot()
Definition: sparsmat.cc:642
void smNewBareiss(int, int)
Definition: sparsmat.cc:549
void smSparseHomog()
int normalize
Definition: sparsmat.cc:133
void smSelectPR()
Definition: sparsmat.cc:1017
void smSign()
Definition: sparsmat.cc:1520
void smRowToCol()
Definition: sparsmat.cc:1101
void smHElim()
Definition: sparsmat.cc:878
void smMultCol()
Definition: sparsmat.cc:1359
void smZeroElim()
Definition: sparsmat.cc:1134
ideal smRes2Mod()
Definition: sparsmat.cc:448
float wpoints
Definition: sparsmat.cc:135
void smNewWeights()
Definition: sparsmat.cc:699
int * perm
Definition: sparsmat.cc:134
void smNewPivot()
Definition: sparsmat.cc:737
smnumber * m_act
Definition: sparsmat.cc:2286
ideal smRes2Ideal()
Definition: sparsmat.cc:2506
smnumber * m_res
Definition: sparsmat.cc:2287
void smZeroToredElim()
Definition: sparsmat.cc:2773
smnumber * m_row
Definition: sparsmat.cc:2288
sparse_number_mat(ideal, const ring)
Definition: sparsmat.cc:2351
static FORCE_INLINE number n_Mult(number a, number b, const coeffs r)
return the product of 'a' and 'b', i.e., a*b
Definition: coeffs.h:633
static FORCE_INLINE number n_Copy(number n, const coeffs r)
return a copy of 'n'
Definition: coeffs.h:448
static FORCE_INLINE number n_Add(number a, number b, const coeffs r)
return the sum of 'a' and 'b', i.e., a+b
Definition: coeffs.h:647
static FORCE_INLINE number n_GetDenom(number &n, const coeffs r)
return the denominator of n (if elements of r are by nature not fractional, result is 1)
Definition: coeffs.h:600
static FORCE_INLINE number n_Invers(number a, const coeffs r)
return the multiplicative inverse of 'a'; raise an error if 'a' is not invertible
Definition: coeffs.h:561
static FORCE_INLINE BOOLEAN n_GreaterZero(number n, const coeffs r)
ordered fields: TRUE iff 'n' is positive; in Z/pZ: TRUE iff 0 < m <= roundedBelow(p/2),...
Definition: coeffs.h:491
static FORCE_INLINE number n_InpNeg(number n, const coeffs r)
in-place negation of n MUST BE USED: n = n_InpNeg(n) (no copy is returned)
Definition: coeffs.h:554
static FORCE_INLINE number n_Div(number a, number b, const coeffs r)
return the quotient of 'a' and 'b', i.e., a/b; raises an error if 'b' is not invertible in r exceptio...
Definition: coeffs.h:612
static FORCE_INLINE BOOLEAN n_Greater(number a, number b, const coeffs r)
ordered fields: TRUE iff 'a' is larger than 'b'; in Z/pZ: TRUE iff la > lb, where la and lb are the l...
Definition: coeffs.h:508
static FORCE_INLINE BOOLEAN n_IsZero(number n, const coeffs r)
TRUE iff 'n' represents the zero element.
Definition: coeffs.h:461
static FORCE_INLINE int n_Size(number n, const coeffs r)
return a non-negative measure for the complexity of n; return 0 only when n represents zero; (used fo...
Definition: coeffs.h:567
static FORCE_INLINE number n_Sub(number a, number b, const coeffs r)
return the difference of 'a' and 'b', i.e., a-b
Definition: coeffs.h:652
static FORCE_INLINE void n_Delete(number *p, const coeffs r)
delete 'p'
Definition: coeffs.h:452
static FORCE_INLINE number n_Init(long i, const coeffs r)
a number representing i in the given coeff field/ring r
Definition: coeffs.h:535
static FORCE_INLINE BOOLEAN n_Equal(number a, number b, const coeffs r)
TRUE iff 'a' and 'b' represent the same number; they may have different representations.
Definition: coeffs.h:457
static FORCE_INLINE void n_Normalize(number &n, const coeffs r)
inplace-normalization of n; produces some canonical representation of n;
Definition: coeffs.h:575
static FORCE_INLINE BOOLEAN n_IsOne(number n, const coeffs r)
TRUE iff 'n' represents the one element.
Definition: coeffs.h:465
static BOOLEAN pa(leftv res, leftv args)
Definition: cohomo.cc:3723
static BOOLEAN pb(leftv res, leftv args)
Definition: cohomo.cc:3747
#define Print
Definition: emacs.cc:80
const CanonicalForm int s
Definition: facAbsFact.cc:51
const CanonicalForm int const CFList const Variable & y
Definition: facAbsFact.cc:53
CanonicalForm res
Definition: facAbsFact.cc:60
const CanonicalForm & w
Definition: facAbsFact.cc:51
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:39
int j
Definition: facHensel.cc:110
void WerrorS(const char *s)
Definition: feFopen.cc:24
#define STATIC_VAR
Definition: globaldefs.h:7
#define VAR
Definition: globaldefs.h:5
STATIC_VAR Poly * h
Definition: janet.cc:971
void kBucketDestroy(kBucket_pt *bucket_pt)
Definition: kbuckets.cc:216
void kBucketInit(kBucket_pt bucket, poly lm, int length)
Definition: kbuckets.cc:493
poly kBucketExtractLm(kBucket_pt bucket)
Definition: kbuckets.cc:511
kBucket_pt kBucketCreate(const ring bucket_ring)
Creation/Destruction of buckets.
Definition: kbuckets.cc:209
void kBucket_Add_q(kBucket_pt bucket, poly q, int *l)
Add to Bucket a poly ,i.e. Bpoly == q+Bpoly.
Definition: kbuckets.cc:660
#define p_GetComp(p, r)
Definition: monomials.h:64
#define pIter(p)
Definition: monomials.h:37
#define pNext(p)
Definition: monomials.h:36
#define pSetCoeff0(p, n)
Definition: monomials.h:59
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition: monomials.h:44
Definition: ap.h:40
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
#define omAlloc(size)
Definition: omAllocDecl.h:210
#define omAllocBin(bin)
Definition: omAllocDecl.h:205
#define omFree(addr)
Definition: omAllocDecl.h:261
#define omAlloc0(size)
Definition: omAllocDecl.h:211
#define omFreeBin(addr, bin)
Definition: omAllocDecl.h:259
#define omGetSpecBin(size)
Definition: omBin.h:11
#define NULL
Definition: omList.c:12
omBin_t * omBin
Definition: omStructs.h:12
#define TEST_OPT_PROT
Definition: options.h:104
#define TEST_OPT_NOT_BUCKETS
Definition: options.h:106
void p_Normalize(poly p, const ring r)
Definition: p_polys.cc:3813
poly p_Cleardenom(poly p, const ring r)
Definition: p_polys.cc:2845
static poly p_Neg(poly p, const ring r)
Definition: p_polys.h:1105
static int pLength(poly a)
Definition: p_polys.h:188
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:934
static void p_LmDelete(poly p, const ring r)
Definition: p_polys.h:721
static void p_ExpVectorAdd(poly p1, poly p2, const ring r)
Definition: p_polys.h:1409
static long p_SubExp(poly p, int v, long ee, ring r)
Definition: p_polys.h:611
static poly pp_Mult_mm(poly p, poly m, const ring r)
Definition: p_polys.h:1029
static unsigned long p_SetExp(poly p, const unsigned long e, const unsigned long iBitmask, const int VarOffset)
set a single variable exponent @Note: VarOffset encodes the position in p->exp
Definition: p_polys.h:486
static void p_ExpVectorDiff(poly pr, poly p1, poly p2, const ring r)
Definition: p_polys.h:1472
static unsigned long p_SetComp(poly p, unsigned long c, ring r)
Definition: p_polys.h:245
static void p_Setm(poly p, const ring r)
Definition: p_polys.h:231
static number p_SetCoeff(poly p, number n, ring r)
Definition: p_polys.h:410
static BOOLEAN p_LmIsConstantComp(const poly p, const ring r)
Definition: p_polys.h:1004
static int p_LmCmp(poly p, poly q, const ring r)
Definition: p_polys.h:1578
static long p_GetExp(const poly p, const unsigned long iBitmask, const int VarOffset)
get a single variable exponent @Note: the integer VarOffset encodes:
Definition: p_polys.h:467
static BOOLEAN p_LmDivisibleByNoComp(poly a, poly b, const ring r)
Definition: p_polys.h:1875
static poly p_New(const ring, omBin bin)
Definition: p_polys.h:662
static long p_MaxComp(poly p, ring lmRing, ring tailRing)
Definition: p_polys.h:290
static poly p_Mult_nn(poly p, number n, const ring r)
Definition: p_polys.h:956
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:899
static poly pp_Mult_qq(poly p, poly q, const ring r)
Definition: p_polys.h:1149
static void p_LmFree(poly p, ring)
Definition: p_polys.h:681
static poly p_Init(const ring r, omBin bin)
Definition: p_polys.h:1318
static poly pp_Mult_Coeff_mm_DivSelect(poly p, const poly m, const ring r)
Definition: p_polys.h:1088
#define p_LmTest(p, r)
Definition: p_polys.h:160
#define p_Test(p, r)
Definition: p_polys.h:159
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
poly prMoveR(poly &p, ring src_r, ring dest_r)
Definition: prCopy.cc:90
ideal idrMoveR(ideal &id, ring src_r, ring dest_r)
Definition: prCopy.cc:248
ideal idrCopyR(ideal id, ring src_r, ring dest_r)
Definition: prCopy.cc:192
ideal idrCopyR_NoSort(ideal id, ring src_r, ring dest_r)
Definition: prCopy.cc:205
void PrintS(const char *s)
Definition: reporter.cc:284
void Werror(const char *fmt,...)
Definition: reporter.cc:189
BOOLEAN rComplete(ring r, int force)
this needs to be called whenever a new ring is created: new fields in ring are created (like VarOffse...
Definition: ring.cc:3450
void rKillModifiedRing(ring r)
Definition: ring.cc:3059
ring rCopy0(const ring r, BOOLEAN copy_qideal, BOOLEAN copy_ordering)
Definition: ring.cc:1421
static BOOLEAN rOrd_is_Comp_dp(const ring r)
Definition: ring.h:772
rRingOrder_t
order stuff
Definition: ring.h:68
@ ringorder_dp
Definition: ring.h:78
@ ringorder_c
Definition: ring.h:72
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:592
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:35
void id_Delete(ideal *h, ring r)
deletes an ideal/module/matrix
BOOLEAN id_IsConstant(ideal id, const ring r)
test if the ideal has only constant polynomials NOTE: zero ideal/module is also constant
long id_RankFreeModule(ideal s, ring lmRing, ring tailRing)
return the maximal component number found in any polynomial in s
VAR omBin sip_sideal_bin
Definition: simpleideals.cc:27
#define IDELEMS(i)
Definition: simpleideals.h:23
#define R
Definition: sirandom.c:27
#define M
Definition: sirandom.c:25
smnumber n
Definition: sparsmat.cc:2259
int e
Definition: sparsmat.cc:51
STATIC_VAR omBin smnrec_bin
Definition: sparsmat.cc:2264
float f
Definition: sparsmat.cc:53
static void sm_ElemDelete(smpoly *, const ring)
Definition: sparsmat.cc:2076
static smnumber sm_Poly2Smnumber(poly, const ring)
Definition: sparsmat.cc:2845
sm_nrec * smnumber
Definition: sparsmat.cc:2257
static void smMinSelect(long *, int, int)
Definition: sparsmat.cc:236
poly sm_MultDiv(poly a, poly b, const poly c, const ring R)
Definition: sparsmat.cc:1759
static BOOLEAN sm_HaveDenom(poly, const ring)
Definition: sparsmat.cc:2204
#define SM_MIN_LENGTH_BUCKET
Definition: sparsmat.cc:40
long sm_ExpBound(ideal m, int di, int ra, int t, const ring currRing)
Definition: sparsmat.cc:188
static void sm_CombineChain(poly *, poly, const ring)
Definition: sparsmat.cc:1973
poly sm_CallDet(ideal I, const ring R)
Definition: sparsmat.cc:302
static poly pp_Mult_Coeff_mm_DivSelect_MultDiv(poly p, int &lp, poly m, poly a, poly b, const ring currRing)
Definition: sparsmat.cc:76
static void sm_NumberDelete(smnumber *, const ring R)
Definition: sparsmat.cc:2825
int pos
Definition: sparsmat.cc:50
void sm_CallBareiss(ideal I, int x, int y, ideal &M, intvec **iv, const ring R)
Definition: sparsmat.cc:347
static float sm_PolyWeight(smpoly, const ring)
Definition: sparsmat.cc:2175
smpoly n
Definition: sparsmat.cc:49
static BOOLEAN smSmaller(poly, poly)
Definition: sparsmat.cc:1962
ideal sm_CallSolv(ideal I, const ring R)
Definition: sparsmat.cc:2316
static smnumber smNumberCopy(smnumber)
Definition: sparsmat.cc:2834
int pos
Definition: sparsmat.cc:2260
VAR omBin smprec_bin
Definition: sparsmat.cc:74
number m
Definition: sparsmat.cc:2261
static poly sm_Smpoly2Poly(smpoly, const ring)
Definition: sparsmat.cc:2134
sm_prec * smpoly
Definition: sparsmat.cc:46
static smpoly smElemCopy(smpoly)
Definition: sparsmat.cc:2085
static number sm_Cleardenom(ideal, const ring)
Definition: sparsmat.cc:2223
static poly sm_SelectCopy_ExpMultDiv(poly p, poly m, poly a, poly b, const ring currRing)
Definition: sparsmat.cc:98
poly m
Definition: sparsmat.cc:52
static void sm_ExpMultDiv(poly, const poly, const poly, const ring)
Definition: sparsmat.cc:1931
ring sm_RingChange(const ring origR, long bound)
Definition: sparsmat.cc:258
static BOOLEAN smCheckSolv(ideal)
Definition: sparsmat.cc:2889
void sm_SpecialPolyDiv(poly a, poly b, const ring R)
Definition: sparsmat.cc:1840
static BOOLEAN sm_IsNegQuot(poly, const poly, const poly, const ring)
Definition: sparsmat.cc:1906
static void sm_ExactPolyDiv(poly, poly, const ring)
Definition: sparsmat.cc:1851
static void sm_FindRef(poly *, poly *, poly, const ring)
Definition: sparsmat.cc:2020
void sm_KillModifiedRing(ring r)
Definition: sparsmat.cc:289
static smpoly sm_Poly2Smpoly(poly, const ring)
Definition: sparsmat.cc:2097
static void sm_PolyDivN(poly, const number, const ring)
Definition: sparsmat.cc:1949
static poly sm_Smnumber2Poly(number, const ring)
Definition: sparsmat.cc:2876
#define SM_DIV
Definition: sparsmat.h:24
#define SM_MULT
Definition: sparsmat.h:23
#define loop
Definition: structs.h:75