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