source: git/Singular/polys.inc @ 949d0d

fieker-DuValspielwiese
Last change on this file since 949d0d was 311499, checked in by Hans Schönemann <hannes@…>, 26 years ago
* hannes: DRING-changes git-svn-id: file:///usr/local/Singular/svn/trunk@2007 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 13.5 KB
Line 
1#ifdef DRING
2// D=k[x,d,y] is the Weyl-Algebra [y], y commuting with all others
3// M=k[x,x^(-1),y] is a D-module
4// all x(1..n),d,x(1..n)^(-1),y(1..k) are considered as "ring variables" v(1..N)
5// the map from x(i) to v:
6//#define pdX(i)  (i)
7// d(i)
8//#define pdDX(i) (pdN+i)
9// x(i)^(-1)
10//#define pdIX(i) (pdN+i)
11// y(i)
12//#define pdY(i)  (pdN*2+i+1)
13// a monomial m belongs to a D-module M iff pdDFlag(m)==0
14// a monomial m belongs to an ideal in the Weyl-Algebra D iff pdDFlag(m)==1
15//#define pdDFlag(m) ((m)->exp[pdN*2+1])
16
17int pdN; /* the number of x(i) / d(i) / x(i)^(-1) */
18int pdK; /* the number of y(i) */
19/* pVariables=pdN*2+1+pdK */
20
21void pdSetDFlagP(poly p, int i)
22{
23  while (p!=NULL)
24  {
25    pdSetDFlag(p,1);
26    pSetm(p);
27    pIter(p);
28  }
29}
30
31extern int binom(int n, int r);
32#define c(A,B) binom(A+B-1,A)
33/*2
34* the commutator relation of var number n and exps d and x
35*/
36poly comm(int n, short d, short x)
37{
38  poly e1=pOne();
39  poly e2=pOne();
40
41  pdSetDFlag(e1,1);
42  pdSetDFlag(e2,1);
43  number t;
44
45  if (x==1)
46  {
47    pSetExp(e1,pdX(n),1);
48    pSetExp(e1,pdDX(n),d);
49    pSetm(e1);
50    pNext(e1)=e2;
51    pSetExp(e2,pdDX(n),d-1);
52    t = nInit(d);
53    pSetCoeff(e2,t);
54    pSetm(e2);
55  }
56  else if (d==1)
57  {
58    pSetExp(e1,pdX(n),x);
59    pSetExp(e1,pdDX(n),1);
60    pSetm(e1);
61    pNext(e1)=e2;
62    pSetExp(e2,n,x-1);
63    t = nInit(x);
64    pSetCoeff(e2,t);
65    pSetm(e2);
66  }
67  else
68  {
69    int i,j,k;
70    int p;
71    int tp;
72    poly h=NULL;
73
74    pSetExp(e1,pdX(n),x);
75    pSetExp(e1,pdDX(n),d);
76    pSetm(e1);
77    for (j=1;j<=min(x,d);j++)
78    {
79      p=1;
80      for(k=0;k<j;k++)
81      {
82        p *= (d-k);
83      }
84      pSetExp(e2,pdX(n),x-j);
85      pSetExp(e2,pdDX(n),d-j);
86      pSetm(e2);
87      tp=0;
88      for(i=0;i<=x-j;i++)
89      {
90        tp+=c(i,j);
91      }
92      t = nInit(p*tp);
93      pSetCoeff(e2,t);
94      h=pAdd(h,e2);
95      e2=pOne();
96    }
97    pNext(e1)=h;
98    pFree1(e2);
99  }
100  return e1;
101}
102#undef c
103
104/*2
105* multiply 2 monomials (assume coeff of b == 1)
106* DRING-case : a and b are operators, pdDFlag==1
107*/
108poly pMultDD(poly a, poly b)
109{
110  int i;
111  poly resl=pOne();
112  poly resr=pOne();
113  poly multiply=NULL;
114  short c; /* the component number of the result*/
115
116  if((c=pGetComp(a))!=0)
117  {
118#ifdef TEST
119    if (pGetComp(b)!=0)
120    {
121      Werror("mult vector * vector");
122      return NULL;
123    }
124#endif
125  }
126  else
127    c=pGetComp(b);
128
129  pdSetDFlag(resl,1);
130  pdSetDFlag(resr,1);
131  pSetCoeff(resl,nCopy(pGetCoeff(a)));
132  // put all x from a to the left result resl,
133  // put all d from b to the right result resr
134  for (i=1;i<=pdN;i++)
135  {
136    pSetExp(resl,pdX(i),pGetExp(a,pdX(i)));
137    pSetExp(resr,pdDX(i),pGetExp(b,pdDX(i)));
138  }
139  // put all commutative vars y to the right result resr
140  for (i=1;i<=pdK;i++)
141  {
142    pSetExp(resr,pdY(i),pGetExp(a,pdY(i))+pGetExp(b,pdY(i)));
143  }
144  // set the component number
145  pSetComp(resl,c);
146  for (i=1;i<=pdN;i++)
147  {
148    if ((a->exp[pdDX(i)] !=0) && (b->exp[pdX(i)] !=0))
149    if ((pGetExp(a,pdDX(i)) !=0) && (pGetExp(b,pdX(i)) !=0))
150    {
151      if (multiply!=NULL)
152      {
153        multiply=pMult(multiply,comm(i,pGetExp(a,pdDX(i)),pGetExp(b,pdX(i))));
154      }
155      else
156      {
157        multiply=comm(i,pGetExp(a,pdDX(i)),pGetExp(b,pdX(i)));
158      }
159    }
160    else
161    {
162      pSetExp(resl,pdX(i),pGetExp(resl,pdX(i))+pGetExp(b,pdX(i)));
163      pSetExp(resr,pdDX(i),pGetExp(resr,pdDX(i))+pGetExp(a,pdDX(i)));
164    }
165  }
166  pSetm(resl); /*poly or vector*/
167  pSetm(resr); /*poly*/
168  if (multiply!=NULL)
169  {
170    resl=pMult(resl,multiply);/*(poly or vector) * poly*/
171    resl=pMult(resl,resr);
172  }
173  else
174  {
175    // now resl has only powers of x(i) and y(i), resr has only powers of d(i):
176    for (i=1;i<=pdN;i++)
177    {
178      pSetExp(resl,pdDX(i),pGetExp(resr,pdDX(i)));
179    }
180    pSetm(resl);
181    pFree1(resr);
182  }
183  return resl;
184}
185
186/*2
187* multiply 2 monomials (assume coeff of b == 1)
188* DRING-case : a is an operator, pdDFlag==1, b is in a D-module (DFlag==0)
189*/
190poly pMultDT(poly a, poly b)
191{
192  int i;
193  short c; /* the component number of the result*/
194
195  if((c=pGetComp(a))!=0)
196  {
197#ifdef TEST
198    if (pGetComp(b)!=0)
199    {
200      Werror("mult vector * vector");
201      return NULL;
202    }
203#endif
204  }
205  else
206    c=pGetComp(b);
207
208  // is the product 0 ?
209  for (i=1;i<=pdN;i++)
210  {
211    if ((pGetExp(a,pdDX(i)) > pGetExp(b,pdX(i)))
212    && (pGetExp(b,pdIX(i))==0))
213      return NULL;
214  }
215  poly resl=pOne();
216  pdSetDFlag(resl,0);
217  pSetCoeff(resl,nCopy(pGetCoeff(a)));
218  // put all x from a to the left result resl
219  for (i=1;i<=pdN;i++)
220  {
221    pSetExp(resl,pdX(i),pGetExp(a,pdX(i)));
222  }
223  // put all commutative vars y to the left result resl
224  for (i=1;i<=pdK;i++)
225  {
226    pSetExp(resl,pdY(i),pGetExp(a,pdY(i))+pGetExp(b,pdY(i)));
227  }
228  // set the component number
229  pSetComp(resl,c);
230  int q,p;
231  number n,h1,h2;
232  for (i=1;i<=pdN;i++)
233  {
234    if (((p=pGetExp(a,pdDX(i))) !=0) && ((q=pGetExp(b,pdX(i))) !=0))
235    {
236      // d^p(x^q): q*(q-1)*...*(q-p+1)* x^(q-p)
237      pSetExp(resl,pdX(i),pGetExp(resl,pdX(i))+q-p);
238      n=nInit(q);
239      q--;p--;
240      while (p>0)
241      {
242        h1=nInit(q);
243        h2=nMult(n,h1);
244        nDelete(&h1);
245        nDelete(&n);
246        n=h2;
247        q--;
248        p--;
249      }
250      h1=nMult(pGetCoeff(resl),n);
251      pSetCoeff(resl,h1);
252      nDelete(&n);
253    }
254    else
255    if (((p=pGetExp(a,pdDX(i))) !=0) && ((q=pGetExp(b,pdIX(i))) !=0))
256    {
257      // d^p(x^(-q)): (-1)^p*q*(q+1)*...*(q+p-1)* x^(-(q+p))
258      pSetExp(resl,pdIX(i),pGetExp(resl,pdIX(i))+q+p);
259      if (p & 1) n=nInit(-q);
260      else       n=nInit(q);
261      // last index: p=q+p-1
262      q++;p--;
263      while (p>0)
264      {
265        h1=nInit(q);
266        h2=nMult(n,h1);
267        nDelete(&h1);
268        nDelete(&n);
269        n=h2;
270        q++;p--;
271      }
272      h1=nMult(pGetCoeff(resl),n);
273      pSetCoeff(resl,h1);
274      nDelete(&n);
275    }
276    else
277    {
278      pSetExp(resl,pdX(i),pGetExp(resl,pdX(i))+pGetExp(b,pdX(i)));
279      pSetExp(resl,pdIX(i),pGetExp(resl,pdIX(i))+pGetExp(b,pdIX(i)));
280    }
281  }
282  for(i=1;i<=pdN;i++)
283  {
284    if (pGetExp(resl,pdX(i))>=pGetExp(resl,pdIX(i)))
285    {
286      pSetExp(resl,pdX(i),pGetExp(resl,pdX(i))-pGetExp(resl,pdIX(i)));
287      pSetExp(resl,pdIX(i),0);
288    }
289    else
290    {
291      pSetExp(resl,pdIX(i),pGetExp(resl,pdIX(i))-pGetExp(resl,pdX(i)));
292      pSetExp(resl,pdX(i),0);
293    }
294  }
295  pSetm(resl); /*poly or vector*/
296  return resl;
297}
298
299/*2
300* multiply 2 monomials (assume coeff of b == 1)
301* DRING-case : a and b is in a D-module (DFlag==0)
302*/
303poly pMultTT(poly a, poly b)
304{
305  int i;
306  short c; /* the component number of the result*/
307
308  if((c=pGetComp(a))!=0)
309  {
310#ifdef TEST
311    if (pGetComp(b)!=0)
312    {
313      Werror("mult vector * vector");
314      return NULL;
315    }
316#endif
317  }
318  else
319    c=pGetComp(b);
320
321  poly resl=pOne();
322  int t;
323  pdSetDFlag(resl,0);
324  pSetCoeff(resl,nCopy(pGetCoeff(a)));
325  for (i=1;i<=pdN;i++)
326  {
327    t=pGetExp(a,pdX(i))+pGetExp(b,pdX(i))-pGetExp(a,pdIX(i))-pGetExp(b,pdIX(i));
328    if (t>=0)  pSetExp(resl,pdX(i),t);
329    else       pSetExp(resl,pdIX(i),-t);
330  }
331  // put all commutative vars y to the result resl
332  for (i=1;i<=pdK;i++)
333  {
334    pSetExp(resl,pdY(i),pGetExp(a,pdY(i))+pGetExp(b,pdY(i)));
335  }
336  // set the component number
337  pSetComp(resl,c);
338  pSetm(resl); /*poly or vector*/
339  return resl;
340}
341#endif
342
343#ifdef SRING
344int psFirst(poly p, int i)
345{
346  loop
347  {
348    if (i>pVariables) return pVariables;
349    if (pGetExp(p,i) != 0) return i;
350    i++;
351  }
352}
353#endif
354
355#ifdef SRING
356int psLast(poly p, int i)
357{
358  loop
359  {
360    if (i<pAltVars) return pAltVars;
361    if (pGetExp(p,i) != 0) return i;
362    i--;
363  }
364}
365#endif
366
367#ifdef SRING
368/*2
369* returns TRUE, if the product of the monomials p1
370* and p2 is not zero
371*/
372BOOLEAN psMultTest(poly p1, poly p2)
373{
374  int i1,i2,j;
375
376  if ((pGetComp(p1) !=0) && (pGetComp(p2) != 0))
377    return FALSE;
378  j=pAltVars;
379  loop
380  {
381    if (j>pVariables) return TRUE;
382    i1=pGetExp(p1,j);
383    i2=pGetExp(p2,j);
384    if ((i1==1) && (i2==1)) return FALSE;
385    if ((i1>1) || (i2>1))
386    {
387      Werror("internal error in psMultTest");
388      HALT();
389    }
390    j++;
391  }
392}
393#endif
394
395#ifdef SRING
396/*2
397* multiply  monomials a and b (without coeffs)
398* take the coeff from a
399*/
400poly psMultM(poly a, poly b)
401{
402  int i, j;
403  BOOLEAN positiv=TRUE;
404  poly res=NULL;
405
406  if (psMultTest(a,b))
407  {
408    res=pNew();
409    memset(res,0,pMonomSize);
410    pGetCoeff(res)=nCopy(pGetCoeff(a));
411    i=pAltVars-1;
412    do
413    {
414      i=psFirst(b,i+1);
415      j=pVariables+1;
416      loop
417      {
418        j=psLast(a,j-1);
419        if (j>i)
420          positiv= !positiv;
421        else
422          break;
423      }
424    } while (i != pVariables);
425    for (i=0;i<=pVariables; i++)
426    {
427      pSetExp(res,i,pGetExp(a,i)+pGetExp(b,i));
428    }
429    pSetm(res);
430    if (!positiv)
431      pGetCoeff(res)=nNeg(pGetCoeff(res));
432  }
433  return res;
434}
435#endif
436
437#ifdef SDRING
438/*2
439*puts a poly into a polyset,
440*returns FALSE, if there is a trivial multiple in the set
441*/
442BOOLEAN psEnter(poly p,polyset *s, int *l, int *m)
443{
444  int i=*l;
445
446  /*is there already a multiple in s ?*/
447  pTest(p);
448  //Print("psEnter:  ");pWrite0(p);
449  for(;i > 0;i--)
450  {
451   pTest((*s)[i]);
452   if (pComparePolys(p,(*s)[i]))
453   {
454     //Print("-- ist vielfaches von ");
455     //pWrite((*s)[i]);
456     pDelete(&p);
457     return FALSE;
458   }
459  }
460  /*is there enough space in s ?*/
461  if (*l==(*m-1))
462  {
463    pEnlargeSet(s,*m,16);
464    *m+=16;
465  }
466  (*l)++;
467  //Print("++an pos %d\n",*l);
468  (*s)[*l]=p;
469  pTest(p);
470  return TRUE;
471}
472#endif
473
474#ifdef SRING
475/*2
476* create the augmentation set of p
477* p and done should be destroyed
478* done is the product of all vars already multiplied with
479*/
480//int auglev=0;
481void psAug(poly p, poly done, polyset *s, int *l, int *m)
482{
483  int an;
484  poly doneCopy, q;
485
486  if (p==NULL)
487  {
488    pDelete(&done);
489    return;
490  }
491//  auglev++;
492//  Print("lev %d, aug of ",auglev); wrp(p); Print("\n");
493  if (pSRING)
494  for (an=pAltVars; an<=pVariables; an++)
495  {
496    if ((pGetExp(p,an)==1) && (pGetExp(done,an)==0))
497    {
498      q=pOne();
499      doneCopy=pCopy(done);
500      pSetExp(doneCopy,an,1);
501      pSetExp(q,an,1);
502      pSetm(q);
503      q=pMult(q,pCopy(p));
504      if (q!=NULL)
505      {
506//        Print("lev %d, to ",auglev); wrp(q); Print("\n");
507        if (psEnter(q,s,l,m))
508        {
509          psAug(pCopy(q),doneCopy,s,l,m);
510        }
511        else
512        {
513          q=NULL;
514          pDelete(&doneCopy);
515        }
516      }
517    }
518  }
519//  Print("end aug");
520//  auglev--;
521  pDelete(&done);
522  pDelete(&p);
523}
524#endif
525#ifdef DRING
526/*2
527* create the augmentation set of p
528* p and done should be destroyed
529* done is the product of all vars already multiplied with
530* x^i --> dx^(i+1) und x^-i --> dx*x^i
531*/
532void pdAug(poly p, polyset *s, int *l, int *m)
533{
534  int an;
535  poly q, dif;
536
537  if ((p==NULL) || pdDFlag(p)==1)
538  {
539    return;
540  }
541  //Print(" aug of ");  wrp(p); Print("\n");
542  if (pDRING)
543  for (an=1; an<=pdN; an++)
544  {
545    if(pGetExp(p,pdIX(an))==0)
546    {
547      q=pOne();
548      pSetExp(q,pdDX(an),pGetExp(p,pdX(an))+1);
549      pdSetDFlag(q,1);
550      pSetm(q);
551      //Print("1:  q ");wrp(q); Print("\n");
552      q=pMult(q,pCopy(p));
553      if (q!=NULL)
554      {
555          //Print("1:  to ");  wrp(q); Print("\n");
556        if (psEnter(q,s,l,m))
557        {
558          pdAug(pCopy(q),s,l,m);
559        }
560        else
561        {
562          pDelete(&q);
563        }
564      }
565    }
566  }
567  for (an=1; an<=pdN; an++)
568  {
569    if (pGetExp(p,pdIX(an))!=0)
570    {
571      q=pOne();
572      dif=pOne();
573      pSetExp(q,pdX(an),pGetExp(p,pdIX(an)));
574      pSetExp(dif,pdDX(an),1);
575      pdSetDFlag(q,1);
576      pdSetDFlag(dif,1);
577      pSetm(q);
578      //Print("2: q ");wrp(q); Print("\n");
579      pSetm(dif);
580      //Print("2: dif ");wrp(dif);  Print("\n");
581      q=pMult(q,pCopy(p));
582      //wrp(q); Print("\n");
583      q=pMult(dif,q);
584      if (q!=NULL)
585      {
586        //Print("2: to ");wrp(q); Print("\n");
587        if (psEnter(q,s,l,m))
588        {
589          pdAug(pCopy(q),s,l,m);
590        }
591        else
592        {
593          pDelete(&q);
594        }
595      }
596    }
597  }
598//Print("end aug");Print("\n");
599  pDelete(&p);
600}
601/*2
602* returns the differential LCM of the head terms of a and b in *m
603*/
604void pdLcm(poly a, poly b, poly m)
605{
606  int i;
607  for (i=2*pdN; i>=0; i--)
608  {
609     pSetExp(m,i, min(pGetExp(a,i),pGetExp(b,i)));
610  }
611  for (i=pdY(1); i<=pdK; i++)
612  {
613     pSetExp(m,i, max(pGetExp(a,i),pGetExp(b,i)));
614  }
615}
616
617BOOLEAN pdIsConstantComp(poly p)
618{
619  if (!pDRING) return FALSE;
620
621  int i;
622
623  for (i=pdN;i>0;i--)
624  {
625    if (pGetExp(p,pdX(i))!=0) return FALSE;
626    if (pGetExp(p,pdDX(i))!=0) return FALSE;
627  }
628  for (i=pdK;i>0;i--)
629  {
630    if (pGetExp(p,pdY(i))!=0) return FALSE;
631  }
632  return TRUE;
633}
634
635
636/*2
637* returns the differential Spolynomial of a and b in *m
638*/
639
640poly pdSpolyCreate(poly p1,poly p2)
641{
642  poly m1=pOne();
643  poly m2=pOne();
644  poly pp1=p1,pp2=p2;
645  number n1,n2;
646  int i,j;
647  int co=0;
648  if (pGetComp(p1)!=pGetComp(p2))
649  {
650    if (pGetComp(p1)==0)
651    {
652      co=1;
653      pSetCompP(p1,pGetComp(p2));
654    }
655    else
656    {
657      co=2;
658      pSetCompP(p2,pGetComp(p1));
659    }
660  }
661  for (i=1;i<=pdN;i++)
662  {
663    j=pGetExp(p2,pdX(i))-pGetExp(p1,pdX(i));
664    if (j>0)
665    {
666      pSetExp(m2,pdDX(i),j);
667    }
668    else
669    {
670      pSetExp(m1,pdDX(i),(-j));
671    }
672  }
673  for (i=1;i<=pdN;i++)
674  {
675    j=pGetExp(p2,pdIX(i))-pGetExp(p1,pdIX(i));
676    if (j>0)
677    {
678      pSetExp(m2,pdX(i),j);
679    }
680    else
681    {
682      pSetExp(m1,pdX(i),(-j));
683    }
684  }
685  for (i=1;i<=pdK;i++)
686  {
687    j=pGetExp(p2,pdY(i))-pGetExp(p1,pdY(i));
688    if (j>0)
689    {
690      pSetExp(m1,pdY(i),j);
691    }
692    else
693    {
694      pSetExp(m2,pdY(i),(-j));
695    }
696  }
697  pSetm(m1);
698  pSetm(m2);
699
700  p1=pMult(m1,pCopy(p1));
701  p2=pMult(m2,pCopy(p2));
702
703  n1=nCopy(pGetCoeff(p2));
704  n2=nCopy(pGetCoeff(p1));
705
706  pDelete1(&p1);
707  pDelete1(&p2);
708
709  n1=nNeg(n1);
710  pMultN(p1,n2);
711  nDelete(&n2);
712  pMultN(p2,n1);
713  nDelete(&n1);
714
715  m1=pAdd(p1,p2);
716  if (co==1) spModuleToPoly(pp1);
717  else if (co==2) spModuleToPoly(pp2);
718  return m1;
719}
720#endif
Note: See TracBrowser for help on using the repository browser.