source: git/Singular/spolys.cc @ b95c83

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