source: git/Singular/spolys.cc @ 6f2edc

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