1 | /**************************************** |
---|
2 | * Computer Algebra System SINGULAR * |
---|
3 | ****************************************/ |
---|
4 | /*************************************************************** |
---|
5 | * File: gring.cc |
---|
6 | * Purpose: noncommutative kernel procedures |
---|
7 | * Author: levandov (Viktor Levandovsky) |
---|
8 | * Created: 8/00 - 11/00 |
---|
9 | * Version: $Id: gring.cc,v 1.5 2004-03-25 21:19:12 levandov Exp $ |
---|
10 | *******************************************************************/ |
---|
11 | #include "mod2.h" |
---|
12 | #ifdef HAVE_PLURAL |
---|
13 | #include "gring.h" |
---|
14 | #include "febase.h" |
---|
15 | #include "ring.h" |
---|
16 | #include "polys.h" |
---|
17 | #include "numbers.h" |
---|
18 | #include "ideals.h" |
---|
19 | #include "matpol.h" |
---|
20 | #include "kbuckets.h" |
---|
21 | #include "kstd1.h" |
---|
22 | #include "sbuckets.h" |
---|
23 | #include "prCopy.h" |
---|
24 | #include "p_Mult_q.h" |
---|
25 | |
---|
26 | /* global nc_macros : */ |
---|
27 | #define freeT(A,v) omFreeSize((ADDRESS)A,(v+1)*sizeof(int)) |
---|
28 | #define freeN(A,k) omFreeSize((ADDRESS)A,k*sizeof(number)) |
---|
29 | |
---|
30 | /* poly functions defined in p_Procs : */ |
---|
31 | poly nc_pp_Mult_mm(poly p, poly m, const ring r, poly &last) |
---|
32 | { |
---|
33 | return( nc_p_Mult_mm_Common(p_Copy(p,r), m, 1, r) ); |
---|
34 | } |
---|
35 | |
---|
36 | poly nc_p_Mult_mm(poly p, const poly m, const ring r) |
---|
37 | { |
---|
38 | return( nc_p_Mult_mm_Common(p, m, 1, r) ); |
---|
39 | } |
---|
40 | |
---|
41 | poly nc_mm_Mult_p(const poly m, poly p, const ring r) |
---|
42 | { |
---|
43 | return( nc_p_Mult_mm_Common(p, m, 0, r) ); |
---|
44 | } |
---|
45 | |
---|
46 | /* poly nc_p_Mult_mm(poly p, poly m, const ring r); defined below */ |
---|
47 | poly nc_p_Minus_mm_Mult_qq(poly p, const poly m, poly q, const ring r) |
---|
48 | { |
---|
49 | poly mc=p_Neg(p_Copy(m,r),r); |
---|
50 | poly mmc=nc_mm_Mult_p(mc,p_Copy(q,r),r); |
---|
51 | p_Delete(&mc,r); |
---|
52 | p=p_Add_q(p,mmc,r); |
---|
53 | return(p); |
---|
54 | } |
---|
55 | |
---|
56 | //----------- auxiliary routines-------------------------- |
---|
57 | poly _nc_p_Mult_q(poly p, poly q, const int copy, const ring r) |
---|
58 | /* destroy p,q unless copy=1 */ |
---|
59 | { |
---|
60 | poly res=NULL; |
---|
61 | poly ghost=NULL; |
---|
62 | poly qq,pp; |
---|
63 | if (copy) |
---|
64 | { |
---|
65 | qq=p_Copy(q,r); |
---|
66 | pp=p_Copy(p,r); |
---|
67 | } |
---|
68 | else |
---|
69 | { |
---|
70 | qq=q; |
---|
71 | pp=p; |
---|
72 | } |
---|
73 | while (qq!=NULL) |
---|
74 | { |
---|
75 | res=p_Add_q(res, nc_pp_Mult_mm(pp, p_Head(qq,r), r, ghost), r); |
---|
76 | qq=p_LmDeleteAndNext(qq,r); |
---|
77 | } |
---|
78 | p_Delete(&pp,r); |
---|
79 | return(res); |
---|
80 | } |
---|
81 | |
---|
82 | poly nc_p_Mult_mm_Common(poly p, const poly m, int side, const ring r) |
---|
83 | /* p is poly, m is mono with coeff, destroys p */ |
---|
84 | /* if side==1, computes p_Mult_mm; otherwise, mm_Mult_p */ |
---|
85 | { |
---|
86 | if ((p==NULL) || (m==NULL)) return NULL; |
---|
87 | /* if (pNext(p)==NULL) return(nc_mm_Mult_nn(p,pCopy(m),r)); */ |
---|
88 | /* excluded - the cycle will do it anyway - OK. */ |
---|
89 | if (p_IsConstant(m,r)) return(p_Mult_nn(p,p_GetCoeff(m,r),r)); |
---|
90 | |
---|
91 | #ifdef PDEBUG |
---|
92 | p_Test(p,r); |
---|
93 | p_Test(m,r); |
---|
94 | #endif |
---|
95 | poly v=NULL; |
---|
96 | int rN=r->N; |
---|
97 | int *P=(int *)omAlloc0((rN+1)*sizeof(int)); |
---|
98 | int *M=(int *)omAlloc0((rN+1)*sizeof(int)); |
---|
99 | /* coefficients: */ |
---|
100 | number cP,cM,cOut; |
---|
101 | p_GetExpV(m, M, r); |
---|
102 | cM=p_GetCoeff(m,r); |
---|
103 | /* components:*/ |
---|
104 | const int expM=p_GetComp(m,r); |
---|
105 | int expP=0; |
---|
106 | int expOut=0; |
---|
107 | /* bucket constraints: */ |
---|
108 | int UseBuckets=1; |
---|
109 | if (pLength(p)< MIN_LENGTH_BUCKET || TEST_OPT_NOT_BUCKETS) UseBuckets=0; |
---|
110 | sBucket_pt bu_out; |
---|
111 | poly out=NULL; |
---|
112 | if (UseBuckets) bu_out=sBucketCreate(r); |
---|
113 | |
---|
114 | while (p!=NULL) |
---|
115 | { |
---|
116 | #ifdef PDEBUG |
---|
117 | p_Test(p,r); |
---|
118 | #endif |
---|
119 | expP=p_GetComp(p,r); |
---|
120 | if (expP==0) |
---|
121 | { |
---|
122 | expOut=expM; |
---|
123 | } |
---|
124 | else |
---|
125 | { |
---|
126 | if (expM==0) |
---|
127 | { |
---|
128 | expOut=expP; |
---|
129 | #ifdef PDEBUG |
---|
130 | if (side) |
---|
131 | { |
---|
132 | Print("Multiplication in the left module from the right"); |
---|
133 | } |
---|
134 | #endif |
---|
135 | } |
---|
136 | else |
---|
137 | { |
---|
138 | /* REPORT_ERROR */ |
---|
139 | const char* s; |
---|
140 | if (side==1) s="nc_p_Mult_mm"; |
---|
141 | else s="nc_mm_Mult_p"; |
---|
142 | Print("%s: exponent mismatch %d and %d\n",s,expP,expM); |
---|
143 | expOut=0; |
---|
144 | } |
---|
145 | } |
---|
146 | p_GetExpV(p,P,r); |
---|
147 | cP=p_GetCoeff(p,r); |
---|
148 | cOut=n_Mult(cP,cM,r); |
---|
149 | if (side==1) |
---|
150 | { |
---|
151 | v = nc_mm_Mult_nn(P, M, r); |
---|
152 | } |
---|
153 | else |
---|
154 | { |
---|
155 | v = nc_mm_Mult_nn(M, P, r); |
---|
156 | } |
---|
157 | v = p_Mult_nn(v,cOut,r); |
---|
158 | p_SetCompP(v,expOut,r); |
---|
159 | if (UseBuckets) sBucket_Add_p(bu_out,v,pLength(v)); |
---|
160 | else out = p_Add_q(out,v,r); |
---|
161 | p_DeleteLm(&p,r); |
---|
162 | } |
---|
163 | freeT(P,rN); |
---|
164 | freeT(M,rN); |
---|
165 | if (UseBuckets) |
---|
166 | { |
---|
167 | out = NULL; |
---|
168 | int len = pLength(out); |
---|
169 | sBucketDestroyAdd(bu_out, &out, &len); |
---|
170 | } |
---|
171 | #ifdef PDEBUG |
---|
172 | p_Test(out,r); |
---|
173 | #endif |
---|
174 | return(out); |
---|
175 | } |
---|
176 | |
---|
177 | poly nc_mm_Mult_nn(int *F0, int *G0, const ring r) |
---|
178 | /* destroys nothing, no coeffs and exps */ |
---|
179 | { |
---|
180 | poly out=NULL; |
---|
181 | int i,j; |
---|
182 | int iF,jG,iG; |
---|
183 | int rN=r->N; |
---|
184 | int ExpSize=(((rN+1)*sizeof(int)+sizeof(long)-1)/sizeof(long))*sizeof(long); |
---|
185 | |
---|
186 | int *F=(int *)omAlloc0((rN+1)*sizeof(int)); |
---|
187 | int *G=(int *)omAlloc0((rN+1)*sizeof(int)); |
---|
188 | |
---|
189 | memcpy(F, F0,(rN+1)*sizeof(int)); |
---|
190 | // pExpVectorCopy(F,F0); |
---|
191 | memcpy(G, G0,(rN+1)*sizeof(int)); |
---|
192 | // pExpVectorCopy(G,G0); |
---|
193 | F[0]=0; /* important for p_MemAdd */ |
---|
194 | G[0]=0; |
---|
195 | |
---|
196 | iF=rN; |
---|
197 | while ((F[iF]==0)&&(iF>=1)) iF--; /* last exp_num of F */ |
---|
198 | if (iF==0) /* F0 is zero vector */ |
---|
199 | { |
---|
200 | out=pOne(); |
---|
201 | p_SetExpV(out,G0,r); |
---|
202 | p_Setm(out,r); |
---|
203 | freeT(F,rN); |
---|
204 | freeT(G,rN); |
---|
205 | return(out); |
---|
206 | } |
---|
207 | jG=1; |
---|
208 | while ((G[jG]==0)&&(jG<rN)) jG++; /* first exp_num of G */ |
---|
209 | iG=rN; |
---|
210 | while ((G[iG]==0)&&(iG>1)) iG--; /* last exp_num of G */ |
---|
211 | |
---|
212 | out=pOne(); |
---|
213 | |
---|
214 | if (iF<=jG) |
---|
215 | /* i.e. no mixed exp_num , MERGE case */ |
---|
216 | { |
---|
217 | p_MemAdd_LengthGeneral(F, G, ExpSize/sizeof(long)); |
---|
218 | p_SetExpV(out,F,r); |
---|
219 | p_Setm(out,r); |
---|
220 | // omFreeSize((ADDRESS)F,ExpSize); |
---|
221 | freeT(F,rN); |
---|
222 | freeT(G,rN); |
---|
223 | return(out); |
---|
224 | } |
---|
225 | |
---|
226 | number cff=n_Init(1,r); |
---|
227 | number tmp_num=NULL; |
---|
228 | int cpower=0; |
---|
229 | |
---|
230 | if (r->nc->type==nc_skew) |
---|
231 | { |
---|
232 | if (r->nc->IsSkewConstant==1) |
---|
233 | { |
---|
234 | int tpower=0; |
---|
235 | for(j=jG; j<=iG; j++) |
---|
236 | { |
---|
237 | if (G[j]!=0) |
---|
238 | { |
---|
239 | cpower = 0; |
---|
240 | for(i=j+1; i<=iF; i++) |
---|
241 | { |
---|
242 | cpower = cpower + F[i]; |
---|
243 | } |
---|
244 | cpower = cpower*G[j]; |
---|
245 | tpower = tpower + cpower; |
---|
246 | } |
---|
247 | } |
---|
248 | cff = n_Copy(p_GetCoeff(MATELEM(r->nc->COM,1,2),r),r); |
---|
249 | nPower(cff,tpower,&tmp_num); |
---|
250 | n_Delete(&cff,r); |
---|
251 | cff = tmp_num; |
---|
252 | } |
---|
253 | else /* skew commutative with nonequal coeffs */ |
---|
254 | { |
---|
255 | number totcff=n_Init(1,r); |
---|
256 | for(j=jG; j<=iG; j++) |
---|
257 | { |
---|
258 | if (G[j]!=0) |
---|
259 | { |
---|
260 | cpower = 0; |
---|
261 | for(i=j+1; i<=iF; i++) |
---|
262 | { |
---|
263 | if (F[i]!=0) |
---|
264 | { |
---|
265 | cpower = F[i]*G[j]; |
---|
266 | cff = n_Copy(p_GetCoeff(MATELEM(r->nc->COM,j,i),r),r); |
---|
267 | nPower(cff,cpower,&tmp_num); |
---|
268 | cff = nMult(totcff,tmp_num); |
---|
269 | nDelete(&totcff); |
---|
270 | nDelete(&tmp_num); |
---|
271 | totcff = n_Copy(cff,r); |
---|
272 | n_Delete(&cff,r); |
---|
273 | } |
---|
274 | } /* end 2nd for */ |
---|
275 | } |
---|
276 | } |
---|
277 | cff=totcff; |
---|
278 | } |
---|
279 | p_MemAdd_LengthGeneral(F, G, ExpSize/sizeof(long)); |
---|
280 | p_SetExpV(out,F,r); |
---|
281 | p_Setm(out,r); |
---|
282 | p_SetCoeff(out,cff,r); |
---|
283 | // p_MemAdd_NegWeightAdjust(p, r); ??? do we need this? |
---|
284 | freeT(F,rN); |
---|
285 | freeT(G,rN); |
---|
286 | return(out); |
---|
287 | } /* end nc_skew */ |
---|
288 | |
---|
289 | /* now we have to destroy out! */ |
---|
290 | p_Delete(&out,r); |
---|
291 | out = NULL; |
---|
292 | |
---|
293 | if (iG==jG) |
---|
294 | /* g is univariate monomial */ |
---|
295 | { |
---|
296 | /* if (ri->nc->type==nc_skew) -- postpone to TU */ |
---|
297 | out = nc_mm_Mult_uu(F,jG,G[jG],r); |
---|
298 | freeT(F,rN); |
---|
299 | freeT(G,rN); |
---|
300 | return(out); |
---|
301 | } |
---|
302 | |
---|
303 | number n1=n_Init(1,r); |
---|
304 | int *Prv=(int *)omAlloc0((rN+1)*sizeof(int)); |
---|
305 | int *Nxt=(int *)omAlloc0((rN+1)*sizeof(int)); |
---|
306 | |
---|
307 | int *log=(int *)omAlloc0((rN+1)*sizeof(int)); |
---|
308 | int cnt=0; int cnf=0; |
---|
309 | |
---|
310 | /* splitting F wrt jG */ |
---|
311 | for (i=1;i<=jG;i++) |
---|
312 | { |
---|
313 | Prv[i]=F[i]; Nxt[i]=0; /* mult at the very end */ |
---|
314 | if (F[i]!=0) cnf++; |
---|
315 | } |
---|
316 | |
---|
317 | if (cnf==0) freeT(Prv,rN); |
---|
318 | |
---|
319 | for (i=jG+1;i<=rN;i++) |
---|
320 | { |
---|
321 | Nxt[i]=F[i]; |
---|
322 | /* if (cnf!=0) Prv[i]=0; */ |
---|
323 | if (F[i]!=0) |
---|
324 | { |
---|
325 | cnt++; |
---|
326 | } /* effective part for F */ |
---|
327 | } |
---|
328 | freeT(F,rN); |
---|
329 | cnt=0; |
---|
330 | |
---|
331 | for (i=1;i<=rN;i++) |
---|
332 | { |
---|
333 | if (G[i]!=0) |
---|
334 | { |
---|
335 | cnt++; |
---|
336 | log[cnt]=i; |
---|
337 | } /* lG for G */ |
---|
338 | } |
---|
339 | |
---|
340 | /* ---------------------- A C T I O N ------------------------ */ |
---|
341 | poly D=NULL; |
---|
342 | poly Rout=NULL; |
---|
343 | number *c=(number *)omAlloc0((rN+1)*sizeof(number)); |
---|
344 | c[0]=n_Init(1,r); |
---|
345 | |
---|
346 | int *Op=Nxt; |
---|
347 | int *On=G; |
---|
348 | int *U=(int *)omAlloc0((rN+1)*sizeof(int)); |
---|
349 | |
---|
350 | for (i=jG;i<=rN;i++) U[i]=Nxt[i]+G[i]; /* make leadterm */ |
---|
351 | Nxt=NULL; |
---|
352 | G=NULL; |
---|
353 | cnt=1; |
---|
354 | int t=0; |
---|
355 | poly w=NULL; |
---|
356 | poly Pn=pOne(); |
---|
357 | p_SetExpV(Pn,On,r); |
---|
358 | p_Setm(Pn,r); |
---|
359 | |
---|
360 | while (On[iG]!=0) |
---|
361 | { |
---|
362 | t=log[cnt]; |
---|
363 | |
---|
364 | w=nc_mm_Mult_uu(Op,t,On[t],r); |
---|
365 | c[cnt]=n_Mult(c[cnt-1],p_GetCoeff(w,r),r); |
---|
366 | D = pNext(w); /* getting coef and rest D */ |
---|
367 | p_DeleteLm(&w,r); |
---|
368 | w=NULL; |
---|
369 | |
---|
370 | Op[t] += On[t]; /* update exp_vectors */ |
---|
371 | On[t] = 0; |
---|
372 | |
---|
373 | if (t!=iG) /* not the last step */ |
---|
374 | { |
---|
375 | p_SetExpV(Pn,On,r); |
---|
376 | p_Setm(Pn,r); |
---|
377 | #ifdef PDEBUG |
---|
378 | p_Test(Pn,r); |
---|
379 | #endif |
---|
380 | |
---|
381 | // if (pNext(D)==0) |
---|
382 | // is D a monomial? could be postponed higher |
---|
383 | // { |
---|
384 | // Rout=nc_mm_Mult_nn(D,Pn,r); |
---|
385 | // } |
---|
386 | // else |
---|
387 | // { |
---|
388 | Rout=nc_p_Mult_mm(D,Pn,r); |
---|
389 | // } |
---|
390 | } |
---|
391 | else |
---|
392 | { |
---|
393 | Rout=D; |
---|
394 | D=NULL; |
---|
395 | } |
---|
396 | |
---|
397 | if (Rout!=NULL) |
---|
398 | { |
---|
399 | Rout=p_Mult_nn(Rout,c[cnt-1],r); /* Rest is ready */ |
---|
400 | out=p_Add_q(out,Rout,r); |
---|
401 | Rout=NULL; |
---|
402 | } |
---|
403 | cnt++; |
---|
404 | } |
---|
405 | freeT(On,rN); |
---|
406 | freeT(Op,rN); |
---|
407 | p_Delete(&Pn,r); |
---|
408 | omFreeSize((ADDRESS)log,(rN+1)*sizeof(int)); |
---|
409 | |
---|
410 | /* leadterm and Prv-part */ |
---|
411 | |
---|
412 | Rout=pOne(); |
---|
413 | /* U is lead.monomial */ |
---|
414 | U[0]=0; |
---|
415 | p_SetExpV(Rout,U,r); |
---|
416 | p_Setm(Rout,r); /* use again this name Rout */ |
---|
417 | #ifdef PDEBUG |
---|
418 | p_Test(Rout,r); |
---|
419 | #endif |
---|
420 | p_SetCoeff(Rout,c[cnt-1],r); |
---|
421 | out=p_Add_q(out,Rout,r); |
---|
422 | freeT(U,rN); |
---|
423 | freeN(c,rN+1); |
---|
424 | if (cnf!=0) /* Prv is non-zero vector */ |
---|
425 | { |
---|
426 | Rout=pOne(); |
---|
427 | Prv[0]=0; |
---|
428 | p_SetExpV(Rout,Prv,r); |
---|
429 | p_Setm(Rout,r); |
---|
430 | #ifdef PDEBUG |
---|
431 | p_Test(Rout,r); |
---|
432 | #endif |
---|
433 | out=nc_mm_Mult_p(Rout,out,r); /* getting the final result */ |
---|
434 | freeT(Prv,rN); |
---|
435 | p_Delete(&Rout,r); |
---|
436 | } |
---|
437 | return (out); |
---|
438 | } |
---|
439 | |
---|
440 | |
---|
441 | poly nc_mm_Mult_uu(int *F,int jG,int bG, const ring r) |
---|
442 | /* f=mono(F),g=(x_iG)^bG */ |
---|
443 | { |
---|
444 | poly out=NULL; |
---|
445 | int i; |
---|
446 | number num=NULL; |
---|
447 | |
---|
448 | int rN=r->N; |
---|
449 | int iF=r->N; |
---|
450 | while ((F[iF]==0)&&(iF>0)) iF-- ; /* last exponent_num of F */ |
---|
451 | |
---|
452 | if (iF==0) /* F==zero vector in other words */ |
---|
453 | { |
---|
454 | out=pOne(); |
---|
455 | p_SetExp(out,jG,bG,r); |
---|
456 | p_Setm(out,r); |
---|
457 | return(out); |
---|
458 | } |
---|
459 | |
---|
460 | int jF=1; |
---|
461 | while ((F[jF]==0)&&(jF<=rN)) jF++; /* first exp of F */ |
---|
462 | |
---|
463 | if (iF<=jG) /* i.e. no mixed exp_num */ |
---|
464 | { |
---|
465 | out=pOne(); |
---|
466 | F[jG]=F[jG]+bG; |
---|
467 | p_SetExpV(out,F,r); |
---|
468 | p_Setm(out,r); |
---|
469 | return(out); |
---|
470 | } |
---|
471 | |
---|
472 | if (iF==jF) /* uni times uni */ |
---|
473 | { |
---|
474 | out=nc_uu_Mult_ww(iF,F[iF],jG,bG,r); |
---|
475 | return(out); |
---|
476 | } |
---|
477 | |
---|
478 | /* Now: F is mono with >=2 exponents, jG<iF */ |
---|
479 | /* check the quasi-commutative case */ |
---|
480 | // matrix LCOM=r->nc->COM; |
---|
481 | // number rescoef=n_Init(1,r); |
---|
482 | // number tmpcoef=n_Init(1,r); |
---|
483 | // int tmpint; |
---|
484 | // i=iF; |
---|
485 | // while (i>=jG+1) |
---|
486 | // /* all the non-zero exponents */ |
---|
487 | // { |
---|
488 | // if (MATELEM(LCOM,jG,i)!=NULL) |
---|
489 | // { |
---|
490 | // tmpcoef=pGetCoeff(MATELEM(LCOM,jG,i)); |
---|
491 | // tmpint=(int)F[i]; |
---|
492 | // nPower(tmpcoef,F[i],&tmpcoef); |
---|
493 | // rescoef=nMult(rescoef,tmpcoef); |
---|
494 | // i--; |
---|
495 | // } |
---|
496 | // else |
---|
497 | // { |
---|
498 | // if (F[i]!=0) break; |
---|
499 | // } |
---|
500 | // } |
---|
501 | // if (iF==i) |
---|
502 | // /* no action took place*/ |
---|
503 | // { |
---|
504 | |
---|
505 | // } |
---|
506 | // else /* power the result up to bG */ |
---|
507 | // { |
---|
508 | // nPower(rescoef,bG,&rescoef); |
---|
509 | // /* + cleanup, post-processing */ |
---|
510 | // } |
---|
511 | |
---|
512 | int *Prv=(int*)omAlloc0((rN+1)*sizeof(int)); |
---|
513 | int *Nxt=(int*)omAlloc0((rN+1)*sizeof(int)); |
---|
514 | int *lF=(int *)omAlloc0((rN+1)*sizeof(int)); |
---|
515 | int cnt=0; int cnf=0; |
---|
516 | /* splitting F wrt jG */ |
---|
517 | for (i=1;i<=jG;i++) /* mult at the very end */ |
---|
518 | { |
---|
519 | Prv[i]=F[i]; Nxt[i]=0; |
---|
520 | if (F[i]!=0) cnf++; |
---|
521 | } |
---|
522 | if (cnf==0) freeT(Prv,rN); |
---|
523 | for (i=jG+1;i<=rN;i++) |
---|
524 | { |
---|
525 | Nxt[i]=F[i]; |
---|
526 | if (cnf!=0) { Prv[i]=0;} |
---|
527 | if (F[i]!=0) |
---|
528 | { |
---|
529 | cnt++; |
---|
530 | lF[cnt]=i; |
---|
531 | } /* eff_part,lF_for_F */ |
---|
532 | } |
---|
533 | |
---|
534 | if (cnt==1) /* Nxt consists of 1 nonzero el-t only */ |
---|
535 | { |
---|
536 | int q=lF[1]; |
---|
537 | poly Rout=pOne(); |
---|
538 | out=nc_uu_Mult_ww(q,Nxt[q],jG,bG,r); |
---|
539 | freeT(Nxt,rN); |
---|
540 | |
---|
541 | if (cnf!=0) |
---|
542 | { |
---|
543 | Prv[0]=0; |
---|
544 | p_SetExpV(Rout,Prv,r); |
---|
545 | p_Setm(Rout,r); |
---|
546 | #ifdef PDEBUG |
---|
547 | p_Test(Rout,r); |
---|
548 | #endif |
---|
549 | freeT(Prv,rN); |
---|
550 | out=nc_mm_Mult_p(Rout,out,r); /* getting the final result */ |
---|
551 | } |
---|
552 | |
---|
553 | omFreeSize((ADDRESS)lF,(rN+1)*sizeof(int)); |
---|
554 | p_Delete(&Rout,r); |
---|
555 | return (out); |
---|
556 | } |
---|
557 | /* -------------------- MAIN ACTION --------------------- */ |
---|
558 | |
---|
559 | poly D=NULL; |
---|
560 | poly Rout=NULL; |
---|
561 | number *c=(number *)omAlloc0((cnt+2)*sizeof(number)); |
---|
562 | c[cnt+1]=n_Init(1,r); |
---|
563 | i=cnt+2; /* later in freeN */ |
---|
564 | int *Op=Nxt; |
---|
565 | int *On=(int *)omAlloc0((rN+1)*sizeof(int)); |
---|
566 | int *U=(int *)omAlloc0((rN+1)*sizeof(int)); |
---|
567 | |
---|
568 | |
---|
569 | // pExpVectorCopy(U,Nxt); |
---|
570 | memcpy(U, Nxt,(rN+1)*sizeof(int)); |
---|
571 | U[jG] = U[jG] + bG; |
---|
572 | |
---|
573 | /* Op=Nxt and initial On=(0); */ |
---|
574 | Nxt=NULL; |
---|
575 | |
---|
576 | poly Pp; |
---|
577 | poly Pn; |
---|
578 | int t=0; |
---|
579 | int first=lF[1]; |
---|
580 | int nlast=lF[cnt]; |
---|
581 | int kk=0; |
---|
582 | /* cnt--; */ |
---|
583 | /* now lF[cnt] should be <=iF-1 */ |
---|
584 | |
---|
585 | while (Op[first]!=0) |
---|
586 | { |
---|
587 | t=lF[cnt]; /* cnt as it was computed */ |
---|
588 | |
---|
589 | poly w=nc_uu_Mult_ww(t,Op[t],jG,bG,r); |
---|
590 | c[cnt]=n_Copy(p_GetCoeff(w,r),r); |
---|
591 | D = pNext(w); /* getting coef and rest D */ |
---|
592 | p_DeleteLm(&w,r); |
---|
593 | w=NULL; |
---|
594 | |
---|
595 | Op[t]= 0; |
---|
596 | Pp=pOne(); |
---|
597 | p_SetExpV(Pp,Op,r); |
---|
598 | p_Setm(Pp,r); |
---|
599 | |
---|
600 | if (t<nlast) |
---|
601 | { |
---|
602 | kk=lF[cnt+1]; |
---|
603 | On[kk]=F[kk]; |
---|
604 | |
---|
605 | Pn=pOne(); |
---|
606 | p_SetExpV(Pn,On,r); |
---|
607 | p_Setm(Pn,r); |
---|
608 | |
---|
609 | if (t!=first) /* typical expr */ |
---|
610 | { |
---|
611 | w=nc_p_Mult_mm(D,Pn,r); |
---|
612 | Rout=nc_mm_Mult_p(Pp,w,r); |
---|
613 | w=NULL; |
---|
614 | } |
---|
615 | else /* last step */ |
---|
616 | { |
---|
617 | On[t]=0; |
---|
618 | p_SetExpV(Pn,On,r); |
---|
619 | p_Setm(Pn,r); |
---|
620 | Rout=nc_p_Mult_mm(D,Pn,r); |
---|
621 | } |
---|
622 | #ifdef PDEBUG |
---|
623 | p_Test(Pp,r); |
---|
624 | #endif |
---|
625 | p_Delete(&Pn,r); |
---|
626 | } |
---|
627 | else /* first step */ |
---|
628 | { |
---|
629 | Rout=nc_mm_Mult_p(Pp,D,r); |
---|
630 | } |
---|
631 | #ifdef PDEBUG |
---|
632 | p_Test(Pp,r); |
---|
633 | #endif |
---|
634 | p_Delete(&Pp,r); |
---|
635 | num=n_Mult(c[cnt+1],c[cnt],r); |
---|
636 | n_Delete(&c[cnt],r); |
---|
637 | c[cnt]=num; |
---|
638 | Rout=p_Mult_nn(Rout,c[cnt+1],r); /* Rest is ready */ |
---|
639 | out=p_Add_q(out,Rout,r); |
---|
640 | Pp=NULL; |
---|
641 | cnt--; |
---|
642 | } |
---|
643 | /* only to feel safe:*/ |
---|
644 | Pn=Pp=NULL; |
---|
645 | freeT(On,rN); |
---|
646 | freeT(Op,rN); |
---|
647 | |
---|
648 | /* leadterm and Prv-part with coef 1 */ |
---|
649 | /* U[0]=exp; */ |
---|
650 | /* U[jG]=U[jG]+bG; */ |
---|
651 | /* make leadterm */ |
---|
652 | /* ??????????? we have done it already :-0 */ |
---|
653 | Rout=pOne(); |
---|
654 | p_SetExpV(Rout,U,r); |
---|
655 | p_Setm(Rout,r); /* use again this name */ |
---|
656 | p_SetCoeff(Rout,c[cnt+1],r); /* last computed coef */ |
---|
657 | out=p_Add_q(out,Rout,r); |
---|
658 | Rout=NULL; |
---|
659 | freeT(U,rN); |
---|
660 | freeN(c,i); |
---|
661 | omFreeSize((ADDRESS)lF,(rN+1)*sizeof(int)); |
---|
662 | |
---|
663 | if (cnf!=0) |
---|
664 | { |
---|
665 | Rout=pOne(); |
---|
666 | p_SetExpV(Rout,Prv,r); |
---|
667 | p_Setm(Rout,r); |
---|
668 | freeT(Prv,rN); |
---|
669 | out=nc_mm_Mult_p(Rout,out,r); /* getting the final result */ |
---|
670 | p_Delete(&Rout,r); |
---|
671 | } |
---|
672 | return (out); |
---|
673 | } |
---|
674 | |
---|
675 | poly nc_uu_Mult_ww_vert (int i, int a, int j, int b, const ring r) |
---|
676 | { |
---|
677 | int k,m; |
---|
678 | int rN=r->N; |
---|
679 | matrix cMT=r->nc->MT[UPMATELEM(j,i,rN)]; /* cMT=current MT */ |
---|
680 | |
---|
681 | poly x=pOne();p_SetExp(x,j,1,r);p_Setm(x,r); |
---|
682 | /* var(j); */ |
---|
683 | poly y=pOne();p_SetExp(y,i,1,r);p_Setm(y,r); |
---|
684 | /*var(i); for convenience */ |
---|
685 | #ifdef PDEBUG |
---|
686 | p_Test(x,r); |
---|
687 | p_Test(y,r); |
---|
688 | #endif |
---|
689 | poly t=NULL; |
---|
690 | /* ------------ Main Cycles ----------------------------*/ |
---|
691 | |
---|
692 | for (k=2;k<=a;k++) |
---|
693 | { |
---|
694 | t = nc_p_CopyGet(MATELEM(cMT,k,1),r); |
---|
695 | |
---|
696 | if (t==NULL) /* not computed yet */ |
---|
697 | { |
---|
698 | t = nc_p_CopyGet(MATELEM(cMT,k-1,1),r); |
---|
699 | // t=p_Copy(MATELEM(cMT,k-1,1),r); |
---|
700 | t = nc_mm_Mult_p(y,t,r); |
---|
701 | MATELEM(cMT,k,1) = nc_p_CopyPut(t,r); |
---|
702 | // omCheckAddr(cMT->m); |
---|
703 | p_Delete(&t,r); |
---|
704 | } |
---|
705 | t=NULL; |
---|
706 | } |
---|
707 | |
---|
708 | for (m=2;m<=b;m++) |
---|
709 | { |
---|
710 | t = nc_p_CopyGet(MATELEM(cMT,a,m),r); |
---|
711 | // t=MATELEM(cMT,a,m); |
---|
712 | if (t==NULL) //not computed yet |
---|
713 | { |
---|
714 | t = nc_p_CopyGet(MATELEM(cMT,a,m-1),r); |
---|
715 | // t=p_Copy(MATELEM(cMT,a,m-1),r); |
---|
716 | t = nc_p_Mult_mm(t,x,r); |
---|
717 | MATELEM(cMT,a,m) = nc_p_CopyPut(t,r); |
---|
718 | // MATELEM(cMT,a,m) = t; |
---|
719 | // omCheckAddr(cMT->m); |
---|
720 | p_Delete(&t,r); |
---|
721 | } |
---|
722 | t=NULL; |
---|
723 | } |
---|
724 | p_Delete(&x,r); |
---|
725 | p_Delete(&y,r); |
---|
726 | // t=MATELEM(cMT,a,b); |
---|
727 | t= nc_p_CopyGet(MATELEM(cMT,a,b),r); |
---|
728 | // return(p_Copy(t,r)); |
---|
729 | /* since the last computed element was cMT[a,b] */ |
---|
730 | return(t); |
---|
731 | } |
---|
732 | |
---|
733 | poly nc_uu_Mult_ww (int i, int a, int j, int b, const ring r) |
---|
734 | /* (x_i)^a times (x_j)^b */ |
---|
735 | /* x_i = y, x_j = x ! */ |
---|
736 | { |
---|
737 | /* Check zero exceptions, (q-)commutativity and is there something to do? */ |
---|
738 | assume(a!=0); |
---|
739 | assume(b!=0); |
---|
740 | poly out=pOne(); |
---|
741 | if (i<=j) |
---|
742 | { |
---|
743 | p_SetExp(out,i,a,r); |
---|
744 | p_AddExp(out,j,b,r); |
---|
745 | p_Setm(out,r); |
---|
746 | return(out); |
---|
747 | }/* zero exeptions and usual case */ |
---|
748 | /* if ((a==0)||(b==0)||(i<=j)) return(out); */ |
---|
749 | |
---|
750 | if (MATELEM(r->nc->COM,j,i)!=NULL) |
---|
751 | /* commutative or quasicommutative case */ |
---|
752 | { |
---|
753 | p_SetExp(out,i,a,r); |
---|
754 | p_AddExp(out,j,b,r); |
---|
755 | p_Setm(out,r); |
---|
756 | if (r->cf->nIsOne(p_GetCoeff(MATELEM(r->nc->COM,j,i),r))) /* commutative case */ |
---|
757 | { |
---|
758 | return(out); |
---|
759 | } |
---|
760 | else |
---|
761 | { |
---|
762 | number tmp_number=p_GetCoeff(MATELEM(r->nc->COM,j,i),r); /* quasicommutative case */ |
---|
763 | nPower(tmp_number,a*b,&tmp_number); |
---|
764 | p_SetCoeff(out,tmp_number,r); |
---|
765 | return(out); |
---|
766 | } |
---|
767 | }/* end_of commutative or quasicommutative case */ |
---|
768 | p_Delete(&out,r); |
---|
769 | |
---|
770 | /* we are here if i>j and variables do not commute or quasicommute */ |
---|
771 | /* in fact, now a>=1 and b>=1; and j<i */ |
---|
772 | /* now check whether the polynomial is already computed */ |
---|
773 | int rN=r->N; |
---|
774 | int vik = UPMATELEM(j,i,rN); |
---|
775 | int cMTsize=r->nc->MTsize[vik]; |
---|
776 | int newcMTsize=0; |
---|
777 | newcMTsize=si_max(a,b); |
---|
778 | |
---|
779 | if (newcMTsize<=cMTsize) |
---|
780 | { |
---|
781 | out = nc_p_CopyGet(MATELEM(r->nc->MT[vik],a,b),r); |
---|
782 | if (out !=NULL) return (out); |
---|
783 | } |
---|
784 | int k,m; |
---|
785 | if (newcMTsize > cMTsize) |
---|
786 | { |
---|
787 | int inM=(((newcMTsize+6)/7)*7); |
---|
788 | assume (inM>=newcMTsize); |
---|
789 | newcMTsize = inM; |
---|
790 | // matrix tmp = (matrix)omAlloc0(inM*inM*sizeof(poly)); |
---|
791 | matrix tmp = mpNew(newcMTsize,newcMTsize); |
---|
792 | |
---|
793 | for (k=1;k<=cMTsize;k++) |
---|
794 | { |
---|
795 | for (m=1;m<=cMTsize;m++) |
---|
796 | { |
---|
797 | out = MATELEM(r->nc->MT[UPMATELEM(j,i,rN)],k,m); |
---|
798 | if ( out != NULL ) |
---|
799 | { |
---|
800 | MATELEM(tmp,k,m) = out;/*MATELEM(r->nc->MT[UPMATELEM(j,i,rN)],k,m)*/ |
---|
801 | // omCheckAddr(tmp->m); |
---|
802 | MATELEM(r->nc->MT[UPMATELEM(j,i,rN)],k,m)=NULL; |
---|
803 | // omCheckAddr(r->nc->MT[UPMATELEM(j,i,rN)]->m); |
---|
804 | } |
---|
805 | } |
---|
806 | } |
---|
807 | id_Delete((ideal *)&(r->nc->MT[UPMATELEM(j,i,rN)]),r); |
---|
808 | r->nc->MT[UPMATELEM(j,i,rN)] = tmp; |
---|
809 | tmp=NULL; |
---|
810 | r->nc->MTsize[UPMATELEM(j,i,rN)] = newcMTsize; |
---|
811 | } |
---|
812 | /* The update of multiplication matrix is finished */ |
---|
813 | pDelete(&out); |
---|
814 | out = nc_uu_Mult_ww_vert(i, a, j, b, r); |
---|
815 | // out = nc_uu_Mult_ww_horvert(i, a, j, b, r); |
---|
816 | return(out); |
---|
817 | } |
---|
818 | |
---|
819 | poly nc_uu_Mult_ww_horvert (int i, int a, int j, int b, const ring r) |
---|
820 | |
---|
821 | { |
---|
822 | int k,m; |
---|
823 | int rN=r->N; |
---|
824 | matrix cMT=r->nc->MT[UPMATELEM(j,i,rN)]; /* cMT=current MT */ |
---|
825 | |
---|
826 | poly x=pOne();p_SetExp(x,j,1,r);p_Setm(x,r);/* var(j); */ |
---|
827 | poly y=pOne();p_SetExp(y,i,1,r);p_Setm(y,r); /*var(i); for convenience */ |
---|
828 | #ifdef PDEBUG |
---|
829 | p_Test(x,r); |
---|
830 | p_Test(y,r); |
---|
831 | #endif |
---|
832 | |
---|
833 | poly t=NULL; |
---|
834 | |
---|
835 | int toXY; |
---|
836 | int toYX; |
---|
837 | |
---|
838 | if (a==1) /* y*x^b, b>=2 */ |
---|
839 | { |
---|
840 | toXY=b-1; |
---|
841 | while ( (MATELEM(cMT,1,toXY)==NULL) && (toXY>=2)) toXY--; |
---|
842 | for (m=toXY+1;m<=b;m++) |
---|
843 | { |
---|
844 | t=MATELEM(cMT,1,m); |
---|
845 | if (t==NULL) /* remove after debug */ |
---|
846 | { |
---|
847 | t = p_Copy(MATELEM(cMT,1,m-1),r); |
---|
848 | t = nc_p_Mult_mm(t,x,r); |
---|
849 | MATELEM(cMT,1,m) = t; |
---|
850 | /* omCheckAddr(cMT->m); */ |
---|
851 | } |
---|
852 | else |
---|
853 | { |
---|
854 | /* Error, should never get there */ |
---|
855 | WarnS("Error: a=1; MATELEM!=0"); |
---|
856 | } |
---|
857 | t=NULL; |
---|
858 | } |
---|
859 | return(p_Copy(MATELEM(cMT,1,b),r)); |
---|
860 | } |
---|
861 | |
---|
862 | if (b==1) /* y^a*x, a>=2 */ |
---|
863 | { |
---|
864 | toYX=a-1; |
---|
865 | while ( (MATELEM(cMT,toYX,1)==NULL) && (toYX>=2)) toYX--; |
---|
866 | for (m=toYX+1;m<=a;m++) |
---|
867 | { |
---|
868 | t=MATELEM(cMT,m,1); |
---|
869 | if (t==NULL) /* remove after debug */ |
---|
870 | { |
---|
871 | t = p_Copy(MATELEM(cMT,m-1,1),r); |
---|
872 | t = nc_mm_Mult_p(y,t,r); |
---|
873 | MATELEM(cMT,m,1) = t; |
---|
874 | /* omCheckAddr(cMT->m); */ |
---|
875 | } |
---|
876 | else |
---|
877 | { |
---|
878 | /* Error, should never get there */ |
---|
879 | WarnS("Error: b=1, MATELEM!=0"); |
---|
880 | } |
---|
881 | t=NULL; |
---|
882 | } |
---|
883 | return(p_Copy(MATELEM(cMT,a,1),r)); |
---|
884 | } |
---|
885 | |
---|
886 | /* ------------ Main Cycles ----------------------------*/ |
---|
887 | /* a>1, b>1 */ |
---|
888 | |
---|
889 | int dXY=0; int dYX=0; |
---|
890 | /* dXY = distance for computing x-mult, then y-mult */ |
---|
891 | /* dYX = distance for computing y-mult, then x-mult */ |
---|
892 | int toX=a-1; int toY=b-1; /* toX = to axe X, toY = to axe Y */ |
---|
893 | toXY=b-1; toYX=a-1; |
---|
894 | /* if toX==0, toXY = dist. to computed y * x^toXY */ |
---|
895 | /* if toY==0, toYX = dist. to computed y^toYX * x */ |
---|
896 | while ( (MATELEM(cMT,toX,b)==NULL) && (toX>=1)) toX--; |
---|
897 | if (toX==0) /* the whole column is not computed yet */ |
---|
898 | { |
---|
899 | while ( (MATELEM(cMT,1,toXY)==NULL) && (toXY>=1)) toXY--; |
---|
900 | /* toXY >=1 */ |
---|
901 | dXY=b-1-toXY; |
---|
902 | } |
---|
903 | dXY=dXY+a-toX; /* the distance to nearest computed y^toX x^b */ |
---|
904 | |
---|
905 | while ( (MATELEM(cMT,a,toY)==NULL) && (toY>=1)) toY--; |
---|
906 | if (toY==0) /* the whole row is not computed yet */ |
---|
907 | { |
---|
908 | while ( (MATELEM(cMT,toYX,1)==NULL) && (toYX>=1)) toYX--; |
---|
909 | /* toYX >=1 */ |
---|
910 | dYX=a-1-toYX; |
---|
911 | } |
---|
912 | dYX=dYX+b-toY; /* the distance to nearest computed y^a x^toY */ |
---|
913 | |
---|
914 | if (dYX>=dXY) |
---|
915 | { |
---|
916 | /* first x, then y */ |
---|
917 | if (toX==0) /* start with the row*/ |
---|
918 | { |
---|
919 | for (m=toXY+1;m<=b;m++) |
---|
920 | { |
---|
921 | t=MATELEM(cMT,1,m); |
---|
922 | if (t==NULL) /* remove after debug */ |
---|
923 | { |
---|
924 | t = p_Copy(MATELEM(cMT,1,m-1),r); |
---|
925 | t = nc_p_Mult_mm(t,x,r); |
---|
926 | MATELEM(cMT,1,m) = t; |
---|
927 | /* omCheckAddr(cMT->m); */ |
---|
928 | } |
---|
929 | else |
---|
930 | { |
---|
931 | /* Error, should never get there */ |
---|
932 | WarnS("dYX>=dXY,toXY; MATELEM==0"); |
---|
933 | } |
---|
934 | t=NULL; |
---|
935 | } |
---|
936 | toX=1; /* y*x^b is computed */ |
---|
937 | } |
---|
938 | /* Now toX>=1 */ |
---|
939 | for (k=toX+1;k<=a;k++) |
---|
940 | { |
---|
941 | t=MATELEM(cMT,k,b); |
---|
942 | if (t==NULL) /* remove after debug */ |
---|
943 | { |
---|
944 | t = p_Copy(MATELEM(cMT,k-1,b),r); |
---|
945 | t = nc_mm_Mult_p(y,t,r); |
---|
946 | MATELEM(cMT,k,b) = t; |
---|
947 | /* omCheckAddr(cMT->m); */ |
---|
948 | } |
---|
949 | else |
---|
950 | { |
---|
951 | /* Error, should never get there */ |
---|
952 | WarnS("dYX>=dXY,toX; MATELEM==0"); |
---|
953 | } |
---|
954 | t=NULL; |
---|
955 | } |
---|
956 | } /* endif (dYX>=dXY) */ |
---|
957 | |
---|
958 | |
---|
959 | if (dYX<dXY) |
---|
960 | { |
---|
961 | /* first y, then x */ |
---|
962 | if (toY==0) /* start with the column*/ |
---|
963 | { |
---|
964 | for (m=toYX+1;m<=a;m++) |
---|
965 | { |
---|
966 | t=MATELEM(cMT,m,1); |
---|
967 | if (t==NULL) /* remove after debug */ |
---|
968 | { |
---|
969 | t = p_Copy(MATELEM(cMT,m-1,1),r); |
---|
970 | t = nc_mm_Mult_p(y,t,r); |
---|
971 | MATELEM(cMT,m,1) = t; |
---|
972 | /* omCheckAddr(cMT->m); */ |
---|
973 | } |
---|
974 | else |
---|
975 | { |
---|
976 | /* Error, should never get there */ |
---|
977 | WarnS("dYX<dXY,toYX; MATELEM==0"); |
---|
978 | } |
---|
979 | t=NULL; |
---|
980 | } |
---|
981 | toY=1; /* y^a*x is computed */ |
---|
982 | } |
---|
983 | /* Now toY>=1 */ |
---|
984 | for (k=toY+1;k<=b;k++) |
---|
985 | { |
---|
986 | t=MATELEM(cMT,a,k); |
---|
987 | if (t==NULL) /* remove after debug */ |
---|
988 | { |
---|
989 | t = p_Copy(MATELEM(cMT,a,k-1),r); |
---|
990 | t = nc_p_Mult_mm(t,x,r); |
---|
991 | MATELEM(cMT,a,k) = t; |
---|
992 | /* omCheckAddr(cMT->m); */ |
---|
993 | } |
---|
994 | else |
---|
995 | { |
---|
996 | /* Error, should never get there */ |
---|
997 | WarnS("dYX<dXY,toY; MATELEM==0"); |
---|
998 | } |
---|
999 | t=NULL; |
---|
1000 | } |
---|
1001 | } /* endif (dYX<dXY) */ |
---|
1002 | |
---|
1003 | p_Delete(&x,r); |
---|
1004 | p_Delete(&y,r); |
---|
1005 | t=p_Copy(MATELEM(cMT,a,b),r); |
---|
1006 | return(t); /* since the last computed element was cMT[a,b] */ |
---|
1007 | } |
---|
1008 | |
---|
1009 | |
---|
1010 | /* ----------------------------- Syzygies ---------------------- */ |
---|
1011 | |
---|
1012 | /*2 |
---|
1013 | * reduction of p2 with p1 |
---|
1014 | * do not destroy p1, but p2 |
---|
1015 | * p1 divides p2 -> for use in NF algorithm |
---|
1016 | */ |
---|
1017 | |
---|
1018 | poly nc_ReduceSpoly(poly p1, poly p2,poly spNoether, const ring r) |
---|
1019 | { |
---|
1020 | if (p_GetComp(p1,r)!=p_GetComp(p2,r) |
---|
1021 | && (p_GetComp(p1,r)!=0) |
---|
1022 | && (p_GetComp(p2,r)!=0)) |
---|
1023 | { |
---|
1024 | #ifdef PDEBUG |
---|
1025 | Print("nc_ReduceSpoly: different components"); |
---|
1026 | #endif |
---|
1027 | return(NULL); |
---|
1028 | } |
---|
1029 | poly m=pOne(); |
---|
1030 | p_ExpVectorDiff(m,p2,p1,r); |
---|
1031 | p_Setm(m,r); |
---|
1032 | #ifdef PDEBUG |
---|
1033 | p_Test(m,r); |
---|
1034 | #endif |
---|
1035 | /* pSetComp(m,r)=0? */ |
---|
1036 | poly N=nc_mm_Mult_p(m,p_Head(p1,r),r); |
---|
1037 | number C=n_Copy(p_GetCoeff(N,r),r); |
---|
1038 | number cF=n_Copy(p_GetCoeff(p2,r),r); |
---|
1039 | /* GCD stuff */ |
---|
1040 | number cG = nGcd(C,cF,r); |
---|
1041 | if (!nEqual(cG,n_Init(1,r))) |
---|
1042 | { |
---|
1043 | cF = nDiv(cF,cG); |
---|
1044 | C = nDiv(C,cG); |
---|
1045 | } |
---|
1046 | p2=p_Mult_nn(p2,C,r); |
---|
1047 | poly out = nc_mm_Mult_p(m, p_Copy(pNext(p1),r), r); |
---|
1048 | N=p_Add_q(N,out,r); |
---|
1049 | number MinusOne=n_Init(-1,r); |
---|
1050 | if (!n_Equal(cF,MinusOne,r)) |
---|
1051 | { |
---|
1052 | cF=n_Neg(cF,r); |
---|
1053 | N=p_Mult_nn(N,cF,r); |
---|
1054 | } |
---|
1055 | out=p_Add_q(p2,N,r); |
---|
1056 | if (out!=NULL) pContent(out); |
---|
1057 | p_Delete(&m,r); |
---|
1058 | n_Delete(&cF,r); |
---|
1059 | n_Delete(&C,r); |
---|
1060 | n_Delete(&MinusOne,r); |
---|
1061 | return(out); |
---|
1062 | } |
---|
1063 | |
---|
1064 | |
---|
1065 | /*3 |
---|
1066 | * reduction of p2 with p1 |
---|
1067 | * do not destroy p1 and p2 |
---|
1068 | * p1 divides p2 -> for use in NF algorithm |
---|
1069 | */ |
---|
1070 | poly nc_ReduceSpolyNew(poly p1, poly p2,poly spNoether, const ring r) |
---|
1071 | { |
---|
1072 | return(nc_ReduceSpoly(p1,p_Copy(p2,r),spNoether,r)); |
---|
1073 | } |
---|
1074 | |
---|
1075 | /*4 |
---|
1076 | * creates the S-polynomial of p1 and p2 |
---|
1077 | * do not destroy p1 and p2 |
---|
1078 | */ |
---|
1079 | poly nc_CreateSpoly(poly p1, poly p2,poly spNoether, const ring r) |
---|
1080 | { |
---|
1081 | if ((p_GetComp(p1,r)!=p_GetComp(p2,r)) |
---|
1082 | && (p_GetComp(p1,r)!=0) |
---|
1083 | && (p_GetComp(p2,r)!=0)) |
---|
1084 | { |
---|
1085 | #ifdef PDEBUG |
---|
1086 | Print("nc_CreateSpoly : different components!"); |
---|
1087 | #endif |
---|
1088 | return(NULL); |
---|
1089 | } |
---|
1090 | if ((r->nc->type==nc_lie) && pHasNotCF(p1,p2)) /* prod crit */ |
---|
1091 | { |
---|
1092 | return(nc_p_Bracket_qq(pCopy(p2),p1)); |
---|
1093 | } |
---|
1094 | poly pL=pOne(); |
---|
1095 | poly m1=pOne(); |
---|
1096 | poly m2=pOne(); |
---|
1097 | pLcm(p1,p2,pL); |
---|
1098 | p_Setm(pL,r); |
---|
1099 | #ifdef PDEBUG |
---|
1100 | p_Test(pL,r); |
---|
1101 | #endif |
---|
1102 | p_ExpVectorDiff(m1,pL,p1,r); |
---|
1103 | //p_SetComp(m1,0,r); |
---|
1104 | p_Setm(m1,r); |
---|
1105 | #ifdef PDEBUG |
---|
1106 | p_Test(m1,r); |
---|
1107 | #endif |
---|
1108 | p_ExpVectorDiff(m2,pL,p2,r); |
---|
1109 | //p_SetComp(m2,0,r); |
---|
1110 | p_Setm(m2,r); |
---|
1111 | #ifdef PDEBUG |
---|
1112 | p_Test(m2,r); |
---|
1113 | #endif |
---|
1114 | p_Delete(&pL,r); |
---|
1115 | /* zero exponents ! */ |
---|
1116 | poly M1 = nc_mm_Mult_p(m1,p_Head(p1,r),r); |
---|
1117 | number C1 = n_Copy(p_GetCoeff(M1,r),r); |
---|
1118 | poly M2 = nc_mm_Mult_p(m2,p_Head(p2,r),r); |
---|
1119 | number C2 = n_Copy(p_GetCoeff(M2,r),r); |
---|
1120 | /* GCD stuff */ |
---|
1121 | number C = nGcd(C1,C2,r); |
---|
1122 | if (!nEqual(C,n_Init(1,r))) |
---|
1123 | { |
---|
1124 | C1=nDiv(C1,C); |
---|
1125 | C2=nDiv(C2,C); |
---|
1126 | } |
---|
1127 | M1=p_Mult_nn(M1,C2,r); |
---|
1128 | p_SetCoeff(m1,C2,r); |
---|
1129 | number MinusOne=n_Init(-1,r); |
---|
1130 | if (n_Equal(C1,MinusOne,r)) |
---|
1131 | { |
---|
1132 | M2=p_Add_q(M1,M2,r); |
---|
1133 | } |
---|
1134 | else |
---|
1135 | { |
---|
1136 | C1=n_Neg(C1,r); |
---|
1137 | M2=p_Mult_nn(M2,C1,r); |
---|
1138 | M2=p_Add_q(M1,M2,r); |
---|
1139 | p_SetCoeff(m2,C1,r); |
---|
1140 | } |
---|
1141 | /* M1 is killed, M2=res = C2 M1 - C1 M2 */ |
---|
1142 | poly tmp=p_Copy(p1,r); |
---|
1143 | tmp=p_LmDeleteAndNext(tmp,r); |
---|
1144 | M1=nc_mm_Mult_p(m1,tmp,r); |
---|
1145 | tmp=p_Copy(p2,r); |
---|
1146 | tmp=p_LmDeleteAndNext(tmp,r); |
---|
1147 | M2=p_Add_q(M2,M1,r); |
---|
1148 | M1=nc_mm_Mult_p(m2,tmp,r); |
---|
1149 | M2=p_Add_q(M2,M1,r); |
---|
1150 | p_Delete(&m1,r); |
---|
1151 | p_Delete(&m2,r); |
---|
1152 | // n_Delete(&C1,r); |
---|
1153 | // n_Delete(&C2,r); |
---|
1154 | n_Delete(&MinusOne,r); |
---|
1155 | #ifdef PDEBUG |
---|
1156 | p_Test(M2,r); |
---|
1157 | #endif |
---|
1158 | if (M2!=NULL) pContent(M2); |
---|
1159 | return(M2); |
---|
1160 | } |
---|
1161 | |
---|
1162 | /*5 |
---|
1163 | * reduction of tail(q) with p1 |
---|
1164 | * lead(p1) divides lead(pNext(q2)) and pNext(q2) is reduced |
---|
1165 | * do not destroy p1, but tail(q) |
---|
1166 | */ |
---|
1167 | void nc_ReduceSpolyTail(poly p1, poly q, poly q2, poly spNoether, const ring r) |
---|
1168 | { |
---|
1169 | poly a1=p_Head(p1,r); |
---|
1170 | poly Q=pNext(q2); |
---|
1171 | number cQ=p_GetCoeff(Q,r); |
---|
1172 | poly m=pOne(); |
---|
1173 | p_ExpVectorDiff(m,Q,p1,r); |
---|
1174 | // p_SetComp(m,0,r); |
---|
1175 | p_Setm(m,r); |
---|
1176 | #ifdef PDEBUG |
---|
1177 | p_Test(m,r); |
---|
1178 | #endif |
---|
1179 | /* pSetComp(m,r)=0? */ |
---|
1180 | poly M=nc_mm_Mult_p(m,p_Copy(p1,r),r); |
---|
1181 | number C=p_GetCoeff(M,r); |
---|
1182 | M=p_Add_q(M,nc_mm_Mult_p(m,p_LmDeleteAndNext(p_Copy(p1,r),r),r),r); |
---|
1183 | q=p_Mult_nn(q,C,r); |
---|
1184 | number MinusOne=n_Init(-1,r); |
---|
1185 | if (!n_Equal(cQ,MinusOne,r)) |
---|
1186 | { |
---|
1187 | cQ=nNeg(cQ); |
---|
1188 | M=p_Mult_nn(M,cQ,r); |
---|
1189 | } |
---|
1190 | Q=p_Add_q(Q,M,r); |
---|
1191 | pNext(q2)=Q; |
---|
1192 | |
---|
1193 | p_Delete(&m,r); |
---|
1194 | n_Delete(&C,r); |
---|
1195 | n_Delete(&cQ,r); |
---|
1196 | n_Delete(&MinusOne,r); |
---|
1197 | /* return(q); */ |
---|
1198 | } |
---|
1199 | |
---|
1200 | /*6 |
---|
1201 | * creates the commutative lcm(lm(p1),lm(p2)) |
---|
1202 | * do not destroy p1 and p2 |
---|
1203 | */ |
---|
1204 | poly nc_CreateShortSpoly(poly p1, poly p2, const ring r) |
---|
1205 | { |
---|
1206 | if (p_GetComp(p1,r)!=p_GetComp(p2,r)) |
---|
1207 | { |
---|
1208 | Print("spShort:exponent mismatch!"); |
---|
1209 | return(NULL); |
---|
1210 | } |
---|
1211 | poly m=pOne(); |
---|
1212 | pLcm(p1,p2,m); |
---|
1213 | p_Setm(m,r); |
---|
1214 | #ifdef PDEBUG |
---|
1215 | p_Test(m,r); |
---|
1216 | #endif |
---|
1217 | return(m); |
---|
1218 | } |
---|
1219 | |
---|
1220 | void nc_kBucketPolyRed(kBucket_pt b, poly p, number *c) |
---|
1221 | { |
---|
1222 | // b will not by multiplied by any constant in this impl. |
---|
1223 | // ==> *c=1 |
---|
1224 | *c=nInit(1); |
---|
1225 | poly m=pOne(); |
---|
1226 | pExpVectorDiff(m,kBucketGetLm(b),p); |
---|
1227 | pSetm(m); |
---|
1228 | #ifdef PDEBUG |
---|
1229 | pTest(m); |
---|
1230 | #endif |
---|
1231 | poly pp=nc_mm_Mult_p(m,pCopy(p),currRing); |
---|
1232 | pDelete(&m); |
---|
1233 | number n=nCopy(pGetCoeff(pp)); |
---|
1234 | number MinusOne=nInit(-1); |
---|
1235 | number nn; |
---|
1236 | if (!nEqual(n,MinusOne)) |
---|
1237 | { |
---|
1238 | nn=nNeg(nInvers(n)); |
---|
1239 | } |
---|
1240 | else nn=nInit(1); |
---|
1241 | nDelete(&n); |
---|
1242 | n=nMult(nn,pGetCoeff(kBucketGetLm(b))); |
---|
1243 | nDelete(&nn); |
---|
1244 | pp=p_Mult_nn(pp,n,currRing); |
---|
1245 | nDelete(&n); |
---|
1246 | nDelete(&MinusOne); |
---|
1247 | int l=pLength(pp); |
---|
1248 | kBucket_Add_q(b,pp,&l); |
---|
1249 | } |
---|
1250 | |
---|
1251 | void nc_PolyPolyRed(poly &b, poly p, number *c) |
---|
1252 | // reduces b with p, do not delete both |
---|
1253 | { |
---|
1254 | // b will not by multiplied by any constant in this impl. |
---|
1255 | // ==> *c=1 |
---|
1256 | *c=nInit(1); |
---|
1257 | poly m=pOne(); |
---|
1258 | pExpVectorDiff(m,pHead(b),p); |
---|
1259 | pSetm(m); |
---|
1260 | #ifdef PDEBUG |
---|
1261 | pTest(m); |
---|
1262 | #endif |
---|
1263 | poly pp=nc_mm_Mult_p(m,pCopy(p),currRing); |
---|
1264 | pDelete(&m); |
---|
1265 | number n=nCopy(pGetCoeff(pp)); |
---|
1266 | number MinusOne=nInit(-1); |
---|
1267 | number nn; |
---|
1268 | if (!nEqual(n,MinusOne)) |
---|
1269 | { |
---|
1270 | nn=nNeg(nInvers(n)); |
---|
1271 | } |
---|
1272 | else nn=nInit(1); |
---|
1273 | nDelete(&n); |
---|
1274 | n=nMult(nn,pGetCoeff(b)); |
---|
1275 | nDelete(&nn); |
---|
1276 | pp=p_Mult_nn(pp,n,currRing); |
---|
1277 | nDelete(&n); |
---|
1278 | nDelete(&MinusOne); |
---|
1279 | b=p_Add_q(b,pp,currRing); |
---|
1280 | } |
---|
1281 | |
---|
1282 | poly nc_p_Bracket_qq(poly p, poly q) |
---|
1283 | /* returns [p,q], destroys p */ |
---|
1284 | { |
---|
1285 | if (!rIsPluralRing(currRing)) return(NULL); |
---|
1286 | if (pComparePolys(p,q)) return(NULL); |
---|
1287 | /* Components !? */ |
---|
1288 | poly Q=NULL; |
---|
1289 | number coef=NULL; |
---|
1290 | poly res=NULL; |
---|
1291 | poly pres=NULL; |
---|
1292 | int UseBuckets=1; |
---|
1293 | if ((pLength(p)< MIN_LENGTH_BUCKET/2) && (pLength(q)< MIN_LENGTH_BUCKET/2) || TEST_OPT_NOT_BUCKETS) UseBuckets=0; |
---|
1294 | sBucket_pt bu_out; |
---|
1295 | if (UseBuckets) bu_out=sBucketCreate(currRing); |
---|
1296 | while (p!=NULL) |
---|
1297 | { |
---|
1298 | Q=q; |
---|
1299 | while(Q!=NULL) |
---|
1300 | { |
---|
1301 | pres=nc_mm_Bracket_nn(p,Q); /* since no coeffs are taken into account there */ |
---|
1302 | if (pres!=NULL) |
---|
1303 | { |
---|
1304 | coef=nMult(pGetCoeff(p),pGetCoeff(Q)); |
---|
1305 | pres=p_Mult_nn(pres,coef,currRing); |
---|
1306 | if (UseBuckets) sBucket_Add_p(bu_out,pres,pLength(pres)); |
---|
1307 | else res=p_Add_q(res,pres,currRing); |
---|
1308 | nDelete(&coef); |
---|
1309 | } |
---|
1310 | pIter(Q); |
---|
1311 | } |
---|
1312 | p=pLmDeleteAndNext(p); |
---|
1313 | } |
---|
1314 | if (UseBuckets) |
---|
1315 | { |
---|
1316 | res = NULL; |
---|
1317 | int len = pLength(res); |
---|
1318 | sBucketDestroyAdd(bu_out, &res, &len); |
---|
1319 | } |
---|
1320 | return(res); |
---|
1321 | } |
---|
1322 | |
---|
1323 | poly nc_mm_Bracket_nn(poly m1, poly m2) |
---|
1324 | /*returns [m1,m2] for two monoms, destroys nothing */ |
---|
1325 | /* without coeffs */ |
---|
1326 | { |
---|
1327 | if (pLmIsConstant(m1) || pLmIsConstant(m1)) return(NULL); |
---|
1328 | if (pLmCmp(m1,m2)==0) return(NULL); |
---|
1329 | int rN=currRing->N; |
---|
1330 | int *M1=(int *)omAlloc0((rN+1)*sizeof(int)); |
---|
1331 | int *M2=(int *)omAlloc0((rN+1)*sizeof(int)); |
---|
1332 | int *PREFIX=(int *)omAlloc0((rN+1)*sizeof(int)); |
---|
1333 | int *SUFFIX=(int *)omAlloc0((rN+1)*sizeof(int)); |
---|
1334 | pGetExpV(m1,M1); |
---|
1335 | pGetExpV(m2,M2); |
---|
1336 | poly res=NULL; |
---|
1337 | poly ares=NULL; |
---|
1338 | poly bres=NULL; |
---|
1339 | poly prefix=NULL; |
---|
1340 | poly suffix=NULL; |
---|
1341 | int nMin,nMax; |
---|
1342 | number nTmp=NULL; |
---|
1343 | int i,j,k; |
---|
1344 | for (i=1;i<=rN;i++) |
---|
1345 | { |
---|
1346 | if (M2[i]!=0) |
---|
1347 | { |
---|
1348 | ares=NULL; |
---|
1349 | for (j=1;j<=rN;j++) |
---|
1350 | { |
---|
1351 | if (M1[j]!=0) |
---|
1352 | { |
---|
1353 | bres=NULL; |
---|
1354 | /* compute [ x_j^M1[j],x_i^M2[i] ] */ |
---|
1355 | if (i<j) {nMax=j; nMin=i;} else {nMax=i; nMin=j;} |
---|
1356 | if ( (i==j) || ((MATELEM(currRing->nc->COM,nMin,nMax)!=NULL) && nIsOne(pGetCoeff(MATELEM(currRing->nc->C,nMin,nMax))) )) /* not (the same exp. or commuting exps)*/ |
---|
1357 | { bres=NULL; } |
---|
1358 | else |
---|
1359 | { |
---|
1360 | if (i<j) { bres=nc_uu_Mult_ww(j,M1[j],i,M2[i],currRing); } |
---|
1361 | else bres=nc_uu_Mult_ww(i,M2[i],j,M1[j],currRing); |
---|
1362 | if (nIsOne(pGetCoeff(bres))) |
---|
1363 | { |
---|
1364 | bres=pLmDeleteAndNext(bres); |
---|
1365 | } |
---|
1366 | else |
---|
1367 | { |
---|
1368 | nTmp=nSub(pGetCoeff(bres),nInit(1)); |
---|
1369 | pSetCoeff(bres,nTmp); /* only lc ! */ |
---|
1370 | } |
---|
1371 | #ifdef PDEBUG |
---|
1372 | pTest(bres); |
---|
1373 | #endif |
---|
1374 | if (i>j) bres=p_Neg(bres, currRing); |
---|
1375 | } |
---|
1376 | if (bres!=NULL) |
---|
1377 | { |
---|
1378 | /* now mult (prefix, bres, suffix) */ |
---|
1379 | memcpy(SUFFIX, M1,(rN+1)*sizeof(int)); |
---|
1380 | memcpy(PREFIX, M1,(rN+1)*sizeof(int)); |
---|
1381 | for (k=1;k<=j;k++) SUFFIX[k]=0; |
---|
1382 | for (k=j;k<=rN;k++) PREFIX[k]=0; |
---|
1383 | SUFFIX[0]=0; |
---|
1384 | PREFIX[0]=0; |
---|
1385 | prefix=pOne(); |
---|
1386 | suffix=pOne(); |
---|
1387 | pSetExpV(prefix,PREFIX); |
---|
1388 | pSetm(prefix); |
---|
1389 | pSetExpV(suffix,SUFFIX); |
---|
1390 | pSetm(suffix); |
---|
1391 | if (!pLmIsConstant(prefix)) bres = nc_mm_Mult_p(prefix, bres,currRing); |
---|
1392 | if (!pLmIsConstant(suffix)) bres = nc_p_Mult_mm(bres, suffix,currRing); |
---|
1393 | ares=p_Add_q(ares, bres,currRing); |
---|
1394 | /* What to give free? */ |
---|
1395 | /* Do we have to free PREFIX/SUFFIX? it seems so */ |
---|
1396 | pDelete(&prefix); |
---|
1397 | pDelete(&suffix); |
---|
1398 | } |
---|
1399 | } |
---|
1400 | } |
---|
1401 | if (ares!=NULL) |
---|
1402 | { |
---|
1403 | /* now mult (prefix, bres, suffix) */ |
---|
1404 | memcpy(SUFFIX, M2,(rN+1)*sizeof(int)); |
---|
1405 | memcpy(PREFIX, M2,(rN+1)*sizeof(int)); |
---|
1406 | for (k=1;k<=i;k++) SUFFIX[k]=0; |
---|
1407 | for (k=i;k<=rN;k++) PREFIX[k]=0; |
---|
1408 | SUFFIX[0]=0; |
---|
1409 | PREFIX[0]=0; |
---|
1410 | prefix=pOne(); |
---|
1411 | suffix=pOne(); |
---|
1412 | pSetExpV(prefix,PREFIX); |
---|
1413 | pSetm(prefix); |
---|
1414 | pSetExpV(suffix,SUFFIX); |
---|
1415 | pSetm(suffix); |
---|
1416 | bres=ares; |
---|
1417 | if (!pLmIsConstant(prefix)) bres = nc_mm_Mult_p(prefix, bres,currRing); |
---|
1418 | if (!pLmIsConstant(suffix)) bres = nc_p_Mult_mm(bres, suffix,currRing); |
---|
1419 | res=p_Add_q(res, bres,currRing); |
---|
1420 | pDelete(&prefix); |
---|
1421 | pDelete(&suffix); |
---|
1422 | } |
---|
1423 | } |
---|
1424 | } |
---|
1425 | freeT(M1, rN); |
---|
1426 | freeT(M2, rN); |
---|
1427 | freeT(PREFIX, rN); |
---|
1428 | freeT(SUFFIX, rN); |
---|
1429 | return(res); |
---|
1430 | } |
---|
1431 | |
---|
1432 | ideal twostd(ideal I) |
---|
1433 | { |
---|
1434 | int i; |
---|
1435 | int j; |
---|
1436 | int s; |
---|
1437 | int flag; |
---|
1438 | poly p=NULL; |
---|
1439 | poly q=NULL; |
---|
1440 | ideal J=kStd(I, currQuotient,testHomog,NULL,NULL,0,0,NULL); |
---|
1441 | idSkipZeroes(J); |
---|
1442 | ideal K=NULL; |
---|
1443 | poly varj=NULL; |
---|
1444 | ideal Q=NULL; |
---|
1445 | ideal id_tmp=NULL; |
---|
1446 | int rN=currRing->N; |
---|
1447 | int iSize=0; |
---|
1448 | loop |
---|
1449 | { |
---|
1450 | flag=0; |
---|
1451 | K=NULL; |
---|
1452 | s=idElem(J); |
---|
1453 | for (i=0;i<=s-1;i++) |
---|
1454 | { |
---|
1455 | p=J->m[i]; |
---|
1456 | for (j=1;j<=rN;j++) |
---|
1457 | { |
---|
1458 | varj = pOne(); |
---|
1459 | pSetExp(varj,j,1); |
---|
1460 | pSetm(varj); |
---|
1461 | q = nc_p_Mult_mm(pCopy(p),varj,currRing); |
---|
1462 | pDelete(&varj); |
---|
1463 | q = nc_ReduceSpoly(p,q,NULL,currRing); |
---|
1464 | q = kNF(J,currQuotient,q,0,0); |
---|
1465 | if (q!=NULL) |
---|
1466 | { |
---|
1467 | if (pIsConstant(q)) |
---|
1468 | { |
---|
1469 | Q=idInit(1,1); |
---|
1470 | Q->m[0]=pOne(); |
---|
1471 | idDelete(&J); |
---|
1472 | pDelete(&q); |
---|
1473 | if (K!=NULL) idDelete(&K); |
---|
1474 | return(Q); |
---|
1475 | } |
---|
1476 | flag=1; |
---|
1477 | Q=idInit(1,1); |
---|
1478 | Q->m[0]=q; |
---|
1479 | id_tmp=idSimpleAdd(K,Q); |
---|
1480 | idDelete(&K); |
---|
1481 | K=id_tmp; |
---|
1482 | idDelete(&Q); |
---|
1483 | } |
---|
1484 | } |
---|
1485 | } |
---|
1486 | if (flag==0) |
---|
1487 | /* i.e. all elements are two-sided */ |
---|
1488 | { |
---|
1489 | idDelete(&K); |
---|
1490 | return(J); |
---|
1491 | } |
---|
1492 | /* now we update GrBasis J with K */ |
---|
1493 | iSize=IDELEMS(J); |
---|
1494 | id_tmp=idSimpleAdd(J,K); |
---|
1495 | idDelete(&K); |
---|
1496 | idDelete(&J); |
---|
1497 | BITSET save_test=test; |
---|
1498 | test|=Sy_bit(OPT_SB_1); |
---|
1499 | J=kStd(id_tmp, currQuotient, testHomog,NULL,NULL,0,iSize); |
---|
1500 | test=save_test; |
---|
1501 | idSkipZeroes(J); |
---|
1502 | } |
---|
1503 | } |
---|
1504 | |
---|
1505 | matrix nc_PrintMat(int a, int b, ring r, int metric) |
---|
1506 | /* returns matrix with the info on noncomm multiplication */ |
---|
1507 | { |
---|
1508 | |
---|
1509 | if ( (a==b) || !rIsPluralRing(r) ) return(NULL); |
---|
1510 | int i; |
---|
1511 | int j; |
---|
1512 | if (a>b) {j=b; i=a;} |
---|
1513 | else {j=a; i=b;} |
---|
1514 | /* i<j */ |
---|
1515 | int rN=r->N; |
---|
1516 | int size=r->nc->MTsize[UPMATELEM(i,j,rN)]; |
---|
1517 | matrix M = r->nc->MT[UPMATELEM(i,j,rN)]; |
---|
1518 | /* return(M); */ |
---|
1519 | int sizeofres; |
---|
1520 | if (metric==0) |
---|
1521 | { |
---|
1522 | sizeofres=sizeof(int); |
---|
1523 | } |
---|
1524 | if (metric==1) |
---|
1525 | { |
---|
1526 | sizeofres=sizeof(number); |
---|
1527 | } |
---|
1528 | matrix res=mpNew(size,size); |
---|
1529 | int s; |
---|
1530 | int t; |
---|
1531 | int length; |
---|
1532 | long totdeg; |
---|
1533 | poly p; |
---|
1534 | for(s=1;s<=size;s++) |
---|
1535 | { |
---|
1536 | for(t=1;t<=size;t++) |
---|
1537 | { |
---|
1538 | p=MATELEM(M,s,t); |
---|
1539 | if (p==NULL) |
---|
1540 | { |
---|
1541 | MATELEM(res,s,t)=0; |
---|
1542 | } |
---|
1543 | else |
---|
1544 | { |
---|
1545 | length = pLength(p); |
---|
1546 | if (metric==0) /* length */ |
---|
1547 | { |
---|
1548 | MATELEM(res,s,t)= p_ISet(length,r); |
---|
1549 | } |
---|
1550 | else if (metric==1) /* sum of deg divided by the length */ |
---|
1551 | { |
---|
1552 | totdeg=0; |
---|
1553 | while (p!=NULL) |
---|
1554 | { |
---|
1555 | totdeg=totdeg+pDeg(p,r); |
---|
1556 | pIter(p); |
---|
1557 | } |
---|
1558 | number ntd = nInit(totdeg); |
---|
1559 | number nln = nInit(length); |
---|
1560 | number nres=nDiv(ntd,nln); |
---|
1561 | nDelete(&ntd); |
---|
1562 | nDelete(&nln); |
---|
1563 | MATELEM(res,s,t)=p_NSet(nres,r); |
---|
1564 | } |
---|
1565 | } |
---|
1566 | } |
---|
1567 | } |
---|
1568 | return(res); |
---|
1569 | } |
---|
1570 | |
---|
1571 | void ncKill(ring r) |
---|
1572 | /* kills the nc extension of ring r */ |
---|
1573 | { |
---|
1574 | int i,j; |
---|
1575 | int rN=r->N; |
---|
1576 | for(i=1;i<rN;i++) |
---|
1577 | { |
---|
1578 | for(j=i+1;j<=rN;j++) |
---|
1579 | { |
---|
1580 | id_Delete((ideal *)&(r->nc->MT[UPMATELEM(i,j,rN)]),r->nc->basering); |
---|
1581 | } |
---|
1582 | } |
---|
1583 | omFreeSize((ADDRESS)r->nc->MT,rN*(rN-1)/2*sizeof(matrix)); |
---|
1584 | omFreeSize((ADDRESS)r->nc->MTsize,rN*(rN-1)/2*sizeof(int)); |
---|
1585 | id_Delete((ideal *)&(r->nc->C),r->nc->basering); |
---|
1586 | id_Delete((ideal *)&(r->nc->D),r->nc->basering); |
---|
1587 | id_Delete((ideal *)&(r->nc->COM),r->nc->basering); |
---|
1588 | omFreeSize((ADDRESS)r->nc,sizeof(nc_struct)); |
---|
1589 | r->nc=NULL; |
---|
1590 | } |
---|
1591 | |
---|
1592 | void ncCleanUp(ring r) |
---|
1593 | { |
---|
1594 | /* small CleanUp of r->nc */ |
---|
1595 | omFreeSize((ADDRESS)r->nc,sizeof(nc_struct)); |
---|
1596 | r->nc = NULL; |
---|
1597 | } |
---|
1598 | |
---|
1599 | poly nc_p_CopyGet(poly a, ring r) |
---|
1600 | /* for use in getting the mult. martix elements*/ |
---|
1601 | { |
---|
1602 | if (!rIsPluralRing(r)) return(p_Copy(a,r)); |
---|
1603 | if (r==r->nc->basering) return(p_Copy(a,r)); |
---|
1604 | else |
---|
1605 | { |
---|
1606 | return(prCopyR_NoSort(a,r->nc->basering,r)); |
---|
1607 | } |
---|
1608 | } |
---|
1609 | |
---|
1610 | poly nc_p_CopyPut(poly a, ring r) |
---|
1611 | /* for use in defining the mult. martix elements*/ |
---|
1612 | { |
---|
1613 | if (!rIsPluralRing(r)) return(p_Copy(a,r)); |
---|
1614 | if (r==r->nc->basering) return(p_Copy(a,r)); |
---|
1615 | else |
---|
1616 | { |
---|
1617 | return(prCopyR_NoSort(a,r,r->nc->basering)); |
---|
1618 | } |
---|
1619 | } |
---|
1620 | |
---|
1621 | int nc_CheckSubalgebra(poly PolyVar, ring r) |
---|
1622 | /* returns TRUE if product of vars from PolyVar defines */ |
---|
1623 | /* an admissible subalgebra of r */ |
---|
1624 | { |
---|
1625 | int rN=r->N; |
---|
1626 | int *ExpVar=(int*)omAlloc0((rN+1)*sizeof(int)); |
---|
1627 | int *ExpTmp=(int*)omAlloc0((rN+1)*sizeof(int)); |
---|
1628 | p_GetExpV(PolyVar, ExpVar, r); |
---|
1629 | int i; int j; int k; |
---|
1630 | poly test=NULL; |
---|
1631 | int OK=1; |
---|
1632 | for (i=1;i<rN;i++) |
---|
1633 | { |
---|
1634 | if (ExpVar[i]==0) /* i.e. not in PolyVar */ |
---|
1635 | { |
---|
1636 | for (j=i+1;j<=rN;j++) |
---|
1637 | { |
---|
1638 | if (ExpVar[j]==0) |
---|
1639 | { |
---|
1640 | test=nc_p_CopyGet(MATELEM(r->nc->D,i,j),r); |
---|
1641 | while (test!=NULL) |
---|
1642 | { |
---|
1643 | p_GetExpV(test, ExpTmp, r); |
---|
1644 | OK=1; |
---|
1645 | for (k=1;k<=rN;k++) |
---|
1646 | { |
---|
1647 | if (ExpTmp[k]!=0) |
---|
1648 | { |
---|
1649 | if (ExpVar[k]!=0) OK=0; |
---|
1650 | } |
---|
1651 | } |
---|
1652 | if (!OK) return(FALSE); |
---|
1653 | pIter(test); |
---|
1654 | } |
---|
1655 | } |
---|
1656 | } |
---|
1657 | } |
---|
1658 | } |
---|
1659 | p_Delete(&test,r); |
---|
1660 | freeT(ExpVar,rN); |
---|
1661 | freeT(ExpTmp,rN); |
---|
1662 | return(TRUE); |
---|
1663 | } |
---|
1664 | |
---|
1665 | // int Commutative_Context(ring r, leftv expression) |
---|
1666 | // /* returns 1 if expression consists */ |
---|
1667 | // /* of commutative elements */ |
---|
1668 | // { |
---|
1669 | // /* crucial: poly -> ideal, module, matrix */ |
---|
1670 | |
---|
1671 | // } |
---|
1672 | |
---|
1673 | // int Comm_Context_Poly(ring r, poly p) |
---|
1674 | // { |
---|
1675 | // poly COMM=r->nc->COMM; |
---|
1676 | // poly pp=pOne(); |
---|
1677 | // memset(pp->exp,0,r->ExpL_Size*sizeof(long)); |
---|
1678 | // while (p!=NULL) |
---|
1679 | // { |
---|
1680 | // for (i=0;i<=r->ExpL_Size;i++) |
---|
1681 | // { |
---|
1682 | // if ((p->exp[i]) && (pp->exp[i])) return(FALSE); |
---|
1683 | // /* nonzero exponent of non-comm variable */ |
---|
1684 | // } |
---|
1685 | // pIter(p); |
---|
1686 | // } |
---|
1687 | // return(TRUE); |
---|
1688 | // } |
---|
1689 | |
---|
1690 | BOOLEAN nc_CallPlural(matrix CCC, matrix DDD, poly CCN, poly DDN, ring r) |
---|
1691 | /* returns TRUE if there were errors */ |
---|
1692 | /* analyze inputs, check them for consistency */ |
---|
1693 | /* detect nc_type, DO NOT initialize multiplication */ |
---|
1694 | /* check the ordering condition and evtl. NDC */ |
---|
1695 | { |
---|
1696 | matrix CC = NULL; |
---|
1697 | matrix DD = NULL; |
---|
1698 | poly CN = NULL; |
---|
1699 | poly DN = NULL; |
---|
1700 | matrix C; |
---|
1701 | matrix D; |
---|
1702 | number nN,pN,qN; |
---|
1703 | int tmpIsSkewConstant; |
---|
1704 | int i,j; |
---|
1705 | if (r->nc != NULL) |
---|
1706 | { |
---|
1707 | WarnS("redefining algebra structure"); |
---|
1708 | if (r->nc->ref>1) /* in use by somebody else */ |
---|
1709 | { |
---|
1710 | r->nc->ref--; |
---|
1711 | } |
---|
1712 | else /* kill the previous nc data */ |
---|
1713 | { |
---|
1714 | ncKill(r); |
---|
1715 | } |
---|
1716 | } |
---|
1717 | r->nc = (nc_struct *)omAlloc0(sizeof(nc_struct)); |
---|
1718 | r->nc->ref = 1; |
---|
1719 | r->nc->basering = r; |
---|
1720 | r->nc->type = nc_undef; |
---|
1721 | |
---|
1722 | /* initialition of the matrix C */ |
---|
1723 | /* check the correctness of arguments */ |
---|
1724 | |
---|
1725 | if ((CCC != NULL) && ( (MATCOLS(CCC)==1) || MATROWS(CCC)==1 ) ) |
---|
1726 | { |
---|
1727 | CN = MATELEM(CCC,1,1); |
---|
1728 | } |
---|
1729 | else |
---|
1730 | { |
---|
1731 | if ((CCC != NULL) && ( (MATCOLS(CCC)!=r->N) || (MATROWS(CCC)!=r->N) )) |
---|
1732 | { |
---|
1733 | Werror("Square %d x %d matrix expected",r->N,r->N); |
---|
1734 | ncCleanUp(r); |
---|
1735 | return TRUE; |
---|
1736 | } |
---|
1737 | } |
---|
1738 | if (( CCC != NULL) && (CC == NULL)) CC = mpCopy(CCC); |
---|
1739 | if (( CCN != NULL) && (CN == NULL)) CN = CCN; |
---|
1740 | |
---|
1741 | /* initialition of the matrix D */ |
---|
1742 | /* check the correctness of arguments */ |
---|
1743 | |
---|
1744 | if ((DDD != NULL) && ( (MATCOLS(DDD)==1) || MATROWS(DDD)==1 ) ) |
---|
1745 | { |
---|
1746 | DN = MATELEM(DDD,1,1); |
---|
1747 | } |
---|
1748 | else |
---|
1749 | { |
---|
1750 | if ((DDD != NULL) && ( (MATCOLS(DDD)!=r->N) || (MATROWS(DDD)!=r->N) )) |
---|
1751 | { |
---|
1752 | Werror("Square %d x %d matrix expected",r->N,r->N); |
---|
1753 | ncCleanUp(r); |
---|
1754 | return TRUE; |
---|
1755 | } |
---|
1756 | } |
---|
1757 | if (( DDD != NULL) && (DD == NULL)) DD = mpCopy(DDD); |
---|
1758 | if (( DDN != NULL) && (DN == NULL)) DN = DDN; |
---|
1759 | |
---|
1760 | /* further checks */ |
---|
1761 | |
---|
1762 | if (CN != NULL) /* create matrix C = CN * Id */ |
---|
1763 | { |
---|
1764 | nN = p_GetCoeff(CN,r); |
---|
1765 | if (n_IsZero(nN,r)) |
---|
1766 | { |
---|
1767 | Werror("Incorrect input : zero coefficients are not allowed"); |
---|
1768 | ncCleanUp(r); |
---|
1769 | return TRUE; |
---|
1770 | } |
---|
1771 | if (nIsOne(nN)) |
---|
1772 | { |
---|
1773 | r->nc->type = nc_lie; |
---|
1774 | } |
---|
1775 | else |
---|
1776 | { |
---|
1777 | r->nc->type = nc_general; |
---|
1778 | } |
---|
1779 | r->nc->IsSkewConstant = 1; |
---|
1780 | C = mpNew(r->N,r->N); |
---|
1781 | for(i=1; i<r->N; i++) |
---|
1782 | { |
---|
1783 | for(j=i+1; j<=r->N; j++) |
---|
1784 | { |
---|
1785 | MATELEM(C,i,j) = nc_p_CopyPut(CN,r); |
---|
1786 | } |
---|
1787 | } |
---|
1788 | } |
---|
1789 | if ( (CN == NULL) && (CC != NULL) ) /* copy matrix C */ |
---|
1790 | { |
---|
1791 | C = mpCopy(CC); |
---|
1792 | /* analyze C */ |
---|
1793 | pN = p_GetCoeff(MATELEM(C,1,2),r); |
---|
1794 | tmpIsSkewConstant = 1; |
---|
1795 | for(i=1; i<r->N; i++) |
---|
1796 | { |
---|
1797 | for(j=i+1; j<=r->N; j++) |
---|
1798 | { |
---|
1799 | qN = p_GetCoeff(MATELEM(C,i,j),r); |
---|
1800 | if ( qN == NULL ) /* check the consistency: Cij!=0 */ |
---|
1801 | { |
---|
1802 | Werror("Incorrect input : matrix of coefficients contains zeros in the upper triangle"); |
---|
1803 | ncCleanUp(r); |
---|
1804 | return TRUE; |
---|
1805 | } |
---|
1806 | if (!nEqual(pN,qN)) tmpIsSkewConstant = 0; |
---|
1807 | } |
---|
1808 | } |
---|
1809 | r->nc->IsSkewConstant=tmpIsSkewConstant; |
---|
1810 | if ( (tmpIsSkewConstant) && (nIsOne(pN)) ) |
---|
1811 | { |
---|
1812 | r->nc->type = nc_lie; |
---|
1813 | } |
---|
1814 | else |
---|
1815 | { |
---|
1816 | r->nc->type = nc_general; |
---|
1817 | } |
---|
1818 | } |
---|
1819 | |
---|
1820 | /* initialition of the matrix D */ |
---|
1821 | if ( DD == NULL ) |
---|
1822 | /* we treat DN only (it could also be NULL) */ |
---|
1823 | { |
---|
1824 | D = mpNew(r->N,r->N); |
---|
1825 | if (DN == NULL) |
---|
1826 | { |
---|
1827 | if ( (currRing->nc->type == nc_lie) || (currRing->nc->type == nc_undef) ) |
---|
1828 | { |
---|
1829 | currRing->nc->type = nc_comm; /* it was nc_skew earlier */ |
---|
1830 | } |
---|
1831 | else /* nc_general, nc_skew */ |
---|
1832 | { |
---|
1833 | currRing->nc->type = nc_skew; |
---|
1834 | } |
---|
1835 | } |
---|
1836 | else /* DN != NULL */ |
---|
1837 | { |
---|
1838 | for(i=1; i<r->N; i++) |
---|
1839 | { |
---|
1840 | for(j=i+1; j<=r->N; j++) |
---|
1841 | { |
---|
1842 | MATELEM(D,i,j) = nc_p_CopyPut(DN,r); |
---|
1843 | } |
---|
1844 | } |
---|
1845 | } |
---|
1846 | } |
---|
1847 | else /* DD != NULL */ |
---|
1848 | { |
---|
1849 | D = mpCopy(DD); |
---|
1850 | } |
---|
1851 | /* analyze D */ |
---|
1852 | /* check the ordering condition for D (both matrix and poly cases) */ |
---|
1853 | poly p,q; |
---|
1854 | int report = 1; |
---|
1855 | for(i=1; i<r->N; i++) |
---|
1856 | { |
---|
1857 | for(j=i+1; j<=r->N; j++) |
---|
1858 | { |
---|
1859 | p = MATELEM(D,i,j); |
---|
1860 | if ( p != NULL) |
---|
1861 | { |
---|
1862 | q = pOne(); |
---|
1863 | p_SetExp(q,i,1,r); |
---|
1864 | p_SetExp(q,j,1,r); |
---|
1865 | p_Setm(q,r); |
---|
1866 | if (p_LmCmp(q,p,r) != 1) /* i.e. lm(p)<=lm(q) */ |
---|
1867 | { |
---|
1868 | Print("Bad ordering at %d,%d",i,j); |
---|
1869 | report = 0; |
---|
1870 | } |
---|
1871 | p_Delete(&q,r); |
---|
1872 | p = NULL; |
---|
1873 | } |
---|
1874 | } |
---|
1875 | } |
---|
1876 | if (!report) |
---|
1877 | { |
---|
1878 | Werror("Matrix of polynomials violates the ordering condition"); |
---|
1879 | ncCleanUp(r); |
---|
1880 | return TRUE; |
---|
1881 | } |
---|
1882 | r->nc->C = C; |
---|
1883 | r->nc->D = D; |
---|
1884 | return nc_InitMultiplication(r); |
---|
1885 | } |
---|
1886 | |
---|
1887 | BOOLEAN nc_InitMultiplication(ring r) |
---|
1888 | { |
---|
1889 | /* returns TRUE if there were errors */ |
---|
1890 | /* initialize the multiplication */ |
---|
1891 | int i,j; |
---|
1892 | matrix COM; |
---|
1893 | r->nc->MT = (matrix *)omAlloc0(r->N*(r->N-1)/2*sizeof(matrix)); |
---|
1894 | r->nc->MTsize = (int *)omAlloc0(r->N*(r->N-1)/2*sizeof(int)); |
---|
1895 | COM = mpCopy(r->nc->C); |
---|
1896 | poly p; |
---|
1897 | short DefMTsize=7; |
---|
1898 | int IsNonComm=0; |
---|
1899 | int tmpIsSkewConstant; |
---|
1900 | |
---|
1901 | for(i=1; i<r->N; i++) |
---|
1902 | { |
---|
1903 | for(j=i+1; j<=r->N; j++) |
---|
1904 | { |
---|
1905 | if ( MATELEM(r->nc->D,i,j) == NULL ) /* quasicommutative case */ |
---|
1906 | { |
---|
1907 | /* 1x1 mult.matrix */ |
---|
1908 | r->nc->MTsize[UPMATELEM(i,j,r->N)] = 1; |
---|
1909 | r->nc->MT[UPMATELEM(i,j,r->N)] = mpNew(1,1); |
---|
1910 | } |
---|
1911 | else /* pure noncommutative case */ |
---|
1912 | { |
---|
1913 | /* TODO check the special multiplication properties */ |
---|
1914 | IsNonComm = 1; |
---|
1915 | MATELEM(COM,i,j) = NULL; |
---|
1916 | r->nc->MTsize[UPMATELEM(i,j,r->N)] = DefMTsize; /* default sizes */ |
---|
1917 | r->nc->MT[UPMATELEM(i,j,r->N)] = mpNew(DefMTsize, DefMTsize); |
---|
1918 | } |
---|
1919 | /* set MT[i,j,1,1] to c_i_j*x_i*x_j + D_i_j */ |
---|
1920 | p = pOne(); |
---|
1921 | p_SetCoeff(p,nCopy(pGetCoeff(MATELEM(r->nc->C,i,j))),r); |
---|
1922 | p_SetExp(p,i,1,r); |
---|
1923 | p_SetExp(p,j,1,r); |
---|
1924 | p_Setm(p,r); |
---|
1925 | p = p_Add_q(p, nc_p_CopyGet(MATELEM(r->nc->D,i,j),r),r); |
---|
1926 | MATELEM(r->nc->MT[UPMATELEM(i,j,r->N)],1,1) = nc_p_CopyPut(p,r); |
---|
1927 | pDelete(&p); |
---|
1928 | p = NULL; |
---|
1929 | } |
---|
1930 | } |
---|
1931 | if (r->nc->type==nc_undef) |
---|
1932 | { |
---|
1933 | if (IsNonComm==1) |
---|
1934 | { |
---|
1935 | // assume(pN!=NULL); |
---|
1936 | // if ((tmpIsSkewConstant==1) && (nIsOne(pGetCoeff(pN)))) r->nc->type=nc_lie; |
---|
1937 | // else r->nc->type=nc_general; |
---|
1938 | } |
---|
1939 | if (IsNonComm==0) |
---|
1940 | { |
---|
1941 | r->nc->type=nc_skew; /* TODO: check whether it is commutative */ |
---|
1942 | r->nc->IsSkewConstant=tmpIsSkewConstant; |
---|
1943 | } |
---|
1944 | } |
---|
1945 | r->nc->COM=COM; |
---|
1946 | return FALSE; |
---|
1947 | } |
---|
1948 | |
---|
1949 | /* substitute the n-th variable by e in p |
---|
1950 | * destroy p |
---|
1951 | * e is not a constant |
---|
1952 | */ |
---|
1953 | poly nc_pSubst(poly p, int n, poly e) |
---|
1954 | { |
---|
1955 | int rN=currRing->N; |
---|
1956 | int *PRE = (int *)omAlloc0((rN+1)*sizeof(int)); |
---|
1957 | int *SUF = (int *)omAlloc0((rN+1)*sizeof(int)); |
---|
1958 | int i,j,pow; |
---|
1959 | poly suf,pre; |
---|
1960 | poly res = NULL; |
---|
1961 | poly out = NULL; |
---|
1962 | while ( p!= NULL ) |
---|
1963 | { |
---|
1964 | pGetExpV(p, PRE); /* faster splitting? */ |
---|
1965 | pow = PRE[n]; PRE[n]=0; |
---|
1966 | res = NULL; |
---|
1967 | if (pow!=0) |
---|
1968 | { |
---|
1969 | for (i=n+1; i<=rN; i++) |
---|
1970 | { |
---|
1971 | SUF[i] = PRE[i]; |
---|
1972 | PRE[i] = 0; |
---|
1973 | } |
---|
1974 | res = pPower(pCopy(e),pow); |
---|
1975 | /* multiply with prefix */ |
---|
1976 | pre = pOne(); |
---|
1977 | pSetExpV(pre,PRE); |
---|
1978 | pSetm(pre); |
---|
1979 | pSetComp(pre,PRE[0]); |
---|
1980 | res = nc_mm_Mult_p(pre,res,currRing); |
---|
1981 | /* multiply with suffix */ |
---|
1982 | suf = pOne(); |
---|
1983 | pSetExpV(suf,SUF); |
---|
1984 | pSetm(suf); |
---|
1985 | pSetComp(suf,PRE[0]); |
---|
1986 | res = nc_p_Mult_mm(res,suf,currRing); |
---|
1987 | } |
---|
1988 | else /* pow==0 */ |
---|
1989 | { |
---|
1990 | res = pHead(p); |
---|
1991 | } |
---|
1992 | p = pLmDeleteAndNext(p); |
---|
1993 | out = pAdd(out,res); |
---|
1994 | } |
---|
1995 | freeT(PRE,rN); |
---|
1996 | freeT(SUF,rN); |
---|
1997 | return(out); |
---|
1998 | } |
---|
1999 | |
---|
2000 | #endif |
---|