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 | |
---|
17 | int pdN; /* the number of x(i) / d(i) / x(i)^(-1) */ |
---|
18 | int pdK; /* the number of y(i) */ |
---|
19 | /* pVariables=pdN*2+1+pdK */ |
---|
20 | |
---|
21 | void pdSetDFlagP(poly p, int i) |
---|
22 | { |
---|
23 | while (p!=NULL) |
---|
24 | { |
---|
25 | pdSetDFlag(p,1); |
---|
26 | pSetm(p); |
---|
27 | pIter(p); |
---|
28 | } |
---|
29 | } |
---|
30 | |
---|
31 | extern 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 | */ |
---|
36 | poly 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 | */ |
---|
108 | poly 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 | */ |
---|
190 | poly 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 | */ |
---|
303 | poly 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 |
---|
344 | int 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 |
---|
356 | int 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 | */ |
---|
372 | BOOLEAN 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 | */ |
---|
400 | poly 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 | */ |
---|
442 | BOOLEAN 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; |
---|
481 | void 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 | */ |
---|
532 | void 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 | */ |
---|
604 | void 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 | |
---|
617 | BOOLEAN 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 | |
---|
640 | poly 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 |
---|