source: git/Singular/spolys.cc @ d18870

fieker-DuValspielwiese
Last change on this file since d18870 was a9a7be, checked in by Olaf Bachmann <obachman@…>, 25 years ago
* walk stuff git-svn-id: file:///usr/local/Singular/svn/trunk@3682 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 15.2 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: spolys.cc,v 1.22 1999-09-27 15:05:32 obachman 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  pTest(p1);
263  pTest(p2);
264  poly a1 = pNext(p1), a2 = pNext(p2);
265  if(a1==NULL)
266  {
267    pFree1(p2);
268    return a2;
269  }
270  BOOLEAN reset_vec=FALSE;
271  if (pGetComp(p1) != pGetComp(p2))
272  {
273    pSetCompP(a1,pGetComp(p2));
274    reset_vec=TRUE;
275  }
276  pTest(p1);
277  pTest(p2);
278  spMonSub(p2,p1);
279  if (1!=(int)pGetCoeff(p2))
280  {
281    if (SpolyLoop != NULL)
282      SpolyLoop(a1, a2, p2,spNoether);
283    else
284      spSpolyLoop_General(a1, a2, p2,spNoether);
285  }
286  else
287    spSpolyLoop1(a1, a2, p2,spNoether);
288  a2 = pNext(p2);
289  if (reset_vec)
290    spModuleToPoly(a1);
291  pFree1(p2);
292  return a2;
293}
294
295/*
296* reduction of tail(q) with p1
297* lead(p1) divides lead(pNext(q2)) and pNext(q2) is reduced
298* do not destroy p1, but tail(q)
299*/
300static void spPSpolyTail(poly p1, poly q, poly q2, poly spNoether,
301                         spSpolyLoopProc spSpolyLoop)
302{
303  poly h = pNext(q2);
304  if (p1 != q)
305  {
306    pNext(q2) = spSpolyRed(p1, h, spNoether, spSpolyLoop);
307  }
308  else
309  {
310    pNext(q2) = spSpolyRedNew(p1, h, spNoether, spSpolyLoop);
311    pDelete(&h);
312  }
313}
314
315/*2
316* reduction of p2 with p1
317* do not destroy p1 and p2
318*/
319static poly spPSpolyRedNew(poly p1, poly p2,poly spNoether,
320                           spSpolyLoopProc SpolyLoop)
321{
322  poly a1 = pNext(p1), a2 = pNext(p2);
323  poly m;
324  if (a2!=NULL)
325    a2 = pCopy(a2);
326  if (a1==NULL)
327    return a2;
328  BOOLEAN reset_vec=FALSE;
329  if (pGetComp(p1) != pGetComp(p2))
330  {
331    pSetCompP(a1,pGetComp(p2));
332    reset_vec=TRUE;
333  }
334  m = pNew();
335  spMemcpy(m,p2);
336  spMonSub(m,p1);
337  if (1!=(int)pGetCoeff(m))
338  {
339    if (SpolyLoop != NULL)
340      SpolyLoop(a1, a2, m,spNoether);
341    else
342      spSpolyLoop_General(a1, a2, m,spNoether);
343  }
344  else
345    spSpolyLoop1(a1, a2, m,spNoether);
346  a2 = pNext(m);
347  if (reset_vec)
348    spModuleToPoly(a1);
349  pFree1(m);
350  return a2;
351}
352
353/*2
354* creates the S-polynomial of p1 and p2
355* do not destroy p1 and p2
356*/
357static poly spPSpolyCreate(poly p1, poly p2,poly spNoether,
358                           spSpolyLoopProc SpolyLoop)
359{
360  poly a1 = pNext(p1), a2 = pNext(p2);
361  poly m, b;
362  int co=0;
363  Exponent_t c;
364  if (pGetComp(p1)!=pGetComp(p2))
365  {
366    if (pGetComp(p1)==0)
367    {
368      co=1;
369      pSetCompP(p1,pGetComp(p2));
370    }
371    else
372    {
373      co=2;
374      pSetCompP(p2,pGetComp(p1));
375    }
376  }
377  b = pInit();
378  m = pInit();
379  for (int i = pVariables; i; i--)
380  {
381    c = pGetExpDiff(p1, p2,i);
382    if (c > 0)
383    {
384      pSetExp(b,i,c);
385      pSetExp(m,i,0);
386    }
387    else
388    {
389      pSetExp(m,i,-c);
390      pSetExp(b,i,0);
391    }
392  }
393  pSetm(m);
394  pSetm(b);
395  if (a2!=NULL)
396  {
397    spMultCopy0(a2, b,spNoether);
398    a2 = pNext(b);
399  }
400  pFree1(b);
401  if (a1!=NULL)
402  {
403    spSpolyLoop1(a1, a2, m,spNoether);
404    a2 = pNext(m);
405  }
406  pFree1(m);
407  if (co==1) spModuleToPoly(p1);
408  else if (co==2) spModuleToPoly(p2);
409  return a2;
410}
411
412/*2
413* creates the leading term of the S-polynomial of p1 and p2
414* do not destroy p1 and p2
415* remarks:
416*   1. the coefficient is undefined
417*   2. pNext is undefined
418*/
419//#ifndef TEST_MAC_ORDER
420#if 1
421static poly spPSpolyShortBba(poly p1, poly p2)
422{
423  poly a1 = pNext(p1), a2 = pNext(p2);
424  Exponent_t c1=pGetComp(p1),c2=pGetComp(p2);
425  Exponent_t c;
426  poly m1,m2;
427  int cm,i;
428
429  if (a1==NULL)
430  {
431    if(a2!=NULL)
432    {
433      m2=pInit();
434x2:
435      for (i = pVariables; i; i--)
436      {
437        c = pGetExpDiff(p1, p2,i);
438        if (c>0)
439        {
440          pSetExp(m2,i,(c+pGetExp(a2,i)));
441        }
442        else
443        {
444          pSetExp(m2,i,pGetExp(a2,i));
445        }
446      }
447      if ((c1==c2)||(c2!=0))
448      {
449        pSetComp(m2,pGetComp(a2));
450      }
451      else
452      {
453        pSetComp(m2,c1);
454      }
455      pSetm(m2);
456      return m2;
457    }
458    else
459      return NULL;
460  }
461  else if (a2==NULL)
462  {
463    m1=pInit();
464x1:
465    for (i = pVariables; i; i--)
466    {
467      c = pGetExpDiff(p2, p1,i);
468      if (c>0)
469      {
470        pSetExp(m1,i,(c+pGetExp(a1,i)));
471      }
472      else
473      {
474        pSetExp(m1,i,pGetExp(a1,i));
475      }
476    }
477    if ((c1==c2)||(c1!=0))
478    {
479      pSetComp(m1,pGetComp(a1));
480    }
481    else
482    {
483      pSetComp(m1,c2);
484    }
485    pSetm(m1);
486    return m1;
487  }
488  m1 = pInit();
489  m2 = pInit();
490  loop
491  {
492    for (i = pVariables; i; i--)
493    {
494      c = pGetExpDiff(p1, p2,i);
495      if (c > 0)
496      {
497        pSetExp(m2,i,(c+pGetExp(a2,i)));
498        pSetExp(m1,i,pGetExp(a1,i));
499      }
500      else
501      {
502        pSetExp(m1,i,(pGetExp(a1,i)-c));
503        pSetExp(m2,i,pGetExp(a2,i));
504      }
505    }
506    if(c1==c2)
507    {
508      pSetComp(m1,pGetComp(a1));
509      pSetComp(m2,pGetComp(a2));
510    }
511    else
512    {
513      if(c1!=0)
514      {
515        pSetComp(m1,pGetComp(a1));
516        pSetComp(m2,c1);
517      }
518      else
519      {
520        pSetComp(m2,pGetComp(a2));
521        pSetComp(m1,c2);
522      }
523    }
524    pSetm(m1);
525    pSetm(m2);
526    cm = pComp0(m1, m2);
527    if (cm!=0)
528    {
529      if(cm==1)
530      {
531        pFree1(m2);
532        return m1;
533      }
534      else
535      {
536        pFree1(m1);
537        return m2;
538      }
539    }
540    if (!npEqualM(pGetCoeff(a2),pGetCoeff(a1)))
541    {
542      pFree1(m2);
543      return m1;
544    }
545    pIter(a1);
546    pIter(a2);
547    if (a2==NULL)
548    {
549      pFree1(m2);
550      if (a1==NULL)
551      {
552        pFree1(m1);
553        return NULL;
554      }
555      goto x1;
556    }
557    if (a1==NULL)
558    {
559      pFree1(m1);
560      goto x2;
561    }
562  }
563}
564#else
565static poly spPSpolyShortBba(poly p1, poly p2)
566{
567  poly a1 = pNext(p1), a2 = pNext(p2);
568  int co=0;
569  // adjust components, store change poly-nr in co
570  if (pGetComp(p1)!=pGetComp(p2))
571  {
572    if (pGetComp(p1)==0)
573    {
574      co=1;
575      pSetCompP(p1,pGetComp(p2));
576    }
577    else
578    {
579      co=2;
580      pSetCompP(p2,pGetComp(p1));
581    }
582  }
583  //------------------------------------------------------------
584  poly b = pInit(); /* the result */
585  poly lcm=pInit();
586  pLcm(p1,p2,lcm);
587  pSetComp(b,pGetComp(p1));
588  loop
589  {
590    if(a1==NULL)
591    {
592      if(a2==NULL)
593      {
594        pFree1(b);
595        b=NULL;
596        break;
597      }
598      else
599      {
600        spMemcpy(b,a2);
601        MyspMonAdd0(b,lcm);
602        spMonSub(b,p2);
603        break;
604      }
605    }
606    else /*a1!=NULL*/
607    {
608      if(a2==NULL)
609      {
610        spMemcpy(b,a1);
611        MyspMonAdd0(b,lcm);
612        spMonSub(b,p1);
613        break;
614      }
615    }
616    /* now a1!=NULL and a2!=NULL */
617    {
618      poly n1=pNew();
619      poly n2=pNew();
620      spMemcpy(n1,a1);
621      spMemcpy(n2,a2);
622      spMonAdd(n1,p2);
623      spMonAdd(n2,p1);
624      int c=pComp0(n1,n2);
625      pFree1(n2);
626      pFree1(n1);
627      if(c==1)
628      {
629        spMemcpy(b,a1);
630        MyspMonAdd0(b,lcm);
631        spMonSub(b,p1);
632        break;
633      }
634      else if(c==-1)
635      {
636        spMemcpy(b,a2);
637        MyspMonAdd0(b,lcm);
638        spMonSub(b,p2);
639        break;
640      }
641      else
642      {
643        if (!npEqualM(pGetCoeff(a2),pGetCoeff(a1)))
644        {
645          spMemcpy(b,a1);
646          MyspMonAdd0(b,lcm);
647          spMonSub(b,p1);
648          break;
649        }
650        pIter(a1);
651        pIter(a2);
652      }
653    }
654  }
655  if (co==1) spModuleToPoly(p1);
656  else if (co==2) spModuleToPoly(p2);
657  pFree1(lcm);
658  return b;
659}
660#endif
661
662#ifdef SDRING
663poly spDSpolyRed(poly p1,poly p2,poly spNoether, spSpolyLoopProc dummy)
664{
665  poly m=pOne();
666  poly tp1;
667  number n1,n2;
668
669  int i;
670
671  pTest(p1);
672  for (i=1;i<=pVariables;i++)
673  {
674    pSetExp(m,i,pGetExp(p2,i)-pGetExp(p1,i));
675  }
676#ifdef DRING
677  if (pDRING)
678  {
679    if (pdDFlag(p1)==0) /* && (pdDFlag(p2==0) */
680    {
681      for (i=pdN+1;i<=2*pdN;i++)
682      {
683        if (pGetExp(m,i)!=0)
684        {
685          pDelete1(&m);
686          pTest(p2);
687          return p2;
688        }
689      }
690    }
691    pdSetDFlag(m,1);
692  }
693#endif
694  pSetm(m);
695  pTest(p1);
696  tp1=pMult(m,pCopy(p1));
697  pTest(tp1);
698  if (tp1!=NULL)
699  {
700    n2=nCopy(pGetCoeff(p2));
701    n2=nNeg(n2);
702    n1=pGetCoeff(tp1);
703    pMultN(tp1,n2);
704    pMultN(p2,n1);
705    nDelete(&n2);
706    //pDelete1(&p2);
707    //pDelete1(&tp1);
708    p2=pAdd(tp1,p2);
709  }
710  pTest(p2);
711  return p2;
712}
713
714poly spDSpolyRedNew(poly p1,poly p2,poly spNoether, spSpolyLoopProc dummy)
715{
716  return spDSpolyRed(p1,pCopy(p2),spNoether, dummy);
717}
718
719poly spDSpolyCreate(poly p1,poly p2,poly spNoether, spSpolyLoopProc dummy)
720{
721  poly m1=pOne();
722  poly m2=pOne();
723  poly pp1=p1,pp2=p2;
724  number n1,n2;
725  int i,j;
726  int co=0;
727  if (pGetComp(p1)!=pGetComp(p2))
728  {
729    if (pGetComp(p1)==0)
730    {
731      co=1;
732      pSetCompP(p1,pGetComp(p2));
733    }
734    else
735    {
736      co=2;
737      pSetCompP(p2,pGetComp(p1));
738    }
739  }
740
741  for (i=1;i<=pVariables;i++)
742  {
743    j=pGetExp(p2,i)-pGetExp(p1,i);
744    if (j>0)
745    {
746      pSetExp(m1,i,j);
747    }
748    else
749    {
750      pSetExp(m2,i,-j);
751    }
752  }
753#ifdef DRING
754  if (pDRING)
755  {
756    if (pdDFlag(p1)==0) /* && (pdDFlag(p2==0) */
757    {
758      for (i=pdN+1;i<=2*pdN;i++)
759      {
760        if ((pGetExp(m1,i)!=0) || (pGetExp(m2,i)!=0))
761        {
762         pDelete1(&m1);
763         pDelete1(&m2);
764         return NULL;
765        }
766      }
767    }
768    pdSetDFlag(m1,1);
769    pdSetDFlag(m2,1);
770  }
771#endif
772  pSetm(m1);
773  pSetm(m2);
774  p1=pMult(m1,pCopy(p1));
775  p2=pMult(m2,pCopy(p2));
776
777  n1=nCopy(pGetCoeff(p1));
778  n2=nCopy(pGetCoeff(p2));
779
780  pDelete1(&p1);
781  pDelete1(&p2);
782
783  n1=nNeg(n1);
784  pMultN(p1,n2);
785  nDelete(&n2);
786  pMultN(p2,n1);
787  nDelete(&n1);
788
789  m1=pAdd(p1,p2);
790  if (co==1) spModuleToPoly(pp1);
791  else if (co==2) spModuleToPoly(pp2);
792  pTest(m1);
793  return m1;
794}
795
796#endif
797
798void spSet(ring r)
799{
800  if ((rField_is_Q())
801  || (rField_is_Extension()) /* Q(a), Fp(a) */
802  || (rField_is_numeric()) /* R, long R, long C*/
803#ifdef SRING
804  && (pSRING==0)
805#endif
806#ifdef DRING
807  && (pDRING==0)
808#endif
809  )
810  {
811    spSpolyRed=spGSpolyRed;
812    spSpolyTail=spGSpolyTail;
813    spSpolyRedNew=spGSpolyRedNew;
814    spSpolyCreate=spGSpolyCreate;
815    spSpolyShortBba=spGSpolyShortBba;
816    return;
817  }
818  if (rField_is_GF(r) /* finite fields GF(p,n) */
819#ifdef SRING
820  && (pSRING==0)
821#endif
822#ifdef DRING
823  && (pDRING==0)
824#endif
825  )
826  {
827    spSpolyRed=spGSpolyRed;
828    spSpolyTail=spGSpolyTail;
829    spSpolyRedNew=spGSpolyRedNew;
830    spSpolyCreate=spGSpolyCreate;
831    spSpolyShortBba=spGSpolyShortBba;
832    return;
833  }
834  if (rField_is_Zp(r) /* Fp */
835#ifdef SRING
836  && (pSRING==0)
837#endif
838#ifdef DRING
839  && (pDRING==0)
840#endif
841  )
842  {
843    spSpolyRed=spPSpolyRed;
844    spSpolyTail=spPSpolyTail;
845    spSpolyRedNew=spPSpolyRedNew;
846    spSpolyCreate=spPSpolyCreate;
847    spSpolyShortBba=spPSpolyShortBba;
848    return;
849  }
850#ifdef SRING
851  if (pSRING!=0)
852  {
853    spSpolyRed=spDSpolyRed;
854    spSpolyRedNew=spDSpolyRedNew;
855    spSpolyCreate=spDSpolyCreate;
856    spSpolyShortBba=NULL;/*spGSpolyShortBba;*/
857    test&=~Sy_bit(OPT_REDTAIL); // no redtail
858    return;
859  }
860#endif
861#ifdef DRING
862  if (pDRING!=0)
863  {
864    spSpolyRed=spDSpolyRed;
865    spSpolyRedNew=spDSpolyRedNew;
866    spSpolyCreate=spDSpolyCreate;
867    spSpolyShortBba=NULL;/*spGSpolyShortBba;*/
868    test&=~Sy_bit(OPT_REDTAIL); // no redtail
869    return;
870  }
871#endif
872  WarnS("spoly init err");
873}
Note: See TracBrowser for help on using the repository browser.