source: git/Singular/polys1.cc @ d2b2a7

spielwiese
Last change on this file since d2b2a7 was 49eeee, checked in by Hans Schönemann <hannes@…>, 26 years ago
* hannes: fixed bug in maps.cc: size of monom is not pMonomSize but mmGetSpecSize() changed header in tesths.cc changed handling of ring changes (polys*.cc ring.cc) git-svn-id: file:///usr/local/Singular/svn/trunk@1462 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 21.9 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: polys1.cc,v 1.13 1998-04-24 16:39:25 Singular Exp $ */
5
6/*
7* ABSTRACT - all basic methods to manipulate polynomials:
8* independent of representation
9*/
10
11/* includes */
12#include <string.h>
13#include "mod2.h"
14#include "structs.h"
15#include "tok.h"
16#include "numbers.h"
17#include "mmemory.h"
18#include "febase.h"
19#include "weight.h"
20#include "intvec.h"
21#include "longalg.h"
22#include "ring.h"
23#include "ideals.h"
24#include "polys.h"
25#include "ipid.h"
26#ifdef HAVE_FACTORY
27#include "clapsing.h"
28#endif
29
30/*-------- several access procedures to monomials -------------------- */
31/*2
32*test if the monomial is a constant
33*/
34BOOLEAN   pIsConstant(poly p)
35{
36  if (p!=NULL)
37  {
38    int i;
39    for (i=pVariables;i;i--)
40    {
41      if (pGetExp(p,i)!=0) return FALSE;
42    }
43    if (pGetComp(p)) return FALSE;
44  }
45  return TRUE;
46}
47
48/*2
49*test if the monomial is a constant as a vector component
50*/
51BOOLEAN   pIsConstantComp(poly p)
52{
53  int i;
54
55  for (i=pVariables;i>0;i--)
56  {
57    if (pGetExp(p,i)!=0) return FALSE;
58  }
59  return TRUE;
60}
61
62/*2
63*test if a monomial is a pure power
64*/
65int pIsPurePower(poly p)
66{
67  pTest(p);
68  int i,k=0;
69
70  for (i=pVariables;i>0;i--)
71  {
72    if (pGetExp(p,i)!=0)
73    {
74      if(k!=0) return 0;
75      k=i;
76    }
77  }
78  return k;
79}
80
81/*-----------------------------------------------------------*/
82/*
83* the module weights for std
84*/
85static pFDegProc pOldFDeg;
86static intvec * pModW;
87static BOOLEAN pOldLexOrder;
88
89static int pModDeg(poly p)
90{
91  return pOldFDeg(p)+(*pModW)[pGetComp(p)-1];
92}
93
94void pSetModDeg(intvec *w)
95{
96  if (w!=NULL)
97  {
98    pModW = w;
99    pOldFDeg = pFDeg;
100    pOldLexOrder = pLexOrder;
101    pFDeg = pModDeg;
102    pLexOrder = TRUE;
103  }
104  else
105  {
106    pModW = NULL;
107    pFDeg = pOldFDeg;
108    pLexOrder = pOldLexOrder;
109  }
110}
111
112/*-------------operations on polynomials:------------*/
113/*2
114* add p1 and p2, p1 and p2 are destroyed
115*/
116poly pAdd(poly p1, poly p2)
117{
118  static poly a1, p, a2, a;
119  int c;
120  number t;
121
122  if (p1==NULL) return p2;
123  if (p2==NULL) return p1;
124  a1 = p1;
125  a2 = p2;
126  a = p  = pNew();
127  nNew(&(p->coef));
128  loop
129  {
130    /* a1 and a2 are non-NULL, so we may use pComp0 instead of pComp */
131    c = pComp0(a1, a2);
132    if (c == 1)
133    {
134      a = pNext(a) = a1;
135      pIter(a1);
136      if (a1==NULL)
137      {
138        pNext(a) = a2;
139        break;
140      }
141    }
142    else if (c == -1)
143    {
144      a = pNext(a) = a2;
145      pIter(a2);
146      if (a2==NULL)
147      {
148        pNext(a) = a1;
149        break;
150      }
151    }
152    else
153    {
154      t=nAdd(a1->coef,a2->coef);
155      pDelete1(&a2);
156      if (nIsZero(t))
157      {
158        nDelete(&t);
159        pDelete1(&a1);
160      }
161      else
162      {
163        pSetCoeff(a1,t);
164        a = pNext(a) = a1;
165        pIter(a1);
166      }
167      if (a1==NULL)
168      {
169        pNext(a) = a2;
170        break;
171      }
172      else if (a2==NULL)
173      {
174        pNext(a) = a1;
175        break;
176      }
177    }
178  }
179  pDelete1(&p);
180  return p;
181}
182
183/*2
184* multiply a polynomial by -1
185*/
186poly pNeg(poly p)
187{
188  poly h = p;
189  while (h!=NULL)
190  {
191    pGetCoeff(h)=nNeg(pGetCoeff(h));
192    pIter(h);
193  }
194  return p;
195}
196
197/*2
198* subtract p2 from p1, p1 and p2 are destroyed
199* do not put attention on speed: the procedure is only used in the interpreter
200*/
201poly pSub(poly p1, poly p2)
202{
203  return pAdd(p1, pNeg(p2));
204}
205
206/*3
207*  create binomial coef.
208*/
209static number* pnBin(int exp)
210{
211  int e, i, h;
212  number x, y, *bin=NULL;
213
214  x = nInit(exp);
215  if (nIsZero(x))
216  {
217    nDelete(&x);
218    return bin;
219  }
220  h = (exp >> 1) + 1;
221  bin = (number *)Alloc0(h*sizeof(number));
222  bin[1] = x;
223  if (exp < 4)
224    return bin;
225  i = exp - 1;
226  for (e=2; e<h; e++)
227  {
228//    if(!nIsZero(bin[e-1]))
229//    {
230      x = nInit(i);
231      i--;
232      y = nMult(x,bin[e-1]);
233      nDelete(&x);
234      x = nInit(e);
235      bin[e] = nIntDiv(y,x);
236      nDelete(&x);
237      nDelete(&y);
238//    }
239//    else
240//    {
241//      bin[e] = nInit(binom(exp,e));
242//    }
243  }
244  return bin;
245}
246
247static void pnFreeBin(number *bin, int exp)
248{
249  int e, h = (exp >> 1) + 1;
250
251  if (bin[1] != NULL)
252  {
253    for (e=1; e<h; e++)
254      nDelete(&(bin[e]));
255  }
256  Free((ADDRESS)bin, h*sizeof(number));
257}
258
259/*3
260* compute for a monomial m
261* the power m^exp, exp > 1
262* destroys p
263*/
264static poly pMonPower(poly p, int exp)
265{
266  number x, y;
267  int i;
268
269  y = pGetCoeff(p);
270  nPower(y,exp,&x);
271  nDelete(&y);
272  pSetCoeff0(p,x);
273  for (i=pVariables; i!=0; i--)
274  {
275    pMultExp(p,i, exp);
276  }
277  #ifdef TEST_MAC_ORDER
278  if (bNoAdd)
279    pSetm(p);
280  else
281  #endif
282    p->Order *= exp;
283  return p;
284}
285
286/*3
287* compute for monomials p*q
288* destroys p, keeps q
289*/
290static void pMonMult(poly p, poly q)
291{
292  number x, y;
293  int i;
294
295  y = pGetCoeff(p);
296  x = nMult(y,pGetCoeff(q));
297  nDelete(&y);
298  pSetCoeff0(p,x);
299  for (i=pVariables; i!=0; i--)
300  {
301    pAddExp(p,i, pGetExp(q,i));
302  }
303  p->Order += q->Order;
304}
305
306/*3
307* compute for monomials p*q
308* keeps p, q
309*/
310static poly pMonMultC(poly p, poly q)
311{
312  number x;
313  int i;
314  poly r = pInit();
315
316  x = nMult(pGetCoeff(p),pGetCoeff(q));
317  pSetCoeff0(r,x);
318  pSetComp(r, pGetComp(p));
319  for (i=pVariables; i!=0; i--)
320  {
321    pSetExp(r,i, pGetExp(p,i) + pGetExp(q,i));
322  }
323  r->Order = p->Order + q->Order;
324  return r;
325}
326
327/*
328*  compute for a poly p = head+tail, tail is monomial
329*          (head + tail)^exp, exp > 1
330*          with binomial coef.
331*/
332static poly pTwoMonPower(poly p, int exp)
333{
334  int eh, e, al;
335  poly *a;
336  poly tail, b, res, h;
337  number x;
338  number *bin = pnBin(exp);
339
340  tail = pNext(p);
341  if (bin == NULL)
342  {
343    pMonPower(p,exp);
344    pMonPower(tail,exp);
345#ifdef PDEBUG
346    pTest(p);
347#endif
348    return p;
349  }
350  eh = exp >> 1;
351  al = (exp + 1) * sizeof(poly);
352  a = (poly *)Alloc(al);
353  a[1] = p;
354  for (e=1; e<exp; e++)
355  {
356    a[e+1] = pMonMultC(a[e],p);
357  }
358  res = a[exp];
359  b = pHead(tail);
360  for (e=exp-1; e>eh; e--)
361  {
362    h = a[e];
363    x = nMult(bin[exp-e],pGetCoeff(h));
364    pSetCoeff(h,x);
365    pMonMult(h,b);
366    res = pNext(res) = h;
367    pMonMult(b,tail);
368  }
369  for (e=eh; e!=0; e--)
370  {
371    h = a[e];
372    x = nMult(bin[e],pGetCoeff(h));
373    pSetCoeff(h,x);
374    pMonMult(h,b);
375    res = pNext(res) = h;
376    pMonMult(b,tail);
377  }
378  pDelete1(&tail);
379  pNext(res) = b;
380  pNext(b) = NULL;
381  res = a[exp];
382  Free((ADDRESS)a, al);
383  pnFreeBin(bin, exp);
384//  tail=res;
385// while((tail!=NULL)&&(pNext(tail)!=NULL))
386// {
387//   if(nIsZero(pGetCoeff(pNext(tail))))
388//   {
389//     pDelete1(&pNext(tail));
390//   }
391//   else
392//     pIter(tail);
393// }
394#ifdef PDEBUG
395  pTest(res);
396#endif
397  return res;
398}
399
400static poly pPow(poly p, int i)
401{
402  poly rc = pCopy(p);
403  i -= 2;
404  do
405  {
406    rc = pMult(rc,pCopy(p));
407    i--;
408  }
409  while (i != 0);
410  return pMult(rc,p);
411}
412
413/*2
414* returns the i-th power of p
415* p will be destroyed
416*/
417poly pPower(poly p, int i)
418{
419  if (i==0)
420    return pOne();
421
422  poly rc=NULL;
423
424  if(p!=NULL)
425  {
426    switch (i)
427    {
428      case 0:
429      {
430        rc=pOne();
431#ifdef DRING
432        if ((pDRING) && (pdDFlag(p)==1))
433        {
434          pdDFlag(rc)=1;
435        }
436#endif
437        pDelete(&p);
438        break;
439      }
440      case 1:
441        rc=p;
442        break;
443      case 2:
444        rc=pMult(pCopy(p),p);
445        break;
446      default:
447        if (i < 0)
448        {
449#ifdef DRING
450          if (pDRING)
451          {
452            int j,t;
453            rc=p;
454            while(p!=NULL)
455            {
456              for(j=1;j<=pdK;j++)
457              {
458                if(pGetExp(p,pdY(j))!=0)
459                {
460                  pDelete(&rc);
461                  return NULL;
462                }
463              }
464              for(j=1;j<=pdN;j++)
465              {
466                t=pGetExp(p,pdX(j));
467                pGetExp(p,pdX(j)]=pGetExp(p,pdIX(j));
468                pGetExp(p,pdIX(j))=t;
469              }
470              pIter(p);
471            }
472            return pPower(rc,-i);
473          }
474          else
475#endif
476          {
477            pDelete(&p);
478            return NULL;
479          }
480        }
481#ifdef SDRING
482        if(pSDRING)
483        {
484          return pPow(p,i);
485        }
486        else
487#endif
488        {
489          rc = pNext(p);
490          if (rc != NULL)
491          {
492            int char_p=rChar(currRing);
493            if ((pNext(rc) != NULL)
494            || ((char_p>0) && (i>char_p)))
495            {
496              return pPow(p,i);
497            }
498            return pTwoMonPower(p,i);
499          }
500          return pMonPower(p,i);
501        }
502      /*end default:*/
503    }
504  }
505  return rc;
506}
507
508/*2
509* returns the partial differentiate of a by the k-th variable
510* does not destroy the input
511*/
512poly pDiff(poly a, int k)
513{
514  poly res, f, last;
515  number t;
516
517  last = res = NULL;
518  while (a!=NULL)
519  {
520    if (pGetExp(a,k)!=0)
521    {
522      f = pNew();
523      pCopy2(f,a);
524      pNext(f)=NULL;
525      t = nInit(pGetExp(a,k));
526      pSetCoeff0(f,nMult(t,pGetCoeff(a)));
527      nDelete(&t);
528      if (nIsZero(pGetCoeff(f)))
529        pDelete1(&f);
530      else
531      {
532        pDecrExp(f,k);
533        pSetm(f);
534        if (res==NULL)
535        {
536          res=last=f;
537        }
538        else
539        {
540          pNext(last)=f;
541          last=f;
542        }
543        //res = pAdd(res, f);
544      }
545    }
546    pIter(a);
547  }
548  return res;
549}
550
551static poly pDiffOpM(poly a, poly b,BOOLEAN multiply)
552{
553  int i,j,s;
554  number n,h,hh;
555  poly p=pOne();
556  n=nMult(pGetCoeff(a),pGetCoeff(b));
557  for(i=pVariables;i>0;i--)
558  {
559    s=pGetExp(b,i);
560    if (s<pGetExp(a,i))
561    {
562      nDelete(&n);
563      pDelete1(&p);
564      return NULL;
565    }
566    if (multiply)
567    {
568      for(j=pGetExp(a,i); j>0;j--)
569      {
570        h = nInit(s);
571        hh=nMult(n,h);
572        nDelete(&h);
573        nDelete(&n);
574        n=hh;
575        s--;
576      }
577      pSetExp(p,i,s);
578    }
579    else
580    {
581      pSetExp(p,i,s-pGetExp(a,i));
582    }
583  }
584  pSetm(p);
585  /*if (multiply)*/ pSetCoeff(p,n);
586  return p;
587}
588
589poly pDiffOp(poly a, poly b,BOOLEAN multiply)
590{
591  poly result=NULL;
592  poly h;
593  for(;a!=NULL;pIter(a))
594  {
595    for(h=b;h!=NULL;pIter(h))
596    {
597      result=pAdd(result,pDiffOpM(a,h,multiply));
598    }
599  }
600  return result;
601}
602
603/*2
604* returns the length of a (numbers of monomials)
605*/
606int pLength(poly a)
607{
608  int l = 0;
609
610  while (a!=NULL)
611  {
612    pIter(a);
613    l++;
614  }
615  return l;
616}
617
618
619void pSplit(poly p, poly *h)
620{
621  *h=pNext(p);
622  pNext(p)=NULL;
623}
624
625/*2
626* returns maximal column number in the modul element a (or 0)
627*/
628int pMaxComp(poly p)
629{
630  int result,i;
631
632  if(p==NULL) return 0;
633  result = pGetComp(p);
634  while (pNext(p)!=NULL)
635  {
636    pIter(p);
637    i = pGetComp(p);
638    if (i>result) result = i;
639  }
640  return result;
641}
642
643/*2
644* returns minimal column number in the modul element a (or 0)
645*/
646int pMinComp(poly p)
647{
648  int result,i;
649
650  if(p==NULL) return 0;
651  result = pGetComp(p);
652  while (pNext(p)!=NULL)
653  {
654    pIter(p);
655    i = pGetComp(p);
656    if (i<result) result = i;
657  }
658  return result;
659}
660
661/*2
662* multiplies the polynomial a by the column generator with number i
663*/
664void pSetCompP(poly p, int i)
665{
666
667  if ((p!=NULL) && (pGetComp(p)==0))
668  {
669    do
670    {
671      pSetComp(p, (Exponent_t)i);
672      pIter(p);
673    } while (p!=NULL);
674  }
675}
676
677/*2
678* handle memory request for sets of polynomials (ideals)
679* l is the length of *p, increment is the difference (may be negative)
680*/
681void pEnlargeSet(polyset *p, int l, int increment)
682{
683  int i;
684  polyset h;
685
686  h=(polyset)ReAlloc((poly*)*p,l*sizeof(poly),(l+increment)*sizeof(poly));
687  if (increment>0)
688  {
689    for (i=l; i<l+increment; i++)
690      h[i]=NULL;
691  }
692  *p=h;
693}
694
695/*2
696* returns a polynomial representing the integer i
697*/
698poly pISet(int i)
699{
700  poly rc = NULL;
701  if (i!=0)
702  {
703    rc = pInit();
704    pSetCoeff0(rc,nInit(i));
705    if (nIsZero(pGetCoeff(rc)))
706      pDelete1(&rc);
707#ifdef TEST_MAC_ORDER
708    else if (bNoAdd)
709      pSetm(rc);
710#endif
711  }
712  return rc;
713}
714
715// /*2
716// * create a new polynomial and init it as 1
717// */
718// poly pOne(void)
719// {
720//   poly p=pInit();
721//   pSetCoeff0(p,nInit(1));
722// #ifdef TEST_MAC_ORDER
723//   pSetm(p);
724// #endif
725//   return p;
726// }
727
728void pContent(poly ph)
729{
730  number h,d;
731  poly p;
732
733  p = ph;
734  if(pNext(p)==NULL)
735  {
736    pSetCoeff(p,nInit(1));
737  }
738  else
739  {
740#ifdef PDEBUG
741    if (!pTest(p)) return;
742#endif
743    nNormalize(pGetCoeff(p));
744    if(!nGreaterZero(pGetCoeff(ph)))
745      ph = pNeg(ph);
746    h=nCopy(pGetCoeff(p));
747    pIter(p);
748    while (p!=NULL)
749    {
750      nNormalize(pGetCoeff(p));
751      d=nGcd(h,pGetCoeff(p));
752      nDelete(&h);
753      h = d;
754      if(nIsOne(h))
755      {
756        break;
757      }
758      pIter(p);
759    }
760    p = ph;
761    //number tmp;
762    if(!nIsOne(h))
763    {
764      while (p!=NULL)
765      {
766        //d = nDiv(pGetCoeff(p),h);
767        //tmp = nIntDiv(pGetCoeff(p),h);
768        //if (!nEqual(d,tmp))
769        //{
770        //  StringSetS("** div0:");nWrite(pGetCoeff(p));StringAppendS("/");
771        //  nWrite(h);StringAppendS("=");nWrite(d);StringAppendS(" int:");
772        //  nWrite(tmp);Print(StringAppendS("\n"));
773        //}
774        //nDelete(&tmp);
775        d = nIntDiv(pGetCoeff(p),h);
776        pSetCoeff(p,d);
777        pIter(p);
778      }
779    }
780    nDelete(&h);
781#ifdef HAVE_FACTORY
782    if ( (nGetChar() == 1) || (nGetChar() < 0) ) /* Q[a],Q(a),Zp[a],Z/p(a) */
783    {
784      pTest(ph);
785      singclap_divide_content(ph);
786      pTest(ph);
787    }
788#endif
789  }
790}
791
792void pCleardenom(poly ph)
793{
794  number d, h;
795  poly p;
796
797  p = ph;
798  if(pNext(p)==NULL)
799  {
800    pSetCoeff(p,nInit(1));
801  }
802  else
803  {
804    h = nInit(1);
805    while (p!=NULL)
806    {
807      nNormalize(pGetCoeff(p));
808      d=nLcm(h,pGetCoeff(p));
809      nDelete(&h);
810      h=d;
811      pIter(p);
812    }
813    /* contains the 1/lcm of all denominators */
814    if(!nIsOne(h))
815    {
816      p = ph;
817      while (p!=NULL)
818      {
819        /* should be:
820        * number hh;
821        * nGetDenom(p->coef,&hh);
822        * nMult(&h,&hh,&d);
823        * nNormalize(d);
824        * nDelete(&hh);
825        * nMult(d,p->coef,&hh);
826        * nDelete(&d);
827        * nDelete(&(p->coef));
828        * p->coef =hh;
829        */
830        d=nMult(h,pGetCoeff(p));
831        nNormalize(d);
832        pSetCoeff(p,d);
833        pIter(p);
834      }
835      nDelete(&h);
836      if (nGetChar()==1)
837      {
838        h = nInit(1);
839        p=ph;
840        while (p!=NULL)
841        {
842          d=nLcm(h,pGetCoeff(p));
843          nDelete(&h);
844          h=d;
845          pIter(p);
846        }
847        /* contains the 1/lcm of all denominators */
848        if(!nIsOne(h))
849        {
850          p = ph;
851          while (p!=NULL)
852          {
853            /* should be:
854            * number hh;
855            * nGetDenom(p->coef,&hh);
856            * nMult(&h,&hh,&d);
857            * nNormalize(d);
858            * nDelete(&hh);
859            * nMult(d,p->coef,&hh);
860            * nDelete(&d);
861            * nDelete(&(p->coef));
862            * p->coef =hh;
863            */
864            d=nMult(h,pGetCoeff(p));
865            nNormalize(d);
866            pSetCoeff(p,d);
867            pIter(p);
868          }
869          nDelete(&h);
870        }
871      }
872    }
873    pContent(ph);
874  }
875}
876
877/*2
878*tests if p is homogeneous with respect to the actual weigths
879*/
880BOOLEAN pIsHomogeneous (poly p)
881{
882  poly qp=p;
883  int o;
884
885  if ((p == NULL) || (pNext(p) == NULL)) return TRUE;
886  pFDegProc d=(pLexOrder ? pTotaldegree : pFDeg );
887  o = d(p);
888  do
889  {
890    if (d(qp) != o) return FALSE;
891    pIter(qp);
892  }
893  while (qp != NULL);
894  return TRUE;
895}
896
897// orders monoms of poly using merge sort (ususally faster than
898// insertion sort). ASSUMES that pSetm was performed on monoms
899// (i.e. that Order field is set correctly)
900poly pOrdPolyMerge(poly p)
901{
902  poly qq,pp,result=NULL;
903
904  if (p == NULL) return NULL;
905
906  for (;;)
907  {
908    qq = p;
909    for (;;)
910    {
911      if (pNext(p) == NULL) return pAdd(result, qq);
912      if (pComp(p,pNext(p)) != 1)
913      {
914        pp = p;
915        pIter(p);
916        pNext(pp) = NULL;
917        result = pAdd(result, qq);
918        break;
919      }
920      pIter(p);
921    }
922  }
923}
924
925// orders monoms of poly using insertion sort, performs pSetm on each
926// monom (i.e. sets Order field)
927poly pOrdPolyInsertSetm(poly p)
928{
929  poly qq,result = NULL;
930
931  while (p != NULL)
932  {
933    qq = p;
934    pIter(p);
935    qq->next = NULL;
936    pSetm(qq);
937    result = pAdd(result,qq);
938    pTest(result);
939  }
940  return result;
941}
942
943/*2
944*returns a re-ordered copy of a polynomial, with permutation of the variables
945*/
946poly pPermPoly (poly p, int * perm, ring oldRing,
947   int *par_perm, int OldPar)
948{
949  int OldpVariables = oldRing->N;
950  poly result = NULL;
951  poly aq=NULL;
952  poly qq;
953  int i;
954
955  while (p != NULL)
956  {
957    if (OldPar==0)
958    {
959      qq = pInit();
960      pGetCoeff(qq)=nMap(pGetCoeff(p));
961    }
962    else
963    {
964      qq=pOne();
965      aq=naPermNumber(pGetCoeff(p),par_perm,OldPar);
966      pTest(aq);
967    }
968    pSetComp(qq, pRingGetComp(oldRing,p));
969    if (nIsZero(pGetCoeff(qq)))
970    {
971      pDelete1(&qq);
972    }
973    else
974    {
975      for(i=1; i<=OldpVariables; i++)
976      {
977        if (pRingGetExp(oldRing,p,i)!=0)
978        {
979          if (perm==NULL)
980          {
981            pSetExp(qq,i, pRingGetExp(oldRing,p,i));
982          }
983          else if (perm[i]>0)
984            pAddExp(qq,perm[i], pRingGetExp(oldRing, p,i));
985          else if (perm[i]<0)
986          {
987            lnumber c=(lnumber)pGetCoeff(qq);
988            c->z->e[-perm[i]-1]+=pRingGetExp(oldRing, p,i);
989          }
990          else
991          {
992            /* this variable maps to 0 !*/
993            pDelete1(&qq);
994            break;
995          }
996        }
997      }
998    }
999    pIter(p);
1000    if (qq!=NULL)
1001    {
1002      pSetm(qq);
1003      pTest(qq);
1004      pTest(aq);
1005      if (aq!=NULL)
1006      {
1007        qq=pMult(aq,qq);
1008        aq=NULL;
1009      }
1010      result = pAdd(result,qq);
1011      pTest(result);
1012    }
1013    else if (aq!=NULL)
1014    {
1015      pDelete(&aq);
1016    }
1017  }
1018  return result;
1019}
1020
1021poly pJet(poly p, int m)
1022{
1023  poly r=NULL;
1024  poly t=NULL;
1025
1026  while (p!=NULL)
1027  {
1028    if (pTotaldegree(p)<=m)
1029    {
1030      if (r==NULL)
1031        r=pHead(p);
1032      else
1033      if (t==NULL)
1034      {
1035        pNext(r)=pHead(p);
1036        t=pNext(r);
1037      }
1038      else
1039      {
1040        pNext(t)=pHead(p);
1041        pIter(t);
1042      }
1043    }
1044    pIter(p);
1045  }
1046  return r;
1047}
1048
1049poly pJetW(poly p, int m, short *w)
1050{
1051  poly r=NULL;
1052  poly t=NULL;
1053  short *wsave=ecartWeights;
1054
1055  ecartWeights=w;
1056
1057  while (p!=NULL)
1058  {
1059    if (totaldegreeWecart(p)<=m)
1060    {
1061      if (r==NULL)
1062        r=pHead(p);
1063      else
1064      if (t==NULL)
1065      {
1066        pNext(r)=pHead(p);
1067        t=pNext(r);
1068      }
1069      else
1070      {
1071        pNext(t)=pHead(p);
1072        pIter(t);
1073      }
1074    }
1075    pIter(p);
1076  }
1077  ecartWeights=wsave;
1078  return r;
1079}
1080
1081int pDegW(poly p, short *w)
1082{
1083  int r=0;
1084  short *wsave=ecartWeights;
1085
1086  ecartWeights=w;
1087
1088  while (p!=NULL)
1089  {
1090    r=max(r, totaldegreeWecart(p));
1091    pIter(p);
1092  }
1093  ecartWeights=wsave;
1094  return r;
1095}
1096
1097/*-----------type conversions ----------------------------*/
1098/*2
1099* input: a set of polys (len elements: p[0]..p[len-1])
1100* output: a vector
1101* p will not be changed
1102*/
1103poly  pPolys2Vec(polyset p, int len)
1104{
1105  poly v=NULL;
1106  poly h;
1107  int i;
1108
1109  for (i=len-1; i>=0; i--)
1110  {
1111    if (p[i])
1112    {
1113      h=pCopy(p[i]);
1114      pSetCompP(h,i+1);
1115      v=pAdd(v,h);
1116    }
1117  }
1118 return v;
1119}
1120
1121/*2
1122* convert a vector to a set of polys,
1123* allocates the polyset, (entries 0..(*len)-1)
1124* the vector will not be changed
1125*/
1126void  pVec2Polys(poly v, polyset *p, int *len)
1127{
1128  poly h;
1129  int k;
1130
1131  *len=pMaxComp(v);
1132  if (*len==0) *len=1;
1133  *p=(polyset)Alloc0((*len)*sizeof(poly));
1134  while (v!=NULL)
1135  {
1136    h=pHead(v);
1137    k=pGetComp(h);
1138    pSetComp(h,0);
1139    (*p)[k-1]=pAdd((*p)[k-1],h);
1140    pIter(v);
1141  }
1142}
1143
1144int pVar(poly m)
1145{
1146  if (m==NULL) return 0;
1147  if (pNext(m)!=NULL) return 0;
1148  int i,e=0;
1149  for (i=pVariables; i>0; i--)
1150  {
1151    if (pGetExp(m,i)==1)
1152    {
1153      if (e==0) e=i;
1154      else return 0;
1155    }
1156  }
1157  return e;
1158}
1159
1160/*----------utilities for syzygies--------------*/
1161//BOOLEAN   pVectorHasUnitM(poly p, int * k)
1162//{
1163//  while (p!=NULL)
1164//  {
1165//    if (pIsConstantComp(p))
1166//    {
1167//      *k = pGetComp(p);
1168//      return TRUE;
1169//    }
1170//    else pIter(p);
1171//  }
1172//  return FALSE;
1173//}
1174
1175BOOLEAN   pVectorHasUnitB(poly p, int * k)
1176{
1177  poly q=p,qq;
1178  int i;
1179
1180  while (q!=NULL)
1181  {
1182    if (pIsConstantComp(q))
1183    {
1184      i = pGetComp(q);
1185      qq = p;
1186      while ((qq != q) && (pGetComp(qq) != i)) pIter(qq);
1187      if (qq == q)
1188      {
1189        *k = i;
1190        return TRUE;
1191      }
1192      else
1193        pIter(q);
1194    }
1195    else pIter(q);
1196  }
1197  return FALSE;
1198}
1199
1200void   pVectorHasUnit(poly p, int * k, int * len)
1201{
1202  poly q=p,qq;
1203  int i,j=0;
1204
1205  *len = 0;
1206  while (q!=NULL)
1207  {
1208    if (pIsConstantComp(q))
1209    {
1210      i = pGetComp(q);
1211      qq = p;
1212      while ((qq != q) && (pGetComp(qq) != i)) pIter(qq);
1213      if (qq == q)
1214      {
1215       j = 0;
1216       while (qq!=NULL)
1217       {
1218         if (pGetComp(qq)==i) j++;
1219        pIter(qq);
1220       }
1221       if ((*len == 0) || (j<*len))
1222       {
1223         *len = j;
1224         *k = i;
1225       }
1226      }
1227    }
1228    pIter(q);
1229  }
1230}
1231
1232/*2
1233* returns TRUE if p1 = p2
1234*/
1235BOOLEAN pEqualPolys(poly p1,poly p2)
1236{
1237  while ((p1 != NULL) && (p2 != NULL))
1238  {
1239    /* p1 and p2 are non-NULL, so we may use pComp0 instead of pComp */
1240    if (pComp0(p1,p2) != 0)
1241    {
1242       return FALSE;
1243    }
1244    if (nEqual(pGetCoeff(p1),pGetCoeff(p2)) == FALSE)
1245    {
1246      return FALSE;
1247    }
1248    pIter(p1);
1249    pIter(p2);
1250  }
1251  return (p1==p2);
1252}
1253
1254/*2
1255*returns TRUE if p1 is a skalar multiple of p2
1256*assume p1 != NULL and p2 != NULL
1257*/
1258BOOLEAN pComparePolys(poly p1,poly p2)
1259{
1260  number n,nn;
1261  int i;
1262
1263  if (pLength(p1) != pLength(p2))
1264    return FALSE;
1265  n=nDiv(pGetCoeff(p1),pGetCoeff(p2));
1266  while ((p1 != NULL) /*&& (p2 != NULL)*/)
1267  {
1268    for (i=1; i<=pVariables; i++)
1269    {
1270      if (pGetExp(p1,i)!=pGetExp(p2,i))
1271      {
1272        nDelete(&n);
1273        return FALSE;
1274      }
1275    }
1276    if (pGetComp(p1) != pGetComp(p2))
1277    {
1278      nDelete(&n);
1279      return FALSE;
1280    }
1281    if (!nEqual(pGetCoeff(p1),nn=nMult(pGetCoeff(p2),n)))
1282    {
1283      nDelete(&n);
1284      nDelete(&nn);
1285      return FALSE;
1286    }
1287    nDelete(&nn);
1288    pIter(p1);
1289    pIter(p2);
1290  }
1291  nDelete(&n);
1292  return TRUE;
1293}
1294//{
1295//  number m=NULL,n=NULL;
1296//
1297//  while ((p1 != NULL) && (p2 != NULL))
1298//  {
1299//    if (pComp(p1,p2) != 0)
1300//    {
1301//      if (n!=NULL) nDelete(&n);
1302//      return FALSE;
1303//    }
1304//    if (n == NULL)
1305//    {
1306//      n=nDiv(pGetCoeff(p1),pGetCoeff(p2));
1307//    }
1308//    else
1309//    {
1310//      m=nDiv(pGetCoeff(p1),pGetCoeff(p2));
1311//      if (! nEqual(m,n))
1312//      {
1313//        nDelete(&n);
1314//        nDelete(&m);
1315//        return FALSE;
1316//      }
1317//      nDelete(&m);
1318//    }
1319//    p1=pNext(p1);
1320//    p2=pNext(p2);
1321//  }
1322//  nDelete(&n);
1323//  //if (p1 != p2)
1324//  //  return FALSE;
1325//  //else
1326//  //  return TRUE;
1327//  return (p1==p2);
1328//}
Note: See TracBrowser for help on using the repository browser.