source: git/Singular/spolys0.cc @ 63be42

spielwiese
Last change on this file since 63be42 was 7a5e7c, checked in by Olaf Bachmann <obachman@…>, 26 years ago
* polys.cc: pComp for simple orderings based on only three routines (no variation for # of variables) * spolys0.cc (spGSpolyLoop): globally replaced by spSpolyLoop_General git-svn-id: file:///usr/local/Singular/svn/trunk@2135 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 14.1 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: spolys0.cc,v 1.12 1998-06-12 17:41:39 obachman Exp $ */
5
6/*
7* ABSTRACT - s-polynomials and reduction in general
8*/
9
10#include <string.h>
11#include "mod2.h"
12#include "tok.h"
13#include "mmemory.h"
14#include "numbers.h"
15#include "polys.h"
16#include "spolys0.h"
17
18/*0 implementation*/
19
20/*
21* input - output: a, b
22* returns:
23*   a := a/gcd(a,b), b := b/gcd(a,b)
24*   and return value
25*       0  ->  a != 1,  b != 1
26*       1  ->  a == 1,  b != 1
27*       2  ->  a != 1,  b == 1
28*       3  ->  a == 1,  b == 1
29*   this value is used to control the spolys
30*/
31static int spCheckCoeff(number *a, number *b)
32{
33  int c = 0;
34  number an = *a, bn = *b;
35  nTest(an);
36  nTest(bn);
37  number cn = nGcd(an, bn);
38
39  if(nIsOne(cn))
40  {
41    an = nCopy(an);
42    bn = nCopy(bn);
43  }
44  else
45  {
46    an = nIntDiv(an, cn);
47    bn = nIntDiv(bn, cn);
48  }
49  nDelete(&cn);
50  if (nIsOne(an))
51  {
52    c = 1;
53  }
54  if (nIsOne(bn))
55  {
56    c += 2;
57  }
58  *a = an;
59  *b = bn;
60  return c;
61}
62
63/*2
64* transforms a1 to a polynomial
65*/
66void spModuleToPoly(poly a1)
67{
68#ifdef PDEBUG
69  poly t=a1;
70  pTest(t);
71#endif
72  do
73  {
74    pSetComp(a1,0);
75    pIter(a1);
76  } while (a1!=NULL);
77#ifdef PDEBUG
78  pTest(t);
79#endif
80}
81
82/*2
83* assume m = L(m) and Lc(m) = 1
84* pNext(m) = result = p*m
85* do not destroy p
86*/
87static void spGMultCopy0(poly p, poly m, poly spNoether)
88{
89  poly a, b;
90
91  a = m;
92  if (spNoether==NULL)
93  {
94    do
95    {
96      a = pNext(a) = pNew();
97      pCopyAddFast(a,p,m);
98      pSetCoeff0(a,nCopy(pGetCoeff(a)));
99      pIter(p);
100    }
101    while (p!=NULL);
102    pNext(a) = NULL;
103    return;
104  }
105  else
106  {
107    do
108    {
109      b = pNext(a) = pNew();
110      spMemcpy(b,p);
111      spMonAdd(b,m);
112      if (pComp0(b, spNoether) == -1)
113      {
114        pFree1(b);
115        pNext(a) = NULL;
116        return;
117      }
118      a = b;
119      pSetCoeff0(a,nCopy(pGetCoeff(a)));
120      pIter(p);
121    } while (p!=NULL);
122    pNext(a) = NULL;
123  }
124}
125
126/*2
127* assume m = L(m) and Lc(m) = -1
128* pNext(n) = result = p*m
129* do not destroy p
130*/
131static void spGMultCopy1(poly p, poly m, poly n,poly spNoether)
132{
133  poly a, b;
134
135  a = n;
136  if (spNoether==NULL)
137  {
138    do
139    {
140      number tmp;
141      a = pNext(a) = pNew();
142      spMemcpy(a,p);
143      spMonAdd(a,m);
144      tmp=nCopy(pGetCoeff(a));
145      tmp=nNeg(tmp);
146      pSetCoeff0(a,tmp);
147      pIter(p);
148    }
149    while (p!=NULL);
150    pNext(a) = NULL;
151    return;
152  }
153  else
154  {
155    do
156    {
157      b = pNext(a) = pNew();
158      spMemcpy(b,p);
159      spMonAdd(b,m);
160      if (pComp(b, spNoether) == -1)
161      {
162        pFree1(b);
163        pNext(a) = NULL;
164        return;
165      }
166      a = b;
167      number tmp;
168      tmp=nCopy(pGetCoeff(a));
169      tmp=nNeg(tmp);
170      pSetCoeff0(a,tmp);
171      pIter(p);
172    }
173    while (p!=NULL);
174    pNext(a) = NULL;
175  }
176}
177
178/*2
179* assume m = L(m) and Lc(m) = 1
180* pNext(m) = result = a2-a1*m
181* do not destroy a1, but a2
182*/
183static void spGSpolyLoop1(poly a1, poly a2, poly m, poly spNoether)
184{
185  poly a, b, s;
186  number tb;
187  int c;
188
189  if (a2==NULL)
190  {
191    spGMultCopy1(a1, m, m,spNoether);
192    return;
193  }
194  a = m;
195  b = pNew();
196  spMemcpy(b,a1);
197  spMonAdd(b,m);
198  loop
199  {
200    c = pComp0(b, a2);
201    if (c == 1)
202    {
203      number tmp=nCopy(pGetCoeff(b));
204      tmp=nNeg(tmp);
205      pSetCoeff0(b,tmp);
206      a = pNext(a) = b;
207      pIter(a1);
208      if (a1!=NULL)
209      {
210        b = pNew();
211        spMemcpy(b,a1);
212        spMonAdd(b,m);
213      }
214      else
215      {
216        pNext(a) = a2;
217        return;
218      }
219    }
220    else if (c == -1)
221    {
222      a = pNext(a) = a2;
223      pIter(a2);
224      if (a2==NULL)
225      {
226        pFree1(b);
227        spGMultCopy1(a1, m, a,spNoether);
228        return;
229      }
230    }
231    else
232    {
233      tb = pGetCoeff(a2);
234      if (!nEqual(tb,pGetCoeff(a1)))
235      {
236        tb = nSub(tb,pGetCoeff(a1));
237        pSetCoeff(a2,tb);
238        a = pNext(a) = a2;
239        pIter(a2);
240      }
241      else
242      {
243        s = a2;
244        pIter(a2);
245        nDelete(&tb);
246        pFree1(s);
247      }
248      pIter(a1);
249      if (a1==NULL)
250      {
251        pFree1(b);
252        pNext(a) = a2;
253        return;
254      }
255      if (a2==NULL)
256      {
257        pFree1(b);
258        spGMultCopy1(a1, m, a,spNoether);
259        return;
260      }
261      spMemcpy(b,a1);
262      spMonAdd(b,m);
263    }
264  }
265}
266
267/*2
268* assume m = L(m) and Lc(m) = exp
269* pNext(n) = result = p*m
270* do not destroy p
271*/
272void spGMultCopyX(poly p, poly m, poly n, number exp, poly spNoether)
273{
274  poly a, b;
275
276  a = n;
277  if (spNoether==NULL)
278  {
279    do
280    {
281      a = pNext(a) = pNew();
282      spMemcpy(a,p);
283      spMonAdd(a,m);
284      pSetCoeff0(a,nMult(pGetCoeff(p),exp));
285      pIter(p);
286    } while (p!=NULL);
287    pNext(a) = NULL;
288    return;
289  }
290  else
291  {
292    do
293    {
294      b = pNext(a) = pNew();
295      spMemcpy(b,p);
296      spMonAdd(b,m);
297      if (pComp(b, spNoether) == -1)
298      {
299        pFree1(b);
300        pNext(a) = NULL;
301        return;
302      }
303      a = b;
304      pSetCoeff0(a,nMult(pGetCoeff(p),exp));
305      pIter(p);
306    } while (p!=NULL);
307    pNext(a) = NULL;
308  }
309}
310
311/*2
312* assume m = L(m)
313* pNext(m) = result = a2-a1*m
314* do not destroy a1, but a2
315*/
316#if 0
317void spGSpolyLoop(poly a1, poly a2, poly m,poly spNoether)
318{
319  poly a, b, s;
320  number tm = pGetCoeff(m);
321  number tneg,tb,t;
322  int c;
323
324  tneg = nCopy(tm);
325  tneg = nNeg(tneg);
326  if (a2==NULL)
327  {
328    spGMultCopyX(a1, m, m, tneg,spNoether);
329    nDelete(&tneg);
330    return;
331  }
332  a = m;
333  b = pNew();
334  spMemcpy(b,a1);
335  spMonAdd(b,m);
336  loop
337  {
338    c = pComp0(b, a2);
339    if (c == 1)
340    {
341      pSetCoeff0(b,nMult(pGetCoeff(a1),tneg));
342      a = pNext(a) = b;
343      pIter(a1);
344      if (a1!=NULL)
345      {
346        b = pNew();
347        spMemcpy(b,a1);
348        spMonAdd(b,m);
349      }
350      else
351      {
352        pNext(a) = a2;
353        nDelete(&tneg);
354        return;
355      }
356    }
357    else if (c == -1)
358    {
359      a = pNext(a) = a2;
360      pIter(a2);
361      if (a2==NULL)
362      {
363        pFree1(b);
364        spGMultCopyX(a1, m, a, tneg,spNoether);
365        nDelete(&tneg);
366        return;
367      }
368    }
369    else
370    {
371      t = pGetCoeff(a2);
372      tb = nMult(pGetCoeff(a1),tm);
373      if (!nEqual(t,tb))
374      {
375        t = nSub(t,tb);
376        pSetCoeff(a2,t);
377        a = pNext(a) = a2;
378        pIter(a2);
379      }
380      else
381      {
382        s = a2;
383        pIter(a2);
384        nDelete(&t);
385        pFree1(s);
386      }
387      nDelete(&tb);
388      pIter(a1);
389      if (a1==NULL)
390      {
391        pFree1(b);
392        pNext(a) = a2;
393        nDelete(&tneg);
394        return;
395      }
396      if (a2==NULL)
397      {
398        pFree1(b);
399        spGMultCopyX(a1, m, a, tneg,spNoether);
400        nDelete(&tneg);
401        return;
402      }
403      spMemcpy(b,a1);
404      spMonAdd(b,m);
405    }
406  }
407}
408#endif
409
410/*2
411* reduction of p2 with p1
412* do not destroy p1, but p2
413*/
414poly spGSpolyRed(poly p1, poly p2,poly spNoether, spSpolyLoopProc spSpolyLoop)
415{
416  poly a1 = pNext(p1), a2 = pNext(p2);
417  number an = pGetCoeff(p1), bn = pGetCoeff(p2);
418  int ct = spCheckCoeff(&an, &bn);
419
420  pSetCoeff(p2,bn);
421  if(a1==NULL)
422  {
423    nDelete(&an);
424    nDelete(&bn);
425    pFree1(p2);
426    return a2;
427  }
428  if ((ct == 0) || (ct == 2))
429  {
430    pMultN(a2, an);
431  }
432  nDelete(&an);
433  BOOLEAN reset_vec=FALSE;
434  if (pGetComp(p1) != pGetComp(p2))
435  {
436    pSetCompP(a1,pGetComp(p2));
437    reset_vec=TRUE;
438  }
439  spMonSub(p2,p1);
440  if (ct < 2)
441  {
442    if (spSpolyLoop != NULL)
443      spSpolyLoop(a1, a2, p2,spNoether);
444    else
445      spSpolyLoop_General(a1, a2, p2,spNoether);
446  }
447  else
448  {
449    spGSpolyLoop1(a1, a2, p2,spNoether);
450  }
451  a2 = pNext(p2);
452  if (reset_vec)
453  {
454    spModuleToPoly(a1);
455  }
456  nDelete(&bn);
457  pFree1(p2);
458  return a2;
459}
460
461/*2
462* reduction of tail(q) with p1
463* lead(p1) divides lead(pNext(q2)) and pNext(q2) is reduced
464* do not destroy p1, but tail(q)
465*/
466void spGSpolyTail(poly p1, poly q, poly q2, poly spNoether,
467                  spSpolyLoopProc spSpolyLoop)
468{
469  number t;
470  poly m, h;
471  poly a1 = pNext(p1), p2 = pNext(q2), a2 = pNext(p2);
472  number an = pGetCoeff(p1), bn = pGetCoeff(p2);
473  int ct = spCheckCoeff(&an, &bn);
474  BOOLEAN reset_vec=FALSE;
475
476  if(a1==NULL)
477  {
478    nDelete(&an);
479    nDelete(&bn);
480    nDelete(&pGetCoeff(p2));
481    pFree1(p2);
482    pNext(q2) = a2;
483    return;
484  }
485  if (p1 != q)
486  {
487    m = p2;
488    pSetCoeff(m,bn);
489  }
490  else
491  {
492    m = pNew();
493    spMemcpy(m,p2);
494    pSetCoeff0(m,bn);
495    a2 = pCopy(a2);
496  }
497  if ((ct == 0) || (ct == 2))
498  {
499    pMultN(a2, an);
500  }
501  if ((pGetComp(p1) != pGetComp(p2))
502  && (pGetComp(p1)==0))
503  {
504    pSetCompP(a1,pGetComp(p2));
505    reset_vec=TRUE;
506  }
507  spMonSub(m,p1);
508  if (ct < 2)
509  {
510    if (spSpolyLoop != NULL)
511      spSpolyLoop(a1, a2, m,spNoether);
512    else
513      spSpolyLoop_General(a1, a2, m,spNoether);
514  }
515  else
516  {
517    spGSpolyLoop1(a1, a2, m,spNoether);
518  }
519  if ((ct == 0) || (ct == 2))
520  {
521    h = q;
522    do
523    {
524      t = nMult(pGetCoeff(h),an);
525      pSetCoeff(h,t);
526      pIter(h);
527    }
528    while (h != p2);
529  }
530  h = pNext(m);
531  nDelete(&an);
532  nDelete(&bn);
533  pFree1(m);
534  pNext(q2) = h;
535  if (reset_vec)
536    spModuleToPoly(a1);
537  if (p1 == q)
538  {
539    pDelete(&p2);
540  }
541}
542
543/*2
544* reduction of p2 with p1
545* do not destroy p1 and p2
546*/
547poly spGSpolyRedNew(poly p1, poly p2,poly spNoether, 
548                    spSpolyLoopProc spSpolyLoop)
549{
550  poly m;
551  poly a1 = pNext(p1), a2 = pNext(p2);
552  number an = pGetCoeff(p1), bn = pGetCoeff(p2);
553  int ct = spCheckCoeff(&an, &bn);
554
555  if(a1==NULL)
556  {
557    nDelete(&an);
558    nDelete(&bn);
559    return pCopy(a2);
560  }
561  if ((ct == 0) || (ct == 2))
562  {
563    a2 = pMultCopyN(a2,an);
564  }
565  else
566  {
567    a2 = pCopy(a2);
568  }
569  pTest(a2);
570  pTest(p2);
571  nDelete(&an);
572  BOOLEAN reset_vec=FALSE;
573  if (pGetComp(p1) != pGetComp(p2))
574  {
575    pSetCompP(a1,pGetComp(p2));
576    reset_vec=TRUE;
577  }
578  m = pNew();
579  spMemcpy(m,p2);
580  spMonSub(m,p1);
581  if (ct < 2)
582  {
583    pSetCoeff0(m,bn);
584    if (spSpolyLoop != NULL)
585      spSpolyLoop(a1, a2, m,spNoether);
586    else
587      spSpolyLoop_General(a1, a2, m,spNoether);
588  }
589  else
590  {
591    spGSpolyLoop1(a1, a2, m,spNoether);
592  }
593  a2 = pNext(m);
594  pTest(a2);
595  pTest(p2);
596  if (reset_vec)
597  {
598    spModuleToPoly(a1);
599  }
600  nDelete(&bn);
601  pFree1(m);
602  pTest(a2);
603  pTest(p2);
604  return a2;
605}
606
607/*2
608* creates the S-polynomial of p1 and p2
609* do not destroy p1 and p2
610*/
611poly spGSpolyCreate(poly p1, poly p2,poly spNoether, 
612                    spSpolyLoopProc spSpolyLoop)
613{
614  Exponent_t x;
615  poly m, b;
616  poly a1 = pNext(p1), a2 = pNext(p2);
617  number an = pGetCoeff(p1), bn = pGetCoeff(p2);
618  int co=0, ct = spCheckCoeff(&an, &bn);
619
620  if (pGetComp(p1)!=pGetComp(p2))
621  {
622    if (pGetComp(p1)==0)
623    {
624      co=1;
625      pSetCompP(p1,pGetComp(p2));
626    }
627    else
628    {
629      co=2;
630      pSetCompP(p2,pGetComp(p1));
631    }
632  }
633  b = pInit();
634  m = pInit();
635  for (int i = pVariables; i; i--)
636  {
637    x = pGetExp(p1,i) - pGetExp(p2,i);
638    if (x > 0)
639    {
640      pSetExp(b,i,x);
641      pSetExp(m,i,0);
642    }
643    else
644    {
645      pSetExp(m,i,-x);
646      pSetExp(b,i,0);
647    }
648  }
649  pSetm(m);
650  pSetm(b);
651  if (a2!=NULL)
652  {
653    if ((ct == 0) || (ct == 2))
654    {
655      spGMultCopyX(a2, b, b, an,spNoether);
656    }
657    else
658    {
659      spGMultCopy0(a2, b,spNoether);
660    }
661    a2 = pNext(b);
662  }
663  nDelete(&an);
664  pFree1(b);
665  if (a1!=NULL)
666  {
667    if (ct < 2)
668    {
669      pSetCoeff0(m,bn);
670    if (spSpolyLoop != NULL)
671      spSpolyLoop(a1, a2, m,spNoether);
672    else
673      spSpolyLoop_General(a1, a2, m,spNoether);
674    }
675    else
676    {
677      spGSpolyLoop1(a1, a2, m,spNoether);
678    }
679    a2 = pNext(m);
680  }
681  nDelete(&bn);
682  pFree1(m);
683  if (co != 0)
684  {
685    if (co==1)
686    {
687      spModuleToPoly(p1);
688    }
689    else
690    {
691      spModuleToPoly(p2);
692    }
693  }
694  return a2;
695}
696
697/*2
698* creates the leading term of the S-polynomial of p1 and p2
699* do not destroy p1 and p2
700* remarks:
701*   1. the coefficient is 0 (nNew)
702*   2. pNext is undefined
703*/
704static void bbb() { int i=0; }
705poly spGSpolyShortBba(poly p1, poly p2)
706{
707  poly a1 = pNext(p1), a2 = pNext(p2);
708  Exponent_t c1=pGetComp(p1),c2=pGetComp(p2);
709  Exponent_t c;
710  poly m1,m2;
711  number t1,t2;
712  int cm,i;
713  BOOLEAN equal;
714
715  if (a1==NULL)
716  {
717    if(a2!=NULL)
718    {
719      m2=pInit();
720x2:
721      for (i = pVariables; i; i--)
722      {
723        c = pGetExpDiff(p1, p2,i);
724        if (c>0)
725        {
726          pSetExp(m2,i,(c+pGetExp(a2,i)));
727        }
728        else
729        {
730          pSetExp(m2,i,pGetExp(a2,i));
731        }
732      }
733      if ((c1==c2)||(c2!=0))
734      {
735        pSetComp(m2,pGetComp(a2));
736      }
737      else
738      {
739        pSetComp(m2,c1);
740      }
741      pSetm(m2);
742      nNew(&(pGetCoeff(m2)));
743      return m2;
744    }
745    else
746      return NULL;
747  }
748  if (a2==NULL)
749  {
750    m1=pInit();
751x1:
752    for (i = pVariables; i; i--)
753    {
754      c = pGetExpDiff(p2, p1,i);
755      if (c>0)
756      {
757        pSetExp(m1,i,(c+pGetExp(a1,i)));
758      }
759      else
760      {
761        pSetExp(m1,i,pGetExp(a1,i));
762      }
763    }
764    if ((c1==c2)||(c1!=0))
765    {
766      pSetComp(m1,pGetComp(a1));
767    }
768    else
769    {
770      pSetComp(m1,c2);
771    }
772    pSetm(m1);
773    nNew(&(pGetCoeff(m1)));
774    return m1;
775  }
776  m1 = pInit();
777  m2 = pInit();
778  loop
779  {
780    for (i = pVariables; i; i--)
781    {
782      c = pGetExpDiff(p1, p2,i);
783      if (c > 0)
784      {
785        pSetExp(m2,i,(c+pGetExp(a2,i)));
786        pSetExp(m1,i,pGetExp(a1,i));
787      }
788      else
789      {
790        pSetExp(m1,i,(pGetExp(a1,i)-c));
791        pSetExp(m2,i,pGetExp(a2,i));
792      }
793    }
794    if(c1==c2)
795    {
796      pSetComp(m1,pGetComp(a1));
797      pSetComp(m2,pGetComp(a2));
798    }
799    else
800    {
801      if(c1!=0)
802      {
803        pSetComp(m1,pGetComp(a1));
804        pSetComp(m2,c1);
805      }
806      else
807      {
808        pSetComp(m2,pGetComp(a2));
809        pSetComp(m1,c2);
810      }
811    }
812    pSetm(m1);
813    pSetm(m2);
814    cm = pComp0(m1, m2);
815    if (cm!=0)
816    {
817      if(cm==1)
818      {
819        pFree1(m2);
820        nNew(&(pGetCoeff(m1)));
821        return m1;
822      }
823      else
824      {
825        pFree1(m1);
826        nNew(&(pGetCoeff(m2)));
827        return m2;
828      }
829    }
830    t1 = nMult(pGetCoeff(a2),pGetCoeff(p1));
831    t2 = nMult(pGetCoeff(a1),pGetCoeff(p2));
832    equal = nEqual(t1,t2);
833    nDelete(&t2);
834    nDelete(&t1);
835    if (!equal)
836    {
837      pFree1(m2);
838      nNew(&(pGetCoeff(m1)));
839      return m1;
840    }
841    pIter(a1);
842    pIter(a2);
843    if (a2==NULL)
844    {
845      pFree1(m2);
846      if (a1==NULL)
847      {
848        pFree1(m1);
849        return NULL;
850      }
851      goto x1;
852    }
853    if (a1==NULL)
854    {
855      pFree1(m1);
856      goto x2;
857    }
858  }
859}
860
Note: See TracBrowser for help on using the repository browser.