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 *******************************************************************/
10
11 #define MYTEST 0
12 #define OUTPUT 0
13
14 #if MYTEST
15 #define OM_CHECK 4
16 #define OM_TRACK 5
17 #endif
18
19
20
21
22 #include "misc/auxiliary.h"
23
24 #ifdef HAVE_PLURAL
25
26 # define PLURAL_INTERNAL_DECLARATIONS
27 #include "nc.h"
28 #include "sca.h"
29 #include "gb_hack.h"
30
31 #include "polys/monomials/ring.h"
32
33 #include "coeffs/numbers.h"
34
35 #include "misc/options.h"
36
37 #include "polys/monomials/ring.h"
38 #include "polys/monomials/p_polys.h"
39
40 #include "polys/simpleideals.h"
41 #include "polys/matpol.h"
42
43 #include "polys/kbuckets.h"
44 #include "polys/sbuckets.h"
45
46 #include "polys/prCopy.h"
47
48 #include "polys/operations/p_Mult_q.h"
49
50 #include "summator.h"
51
52 #include "ncSAMult.h" // for CMultiplier etc classes
53 #include "ncSAFormula.h" // for CFormulaPowerMultiplier and enum Enum_ncSAType
54
55 // #ifdef HAVE_RATGRING
56 // #include "polys/ratgring.h"
57 // #endif
58
NF_Proc_Dummy(ideal,ideal,poly,int,int,const ring)59 static poly NF_Proc_Dummy(ideal, ideal, poly, int, int, const ring)
60 { WerrorS("nc_NF not defined"); return NULL; }
BBA_Proc_Dummy(const ideal,const ideal,const intvec *,const intvec *,kStrategy,const ring)61 static ideal BBA_Proc_Dummy (const ideal, const ideal, const intvec *, const intvec *, kStrategy, const ring)
62 { WerrorS("nc_NF not defined"); return NULL; }
63
64 // the following funtion poiters are quasi-static:
65 // they will be set in siInit and never changes afterwards:
66 VAR NF_Proc nc_NF=NF_Proc_Dummy;
67 VAR BBA_Proc gnc_gr_bba=BBA_Proc_Dummy;
68 VAR BBA_Proc gnc_gr_mora=BBA_Proc_Dummy;
69 VAR BBA_Proc sca_bba=BBA_Proc_Dummy;
70 VAR BBA_Proc sca_mora=BBA_Proc_Dummy;
71 VAR BBA_Proc sca_gr_bba=BBA_Proc_Dummy;
72
73 /* copy : */
74 poly nc_p_CopyGet(poly a, const ring r);
75 poly nc_p_CopyPut(poly a, const ring r);
76
77 poly nc_p_Bracket_qq(poly p, const poly q, const ring r);
78
79 // only SCA can be used by default, formulas are off by default
80 VAR int iNCExtensions = SCAMASK | NOFORMULAMASK;
81
getNCExtensions()82 int& getNCExtensions()
83 {
84 return (iNCExtensions);
85 }
86
setNCExtensions(int iMask)87 int setNCExtensions(int iMask)
88 {
89 const int iOld = getNCExtensions();
90 getNCExtensions() = iMask;
91 return (iOld);
92 }
93
ncExtensions(int iMask)94 bool ncExtensions(int iMask) // = 0x0FFFF
95 {
96 return ((getNCExtensions() & iMask) == iMask);
97 }
98
99 /* global nc_macros : */
100
101 #define freeT(A,v) omFreeSize((ADDRESS)A,(v+1)*sizeof(int))
102 #define freeN(A,k) omFreeSize((ADDRESS)A,k*sizeof(number))
103
104
105 // some forward declarations:
106
107
108 // polynomial multiplication functions for p_Procs :
109 poly gnc_pp_Mult_mm(const poly p, const poly m, const ring r, poly &last);
110 poly gnc_p_Mult_mm(poly p, const poly m, const ring r);
111 poly gnc_p_mm_Mult(poly m, const poly p, const ring r);
112 poly gnc_pp_mm_Mult(const poly p, const poly m, const ring r);
113
114
115 /* syzygies : */
116 poly gnc_CreateSpolyOld(const poly p1, const poly p2/*, poly spNoether*/, const ring r);
117 poly gnc_ReduceSpolyOld(const poly p1, poly p2/*, poly spNoether*/, const ring r);
118
119 poly gnc_CreateSpolyNew(const poly p1, const poly p2/*, poly spNoether*/, const ring r);
120 poly gnc_ReduceSpolyNew(const poly p1, poly p2/*, poly spNoether*/, const ring r);
121
122
123
124 void gnc_kBucketPolyRedNew(kBucket_pt b, poly p, number *c);
125 void gnc_kBucketPolyRed_ZNew(kBucket_pt b, poly p, number *c);
126
127 void gnc_kBucketPolyRedOld(kBucket_pt b, poly p, number *c);
128 void gnc_kBucketPolyRed_ZOld(kBucket_pt b, poly p, number *c);
129
130
131 // poly gnc_ReduceSpolyNew(poly p1, poly p2, poly spNoether, const ring r);
132 // void gnc_ReduceSpolyTail(poly p1, poly q, poly q2, poly spNoether, const ring r);
133
134 // void nc_kBucketPolyRed(kBucket_pt b, poly p);
135
136 void nc_CleanUp(nc_struct* p); // just free memory!
137 void nc_rCleanUp(ring r); // smaller than kill: just free mem
138
139
140 #if 0
141 // deprecated functions:
142 // poly gnc_p_Minus_mm_Mult_qq_ign(poly p, const poly m, poly q, int & d1, poly d2, const ring ri, poly &d3);
143 // poly gnc_p_Minus_mm_Mult_qq(poly p, const poly m, poly q, const ring r);
144 // poly nc_p_Minus_mm_Mult_qq(poly p, const poly m, const poly q, int &lp, int lq, const ring r);
145 // poly nc_p_Plus_mm_Mult_qq (poly p, const poly m, const poly q, int &lp, int lq, const ring r);
146 #endif
147
148
149 ///////////////////////////////////////////////////////////////////////////////
nc_p_Minus_mm_Mult_qq(poly p,const poly m,const poly q,int & shorter,const poly,const ring r)150 poly nc_p_Minus_mm_Mult_qq(poly p, const poly m, const poly q, int &shorter,
151 const poly, const ring r)
152 {
153 poly mc = p_Neg( p_Copy(m, r), r );
154 poly mmc = nc_mm_Mult_pp( mc, q, r );
155 p_Delete(&mc, r);
156
157 int org_p=pLength(p);
158 int org_q=pLength(q);
159
160 p = p_Add_q(p, mmc, r);
161
162 shorter = pLength(p)-org_p-org_q; // ring independent!
163
164 return(p);
165 }
166
167 // returns p + m*q destroys p, const: q, m
nc_p_Plus_mm_Mult_qq(poly p,const poly m,const poly q,int & lp,const int,const ring r)168 poly nc_p_Plus_mm_Mult_qq(poly p, const poly m, const poly q, int &lp,
169 const int, const ring r)
170 {
171 p = p_Add_q(p, nc_mm_Mult_pp( m, q, r ), r);
172
173 lp = pLength(p);
174
175 return(p);
176 }
177
178 #if 0
179 poly gnc_p_Minus_mm_Mult_qq_ign(poly p, const poly m, poly q, int & d1, poly d2, const ring r, poly &d3)
180 {
181 poly t;
182 int i;
183
184 return gnc_p_Minus_mm_Mult_qq(p, m, q, d1, i, t, r);
185 }
186 #endif
187
188
189 //----------- auxiliary routines--------------------------
_gnc_p_Mult_q(poly p,poly q,const int copy,const ring r)190 poly _gnc_p_Mult_q(poly p, poly q, const int copy, const ring r) // not used anymore!
191 /* destroy p,q unless copy=1 */
192 {
193 poly res=NULL;
194 poly qq,pp;
195 if (copy)
196 {
197 qq=p_Copy(q,r);
198 pp=p_Copy(p,r);
199 }
200 else
201 {
202 qq=q;
203 pp=p;
204 }
205 while (qq!=NULL)
206 {
207 res=p_Add_q(res, pp_Mult_mm(pp, qq, r), r); // p_Head(qq, r)?
208 qq=p_LmDeleteAndNext(qq,r);
209 }
210 p_Delete(&pp,r);
211 return(res);
212 }
213
214 // return pPolyP * pPolyQ; destroy or reuse pPolyP and pPolyQ
_nc_p_Mult_q(poly pPolyP,poly pPolyQ,const ring rRing)215 poly _nc_p_Mult_q(poly pPolyP, poly pPolyQ, const ring rRing)
216 {
217 assume( rIsNCRing(rRing) );
218 #ifdef PDEBUG
219 p_Test(pPolyP, rRing);
220 p_Test(pPolyQ, rRing);
221 #endif
222 #ifdef RDEBUG
223 rTest(rRing);
224 #endif
225
226 int lp, lq;
227
228 pqLength(pPolyP, pPolyQ, lp, lq, MIN_LENGTH_BUCKET);
229
230 bool bUsePolynomial = TEST_OPT_NOT_BUCKETS || (si_max(lp, lq) < MIN_LENGTH_BUCKET); // ???
231
232 CPolynomialSummator sum(rRing, bUsePolynomial);
233
234 if (lq <= lp) // ?
235 {
236 // always length(q) times "p * q[j]"
237 for( ; pPolyQ!=NULL; pPolyQ = p_LmDeleteAndNext( pPolyQ, rRing ) )
238 sum += pp_Mult_mm( pPolyP, pPolyQ, rRing);
239
240 p_Delete( &pPolyP, rRing );
241 } else
242 {
243 // always length(p) times "p[i] * q"
244 for( ; pPolyP!=NULL; pPolyP = p_LmDeleteAndNext( pPolyP, rRing ) )
245 sum += nc_mm_Mult_pp( pPolyP, pPolyQ, rRing);
246
247 p_Delete( &pPolyQ, rRing );
248 }
249
250 return(sum);
251 }
252
253 // return pPolyP * pPolyQ; preserve pPolyP and pPolyQ
_nc_pp_Mult_qq(const poly pPolyP,const poly pPolyQ,const ring rRing)254 poly _nc_pp_Mult_qq(const poly pPolyP, const poly pPolyQ, const ring rRing)
255 {
256 assume( rIsNCRing(rRing) );
257 #ifdef PDEBUG
258 p_Test(pPolyP, rRing);
259 p_Test(pPolyQ, rRing);
260 #endif
261 #ifdef RDEBUG
262 rTest(rRing);
263 #endif
264
265 int lp, lq;
266
267 pqLength(pPolyP, pPolyQ, lp, lq, MIN_LENGTH_BUCKET);
268
269 bool bUsePolynomial = TEST_OPT_NOT_BUCKETS || (si_max(lp, lq) < MIN_LENGTH_BUCKET); // ???
270
271 CPolynomialSummator sum(rRing, bUsePolynomial);
272
273 if (lq <= lp) // ?
274 {
275 // always length(q) times "p * q[j]"
276 for( poly q = pPolyQ; q !=NULL; q = pNext(q) )
277 sum += pp_Mult_mm(pPolyP, q, rRing);
278 } else
279 {
280 // always length(p) times "p[i] * q"
281 for( poly p = pPolyP; p !=NULL; p = pNext(p) )
282 sum += nc_mm_Mult_pp( p, pPolyQ, rRing);
283 }
284
285 return(sum);
286 }
287
288
289
290 poly gnc_mm_Mult_nn (int *F, int *G, const ring r);
291 poly gnc_mm_Mult_uu (int *F,int jG,int bG, const ring r);
292
293 /* #define nc_uu_Mult_ww nc_uu_Mult_ww_vert */
294 poly gnc_uu_Mult_ww (int i, int a, int j, int b, const ring r);
295 /* poly nc_uu_Mult_ww_vert (int i, int a, int j, int b, const ring r); */
296 /* poly nc_uu_Mult_ww_horvert (int i, int a, int j, int b, const ring r); */
297 /* poly nc_uu_Mult_ww_hvdiag (int i, int a, int j, int b, const ring r); */
298 /* not written yet */
299
300
gnc_p_Mult_mm_Common(poly p,const poly m,int side,const ring r)301 poly gnc_p_Mult_mm_Common(poly p, const poly m, int side, const ring r)
302 /* p is poly, m is mono with coeff, destroys p */
303 /* if side==1, computes p_Mult_mm; otherwise, mm_Mult_p */
304 {
305 if ((p==NULL) || (m==NULL)) return NULL;
306 /* if (pNext(p)==NULL) return(nc_mm_Mult_nn(p,pCopy(m),r)); */
307 /* excluded - the cycle will do it anyway - OK. */
308 if (p_IsConstant(m,r)) return(__p_Mult_nn(p,p_GetCoeff(m,r),r));
309
310 #ifdef PDEBUG
311 p_Test(p,r);
312 p_Test(m,r);
313 #endif
314 poly v=NULL;
315 int rN=r->N;
316 int *P=(int *)omAlloc0((rN+1)*sizeof(int));
317 int *M=(int *)omAlloc0((rN+1)*sizeof(int));
318 /* coefficients: */
319 number cP,cM,cOut;
320 p_GetExpV(m, M, r);
321 cM=p_GetCoeff(m,r);
322 /* components:*/
323 const int expM=p_GetComp(m,r);
324 int expP=0;
325 int expOut=0;
326 /* bucket constraints: */
327 int UseBuckets=1;
328 if (pLength(p)< MIN_LENGTH_BUCKET || TEST_OPT_NOT_BUCKETS) UseBuckets=0;
329
330 CPolynomialSummator sum(r, UseBuckets == 0);
331
332 while (p!=NULL)
333 {
334 #ifdef PDEBUG
335 p_Test(p,r);
336 #endif
337 expP=p_GetComp(p,r);
338 if (expP==0)
339 {
340 expOut=expM;
341 }
342 else
343 {
344 if (expM==0)
345 {
346 expOut=expP;
347 #ifdef PDEBUG
348 // if (side)
349 // {
350 // PrintS("gnc_p_Mult_mm: Multiplication in the left module from the right");
351 // }
352 #endif
353 }
354 else
355 {
356 /* REPORT_ERROR */
357 #ifdef PDEBUG
358 const char* s;
359 if (side==1) s="gnc_p_Mult_mm";
360 else s="gnc_p_mm_Mult";
361 Print("%s: exponent mismatch %d and %d\n",s,expP,expM);
362 #endif
363 expOut=0;
364 }
365 }
366 p_GetExpV(p,P,r);
367 cP=pGetCoeff(p);
368 cOut=n_Mult(cP,cM,r->cf);
369 if (side==1)
370 {
371 v = gnc_mm_Mult_nn(P, M, r);
372 }
373 else
374 {
375 v = gnc_mm_Mult_nn(M, P, r);
376 }
377 v = __p_Mult_nn(v,cOut,r);
378 n_Delete(&cOut,r->cf);
379 p_SetCompP(v,expOut,r);
380
381 sum += v;
382
383 p_LmDelete(&p,r);
384 }
385 freeT(P,rN);
386 freeT(M,rN);
387
388 return(sum);
389 }
390
391 /* poly functions defined in p_Procs : */
gnc_pp_Mult_mm(const poly p,const poly m,const ring r)392 poly gnc_pp_Mult_mm(const poly p, const poly m, const ring r)
393 {
394 return( gnc_p_Mult_mm_Common(p_Copy(p,r), m, 1, r) );
395 }
396
gnc_p_Mult_mm(poly p,const poly m,const ring r)397 poly gnc_p_Mult_mm(poly p, const poly m, const ring r)
398 {
399 return( gnc_p_Mult_mm_Common(p, m, 1, r) );
400 }
401
402 /* m * p */
gnc_p_mm_Mult(poly p,const poly m,const ring r)403 poly gnc_p_mm_Mult(poly p, const poly m, const ring r)
404 {
405 return( gnc_p_Mult_mm_Common(p, m, 0, r) );
406 }
407
gnc_pp_mm_Mult(const poly p,const poly m,const ring r)408 poly gnc_pp_mm_Mult(const poly p, const poly m, const ring r)
409 {
410 return( gnc_p_Mult_mm_Common(p_Copy(p,r), m, 0, r) );
411 }
412
413
414
gnc_mm_Mult_nn(int * F0,int * G0,const ring r)415 poly gnc_mm_Mult_nn(int *F0, int *G0, const ring r)
416 /* destroys nothing, no coeffs and exps */
417 {
418 poly out=NULL;
419 int i,j;
420 int iF,jG,iG;
421 int rN=r->N;
422
423 int *F=(int *)omAlloc0((rN+1)*sizeof(int));
424 int *G=(int *)omAlloc0((rN+1)*sizeof(int));
425
426 memcpy(F, F0,(rN+1)*sizeof(int));
427 // pExpVectorCopy(F,F0);
428 memcpy(G, G0,(rN+1)*sizeof(int));
429 // pExpVectorCopy(G,G0);
430 F[0]=0;
431 G[0]=0;
432
433 iF=rN;
434 while ((F[iF]==0)&&(iF>=1)) iF--; /* last exp_num of F */
435 if (iF==0) /* F0 is zero vector */
436 {
437 out=p_One(r);
438 p_SetExpV(out,G0,r);
439 p_Setm(out,r);
440 freeT(F,rN);
441 freeT(G,rN);
442 return(out);
443 }
444 jG=1;
445 while ((G[jG]==0)&&(jG<rN)) jG++; /* first exp_num of G */
446 iG=rN;
447 while ((G[iG]==0)&&(iG>1)) iG--; /* last exp_num of G */
448
449 out=p_One(r);
450
451 if (iF<=jG)
452 /* i.e. no mixed exp_num , MERGE case */
453 {
454 { for(int ii=rN;ii>0;ii--) F[ii]+=G[ii]; }
455 p_SetExpV(out,F,r);
456 p_Setm(out,r);
457 freeT(F,rN);
458 freeT(G,rN);
459 return(out);
460 }
461
462 number cff=n_Init(1,r->cf);
463 number tmp_num=NULL;
464 int cpower=0;
465
466 if (ncRingType(r)==nc_skew)
467 {
468 if (r->GetNC()->IsSkewConstant==1)
469 {
470 int tpower=0;
471 for(j=jG; j<=iG; j++)
472 {
473 if (G[j]!=0)
474 {
475 cpower = 0;
476 for(i=j+1; i<=iF; i++)
477 {
478 cpower = cpower + F[i];
479 }
480 cpower = cpower*G[j]; // bug! here may happen an arithmetic overflow!!!
481 tpower = tpower + cpower;
482 }
483 }
484 cff = n_Copy(pGetCoeff(MATELEM(r->GetNC()->COM,1,2)),r->cf);
485 n_Power(cff,tpower,&tmp_num, r->cf);
486 n_Delete(&cff,r->cf);
487 cff = tmp_num;
488 }
489 else /* skew commutative with nonequal coeffs */
490 {
491 number totcff=n_Init(1,r->cf);
492 for(j=jG; j<=iG; j++)
493 {
494 if (G[j]!=0)
495 {
496 cpower = 0;
497 for(i=j+1; i<=iF; i++)
498 {
499 if (F[i]!=0)
500 {
501 cpower = F[i]*G[j]; // bug! overflow danger!!!
502 cff = n_Copy(pGetCoeff(MATELEM(r->GetNC()->COM,j,i)),r->cf);
503 n_Power(cff,cpower,&tmp_num, r->cf);
504 cff = n_Mult(totcff,tmp_num, r->cf);
505 n_Delete(&totcff, r->cf);
506 n_Delete(&tmp_num, r->cf);
507 totcff = n_Copy(cff,r->cf);
508 n_Delete(&cff,r->cf);
509 }
510 } /* end 2nd for */
511 }
512 }
513 cff=totcff;
514 }
515 { for(int ii=rN;ii>0;ii--) F[ii]+=G[ii]; }
516 p_SetExpV(out,F,r);
517 p_Setm(out,r);
518 p_SetCoeff(out,cff,r);
519 freeT(F,rN);
520 freeT(G,rN);
521 return(out);
522 } /* end nc_skew */
523
524 /* now we have to destroy out! */
525 p_Delete(&out,r);
526
527 if (iG==jG)
528 /* g is univariate monomial */
529 {
530 /* if (ri->GetNC()->type==nc_skew) -- postpone to TU */
531 out = gnc_mm_Mult_uu(F,jG,G[jG],r);
532 freeT(F,rN);
533 freeT(G,rN);
534 return(out);
535 }
536
537 int *Prv=(int *)omAlloc0((rN+1)*sizeof(int));
538 int *Nxt=(int *)omAlloc0((rN+1)*sizeof(int));
539
540 int *log=(int *)omAlloc0((rN+1)*sizeof(int));
541 int cnt=0; int cnf=0;
542
543 /* splitting F wrt jG */
544 for (i=1;i<=jG;i++)
545 {
546 Prv[i]=F[i]; Nxt[i]=0; /* mult at the very end */
547 if (F[i]!=0) cnf++;
548 }
549
550 if (cnf==0) freeT(Prv,rN);
551
552 for (i=jG+1;i<=rN;i++)
553 {
554 Nxt[i]=F[i];
555 /* if (cnf!=0) Prv[i]=0; */
556 if (F[i]!=0)
557 {
558 cnt++;
559 } /* effective part for F */
560 }
561 freeT(F,rN);
562 cnt=0;
563
564 for (i=1;i<=rN;i++)
565 {
566 if (G[i]!=0)
567 {
568 cnt++;
569 log[cnt]=i;
570 } /* lG for G */
571 }
572
573 /* ---------------------- A C T I O N ------------------------ */
574 poly D=NULL;
575 poly Rout=NULL;
576 number *c=(number *)omAlloc0((rN+1)*sizeof(number));
577 c[0]=n_Init(1,r->cf);
578
579 int *Op=Nxt;
580 int *On=G;
581 int *U=(int *)omAlloc0((rN+1)*sizeof(int));
582
583 for (i=jG;i<=rN;i++) U[i]=Nxt[i]+G[i]; /* make leadterm */
584 Nxt=NULL;
585 G=NULL;
586 cnt=1;
587 int t=0;
588 poly w=NULL;
589 poly Pn=p_One(r);
590 p_SetExpV(Pn,On,r);
591 p_Setm(Pn,r);
592
593 while (On[iG]!=0)
594 {
595 t=log[cnt];
596
597 w=gnc_mm_Mult_uu(Op,t,On[t],r);
598 c[cnt]=n_Mult(c[cnt-1],pGetCoeff(w),r->cf);
599 D = pNext(w); /* getting coef and rest D */
600 p_LmDelete(&w,r);
601 w=NULL;
602
603 Op[t] += On[t]; /* update exp_vectors */
604 On[t] = 0;
605
606 if (t!=iG) /* not the last step */
607 {
608 p_SetExpV(Pn,On,r);
609 p_Setm(Pn,r);
610 #ifdef PDEBUG
611 p_Test(Pn,r);
612 #endif
613
614 // if (pNext(D)==0)
615 // is D a monomial? could be postponed higher
616 // {
617 // Rout=nc_mm_Mult_nn(D,Pn,r);
618 // }
619 // else
620 // {
621 Rout=gnc_p_Mult_mm(D,Pn,r);
622 // }
623 }
624 else
625 {
626 Rout=D;
627 D=NULL;
628 }
629
630 if (Rout!=NULL)
631 {
632 Rout=__p_Mult_nn(Rout,c[cnt-1],r); /* Rest is ready */
633 out=p_Add_q(out,Rout,r);
634 Rout=NULL;
635 }
636 cnt++;
637 }
638 freeT(On,rN);
639 freeT(Op,rN);
640 p_Delete(&Pn,r);
641 omFreeSize((ADDRESS)log,(rN+1)*sizeof(int));
642
643 /* leadterm and Prv-part */
644
645 Rout=p_One(r);
646 /* U is lead.monomial */
647 U[0]=0;
648 p_SetExpV(Rout,U,r);
649 p_Setm(Rout,r); /* use again this name Rout */
650 #ifdef PDEBUG
651 p_Test(Rout,r);
652 #endif
653 p_SetCoeff(Rout,c[cnt-1],r);
654 out=p_Add_q(out,Rout,r);
655 freeT(U,rN);
656 freeN(c,rN+1);
657 if (cnf!=0) /* Prv is non-zero vector */
658 {
659 Rout=p_One(r);
660 Prv[0]=0;
661 p_SetExpV(Rout,Prv,r);
662 p_Setm(Rout,r);
663 #ifdef PDEBUG
664 p_Test(Rout,r);
665 #endif
666 out=gnc_p_mm_Mult(out,Rout,r); /* getting the final result */
667 freeT(Prv,rN);
668 p_Delete(&Rout,r);
669 }
670 return (out);
671 }
672
673
gnc_mm_Mult_uu(int * F,int jG,int bG,const ring r)674 poly gnc_mm_Mult_uu(int *F,int jG,int bG, const ring r)
675 /* f=mono(F),g=(x_iG)^bG */
676 {
677 poly out=NULL;
678 int i;
679 number num=NULL;
680
681 int rN=r->N;
682 int iF=r->N;
683 while ((F[iF]==0)&&(iF>0)) iF-- ; /* last exponent_num of F */
684
685 if (iF==0) /* F==zero vector in other words */
686 {
687 out=p_One(r);
688 p_SetExp(out,jG,bG,r);
689 p_Setm(out,r);
690 return(out);
691 }
692
693 int jF=1;
694 while ((F[jF]==0)&&(jF<=rN)) jF++; /* first exp of F */
695
696 if (iF<=jG) /* i.e. no mixed exp_num */
697 {
698 out=p_One(r);
699 F[jG]=F[jG]+bG;
700 p_SetExpV(out,F,r);
701 p_Setm(out,r);
702 return(out);
703 }
704
705 if (iF==jF) /* uni times uni */
706 {
707 out=gnc_uu_Mult_ww(iF,F[iF],jG,bG,r);
708 return(out);
709 }
710
711 /* Now: F is mono with >=2 exponents, jG<iF */
712 /* check the quasi-commutative case */
713 // matrix LCOM=r->GetNC()->COM;
714 // number rescoef=n_Init(1,r);
715 // number tmpcoef=n_Init(1,r);
716 // int tmpint;
717 // i=iF;
718 // while (i>=jG+1)
719 // /* all the non-zero exponents */
720 // {
721 // if (MATELEM(LCOM,jG,i)!=NULL)
722 // {
723 // tmpcoef=pGetCoeff(MATELEM(LCOM,jG,i));
724 // tmpint=(int)F[i];
725 // nPower(tmpcoef,F[i],&tmpcoef);
726 // rescoef=nMult(rescoef,tmpcoef);
727 // i--;
728 // }
729 // else
730 // {
731 // if (F[i]!=0) break;
732 // }
733 // }
734 // if (iF==i)
735 // /* no action took place*/
736 // {
737
738 // }
739 // else /* power the result up to bG */
740 // {
741 // nPower(rescoef,bG,&rescoef);
742 // /* + cleanup, post-processing */
743 // }
744
745 int *Prv=(int*)omAlloc0((rN+1)*sizeof(int));
746 int *Nxt=(int*)omAlloc0((rN+1)*sizeof(int));
747 int *lF=(int *)omAlloc0((rN+1)*sizeof(int));
748
749 int cnt=0; int cnf=0;
750 /* splitting F wrt jG */
751 for (i=1;i<=jG;i++) /* mult at the very end */
752 {
753 Prv[i]=F[i]; Nxt[i]=0;
754 if (F[i]!=0) cnf++;
755 }
756
757 if (cnf==0)
758 {
759 freeT(Prv,rN); Prv = NULL;
760 }
761
762 for (i=jG+1;i<=rN;i++)
763 {
764 Nxt[i]=F[i];
765 if (cnf!=0) { Prv[i]=0;}
766 if (F[i]!=0)
767 {
768 cnt++;
769 lF[cnt]=i;
770 } /* eff_part,lF_for_F */
771 }
772
773 if (cnt==1) /* Nxt consists of 1 nonzero el-t only */
774 {
775 int q=lF[1];
776 poly Rout=p_One(r);
777 out=gnc_uu_Mult_ww(q,Nxt[q],jG,bG,r);
778
779 freeT(Nxt,rN); Nxt = NULL;
780
781 if (cnf!=0)
782 {
783 Prv[0]=0;
784 p_SetExpV(Rout,Prv,r);
785 p_Setm(Rout,r);
786
787 #ifdef PDEBUG
788 p_Test(Rout,r);
789 #endif
790
791 freeT(Prv,rN);
792 Prv = NULL;
793
794 out=gnc_p_mm_Mult(out,Rout,r); /* getting the final result */
795 }
796
797 freeT(lF,rN);
798 lF = NULL;
799
800 p_Delete(&Rout,r);
801
802 assume(Nxt == NULL);
803 assume(lF == NULL);
804 assume(Prv == NULL);
805
806 return (out);
807 }
808 /* -------------------- MAIN ACTION --------------------- */
809
810 poly D=NULL;
811 poly Rout=NULL;
812 number *c=(number *)omAlloc0((cnt+2)*sizeof(number));
813 c[cnt+1]=n_Init(1,r->cf);
814 i=cnt+2; /* later in freeN */
815 int *Op=Nxt;
816
817 int *On=(int *)omAlloc0((rN+1)*sizeof(int));
818 int *U=(int *)omAlloc0((rN+1)*sizeof(int));
819
820
821 // pExpVectorCopy(U,Nxt);
822 memcpy(U, Nxt,(rN+1)*sizeof(int));
823 U[jG] = U[jG] + bG;
824
825 /* Op=Nxt and initial On=(0); */
826 Nxt=NULL;
827
828 poly Pp;
829 poly Pn;
830 int t=0;
831 int first=lF[1];
832 int nlast=lF[cnt];
833 int kk=0;
834 /* cnt--; */
835 /* now lF[cnt] should be <=iF-1 */
836
837 while (Op[first]!=0)
838 {
839 t=lF[cnt]; /* cnt as it was computed */
840
841 poly w=gnc_uu_Mult_ww(t,Op[t],jG,bG,r);
842 c[cnt]=n_Copy(pGetCoeff(w),r->cf);
843 D = pNext(w); /* getting coef and rest D */
844 p_LmDelete(&w,r);
845 w=NULL;
846
847 Op[t]= 0;
848 Pp=p_One(r);
849 p_SetExpV(Pp,Op,r);
850 p_Setm(Pp,r);
851
852 if (t<nlast)
853 {
854 kk=lF[cnt+1];
855 On[kk]=F[kk];
856
857 Pn=p_One(r);
858 p_SetExpV(Pn,On,r);
859 p_Setm(Pn,r);
860
861 if (t!=first) /* typical expr */
862 {
863 w=gnc_p_Mult_mm(D,Pn,r);
864 Rout=gnc_p_mm_Mult(w,Pp,r);
865 w=NULL;
866 }
867 else /* last step */
868 {
869 On[t]=0;
870 p_SetExpV(Pn,On,r);
871 p_Setm(Pn,r);
872 Rout=gnc_p_Mult_mm(D,Pn,r);
873 }
874 #ifdef PDEBUG
875 p_Test(Pp,r);
876 #endif
877 p_Delete(&Pn,r);
878 }
879 else /* first step */
880 {
881 Rout=gnc_p_mm_Mult(D,Pp,r);
882 }
883 #ifdef PDEBUG
884 p_Test(Pp,r);
885 #endif
886 p_Delete(&Pp,r);
887 num=n_Mult(c[cnt+1],c[cnt],r->cf);
888 n_Delete(&c[cnt],r->cf);
889 c[cnt]=num;
890 Rout=__p_Mult_nn(Rout,c[cnt+1],r); /* Rest is ready */
891 out=p_Add_q(out,Rout,r);
892 Pp=NULL;
893 cnt--;
894 }
895 /* only to feel safe:*/
896 Pn=Pp=NULL;
897 freeT(On,rN);
898 freeT(Op,rN);
899
900 /* leadterm and Prv-part with coef 1 */
901 /* U[0]=exp; */
902 /* U[jG]=U[jG]+bG; */
903 /* make leadterm */
904 /* ??????????? we have done it already :-0 */
905
906 Rout=p_One(r);
907 p_SetExpV(Rout,U,r);
908 p_Setm(Rout,r); /* use again this name */
909 p_SetCoeff(Rout,c[cnt+1],r); /* last computed coef */
910
911 out=p_Add_q(out,Rout,r);
912
913 Rout=NULL;
914
915 freeT(U, rN);
916 freeN(c, i);
917 freeT(lF, rN);
918
919 if (cnf!=0)
920 {
921 Rout=p_One(r);
922 p_SetExpV(Rout,Prv,r);
923 p_Setm(Rout,r);
924 freeT(Prv, rN);
925 out=gnc_p_mm_Mult(out,Rout,r); /* getting the final result */
926 p_Delete(&Rout,r);
927 }
928
929 return (out);
930 }
931
gnc_uu_Mult_ww_vert(int i,int a,int j,int b,const ring r)932 poly gnc_uu_Mult_ww_vert (int i, int a, int j, int b, const ring r)
933 {
934 int k,m;
935 int rN=r->N;
936 const int cMTindex = UPMATELEM(j,i,rN);
937 matrix cMT=r->GetNC()->MT[cMTindex]; /* cMT=current MT */
938
939 poly x=p_One(r);p_SetExp(x,j,1,r);p_Setm(x,r);
940 /* var(j); */
941 poly y=p_One(r);p_SetExp(y,i,1,r);p_Setm(y,r);
942 /*var(i); for convenience */
943 #ifdef PDEBUG
944 p_Test(x,r);
945 p_Test(y,r);
946 #endif
947 poly t=NULL;
948 /* ------------ Main Cycles ----------------------------*/
949
950 for (k=2;k<=a;k++)
951 {
952 t = MATELEM(cMT,k,1);
953
954 if (t==NULL) /* not computed yet */
955 {
956 t = nc_p_CopyGet(MATELEM(cMT,k-1,1),r);
957 // t=p_Copy(MATELEM(cMT,k-1,1),r);
958 t = gnc_p_mm_Mult(t,y,r);
959 cMT=r->GetNC()->MT[cMTindex]; // since multiplication can change the MT table...
960 assume( t != NULL );
961 #ifdef PDEBUG
962 p_Test(t,r);
963 #endif
964 MATELEM(cMT,k,1) = nc_p_CopyPut(t,r);
965 // omCheckAddr(cMT->m);
966 p_Delete(&t,r);
967 }
968 t=NULL;
969 }
970
971 for (m=2;m<=b;m++)
972 {
973 t = MATELEM(cMT,a,m);
974 // t=MATELEM(cMT,a,m);
975 if (t==NULL) //not computed yet
976 {
977 t = nc_p_CopyGet(MATELEM(cMT,a,m-1),r);
978 assume( t != NULL );
979 // t=p_Copy(MATELEM(cMT,a,m-1),r);
980 t = gnc_p_Mult_mm(t,x,r);
981 cMT=r->GetNC()->MT[cMTindex]; // since multiplication can change the MT table...
982 #ifdef PDEBUG
983 p_Test(t,r);
984 #endif
985 MATELEM(cMT,a,m) = nc_p_CopyPut(t,r);
986 // MATELEM(cMT,a,m) = t;
987 // omCheckAddr(cMT->m);
988 p_Delete(&t,r);
989 }
990 t=NULL;
991 }
992 p_Delete(&x,r);
993 p_Delete(&y,r);
994 t=MATELEM(cMT,a,b);
995 assume( t != NULL );
996
997 t= nc_p_CopyGet(t,r);
998 #ifdef PDEBUG
999 p_Test(t,r);
1000 #endif
1001 // return(p_Copy(t,r));
1002 /* since the last computed element was cMT[a,b] */
1003 return(t);
1004 }
1005
1006
gnc_uu_Mult_ww_formula(int i,int a,int j,int b,const ring r)1007 static inline poly gnc_uu_Mult_ww_formula (int i, int a, int j, int b, const ring r)
1008 {
1009 if(ncExtensions(NOFORMULAMASK))
1010 return gnc_uu_Mult_ww_vert(i, a, j, b, r);
1011
1012 CFormulaPowerMultiplier* FormulaMultiplier = GetFormulaPowerMultiplier(r);
1013 Enum_ncSAType PairType = _ncSA_notImplemented;
1014
1015 if( FormulaMultiplier != NULL )
1016 PairType = FormulaMultiplier->GetPair(j, i);
1017
1018
1019 if( PairType == _ncSA_notImplemented )
1020 return gnc_uu_Mult_ww_vert(i, a, j, b, r);
1021
1022
1023 // return FormulaMultiplier->Multiply(j, i, b, a);
1024 poly t = CFormulaPowerMultiplier::Multiply( PairType, j, i, b, a, r);
1025
1026 int rN=r->N;
1027 matrix cMT = r->GetNC()->MT[UPMATELEM(j,i,rN)]; /* cMT=current MT */
1028
1029
1030 MATELEM(cMT, a, b) = nc_p_CopyPut(t,r);
1031
1032 // t=MATELEM(cMT,a,b);
1033 // t= nc_p_CopyGet(MATELEM(cMT,a,b),r);
1034 // return(p_Copy(t,r));
1035 /* since the last computed element was cMT[a,b] */
1036 return(t);
1037 }
1038
1039
gnc_uu_Mult_ww(int i,int a,int j,int b,const ring r)1040 poly gnc_uu_Mult_ww (int i, int a, int j, int b, const ring r)
1041 /* (x_i)^a times (x_j)^b */
1042 /* x_i = y, x_j = x ! */
1043 {
1044 /* Check zero exceptions, (q-)commutativity and is there something to do? */
1045 assume(a!=0);
1046 assume(b!=0);
1047 poly out=p_One(r);
1048 if (i<=j)
1049 {
1050 p_SetExp(out,i,a,r);
1051 p_AddExp(out,j,b,r);
1052 p_Setm(out,r);
1053 return(out);
1054 }/* zero exeptions and usual case */
1055 /* if ((a==0)||(b==0)||(i<=j)) return(out); */
1056
1057 if (MATELEM(r->GetNC()->COM,j,i)!=NULL)
1058 /* commutative or quasicommutative case */
1059 {
1060 p_SetExp(out,i,a,r);
1061 p_AddExp(out,j,b,r);
1062 p_Setm(out,r);
1063 if (n_IsOne(pGetCoeff(MATELEM(r->GetNC()->COM,j,i)),r->cf)) /* commutative case */
1064 {
1065 return(out);
1066 }
1067 else
1068 {
1069 number tmp_number=pGetCoeff(MATELEM(r->GetNC()->COM,j,i)); /* quasicommutative case */
1070 n_Power(tmp_number,a*b,&tmp_number, r->cf); // BUG! ;-(
1071 p_SetCoeff(out,tmp_number,r);
1072 return(out);
1073 }
1074 }/* end_of commutative or quasicommutative case */
1075 p_Delete(&out,r);
1076
1077
1078 if(ncExtensions(NOCACHEMASK) && !ncExtensions(NOFORMULAMASK)) // don't use cache whenever possible!
1079 { // without cache!?
1080 CFormulaPowerMultiplier* FormulaMultiplier = GetFormulaPowerMultiplier(r);
1081 Enum_ncSAType PairType = _ncSA_notImplemented;
1082
1083 if( FormulaMultiplier != NULL )
1084 PairType = FormulaMultiplier->GetPair(j, i);
1085
1086 if( PairType != _ncSA_notImplemented )
1087 // // return FormulaMultiplier->Multiply(j, i, b, a);
1088 return CFormulaPowerMultiplier::Multiply( PairType, j, i, b, a, r);
1089 }
1090
1091
1092 /* we are here if i>j and variables do not commute or quasicommute */
1093 /* in fact, now a>=1 and b>=1; and j<i */
1094 /* now check whether the polynomial is already computed */
1095 int rN=r->N;
1096 int vik = UPMATELEM(j,i,rN);
1097 int cMTsize=r->GetNC()->MTsize[vik];
1098 int newcMTsize=0;
1099 newcMTsize=si_max(a,b);
1100
1101 if (newcMTsize<=cMTsize)
1102 {
1103 out = nc_p_CopyGet(MATELEM(r->GetNC()->MT[vik],a,b),r);
1104 if (out !=NULL) return (out);
1105 }
1106 int k,m;
1107 if (newcMTsize > cMTsize)
1108 {
1109 int inM=(((newcMTsize+6)/7)*7);
1110 assume (inM>=newcMTsize);
1111 newcMTsize = inM;
1112 // matrix tmp = (matrix)omAlloc0(inM*inM*sizeof(poly));
1113 matrix tmp = mpNew(newcMTsize,newcMTsize);
1114
1115 for (k=1;k<=cMTsize;k++)
1116 {
1117 for (m=1;m<=cMTsize;m++)
1118 {
1119 out = MATELEM(r->GetNC()->MT[UPMATELEM(j,i,rN)],k,m);
1120 if ( out != NULL )
1121 {
1122 MATELEM(tmp,k,m) = out;/*MATELEM(r->GetNC()->MT[UPMATELEM(j,i,rN)],k,m)*/
1123 // omCheckAddr(tmp->m);
1124 MATELEM(r->GetNC()->MT[UPMATELEM(j,i,rN)],k,m)=NULL;
1125 // omCheckAddr(r->GetNC()->MT[UPMATELEM(j,i,rN)]->m);
1126 out=NULL;
1127 }
1128 }
1129 }
1130 id_Delete((ideal *)&(r->GetNC()->MT[UPMATELEM(j,i,rN)]),r);
1131 r->GetNC()->MT[UPMATELEM(j,i,rN)] = tmp;
1132 tmp=NULL;
1133 r->GetNC()->MTsize[UPMATELEM(j,i,rN)] = newcMTsize;
1134 }
1135 /* The update of multiplication matrix is finished */
1136
1137
1138 return gnc_uu_Mult_ww_formula(i, a, j, b, r);
1139
1140 out = gnc_uu_Mult_ww_vert(i, a, j, b, r);
1141 // out = nc_uu_Mult_ww_horvert(i, a, j, b, r);
1142 return(out);
1143 }
1144
gnc_uu_Mult_ww_horvert(int i,int a,int j,int b,const ring r)1145 poly gnc_uu_Mult_ww_horvert (int i, int a, int j, int b, const ring r)
1146
1147 {
1148 int k,m;
1149 int rN=r->N;
1150 matrix cMT=r->GetNC()->MT[UPMATELEM(j,i,rN)]; /* cMT=current MT */
1151
1152 poly x=p_One(r);p_SetExp(x,j,1,r);p_Setm(x,r);/* var(j); */
1153 poly y=p_One(r);p_SetExp(y,i,1,r);p_Setm(y,r); /*var(i); for convenience */
1154 #ifdef PDEBUG
1155 p_Test(x,r);
1156 p_Test(y,r);
1157 #endif
1158
1159 poly t=NULL;
1160
1161 int toXY;
1162 int toYX;
1163
1164 if (a==1) /* y*x^b, b>=2 */
1165 {
1166 toXY=b-1;
1167 while ( (MATELEM(cMT,1,toXY)==NULL) && (toXY>=2)) toXY--;
1168 for (m=toXY+1;m<=b;m++)
1169 {
1170 t=MATELEM(cMT,1,m);
1171 if (t==NULL) /* remove after debug */
1172 {
1173 t = p_Copy(MATELEM(cMT,1,m-1),r);
1174 t = gnc_p_Mult_mm(t,x,r);
1175 MATELEM(cMT,1,m) = t;
1176 /* omCheckAddr(cMT->m); */
1177 }
1178 else
1179 {
1180 /* Error, should never get there */
1181 WarnS("Error: a=1; MATELEM!=0");
1182 }
1183 t=NULL;
1184 }
1185 return(p_Copy(MATELEM(cMT,1,b),r));
1186 }
1187
1188 if (b==1) /* y^a*x, a>=2 */
1189 {
1190 toYX=a-1;
1191 while ( (MATELEM(cMT,toYX,1)==NULL) && (toYX>=2)) toYX--;
1192 for (m=toYX+1;m<=a;m++)
1193 {
1194 t=MATELEM(cMT,m,1);
1195 if (t==NULL) /* remove after debug */
1196 {
1197 t = p_Copy(MATELEM(cMT,m-1,1),r);
1198 t = gnc_p_mm_Mult(t,y,r);
1199 MATELEM(cMT,m,1) = t;
1200 /* omCheckAddr(cMT->m); */
1201 }
1202 else
1203 {
1204 /* Error, should never get there */
1205 WarnS("Error: b=1, MATELEM!=0");
1206 }
1207 t=NULL;
1208 }
1209 return(p_Copy(MATELEM(cMT,a,1),r));
1210 }
1211
1212 /* ------------ Main Cycles ----------------------------*/
1213 /* a>1, b>1 */
1214
1215 int dXY=0; int dYX=0;
1216 /* dXY = distance for computing x-mult, then y-mult */
1217 /* dYX = distance for computing y-mult, then x-mult */
1218 int toX=a-1; int toY=b-1; /* toX = to axe X, toY = to axe Y */
1219 toXY=b-1; toYX=a-1;
1220 /* if toX==0, toXY = dist. to computed y * x^toXY */
1221 /* if toY==0, toYX = dist. to computed y^toYX * x */
1222 while ( (MATELEM(cMT,toX,b)==NULL) && (toX>=1)) toX--;
1223 if (toX==0) /* the whole column is not computed yet */
1224 {
1225 while ( (MATELEM(cMT,1,toXY)==NULL) && (toXY>=1)) toXY--;
1226 /* toXY >=1 */
1227 dXY=b-1-toXY;
1228 }
1229 dXY=dXY+a-toX; /* the distance to nearest computed y^toX x^b */
1230
1231 while ( (MATELEM(cMT,a,toY)==NULL) && (toY>=1)) toY--;
1232 if (toY==0) /* the whole row is not computed yet */
1233 {
1234 while ( (MATELEM(cMT,toYX,1)==NULL) && (toYX>=1)) toYX--;
1235 /* toYX >=1 */
1236 dYX=a-1-toYX;
1237 }
1238 dYX=dYX+b-toY; /* the distance to nearest computed y^a x^toY */
1239
1240 if (dYX>=dXY)
1241 {
1242 /* first x, then y */
1243 if (toX==0) /* start with the row*/
1244 {
1245 for (m=toXY+1;m<=b;m++)
1246 {
1247 t=MATELEM(cMT,1,m);
1248 if (t==NULL) /* remove after debug */
1249 {
1250 t = p_Copy(MATELEM(cMT,1,m-1),r);
1251 t = gnc_p_Mult_mm(t,x,r);
1252 MATELEM(cMT,1,m) = t;
1253 /* omCheckAddr(cMT->m); */
1254 }
1255 else
1256 {
1257 /* Error, should never get there */
1258 WarnS("dYX>=dXY,toXY; MATELEM==0");
1259 }
1260 t=NULL;
1261 }
1262 toX=1; /* y*x^b is computed */
1263 }
1264 /* Now toX>=1 */
1265 for (k=toX+1;k<=a;k++)
1266 {
1267 t=MATELEM(cMT,k,b);
1268 if (t==NULL) /* remove after debug */
1269 {
1270 t = p_Copy(MATELEM(cMT,k-1,b),r);
1271 t = gnc_p_mm_Mult(t,y,r);
1272 MATELEM(cMT,k,b) = t;
1273 /* omCheckAddr(cMT->m); */
1274 }
1275 else
1276 {
1277 /* Error, should never get there */
1278 WarnS("dYX>=dXY,toX; MATELEM==0");
1279 }
1280 t=NULL;
1281 }
1282 } /* endif (dYX>=dXY) */
1283
1284
1285 if (dYX<dXY)
1286 {
1287 /* first y, then x */
1288 if (toY==0) /* start with the column*/
1289 {
1290 for (m=toYX+1;m<=a;m++)
1291 {
1292 t=MATELEM(cMT,m,1);
1293 if (t==NULL) /* remove after debug */
1294 {
1295 t = p_Copy(MATELEM(cMT,m-1,1),r);
1296 t = gnc_p_mm_Mult(t,y,r);
1297 MATELEM(cMT,m,1) = t;
1298 /* omCheckAddr(cMT->m); */
1299 }
1300 else
1301 {
1302 /* Error, should never get there */
1303 WarnS("dYX<dXY,toYX; MATELEM==0");
1304 }
1305 t=NULL;
1306 }
1307 toY=1; /* y^a*x is computed */
1308 }
1309 /* Now toY>=1 */
1310 for (k=toY+1;k<=b;k++)
1311 {
1312 t=MATELEM(cMT,a,k);
1313 if (t==NULL) /* remove after debug */
1314 {
1315 t = p_Copy(MATELEM(cMT,a,k-1),r);
1316 t = gnc_p_Mult_mm(t,x,r);
1317 MATELEM(cMT,a,k) = t;
1318 /* omCheckAddr(cMT->m); */
1319 }
1320 else
1321 {
1322 /* Error, should never get there */
1323 WarnS("dYX<dXY,toY; MATELEM==0");
1324 }
1325 t=NULL;
1326 }
1327 } /* endif (dYX<dXY) */
1328
1329 p_Delete(&x,r);
1330 p_Delete(&y,r);
1331 t=p_Copy(MATELEM(cMT,a,b),r);
1332 return(t); /* since the last computed element was cMT[a,b] */
1333 }
1334
1335
1336 /* ----------------------------- Syzygies ---------------------- */
1337
1338 /*2
1339 * reduction of p2 with p1
1340 * do not destroy p1, but p2
1341 * p1 divides p2 -> for use in NF algorithm
1342 */
gnc_ReduceSpolyOld(const poly p1,poly p2,const ring r)1343 poly gnc_ReduceSpolyOld(const poly p1, poly p2/*,poly spNoether*/, const ring r)
1344 {
1345 assume(p_LmDivisibleBy(p1, p2, r));
1346
1347 #ifdef PDEBUG
1348 if (p_GetComp(p1,r)!=p_GetComp(p2,r)
1349 && (p_GetComp(p1,r)!=0)
1350 && (p_GetComp(p2,r)!=0))
1351 {
1352 dReportError("nc_ReduceSpolyOld: different components");
1353 return(NULL);
1354 }
1355 #endif
1356 poly m = p_One(r);
1357 p_ExpVectorDiff(m,p2,p1,r);
1358 //p_Setm(m,r);
1359 #ifdef PDEBUG
1360 p_Test(m,r);
1361 #endif
1362 /* pSetComp(m,r)=0? */
1363 poly N = nc_mm_Mult_p(m, p_Head(p1,r), r);
1364 number C = p_GetCoeff(N, r);
1365 number cF = p_GetCoeff(p2, r);
1366 /* GCD stuff */
1367 number cG = n_SubringGcd(C, cF, r->cf);
1368 if ( !n_IsOne(cG,r->cf) )
1369 {
1370 cF = n_Div(cF, cG, r->cf); n_Normalize(cF, r->cf);
1371 C = n_Div(C, cG, r->cf); n_Normalize(C, r->cf);
1372 }
1373 else
1374 {
1375 cF = n_Copy(cF, r->cf);
1376 C = n_Copy(C, r->cf);
1377 }
1378 n_Delete(&cG,r->cf);
1379 p2 = __p_Mult_nn(p2, C, r);
1380 poly out = nc_mm_Mult_pp(m, pNext(p1), r);
1381 N = p_Add_q(N, out, r);
1382 p_Test(p2,r);
1383 p_Test(N,r);
1384 if (!n_IsMOne(cF,r->cf))
1385 {
1386 cF = n_InpNeg(cF,r->cf);
1387 N = __p_Mult_nn(N, cF, r);
1388 p_Test(N,r);
1389 }
1390 out = p_Add_q(p2,N,r);
1391 p_Test(out,r);
1392 if ( out!=NULL ) p_Cleardenom(out,r);
1393 p_Delete(&m,r);
1394 n_Delete(&cF,r->cf);
1395 n_Delete(&C,r->cf);
1396 return(out);
1397 }
1398
gnc_ReduceSpolyNew(const poly p1,poly p2,const ring r)1399 poly gnc_ReduceSpolyNew(const poly p1, poly p2, const ring r)
1400 {
1401 assume(p_LmDivisibleBy(p1, p2, r));
1402
1403 const long lCompP1 = p_GetComp(p1,r);
1404 const long lCompP2 = p_GetComp(p2,r);
1405
1406 if ((lCompP1!=lCompP2) && (lCompP1!=0) && (lCompP2!=0))
1407 {
1408 #ifdef PDEBUG
1409 WerrorS("gnc_ReduceSpolyNew: different non-zero components!");
1410 #endif
1411 return(NULL);
1412 }
1413
1414 poly m = p_One(r);
1415 p_ExpVectorDiff(m, p2, p1, r);
1416 //p_Setm(m,r);
1417 #ifdef PDEBUG
1418 p_Test(m,r);
1419 #endif
1420
1421 /* pSetComp(m,r)=0? */
1422 poly N = nc_mm_Mult_p(m, p_Head(p1,r), r);
1423
1424 number C = p_GetCoeff(N, r);
1425 number cF = p_GetCoeff(p2, r);
1426
1427 /* GCD stuff */
1428 number cG = n_SubringGcd(C, cF, r->cf);
1429
1430 if (!n_IsOne(cG, r->cf))
1431 {
1432 cF = n_Div(cF, cG, r->cf); n_Normalize(cF, r->cf);
1433 C = n_Div(C, cG, r->cf); n_Normalize(C, r->cf);
1434 }
1435 else
1436 {
1437 cF = n_Copy(cF, r->cf);
1438 C = n_Copy(C, r->cf);
1439 }
1440 n_Delete(&cG,r->cf);
1441
1442 p2 = __p_Mult_nn(p2, C, r); // p2 !!!
1443 p_Test(p2,r);
1444 n_Delete(&C,r->cf);
1445
1446 poly out = nc_mm_Mult_pp(m, pNext(p1), r);
1447 p_Delete(&m,r);
1448
1449 N = p_Add_q(N, out, r);
1450 p_Test(N,r);
1451
1452 if (!n_IsMOne(cF,r->cf)) // ???
1453 {
1454 cF = n_InpNeg(cF,r->cf);
1455 N = __p_Mult_nn(N, cF, r);
1456 p_Test(N,r);
1457 }
1458 n_Delete(&cF,r->cf);
1459
1460 out = p_Add_q(p2,N,r); // delete N, p2
1461 p_Test(out,r);
1462 if ( out!=NULL ) p_Cleardenom(out,r);
1463 return(out);
1464 }
1465
1466
1467 /*4
1468 * creates the S-polynomial of p1 and p2
1469 * do not destroy p1 and p2
1470 */
gnc_CreateSpolyOld(poly p1,poly p2,const ring r)1471 poly gnc_CreateSpolyOld(poly p1, poly p2/*,poly spNoether*/, const ring r)
1472 {
1473 #ifdef PDEBUG
1474 if ((p_GetComp(p1,r)!=p_GetComp(p2,r))
1475 && (p_GetComp(p1,r)!=0)
1476 && (p_GetComp(p2,r)!=0))
1477 {
1478 dReportError("gnc_CreateSpolyOld : different components!");
1479 return(NULL);
1480 }
1481 #endif
1482 if ((ncRingType(r)==nc_lie) && p_HasNotCF(p1,p2, r)) /* prod crit */
1483 {
1484 return(nc_p_Bracket_qq(p_Copy(p2, r),p1, r));
1485 }
1486 poly pL=p_One(r);
1487 poly m1=p_One(r);
1488 poly m2=p_One(r);
1489 pL = p_Lcm(p1,p2,r);
1490 p_Setm(pL,r);
1491 #ifdef PDEBUG
1492 p_Test(pL,r);
1493 #endif
1494 p_ExpVectorDiff(m1,pL,p1,r);
1495 //p_SetComp(m1,0,r);
1496 //p_Setm(m1,r);
1497 #ifdef PDEBUG
1498 p_Test(m1,r);
1499 #endif
1500 p_ExpVectorDiff(m2,pL,p2,r);
1501 //p_SetComp(m2,0,r);
1502 //p_Setm(m2,r);
1503 #ifdef PDEBUG
1504 p_Test(m2,r);
1505 #endif
1506 p_Delete(&pL,r);
1507 /* zero exponents ! */
1508 poly M1 = nc_mm_Mult_p(m1,p_Head(p1,r),r);
1509 number C1 = p_GetCoeff(M1,r);
1510 poly M2 = nc_mm_Mult_p(m2,p_Head(p2,r),r);
1511 number C2 = p_GetCoeff(M2,r);
1512 /* GCD stuff */
1513 number C = n_SubringGcd(C1,C2,r->cf);
1514 if (!n_IsOne(C,r->cf))
1515 {
1516 C1=n_Div(C1,C, r->cf);n_Normalize(C1,r->cf);
1517 C2=n_Div(C2,C, r->cf);n_Normalize(C2,r->cf);
1518 }
1519 else
1520 {
1521 C1=n_Copy(C1, r->cf);
1522 C2=n_Copy(C2, r->cf);
1523 }
1524 n_Delete(&C,r->cf);
1525 M1=__p_Mult_nn(M1,C2,r);
1526 p_SetCoeff(m1,C2,r);
1527 if (n_IsMOne(C1,r->cf))
1528 {
1529 M2=p_Add_q(M1,M2,r);
1530 }
1531 else
1532 {
1533 C1=n_InpNeg(C1,r->cf);
1534 M2=__p_Mult_nn(M2,C1,r);
1535 M2=p_Add_q(M1,M2,r);
1536 p_SetCoeff(m2,C1,r);
1537 }
1538 /* M1 is killed, M2=res = C2 M1 - C1 M2 */
1539 poly tmp=p_Copy(p1,r);
1540 tmp=p_LmDeleteAndNext(tmp,r);
1541 M1=nc_mm_Mult_p(m1,tmp,r);
1542 tmp=p_Copy(p2,r);
1543 tmp=p_LmDeleteAndNext(tmp,r);
1544 M2=p_Add_q(M2,M1,r);
1545 M1=nc_mm_Mult_p(m2,tmp,r);
1546 M2=p_Add_q(M2,M1,r);
1547 p_Delete(&m1,r);
1548 p_Delete(&m2,r);
1549 // n_Delete(&C1,r);
1550 // n_Delete(&C2,r);
1551 #ifdef PDEBUG
1552 p_Test(M2,r);
1553 #endif
1554 if (M2!=NULL) M2=p_Cleardenom(M2,r);
1555 return(M2);
1556 }
1557
gnc_CreateSpolyNew(poly p1,poly p2,const ring r)1558 poly gnc_CreateSpolyNew(poly p1, poly p2/*,poly spNoether*/, const ring r)
1559 {
1560 #ifdef PDEBUG
1561 p_Test(p1, r);
1562 p_Test(p2, r);
1563 #if MYTEST
1564 PrintS("p1: "); p_Write(p1, r);
1565 PrintS("p2: "); p_Write(p2, r);
1566 #endif
1567 #endif
1568
1569 const long lCompP1 = p_GetComp(p1,r);
1570 const long lCompP2 = p_GetComp(p2,r);
1571
1572 if ((lCompP1!=lCompP2) && (lCompP1!=0) && (lCompP2!=0))
1573 {
1574 #ifdef PDEBUG
1575 WerrorS("gnc_CreateSpolyNew: different non-zero components!");
1576 assume(0);
1577 #endif
1578 return(NULL);
1579 }
1580
1581 // if ((r->GetNC()->type==nc_lie) && pHasNotCF(p1,p2)) /* prod crit */
1582 // {
1583 // return(nc_p_Bracket_qq(pCopy(p2),p1));
1584 // }
1585
1586 // poly pL=p_One( r);
1587
1588 poly m1=p_One( r);
1589 poly m2=p_One( r);
1590
1591 poly pL = p_Lcm(p1,p2,r); // pL = lcm( lm(p1), lm(p2) )
1592
1593
1594 #ifdef PDEBUG
1595 // p_Test(pL,r);
1596 #endif
1597
1598 p_ExpVectorDiff(m1, pL, p1, r); // m1 = pL / lm(p1)
1599 //p_SetComp(m1,0,r);
1600 //p_Setm(m1,r);
1601
1602 #ifdef PDEBUG
1603 p_Test(m1,r);
1604 #endif
1605 // assume(p_GetComp(m1,r) == 0);
1606
1607 p_ExpVectorDiff(m2, pL, p2, r); // m2 = pL / lm(p2)
1608
1609 //p_SetComp(m2,0,r);
1610 //p_Setm(m2,r);
1611 #ifdef PDEBUG
1612 p_Test(m2,r);
1613 #endif
1614
1615 #ifdef PDEBUG
1616 #if MYTEST
1617 PrintS("m1: "); pWrite(m1);
1618 PrintS("m2: "); pWrite(m2);
1619 #endif
1620 #endif
1621
1622
1623 // assume(p_GetComp(m2,r) == 0);
1624
1625 #ifdef PDEBUG
1626 #if 0
1627 if( (p_GetComp(m2,r) != 0) || (p_GetComp(m1,r) != 0) )
1628 {
1629 WarnS("gnc_CreateSpolyNew: wrong monomials!");
1630
1631
1632 #ifdef RDEBUG
1633 PrintS("m1 = "); p_Write(m1, r);
1634 p_DebugPrint(m1, r);
1635
1636 PrintS("m2 = "); p_Write(m2, r);
1637 p_DebugPrint(m2, r);
1638
1639 PrintS("p1 = "); p_Write(p1, r);
1640 p_DebugPrint(p1, r);
1641
1642 PrintS("p2 = "); p_Write(p2, r);
1643 p_DebugPrint(p2, r);
1644
1645 PrintS("pL = "); p_Write(pL, r);
1646 p_DebugPrint(pL, r);
1647 #endif
1648
1649 }
1650
1651 #endif
1652 #endif
1653
1654 p_LmFree(&pL,r);
1655
1656 /* zero exponents !? */
1657 poly M1 = nc_mm_Mult_p(m1,p_Head(p1,r),r); // M1 = m1 * lt(p1)
1658 poly M2 = nc_mm_Mult_p(m2,p_Head(p2,r),r); // M2 = m2 * lt(p2)
1659
1660 #ifdef PDEBUG
1661 p_Test(M1,r);
1662 p_Test(M2,r);
1663
1664 #if MYTEST
1665 PrintS("M1: "); pWrite(M1);
1666 PrintS("M2: "); pWrite(M2);
1667 #endif
1668 #endif
1669
1670 if(M1 == NULL || M2 == NULL)
1671 {
1672 #ifdef PDEBUG
1673 PrintS("\np1 = ");
1674 p_Write(p1, r);
1675
1676 PrintS("m1 = ");
1677 p_Write(m1, r);
1678
1679 PrintS("p2 = ");
1680 p_Write(p2, r);
1681
1682 PrintS("m2 = ");
1683 p_Write(m2, r);
1684
1685 WerrorS("ERROR in nc_CreateSpoly: result of multiplication is Zero!\n");
1686 #endif
1687 return(NULL);
1688 }
1689
1690 number C1 = p_GetCoeff(M1,r); // C1 = lc(M1)
1691 number C2 = p_GetCoeff(M2,r); // C2 = lc(M2)
1692
1693 /* GCD stuff */
1694 number C = n_SubringGcd(C1, C2, r->cf); // C = gcd(C1, C2)
1695
1696 if (!n_IsOne(C, r->cf)) // if C != 1
1697 {
1698 C1=n_Div(C1, C, r->cf);n_Normalize(C1,r->cf); // C1 = C1 / C
1699 C2=n_Div(C2, C, r->cf);n_Normalize(C2,r->cf); // C2 = C2 / C
1700 }
1701 else
1702 {
1703 C1=n_Copy(C1,r->cf);
1704 C2=n_Copy(C2,r->cf);
1705 }
1706
1707 n_Delete(&C,r->cf); // destroy the number C
1708
1709 C1=n_InpNeg(C1,r->cf);
1710
1711 // number MinusOne=n_Init(-1,r);
1712 // if (n_Equal(C1,MinusOne,r)) // lc(M1) / gcd( lc(M1), lc(M2)) == -1 ????
1713 // {
1714 // M2=p_Add_q(M1,M2,r); // ?????
1715 // }
1716 // else
1717 // {
1718 M1=__p_Mult_nn(M1,C2,r); // M1 = (C2*lc(p1)) * (lcm(lm(p1),lm(p2)) / lm(p1)) * lm(p1)
1719
1720 #ifdef PDEBUG
1721 p_Test(M1,r);
1722 #endif
1723
1724 M2=__p_Mult_nn(M2,C1,r); // M2 =(-C1*lc(p2)) * (lcm(lm(p1),lm(p2)) / lm(p2)) * lm(p2)
1725
1726
1727
1728 #ifdef PDEBUG
1729 p_Test(M2,r);
1730
1731 #if MYTEST
1732 PrintS("M1: "); pWrite(M1);
1733 PrintS("M2: "); pWrite(M2);
1734 #endif
1735 #endif
1736
1737
1738 M2=p_Add_q(M1,M2,r); // M1 is killed, M2 = spoly(lt(p1), lt(p2)) = C2*M1 - C1*M2
1739
1740 #ifdef PDEBUG
1741 p_Test(M2,r);
1742
1743 #if MYTEST
1744 PrintS("M2: "); pWrite(M2);
1745 #endif
1746
1747 #endif
1748
1749 // M2 == 0 for supercommutative algebras!
1750 // }
1751 // n_Delete(&MinusOne,r);
1752
1753 p_SetCoeff(m1,C2,r); // lc(m1) = C2!!!
1754 p_SetCoeff(m2,C1,r); // lc(m2) = C1!!!
1755
1756 #ifdef PDEBUG
1757 p_Test(m1,r);
1758 p_Test(m2,r);
1759 #endif
1760
1761 // poly tmp = p_Copy(p1,r); // tmp = p1
1762 // tmp=p_LmDeleteAndNext(tmp,r); // tmp = tail(p1)
1763 //#ifdef PDEBUG
1764 // p_Test(tmp,r);
1765 //#endif
1766
1767 M1 = nc_mm_Mult_pp(m1, pNext(p1), r); // M1 = m1 * tail(p1), delete tmp // ???
1768
1769 #ifdef PDEBUG
1770 p_Test(M1,r);
1771
1772 #if MYTEST
1773 PrintS("M1: "); pWrite(M1);
1774 #endif
1775
1776 #endif
1777
1778 M2=p_Add_q(M2,M1,r); // M2 = spoly(lt(p1), lt(p2)) + m1 * tail(p1), delete M1
1779 #ifdef PDEBUG
1780 M1=NULL;
1781 p_Test(M2,r);
1782
1783 #if MYTEST
1784 PrintS("M2: "); pWrite(M2);
1785 #endif
1786
1787 #endif
1788
1789 // tmp=p_Copy(p2,r); // tmp = p2
1790 // tmp=p_LmDeleteAndNext(tmp,r); // tmp = tail(p2)
1791
1792 //#ifdef PDEBUG
1793 // p_Test(tmp,r);
1794 //#endif
1795
1796 M1 = nc_mm_Mult_pp(m2, pNext(p2), r); // M1 = m2 * tail(p2), detele tmp
1797
1798 #ifdef PDEBUG
1799 p_Test(M1,r);
1800
1801 #if MYTEST
1802 PrintS("M1: "); pWrite(M1);
1803 #endif
1804
1805 #endif
1806
1807 M2 = p_Add_q(M2,M1,r); // M2 = spoly(lt(p1), lt(p2)) + m1 * tail(p1) + m2*tail(p2)
1808
1809 #ifdef PDEBUG
1810 M1=NULL;
1811 p_Test(M2,r);
1812
1813 #if MYTEST
1814 PrintS("M2: "); pWrite(M2);
1815 #endif
1816
1817 #endif
1818
1819 p_Delete(&m1,r); // => n_Delete(&C1,r);
1820 p_Delete(&m2,r); // => n_Delete(&C2,r);
1821
1822 #ifdef PDEBUG
1823 p_Test(M2,r);
1824 #endif
1825
1826 if (M2!=NULL) p_Cleardenom(M2,r);
1827
1828 return(M2);
1829 }
1830
1831
1832
1833
1834 #if 0
1835 /*5
1836 * reduction of tail(q) with p1
1837 * lead(p1) divides lead(pNext(q2)) and pNext(q2) is reduced
1838 * do not destroy p1, but tail(q)
1839 */
1840 void gnc_ReduceSpolyTail(poly p1, poly q, poly q2, poly spNoether, const ring r)
1841 {
1842 poly a1=p_Head(p1,r);
1843 poly Q=pNext(q2);
1844 number cQ=p_GetCoeff(Q,r);
1845 poly m=p_One(r);
1846 p_ExpVectorDiff(m,Q,p1,r);
1847 // p_SetComp(m,0,r);
1848 //p_Setm(m,r);
1849 #ifdef PDEBUG
1850 p_Test(m,r);
1851 #endif
1852 /* pSetComp(m,r)=0? */
1853 poly M = nc_mm_Mult_pp(m, p1,r);
1854 number C=p_GetCoeff(M,r);
1855 M=p_Add_q(M,nc_mm_Mult_p(m,p_LmDeleteAndNext(p_Copy(p1,r),r),r),r); // _pp?
1856 q=__p_Mult_nn(q,C,r);
1857 number MinusOne=n_Init(-1,r->cf);
1858 if (!n_Equal(cQ,MinusOne,r->cf))
1859 {
1860 cQ=nInpNeg(cQ);
1861 M=__p_Mult_nn(M,cQ,r);
1862 }
1863 Q=p_Add_q(Q,M,r);
1864 pNext(q2)=Q;
1865
1866 p_Delete(&m,r);
1867 n_Delete(&C,r->cf);
1868 n_Delete(&cQ,r->cf);
1869 n_Delete(&MinusOne,r->cf);
1870 /* return(q); */
1871 }
1872 #endif
1873
1874
1875 /*6
1876 * creates the commutative lcm(lm(p1),lm(p2))
1877 * do not destroy p1 and p2
1878 */
nc_CreateShortSpoly(poly p1,poly p2,const ring r)1879 poly nc_CreateShortSpoly(poly p1, poly p2, const ring r)
1880 {
1881 #ifdef PDEBUG
1882 p_Test(p1, r);
1883 p_Test(p2, r);
1884 #endif
1885
1886 const long lCompP1 = p_GetComp(p1,r);
1887 const long lCompP2 = p_GetComp(p2,r);
1888
1889 if ((lCompP1!=lCompP2) && (lCompP1!=0) && (lCompP2!=0))
1890 {
1891 #ifdef PDEBUG
1892 WerrorS("nc_CreateShortSpoly: wrong module components!"); // !!!!
1893 #endif
1894 return(NULL);
1895 }
1896
1897 poly m;
1898
1899 #ifdef HAVE_RATGRING
1900 if ( rIsRatGRing(r))
1901 {
1902 /* rational version */
1903 m = p_LcmRat(p1, p2, si_max(lCompP1, lCompP2), r);
1904 } else
1905 #endif
1906 {
1907 m = p_Lcm(p1, p2, r);
1908 }
1909
1910 pSetCoeff0(m,NULL);
1911
1912 return(m);
1913 }
1914
gnc_kBucketPolyRedOld(kBucket_pt b,poly p,number * c)1915 void gnc_kBucketPolyRedOld(kBucket_pt b, poly p, number *c)
1916 {
1917 const ring r = b->bucket_ring;
1918 // b will not be multiplied by any constant in this impl.
1919 // ==> *c=1
1920 if (c!=NULL) *c=n_Init(1, r->cf);
1921 poly m=p_One(r);
1922 p_ExpVectorDiff(m,kBucketGetLm(b),p, r);
1923 //pSetm(m);
1924 #ifdef PDEBUG
1925 p_Test(m, r);
1926 #endif
1927 poly pp= nc_mm_Mult_pp(m,p, r);
1928 assume(pp!=NULL);
1929 p_Delete(&m, r);
1930 number n=pGetCoeff(pp);
1931 number nn;
1932 if (!n_IsMOne(n, r->cf))
1933 {
1934 nn=n_InpNeg(n_Invers(n, r->cf), r->cf);
1935 n= n_Mult(nn,pGetCoeff(kBucketGetLm(b)), r->cf);
1936 n_Delete(&nn, r->cf);
1937 pp=__p_Mult_nn(pp,n,r);
1938 n_Delete(&n, r->cf);
1939 }
1940 else
1941 {
1942 pp=__p_Mult_nn(pp,p_GetCoeff(kBucketGetLm(b), r),r);
1943 }
1944 int l=pLength(pp);
1945 kBucket_Add_q(b,pp,&l);
1946 }
1947
gnc_kBucketPolyRedNew(kBucket_pt b,poly p,number * c)1948 void gnc_kBucketPolyRedNew(kBucket_pt b, poly p, number *c)
1949 {
1950 const ring r = b->bucket_ring;
1951 #ifdef PDEBUG
1952 // PrintS(">*");
1953 #endif
1954
1955 #ifdef KDEBUG
1956 if( !kbTest(b) ) WerrorS("nc_kBucketPolyRed: broken bucket!");
1957 #endif
1958
1959 #ifdef PDEBUG
1960 p_Test(p, r);
1961 #if MYTEST
1962 PrintS("p: "); p_Write(p, r);
1963 #endif
1964 #endif
1965
1966 // b will not be multiplied by any constant in this impl.
1967 // ==> *c=1
1968 if (c!=NULL) *c=n_Init(1, r->cf);
1969 poly m = p_One(r);
1970 const poly pLmB = kBucketGetLm(b); // no new copy!
1971
1972 assume( pLmB != NULL );
1973
1974 #ifdef PDEBUG
1975 p_Test(pLmB, r);
1976
1977 #if MYTEST
1978 PrintS("pLmB: "); p_Write(pLmB, r);
1979 #endif
1980 #endif
1981
1982 p_ExpVectorDiff(m, pLmB, p, r);
1983 //pSetm(m);
1984
1985 #ifdef PDEBUG
1986 p_Test(m, r);
1987 #if MYTEST
1988 PrintS("m: "); p_Write(m, r);
1989 #endif
1990 #endif
1991
1992 poly pp = nc_mm_Mult_pp(m, p, r);
1993 p_Delete(&m, r);
1994
1995 assume( pp != NULL );
1996 const number n = pGetCoeff(pp); // bug!
1997
1998 if (!n_IsMOne(n, r->cf) ) // does this improve performance??!? also see below... // TODO: check later on.
1999 // if n == -1 => nn = 1 and -1/n
2000 {
2001 number nn=n_InpNeg(n_Invers(n, r->cf), r->cf);
2002 number t = n_Mult(nn,pGetCoeff(pLmB), r->cf);
2003 n_Delete(&nn, r->cf);
2004 pp = __p_Mult_nn(pp,t,r);
2005 n_Delete(&t, r->cf);
2006 }
2007 else
2008 {
2009 pp = __p_Mult_nn(pp,p_GetCoeff(pLmB, r), r);
2010 }
2011
2012 int l = pLength(pp);
2013
2014 #ifdef PDEBUG
2015 p_Test(pp, r);
2016 // PrintS("PP: "); pWrite(pp);
2017 #endif
2018
2019 kBucket_Add_q(b,pp,&l);
2020
2021
2022 #ifdef PDEBUG
2023 // PrintS("*>");
2024 #endif
2025 }
2026
2027
gnc_kBucketPolyRed_ZOld(kBucket_pt b,poly p,number * c)2028 void gnc_kBucketPolyRed_ZOld(kBucket_pt b, poly p, number *c)
2029 {
2030 const ring r = b->bucket_ring;
2031 // b is multiplied by a constant in this impl.
2032 number ctmp;
2033 poly m=p_One(r);
2034 p_ExpVectorDiff(m,kBucketGetLm(b),p, r);
2035 //pSetm(m);
2036 #ifdef PDEBUG
2037 p_Test(m, r);
2038 #endif
2039 if(p_IsConstant(m,r))
2040 {
2041 p_Delete(&m, r);
2042 ctmp = kBucketPolyRed(b,p,pLength(p),NULL);
2043 }
2044 else
2045 {
2046 poly pp = nc_mm_Mult_pp(m,p,r);
2047 number c2;
2048 p_Cleardenom_n(pp,r,c2);
2049 p_Delete(&m, r);
2050 ctmp = kBucketPolyRed(b,pp,pLength(pp),NULL);
2051 //cc=*c;
2052 //*c=nMult(*c,c2);
2053 n_Delete(&c2, r->cf);
2054 //nDelete(&cc);
2055 p_Delete(&pp, r);
2056 }
2057 if (c!=NULL) *c=ctmp;
2058 else n_Delete(&ctmp, r->cf);
2059 }
2060
gnc_kBucketPolyRed_ZNew(kBucket_pt b,poly p,number * c)2061 void gnc_kBucketPolyRed_ZNew(kBucket_pt b, poly p, number *c)
2062 {
2063 const ring r = b->bucket_ring;
2064 // b is multiplied by a constant in this impl.
2065 number ctmp;
2066 poly m=p_One(r);
2067 p_ExpVectorDiff(m,kBucketGetLm(b),p, r);
2068 //pSetm(m);
2069 #ifdef PDEBUG
2070 p_Test(m, r);
2071 #endif
2072
2073 if(p_IsConstant(m,r))
2074 {
2075 p_Delete(&m, r);
2076 ctmp = kBucketPolyRed(b,p,pLength(p),NULL);
2077 }
2078 else
2079 {
2080 poly pp = nc_mm_Mult_pp(m,p,r);
2081 number c2;
2082 p_Cleardenom_n(pp,r,c2);
2083 p_Delete(&m, r);
2084 ctmp = kBucketPolyRed(b,pp,pLength(pp),NULL);
2085 //cc=*c;
2086 //*c=nMult(*c,c2);
2087 n_Delete(&c2, r->cf);
2088 //nDelete(&cc);
2089 p_Delete(&pp, r);
2090 }
2091 if (c!=NULL) *c=ctmp;
2092 else n_Delete(&ctmp, r->cf);
2093 }
2094
2095
nc_PolyPolyRedOld(poly & b,poly p,number * c,const ring r)2096 inline void nc_PolyPolyRedOld(poly &b, poly p, number *c, const ring r)
2097 // reduces b with p, do not delete both
2098 {
2099 // b will not by multiplied by any constant in this impl.
2100 // ==> *c=1
2101 if (c!=NULL) *c=n_Init(1, r->cf);
2102 poly m=p_One(r);
2103 p_ExpVectorDiff(m,p_Head(b, r),p, r);
2104 //pSetm(m);
2105 #ifdef PDEBUG
2106 p_Test(m, r);
2107 #endif
2108 poly pp=nc_mm_Mult_pp(m,p,r);
2109 assume(pp!=NULL);
2110
2111 p_Delete(&m, r);
2112 number n=pGetCoeff(pp);
2113 number nn;
2114 if (!n_IsMOne(n, r->cf))
2115 {
2116 nn=n_InpNeg(n_Invers(n, r->cf), r->cf);
2117 n =n_Mult(nn,pGetCoeff(b), r->cf);
2118 n_Delete(&nn, r->cf);
2119 pp=__p_Mult_nn(pp,n,r);
2120 n_Delete(&n, r->cf);
2121 }
2122 else
2123 {
2124 pp=__p_Mult_nn(pp,p_GetCoeff(b, r),r);
2125 }
2126 b=p_Add_q(b,pp,r);
2127 }
2128
2129
nc_PolyPolyRedNew(poly & b,poly p,number * c,const ring r)2130 inline void nc_PolyPolyRedNew(poly &b, poly p, number *c, const ring r)
2131 // reduces b with p, do not delete both
2132 {
2133 #ifdef PDEBUG
2134 p_Test(b, r);
2135 p_Test(p, r);
2136 #endif
2137
2138 #if MYTEST
2139 PrintS("nc_PolyPolyRedNew(");
2140 p_Write0(b, r);
2141 PrintS(", ");
2142 p_Write0(p, r);
2143 PrintS(", *c): ");
2144 #endif
2145
2146 // b will not by multiplied by any constant in this impl.
2147 // ==> *c=1
2148 if (c!=NULL) *c=n_Init(1, r->cf);
2149
2150 poly pp = NULL;
2151
2152 // there is a problem when p is a square(=>0!)
2153
2154 while((b != NULL) && (pp == NULL))
2155 {
2156
2157 // poly pLmB = p_Head(b, r);
2158 poly m = p_One(r);
2159 p_ExpVectorDiff(m, b, p, r);
2160 // pDelete(&pLmB);
2161 //pSetm(m);
2162
2163 #ifdef PDEBUG
2164 p_Test(m, r);
2165 p_Test(b, r);
2166 #endif
2167
2168 pp = nc_mm_Mult_pp(m, p, r);
2169
2170 #if MYTEST
2171 PrintS("\n{b': ");
2172 p_Write0(b, r);
2173 PrintS(", m: ");
2174 p_Write0(m, r);
2175 PrintS(", pp: ");
2176 p_Write0(pp, r);
2177 PrintS(" }\n");
2178 #endif
2179
2180 p_Delete(&m, r); // one m for all tries!
2181
2182 // assume( pp != NULL );
2183
2184 if( pp == NULL )
2185 {
2186 b = p_LmDeleteAndNext(b, r);
2187
2188 if( !p_DivisibleBy(p, b, r) )
2189 return;
2190
2191 }
2192 }
2193
2194 #if MYTEST
2195 PrintS("{b': ");
2196 p_Write0(b, r);
2197 PrintS(", pp: ");
2198 p_Write0(pp, r);
2199 PrintS(" }\n");
2200 #endif
2201
2202
2203 if(b == NULL) return;
2204
2205
2206 assume(pp != NULL);
2207
2208 const number n = pGetCoeff(pp); // no new copy
2209
2210 number nn;
2211
2212 if (!n_IsMOne(n, r->cf)) // TODO: as above.
2213 {
2214 nn=n_InpNeg(n_Invers(n, r->cf), r->cf);
2215 number t = n_Mult(nn, pGetCoeff(b), r->cf);
2216 n_Delete(&nn, r->cf);
2217 pp=__p_Mult_nn(pp, t, r);
2218 n_Delete(&t, r->cf);
2219 }
2220 else
2221 {
2222 pp=__p_Mult_nn(pp, pGetCoeff(b), r);
2223 }
2224
2225
2226 b=p_Add_q(b,pp,r);
2227
2228 }
2229
nc_PolyPolyRed(poly & b,poly p,number * c,const ring r)2230 void nc_PolyPolyRed(poly &b, poly p, number *c, const ring r)
2231 {
2232 #if 0
2233 nc_PolyPolyRedOld(b, p, c, r);
2234 #else
2235 nc_PolyPolyRedNew(b, p, c, r);
2236 #endif
2237 }
2238
2239
2240 poly nc_mm_Bracket_nn(poly m1, poly m2, const ring r);
2241
2242 /// returns [p,q], destroys p
nc_p_Bracket_qq(poly p,const poly q,const ring r)2243 poly nc_p_Bracket_qq(poly p, const poly q, const ring r)
2244 {
2245 assume(p != NULL && q!= NULL);
2246
2247 if (!rIsPluralRing(r)) return(NULL);
2248 if (p_ComparePolys(p,q, r)) return(NULL);
2249 /* Components !? */
2250 poly Q=NULL;
2251 number coef=NULL;
2252 poly pres=NULL;
2253 int UseBuckets=1;
2254 if (((pLength(p)< MIN_LENGTH_BUCKET/2) && (pLength(q)< MIN_LENGTH_BUCKET/2))
2255 || TEST_OPT_NOT_BUCKETS)
2256 UseBuckets=0;
2257
2258
2259 CPolynomialSummator sum(r, UseBuckets == 0);
2260
2261 while (p!=NULL)
2262 {
2263 Q=q;
2264 while(Q!=NULL)
2265 {
2266 pres=nc_mm_Bracket_nn(p,Q, r); /* since no coeffs are taken into account there */
2267 if (pres!=NULL)
2268 {
2269 coef = n_Mult(pGetCoeff(p),pGetCoeff(Q), r->cf);
2270 pres = __p_Mult_nn(pres,coef,r);
2271
2272 sum += pres;
2273 n_Delete(&coef, r->cf);
2274 }
2275 pIter(Q);
2276 }
2277 p=p_LmDeleteAndNext(p, r);
2278 }
2279 return(sum);
2280 }
2281
2282 /// returns [m1,m2] for two monoms, destroys nothing
2283 /// without coeffs
nc_mm_Bracket_nn(poly m1,poly m2,const ring r)2284 poly nc_mm_Bracket_nn(poly m1, poly m2, const ring r)
2285 {
2286 if (p_LmIsConstant(m1, r) || p_LmIsConstant(m1, r)) return(NULL);
2287 if (p_LmCmp(m1,m2, r)==0) return(NULL);
2288 int rN=r->N;
2289 int *M1=(int *)omAlloc0((rN+1)*sizeof(int));
2290 int *M2=(int *)omAlloc0((rN+1)*sizeof(int));
2291 int *aPREFIX=(int *)omAlloc0((rN+1)*sizeof(int));
2292 int *aSUFFIX=(int *)omAlloc0((rN+1)*sizeof(int));
2293 p_GetExpV(m1,M1, r);
2294 p_GetExpV(m2,M2, r);
2295 poly res=NULL;
2296 poly ares=NULL;
2297 poly bres=NULL;
2298 poly prefix=NULL;
2299 poly suffix=NULL;
2300 int nMin,nMax;
2301 number nTmp=NULL;
2302 int i,j,k;
2303 for (i=1;i<=rN;i++)
2304 {
2305 if (M2[i]!=0)
2306 {
2307 ares=NULL;
2308 for (j=1;j<=rN;j++)
2309 {
2310 if (M1[j]!=0)
2311 {
2312 bres=NULL;
2313 /* compute [ x_j^M1[j],x_i^M2[i] ] */
2314 if (i<j) {nMax=j; nMin=i;} else {nMax=i; nMin=j;}
2315 if ( (i==j) || ((MATELEM(r->GetNC()->COM,nMin,nMax)!=NULL) && n_IsOne(pGetCoeff(MATELEM(r->GetNC()->C,nMin,nMax)), r->cf) )) /* not (the same exp. or commuting exps)*/
2316 { bres=NULL; }
2317 else
2318 {
2319 if (i<j) { bres=gnc_uu_Mult_ww(j,M1[j],i,M2[i], r); }
2320 else bres=gnc_uu_Mult_ww(i,M2[i],j,M1[j], r);
2321 if (n_IsOne(pGetCoeff(bres), r->cf))
2322 {
2323 bres=p_LmDeleteAndNext(bres, r);
2324 }
2325 else
2326 {
2327 nTmp=n_Sub(pGetCoeff(bres),n_Init(1, r->cf), r->cf);
2328 p_SetCoeff(bres,nTmp, r); /* only lc ! */
2329 }
2330 #ifdef PDEBUG
2331 p_Test(bres, r);
2332 #endif
2333 if (i>j) bres=p_Neg(bres, r);
2334 }
2335 if (bres!=NULL)
2336 {
2337 /* now mult (prefix, bres, suffix) */
2338 memcpy(aSUFFIX, M1,(rN+1)*sizeof(int));
2339 memcpy(aPREFIX, M1,(rN+1)*sizeof(int));
2340 for (k=1;k<=j;k++) aSUFFIX[k]=0;
2341 for (k=j;k<=rN;k++) aPREFIX[k]=0;
2342 aSUFFIX[0]=0;
2343 aPREFIX[0]=0;
2344 prefix=p_One(r);
2345 suffix=p_One(r);
2346 p_SetExpV(prefix,aPREFIX, r);
2347 p_Setm(prefix, r);
2348 p_SetExpV(suffix,aSUFFIX, r);
2349 p_Setm(suffix, r);
2350 if (!p_LmIsConstant(prefix, r)) bres = gnc_p_mm_Mult(bres, prefix, r);
2351 if (!p_LmIsConstant(suffix, r)) bres = gnc_p_Mult_mm(bres, suffix, r);
2352 ares=p_Add_q(ares, bres, r);
2353 /* What to give free? */
2354 /* Do we have to free aPREFIX/aSUFFIX? it seems so */
2355 p_Delete(&prefix, r);
2356 p_Delete(&suffix, r);
2357 }
2358 }
2359 }
2360 if (ares!=NULL)
2361 {
2362 /* now mult (prefix, bres, suffix) */
2363 memcpy(aSUFFIX, M2,(rN+1)*sizeof(int));
2364 memcpy(aPREFIX, M2,(rN+1)*sizeof(int));
2365 for (k=1;k<=i;k++) aSUFFIX[k]=0;
2366 for (k=i;k<=rN;k++) aPREFIX[k]=0;
2367 aSUFFIX[0]=0;
2368 aPREFIX[0]=0;
2369 prefix=p_One(r);
2370 suffix=p_One(r);
2371 p_SetExpV(prefix,aPREFIX, r);
2372 p_Setm(prefix, r);
2373 p_SetExpV(suffix,aSUFFIX, r);
2374 p_Setm(suffix, r);
2375 bres=ares;
2376 if (!p_LmIsConstant(prefix, r)) bres = gnc_p_mm_Mult(bres, prefix, r);
2377 if (!p_LmIsConstant(suffix, r)) bres = gnc_p_Mult_mm(bres, suffix, r);
2378 res=p_Add_q(res, bres, r);
2379 p_Delete(&prefix, r);
2380 p_Delete(&suffix, r);
2381 }
2382 }
2383 }
2384 freeT(M1, rN);
2385 freeT(M2, rN);
2386 freeT(aPREFIX, rN);
2387 freeT(aSUFFIX, rN);
2388 #ifdef PDEBUG
2389 p_Test(res, r);
2390 #endif
2391 return(res);
2392 }
2393 /// returns matrix with the info on noncomm multiplication
nc_PrintMat(int a,int b,ring r,int metric)2394 matrix nc_PrintMat(int a, int b, ring r, int metric)
2395 {
2396
2397 if ( (a==b) || !rIsPluralRing(r) ) return(NULL);
2398 int i;
2399 int j;
2400 if (a>b) {j=b; i=a;}
2401 else {j=a; i=b;}
2402 /* i<j */
2403 int rN=r->N;
2404 int size=r->GetNC()->MTsize[UPMATELEM(i,j,rN)];
2405 matrix M = r->GetNC()->MT[UPMATELEM(i,j,rN)];
2406 /* return(M); */
2407 /*
2408 int sizeofres;
2409 if (metric==0)
2410 {
2411 sizeofres=sizeof(int);
2412 }
2413 if (metric==1)
2414 {
2415 sizeofres=sizeof(number);
2416 }
2417 */
2418 matrix res=mpNew(size,size);
2419 int s;
2420 int t;
2421 int length;
2422 long totdeg;
2423 poly p;
2424 for(s=1;s<=size;s++)
2425 {
2426 for(t=1;t<=size;t++)
2427 {
2428 p=MATELEM(M,s,t);
2429 if (p==NULL)
2430 {
2431 MATELEM(res,s,t)=0;
2432 }
2433 else
2434 {
2435 length = pLength(p);
2436 if (metric==0) /* length */
2437 {
2438 MATELEM(res,s,t)= p_ISet(length,r);
2439 }
2440 else if (metric==1) /* sum of deg divided by the length */
2441 {
2442 totdeg=0;
2443 while (p!=NULL)
2444 {
2445 totdeg=totdeg+p_Deg(p,r);
2446 pIter(p);
2447 }
2448 number ntd = n_Init(totdeg, r->cf);
2449 number nln = n_Init(length, r->cf);
2450 number nres= n_Div(ntd,nln, r->cf);
2451 n_Delete(&ntd, r->cf);
2452 n_Delete(&nln, r->cf);
2453 MATELEM(res,s,t)=p_NSet(nres,r);
2454 }
2455 }
2456 }
2457 }
2458 return(res);
2459 }
2460
nc_CleanUp(nc_struct * p)2461 inline void nc_CleanUp(nc_struct* p)
2462 {
2463 assume(p != NULL);
2464 omFreeSize((ADDRESS)p,sizeof(nc_struct));
2465 }
2466
nc_CleanUp(ring r)2467 inline void nc_CleanUp(ring r)
2468 {
2469 /* small CleanUp of r->GetNC() */
2470 assume(r != NULL);
2471 nc_CleanUp(r->GetNC());
2472 r->GetNC() = NULL;
2473 }
2474
nc_rKill(ring r)2475 void nc_rKill(ring r)
2476 // kills the nc extension of ring r
2477 {
2478 if( r->GetNC()->GetGlobalMultiplier() != NULL )
2479 {
2480 delete r->GetNC()->GetGlobalMultiplier();
2481 r->GetNC()->GetGlobalMultiplier() = NULL;
2482 }
2483
2484 if( r->GetNC()->GetFormulaPowerMultiplier() != NULL )
2485 {
2486 delete r->GetNC()->GetFormulaPowerMultiplier();
2487 r->GetNC()->GetFormulaPowerMultiplier() = NULL;
2488 }
2489
2490
2491 int i,j;
2492 int rN=r->N;
2493 if ( rN > 1 )
2494 {
2495 for(i=1;i<rN;i++)
2496 {
2497 for(j=i+1;j<=rN;j++)
2498 {
2499 id_Delete((ideal *)&(r->GetNC()->MT[UPMATELEM(i,j,rN)]),r);
2500 }
2501 }
2502 omFreeSize((ADDRESS)r->GetNC()->MT,rN*(rN-1)/2*sizeof(matrix));
2503 omFreeSize((ADDRESS)r->GetNC()->MTsize,rN*(rN-1)/2*sizeof(int));
2504 id_Delete((ideal *)&(r->GetNC()->COM),r);
2505 }
2506 id_Delete((ideal *)&(r->GetNC()->C),r);
2507 id_Delete((ideal *)&(r->GetNC()->D),r);
2508
2509 if( rIsSCA(r) && (r->GetNC()->SCAQuotient() != NULL) )
2510 {
2511 id_Delete(&r->GetNC()->SCAQuotient(), r); // Custom SCA destructor!!!
2512 }
2513
2514
2515 nc_CleanUp(r);
2516 }
2517
2518
2519 ////////////////////////////////////////////////////////////////////////////////////////////////
2520
2521 // deprecated:
2522 /* for use in getting the mult. matrix elements*/
2523 /* ring r must be a currRing! */
2524 /* for consistency, copies a poly from the comm. r->GetNC()->basering */
2525 /* to its image in NC ring */
nc_p_CopyGet(poly a,const ring r)2526 poly nc_p_CopyGet(poly a, const ring r)
2527 {
2528 #ifndef PDEBUG
2529 p_Test(a, r);
2530 #endif
2531
2532 // if (r != currRing)
2533 // {
2534 //#ifdef PDEBUF
2535 // WerrorS("nc_p_CopyGet call not in currRing");
2536 //#endif
2537 // return(NULL);
2538 // }
2539 return(p_Copy(a,r));
2540 }
2541
2542 // deprecated:
2543 /* for use in defining the mult. matrix elements*/
2544 /* ring r must be a currRing! */
2545 /* for consistency, puts a polynomial from the NC ring */
2546 /* to its presentation in the comm. r->GetNC()->basering */
nc_p_CopyPut(poly a,const ring r)2547 poly nc_p_CopyPut(poly a, const ring r)
2548 {
2549 #ifndef PDEBUG
2550 p_Test(a, r);
2551 #endif
2552
2553 // if (r != currRing)
2554 // {
2555 //#ifdef PDEBUF
2556 // WerrorS("nc_p_CopyGet call not in currRing");
2557 //#endif
2558 // return(NULL);
2559 // }
2560
2561 return(p_Copy(a,r));
2562 }
2563
2564 /* returns TRUE if there were errors */
2565 /* checks whether product of vars from PolyVar defines */
2566 /* an admissible subalgebra of r */
2567 /* r is indeed currRing and assumed to be noncomm. */
nc_CheckSubalgebra(poly PolyVar,ring r)2568 BOOLEAN nc_CheckSubalgebra(poly PolyVar, ring r)
2569 {
2570 // ring save = currRing;
2571 // int WeChangeRing = 0;
2572 // if (currRing != r)
2573 // rChangeCurrRing(r);
2574 // WeChangeRing = 1;
2575 // }
2576 int rN=r->N;
2577 int *ExpVar=(int*)omAlloc0((rN+1)*sizeof(int));
2578 int *ExpTmp=(int*)omAlloc0((rN+1)*sizeof(int));
2579 p_GetExpV(PolyVar, ExpVar, r);
2580 int i; int j; int k;
2581 poly test=NULL;
2582 int OK=1;
2583 for (i=1; i<rN; i++)
2584 {
2585 if (ExpVar[i]==0) /* i.e. not in PolyVar */
2586 {
2587 for (j=i+1; j<=rN; j++)
2588 {
2589 if (ExpVar[j]==0)
2590 {
2591 test = MATELEM(r->GetNC()->D,i,j);
2592 while (test!=NULL)
2593 {
2594 p_GetExpV(test, ExpTmp, r);
2595 OK=1;
2596 for (k=1;k<=rN;k++)
2597 {
2598 if (ExpTmp[k]!=0)
2599 {
2600 if (ExpVar[k]!=0) OK=0;
2601 }
2602 }
2603 if (!OK)
2604 {
2605 // if ( WeChangeRing )
2606 // rChangeCurrRing(save);
2607 return(TRUE);
2608 }
2609 pIter(test);
2610 }
2611 }
2612 }
2613 }
2614 }
2615 freeT(ExpVar,rN);
2616 freeT(ExpTmp,rN);
2617 // if ( WeChangeRing )
2618 // rChangeCurrRing(save);
2619 return(FALSE);
2620 }
2621
2622
2623 /* returns TRUE if there were errors */
2624 /* checks whether the current ordering */
2625 /* is admissible for r and D == r->GetNC()->D */
2626 /* to be executed in a currRing */
gnc_CheckOrdCondition(matrix D,ring r)2627 BOOLEAN gnc_CheckOrdCondition(matrix D, ring r)
2628 {
2629 /* analyze D: an upper triangular matrix of polys */
2630 /* check the ordering condition for D */
2631 // ring save = currRing;
2632 // int WeChangeRing = 0;
2633 // if (r != currRing)
2634 // {
2635 // rChangeCurrRing(r);
2636 // WeChangeRing = 1;
2637 // }
2638 poly p,q;
2639 int i,j;
2640 int report = 0;
2641 for(i=1; i<r->N; i++)
2642 {
2643 for(j=i+1; j<=r->N; j++)
2644 {
2645 p = nc_p_CopyGet(MATELEM(D,i,j),r);
2646 if ( p != NULL)
2647 {
2648 q = p_One(r);
2649 p_SetExp(q,i,1,r);
2650 p_SetExp(q,j,1,r);
2651 p_Setm(q,r);
2652 if (p_LmCmp(q,p,r) != 1) /* i.e. lm(p)==xy < lm(q)==D_ij */
2653 {
2654 Werror("Bad ordering at %d,%d\n",i,j);
2655 #if 0 /*Singularg should not differ from Singular except in error case*/
2656 p_Write(p,r);
2657 p_Write(q,r);
2658 #endif
2659 report = 1;
2660 }
2661 p_Delete(&q,r);
2662 p_Delete(&p,r);
2663 p = NULL;
2664 }
2665 }
2666 }
2667 // if ( WeChangeRing )
2668 // rChangeCurrRing(save);
2669 return(report);
2670 }
2671
2672
2673
2674 BOOLEAN gnc_InitMultiplication(ring r, bool bSetupQuotient = false); // just for a moment
2675
2676 /// returns TRUE if there were errors
2677 /// analyze inputs, check them for consistency
2678 /// detects nc_type, DO NOT initialize multiplication but call for it at the end
2679 /// checks the ordering condition and evtl. NDC
2680 /// NOTE: all the data belong to the curr,
2681 /// we change r which may be the same ring, and must have the same representation!
nc_CallPlural(matrix CCC,matrix DDD,poly CCN,poly DDN,ring r,bool bSetupQuotient,bool bCopyInput,bool bBeQuiet,ring curr,bool dummy_ring)2682 BOOLEAN nc_CallPlural(matrix CCC, matrix DDD,
2683 poly CCN, poly DDN,
2684 ring r,
2685 bool bSetupQuotient, bool bCopyInput, bool bBeQuiet,
2686 ring curr, bool dummy_ring /*=false*/)
2687 {
2688 assume( r != NULL );
2689 assume( curr != NULL );
2690
2691 if( !bSetupQuotient)
2692 assume( (r->qideal == NULL) ); // The basering must NOT be a qring!??
2693
2694 assume( rSamePolyRep(r, curr) || bCopyInput ); // wrong assumption?
2695
2696
2697 if( r->N == 1 ) // clearly commutative!!!
2698 {
2699 assume(
2700 ( (CCC != NULL) && (MATCOLS(CCC) == 1) && (MATROWS(CCC) == 1) && (MATELEM(CCC,1,1) == NULL) ) ||
2701 ( (CCN == NULL) )
2702 );
2703
2704 assume(
2705 ( (DDD != NULL) && (MATCOLS(DDD) == 1) && (MATROWS(DDD) == 1) && (MATELEM(DDD,1,1) == NULL) ) ||
2706 ( (DDN == NULL) )
2707 );
2708 if(!dummy_ring)
2709 {
2710 WarnS("commutative ring with 1 variable");
2711 return FALSE;
2712 }
2713 }
2714
2715 // there must be:
2716 assume( (CCC != NULL) != (CCN != NULL) ); // exactly one data about coeffs (C).
2717 assume( !((DDD != NULL) && (DDN != NULL)) ); // at most one data about tails (D).
2718
2719 #if OUTPUT
2720 if( CCC != NULL )
2721 {
2722 PrintS("nc_CallPlural(), Input data, CCC: \n");
2723 iiWriteMatrix(CCC, "C", 2, curr, 4);
2724 }
2725 if( DDD != NULL )
2726 {
2727 PrintS("nc_CallPlural(), Input data, DDD: \n");
2728 iiWriteMatrix(DDD, "D", 2, curr, 4);
2729 }
2730 #endif
2731
2732
2733 #ifndef SING_NDEBUG
2734 if (CCC!=NULL) id_Test((ideal)CCC, curr);
2735 if (DDD!=NULL) id_Test((ideal)DDD, curr);
2736 p_Test(CCN, curr);
2737 p_Test(DDN, curr);
2738 #endif
2739
2740 if( (!bBeQuiet) && (r->GetNC() != NULL) )
2741 WarnS("going to redefine the algebra structure");
2742
2743 matrix CC = NULL;
2744 poly CN = NULL;
2745 matrix C; bool bCnew = false;
2746
2747 matrix DD = NULL;
2748 poly DN = NULL;
2749 matrix D; bool bDnew = false;
2750
2751 number nN, pN, qN;
2752
2753 bool IsSkewConstant = false, tmpIsSkewConstant;
2754 int i, j;
2755
2756 nc_type nctype = nc_undef;
2757
2758 //////////////////////////////////////////////////////////////////
2759 // check the correctness of arguments, without any real chagnes!!!
2760
2761
2762
2763 // check C
2764 if ((CCC != NULL) && ( (MATCOLS(CCC)==1) || MATROWS(CCC)==1 ) )
2765 {
2766 CN = MATELEM(CCC,1,1);
2767 }
2768 else
2769 {
2770 if ((CCC != NULL) && ( (MATCOLS(CCC)!=r->N) || (MATROWS(CCC)!=r->N) ))
2771 {
2772 Werror("Square %d x %d matrix expected", r->N, r->N);
2773 return TRUE;
2774 }
2775 }
2776 if (( CCC != NULL) && (CC == NULL)) CC = CCC; // mp_Copy(CCC, ?); // bug!?
2777 if (( CCN != NULL) && (CN == NULL)) CN = CCN;
2778
2779 // check D
2780 if ((DDD != NULL) && ( (MATCOLS(DDD)==1) || MATROWS(DDD)==1 ) )
2781 {
2782 DN = MATELEM(DDD,1,1);
2783 }
2784 else
2785 {
2786 if ((DDD != NULL) && ( (MATCOLS(DDD)!=r->N) || (MATROWS(DDD)!=r->N) ))
2787 {
2788 Werror("Square %d x %d matrix expected",r->N,r->N);
2789 return TRUE;
2790 }
2791 }
2792
2793 if (( DDD != NULL) && (DD == NULL)) DD = DDD; // mp_Copy(DDD, ?); // ???
2794 if (( DDN != NULL) && (DN == NULL)) DN = DDN;
2795
2796 // further checks and some analysis:
2797 // all data in 'curr'!
2798 if (CN != NULL) /* create matrix C = CN * Id */
2799 {
2800 if (!p_IsConstant(CN,curr))
2801 {
2802 WerrorS("Incorrect input : non-constants are not allowed as coefficients (first argument)");
2803 return TRUE;
2804 }
2805 assume(p_IsConstant(CN,curr));
2806
2807 nN = pGetCoeff(CN);
2808 if (n_IsZero(nN, curr->cf))
2809 {
2810 WerrorS("Incorrect input : zero coefficients are not allowed");
2811 return TRUE;
2812 }
2813
2814 if (n_IsOne(nN, curr->cf))
2815 nctype = nc_lie;
2816 else
2817 nctype = nc_general;
2818
2819 IsSkewConstant = true;
2820
2821 C = mpNew(r->N,r->N); // ring independent!
2822 bCnew = true;
2823
2824 for(i=1; i<r->N; i++)
2825 for(j=i+1; j<=r->N; j++)
2826 MATELEM(C,i,j) = prCopyR_NoSort(CN, curr, r); // nc_p_CopyPut(CN, r); // copy CN from curr into r
2827
2828 #ifndef SING_NDEBUG
2829 id_Test((ideal)C, r);
2830 #endif
2831
2832 }
2833 else if ( (CN == NULL) && (CC != NULL) ) /* copy matrix C */
2834 {
2835 /* analyze C */
2836
2837 BOOLEAN pN_set=FALSE;
2838 pN = n_Init(0,curr->cf);
2839
2840 if( r->N > 1 )
2841 if ( MATELEM(CC,1,2) != NULL )
2842 {
2843 if (!pN_set) n_Delete(&pN,curr->cf); // free initial nInit(0)
2844 pN = p_GetCoeff(MATELEM(CC,1,2), curr);
2845 pN_set=TRUE;
2846 }
2847
2848 tmpIsSkewConstant = true;
2849
2850 for(i=1; i<r->N; i++)
2851 for(j=i+1; j<=r->N; j++)
2852 {
2853 if (MATELEM(CC,i,j) == NULL)
2854 qN = NULL;
2855 else
2856 {
2857 if (!p_IsConstant(MATELEM(CC,i,j),curr))
2858 {
2859 Werror("Incorrect input : non-constants are not allowed as coefficients (first argument at [%d, %d])", i, j);
2860 return TRUE;
2861 }
2862 assume(p_IsConstant(MATELEM(CC,i,j),curr));
2863 qN = p_GetCoeff(MATELEM(CC,i,j),curr);
2864 }
2865
2866 if ( qN == NULL ) /* check the consistency: Cij!=0 */
2867 // find also illegal pN
2868 {
2869 WerrorS("Incorrect input : matrix of coefficients contains zeros in the upper triangle");
2870 return TRUE;
2871 }
2872
2873 if (!n_Equal(pN, qN, curr->cf)) tmpIsSkewConstant = false;
2874 }
2875
2876 if( bCopyInput )
2877 {
2878 C = mp_Copy(CC, curr, r); // Copy C into r!!!???
2879 #ifndef SING_NDEBUG
2880 id_Test((ideal)C, r);
2881 #endif
2882 bCnew = true;
2883 }
2884 else
2885
2886 C = CC;
2887
2888 IsSkewConstant = tmpIsSkewConstant;
2889
2890 if ( tmpIsSkewConstant && n_IsOne(pN, curr->cf) )
2891 nctype = nc_lie;
2892 else
2893 nctype = nc_general;
2894 if (!pN_set) n_Delete(&pN,curr->cf); // free initial nInit(0)
2895 }
2896
2897 /* initialition of the matrix D */
2898 if ( DD == NULL ) /* we treat DN only (it could also be NULL) */
2899 {
2900 D = mpNew(r->N,r->N); bDnew = true;
2901
2902 if (DN == NULL)
2903 {
2904 if ( (nctype == nc_lie) || (nctype == nc_undef) )
2905 nctype = nc_comm; /* it was nc_skew earlier */
2906 else /* nc_general, nc_skew */
2907 nctype = nc_skew;
2908 }
2909 else /* DN != NULL */
2910 for(i=1; i<r->N; i++)
2911 for(j=i+1; j<=r->N; j++)
2912 MATELEM(D,i,j) = prCopyR_NoSort(DN, curr, r); // project DN into r->GetNC()->basering!
2913 #ifndef SING_NDEBUG
2914 id_Test((ideal)D, r);
2915 #endif
2916 }
2917 else /* DD != NULL */
2918 {
2919 bool b = true; // DD == null ?
2920
2921 for(int i = 1; (i < r->N) && b; i++)
2922 for(int j = i+1; (j <= r->N) && b; j++)
2923 if (MATELEM(DD, i, j) != NULL)
2924 {
2925 b = false;
2926 break;
2927 }
2928
2929 if (b) // D == NULL!!!
2930 {
2931 if ( (nctype == nc_lie) || (nctype == nc_undef) )
2932 nctype = nc_comm; /* it was nc_skew earlier */
2933 else /* nc_general, nc_skew */
2934 nctype = nc_skew;
2935 }
2936
2937 if( bCopyInput )
2938 {
2939 D = mp_Copy(DD, curr, r); // Copy DD into r!!!
2940 #ifndef SING_NDEBUG
2941 id_Test((ideal)D, r);
2942 #endif
2943 bDnew = true;
2944 }
2945 else
2946 D = DD;
2947 }
2948
2949 assume( C != NULL );
2950 assume( D != NULL );
2951
2952 #if OUTPUT
2953 PrintS("nc_CallPlural(), Computed data, C: \n");
2954 iiWriteMatrix(C, "C", 2, r, 4);
2955
2956 PrintS("nc_CallPlural(), Computed data, D: \n");
2957 iiWriteMatrix(D, "D", 2, r, 4);
2958
2959 Print("\nTemporary: type = %d, IsSkewConstant = %d\n", nctype, IsSkewConstant);
2960 #endif
2961
2962
2963 // check the ordering condition for D (both matrix and poly cases):
2964 if ( gnc_CheckOrdCondition(D, r) )
2965 {
2966 if( bCnew ) mp_Delete( &C, r );
2967 if( bDnew ) mp_Delete( &D, r );
2968
2969 WerrorS("Matrix of polynomials violates the ordering condition");
2970 return TRUE;
2971 }
2972
2973 // okay now we are ready for this!!!
2974
2975 // create new non-commutative structure
2976 nc_struct *nc_new = (nc_struct *)omAlloc0(sizeof(nc_struct));
2977
2978 ncRingType(nc_new, nctype);
2979
2980 nc_new->C = C; // if C and D were given by matrices at the beginning they are in r
2981 nc_new->D = D; // otherwise they should be in r->GetNC()->basering(polynomial * Id_{N})
2982
2983 nc_new->IsSkewConstant = (IsSkewConstant?1:0);
2984
2985 // Setup new NC structure!!!
2986 if (r->GetNC() != NULL)
2987 {
2988 #ifndef SING_NDEBUG
2989 WarnS("Changing the NC-structure of an existing NC-ring!!!");
2990 #endif
2991 nc_rKill(r);
2992 }
2993
2994 r->GetNC() = nc_new;
2995
2996 r->ext_ref=NULL;
2997
2998 return gnc_InitMultiplication(r, bSetupQuotient);
2999 }
3000
3001 //////////////////////////////////////////////////////////////////////////////
3002
nc_rCopy(ring res,const ring r,bool bSetupQuotient)3003 bool nc_rCopy(ring res, const ring r, bool bSetupQuotient)
3004 {
3005 if (nc_CallPlural(r->GetNC()->C, r->GetNC()->D, NULL, NULL, res, bSetupQuotient, true, true, r))
3006 {
3007 WarnS("Error occurred while coping/setuping the NC structure!"); // No reaction!???
3008 return true; // error
3009 }
3010
3011 return false;
3012 }
3013
3014 //////////////////////////////////////////////////////////////////////////////
gnc_InitMultiplication(ring r,bool bSetupQuotient)3015 BOOLEAN gnc_InitMultiplication(ring r, bool bSetupQuotient)
3016 {
3017 /* returns TRUE if there were errors */
3018 /* initialize the multiplication: */
3019 /* r->GetNC()->MTsize, r->GetNC()->MT, r->GetNC()->COM, */
3020 /* and r->GetNC()->IsSkewConstant for the skew case */
3021 if (rVar(r)==1)
3022 {
3023 ncRingType(r, nc_comm);
3024 r->GetNC()->IsSkewConstant=1;
3025 return FALSE;
3026 }
3027
3028 // ring save = currRing;
3029 // int WeChangeRing = 0;
3030
3031 // if (currRing!=r)
3032 // {
3033 // rChangeCurrRing(r);
3034 // WeChangeRing = 1;
3035 // }
3036 // assume( (currRing == r)
3037 // && (currRing->GetNC()!=NULL) ); // otherwise we cannot work with all these matrices!
3038
3039 int i,j;
3040 r->GetNC()->MT = (matrix *)omAlloc0((r->N*(r->N-1))/2*sizeof(matrix));
3041 r->GetNC()->MTsize = (int *)omAlloc0((r->N*(r->N-1))/2*sizeof(int));
3042 id_Test((ideal)r->GetNC()->C, r);
3043 matrix COM = mp_Copy(r->GetNC()->C, r);
3044 poly p,q;
3045 short DefMTsize=7;
3046 int IsNonComm=0;
3047 // bool tmpIsSkewConstant = false;
3048
3049 for(i=1; i<r->N; i++)
3050 {
3051 for(j=i+1; j<=r->N; j++)
3052 {
3053 if ( MATELEM(r->GetNC()->D,i,j) == NULL ) /* quasicommutative case */
3054 {
3055 /* 1x1 mult.matrix */
3056 r->GetNC()->MTsize[UPMATELEM(i,j,r->N)] = 1;
3057 r->GetNC()->MT[UPMATELEM(i,j,r->N)] = mpNew(1,1);
3058 }
3059 else /* pure noncommutative case */
3060 {
3061 /* TODO check the special multiplication properties */
3062 IsNonComm = 1;
3063 p_Delete(&(MATELEM(COM,i,j)),r);
3064 //MATELEM(COM,i,j) = NULL; // done by p_Delete
3065 r->GetNC()->MTsize[UPMATELEM(i,j,r->N)] = DefMTsize; /* default sizes */
3066 r->GetNC()->MT[UPMATELEM(i,j,r->N)] = mpNew(DefMTsize, DefMTsize);
3067 }
3068 /* set MT[i,j,1,1] to c_i_j*x_i*x_j + D_i_j */
3069 p = p_One(r);
3070 if (MATELEM(r->GetNC()->C,i,j)!=NULL)
3071 p_SetCoeff(p,n_Copy(pGetCoeff(MATELEM(r->GetNC()->C,i,j)),r->cf),r);
3072 p_SetExp(p,i,1,r);
3073 p_SetExp(p,j,1,r);
3074 p_Setm(p,r);
3075 p_Test(MATELEM(r->GetNC()->D,i,j),r);
3076 q = nc_p_CopyGet(MATELEM(r->GetNC()->D,i,j),r);
3077 p = p_Add_q(p,q,r);
3078 MATELEM(r->GetNC()->MT[UPMATELEM(i,j,r->N)],1,1) = nc_p_CopyPut(p,r);
3079 p_Delete(&p,r);
3080 // p = NULL;// done by p_Delete
3081 }
3082 }
3083 if (ncRingType(r)==nc_undef)
3084 {
3085 if (IsNonComm==1)
3086 {
3087 // assume(pN!=NULL);
3088 // if ((tmpIsSkewConstant==1) && (nIsOne(pGetCoeff(pN)))) r->GetNC()->type=nc_lie;
3089 // else r->GetNC()->type=nc_general;
3090 }
3091 if (IsNonComm==0)
3092 {
3093 ncRingType(r, nc_skew); // TODO: check whether it is commutative
3094 r->GetNC()->IsSkewConstant = 0; // true; //tmpIsSkewConstant; // BUG???
3095 } else
3096 assume( FALSE );
3097 }
3098 r->GetNC()->COM=COM;
3099
3100 nc_p_ProcsSet(r, r->p_Procs);
3101
3102 if(bSetupQuotient) // Test me!!!
3103 nc_SetupQuotient(r, NULL, false); // no copy!
3104
3105
3106 // if (save != currRing)
3107 // rChangeCurrRing(save);
3108
3109 return FALSE;
3110 }
3111
3112
3113 // set pProcs for r and global variable p_Procs as for general non-commutative algebras.
3114 static inline
gnc_p_ProcsSet(ring rGR,p_Procs_s * p_Procs)3115 void gnc_p_ProcsSet(ring rGR, p_Procs_s* p_Procs)
3116 {
3117 // "commutative"
3118 p_Procs->p_Mult_mm = rGR->p_Procs->p_Mult_mm = gnc_p_Mult_mm;
3119 p_Procs->pp_Mult_mm = rGR->p_Procs->pp_Mult_mm = gnc_pp_Mult_mm;
3120 p_Procs->p_Minus_mm_Mult_qq = rGR->p_Procs->p_Minus_mm_Mult_qq = nc_p_Minus_mm_Mult_qq;
3121
3122 // non-commutaitve multiplication by monomial from the left
3123 p_Procs->p_mm_Mult = gnc_p_mm_Mult;
3124 p_Procs->pp_mm_Mult = gnc_pp_mm_Mult;
3125
3126 #if 0
3127 // Previous Plural's implementation...
3128 rGR->GetNC()->p_Procs.SPoly = gnc_CreateSpolyOld;
3129 rGR->GetNC()->p_Procs.ReduceSPoly = gnc_ReduceSpolyOld;
3130
3131 rGR->GetNC()->p_Procs.BucketPolyRed = gnc_kBucketPolyRedOld;
3132 rGR->GetNC()->p_Procs.BucketPolyRed_Z= gnc_kBucketPolyRed_ZOld;
3133 #else
3134 // A bit cleaned up and somewhat rewritten functions...
3135 rGR->GetNC()->p_Procs.SPoly = gnc_CreateSpolyNew;
3136 rGR->GetNC()->p_Procs.ReduceSPoly = gnc_ReduceSpolyNew;
3137
3138 rGR->GetNC()->p_Procs.BucketPolyRed_NF= gnc_kBucketPolyRedNew;
3139 rGR->GetNC()->p_Procs.BucketPolyRed_Z= gnc_kBucketPolyRed_ZNew;
3140 #endif
3141
3142 // warning: ISO C++ forbids casting between pointer-to-function and pointer-to-object?
3143 if (rHasLocalOrMixedOrdering(rGR))
3144 rGR->GetNC()->p_Procs.GB = cast_A_to_vptr(gnc_gr_mora);
3145 else
3146 rGR->GetNC()->p_Procs.GB = cast_A_to_vptr(gnc_gr_bba);
3147
3148 /////////// rGR->GetNC()->p_Procs.GB = gnc_gr_bba; // bba even for local case!
3149 // old /// r->GetNC()->GB() = gnc_gr_bba;
3150 // rGR->GetNC()->p_Procs.GlobalGB = gnc_gr_bba;
3151 // rGR->GetNC()->p_Procs.LocalGB = gnc_gr_mora;
3152 // const ring save = currRing; if( save != r ) rChangeCurrRing(r);
3153 // ideal res = gnc_gr_bba(F, Q, w, hilb, strat/*, r*/);
3154 // if( save != r ) rChangeCurrRing(save); return (res);
3155
3156
3157 #if 0
3158 // Old Stuff
3159 p_Procs->p_Mult_mm = gnc_p_Mult_mm;
3160 _p_procs->p_Mult_mm = gnc_p_Mult_mm;
3161
3162 p_Procs->pp_Mult_mm = gnc_pp_Mult_mm;
3163 _p_procs->pp_Mult_mm = gnc_pp_Mult_mm;
3164
3165 p_Procs->p_Minus_mm_Mult_qq = NULL; // gnc_p_Minus_mm_Mult_qq_ign;
3166 _p_procs->p_Minus_mm_Mult_qq= NULL; // gnc_p_Minus_mm_Mult_qq_ign;
3167
3168 r->GetNC()->mmMultP() = gnc_mm_Mult_p;
3169 r->GetNC()->mmMultPP() = gnc_mm_Mult_pp;
3170
3171 r->GetNC()->SPoly() = gnc_CreateSpoly;
3172 r->GetNC()->ReduceSPoly() = gnc_ReduceSpoly;
3173
3174 #endif
3175 }
3176
3177
3178 // set pProcs table for rGR and global variable p_Procs
nc_p_ProcsSet(ring rGR,p_Procs_s * p_Procs)3179 void nc_p_ProcsSet(ring rGR, p_Procs_s* p_Procs)
3180 {
3181 assume(rIsPluralRing(rGR));
3182 assume(p_Procs!=NULL);
3183
3184 gnc_p_ProcsSet(rGR, p_Procs);
3185
3186 if(rIsSCA(rGR) && ncExtensions(SCAMASK) )
3187 {
3188 sca_p_ProcsSet(rGR, p_Procs);
3189 }
3190
3191 if( ncExtensions(NOPLURALMASK) )
3192 ncInitSpecialPairMultiplication(rGR);
3193
3194 if(!rIsSCA(rGR) && !ncExtensions(NOFORMULAMASK))
3195 ncInitSpecialPowersMultiplication(rGR);
3196
3197 }
3198
3199
3200 /// substitute the n-th variable by e in p
3201 /// destroy p
3202 /// e is not a constant
nc_pSubst(poly p,int n,poly e,const ring r)3203 poly nc_pSubst(poly p, int n, poly e, const ring r)
3204 {
3205 int rN = r->N;
3206 int *PRE = (int *)omAlloc0((rN+1)*sizeof(int));
3207 int *SUF = (int *)omAlloc0((rN+1)*sizeof(int));
3208 int i,pow;
3209 number C;
3210 poly suf,pre;
3211 poly res = NULL;
3212 poly out = NULL;
3213 while ( p!= NULL )
3214 {
3215 C = p_GetCoeff(p, r);
3216 p_GetExpV(p, PRE, r); /* faster splitting? */
3217 pow = PRE[n]; PRE[n]=0;
3218 res = NULL;
3219 if (pow!=0)
3220 {
3221 for (i=n+1; i<=rN; i++)
3222 {
3223 SUF[i] = PRE[i];
3224 PRE[i] = 0;
3225 }
3226 res = p_Power(p_Copy(e, r),pow, r);
3227 /* multiply with prefix */
3228 pre = p_One(r);
3229 p_SetExpV(pre,PRE, r);
3230 p_Setm(pre, r);
3231 res = nc_mm_Mult_p(pre,res, r);
3232 /* multiply with suffix */
3233 suf = p_One(r);
3234 p_SetExpV(suf,SUF, r);
3235 p_Setm(suf, r);
3236 res = p_Mult_mm(res,suf, r);
3237 res = __p_Mult_nn(res,C, r);
3238 p_SetComp(res,PRE[0], r);
3239 }
3240 else /* pow==0 */
3241 {
3242 res = p_Head(p, r);
3243 }
3244 p = p_LmDeleteAndNext(p, r);
3245 out = p_Add_q(out,res, r);
3246 }
3247 freeT(PRE,rN);
3248 freeT(SUF,rN);
3249 return(out);
3250 }
3251
3252
3253 // creates a commutative nc extension; "converts" comm.ring to a Plural ring
nc_rCreateNCcomm(ring r)3254 ring nc_rCreateNCcomm(ring r)
3255 {
3256 if (rIsPluralRing(r)) return r;
3257
3258 ring rr = rCopy(r);
3259
3260 matrix C = mpNew(rr->N,rr->N); // ring-independent!?!
3261 matrix D = mpNew(rr->N,rr->N);
3262
3263 for(int i=1; i<rr->N; i++)
3264 for(int j=i+1; j<=rr->N; j++)
3265 MATELEM(C,i,j) = p_One(rr);
3266
3267 if (nc_CallPlural(C, D, NULL, NULL, rr, false, true, false, rr, TRUE)) // TODO: what about quotient ideal?
3268 WarnS("Error initializing multiplication!"); // No reaction!???
3269
3270 return rr;
3271 }
3272
3273 /* NOT USED ANYMORE: replaced by maFindPerm in ring.cc */
3274 /* for use with embeddings: currRing is a sum of smaller rings */
3275 /* and srcRing is one of such smaller rings */
3276 /* shift defines the position of a subring in srcRing */
3277 /* par_shift defines the position of a subfield in basefield of CurrRing */
p_CopyEmbed(poly p,ring srcRing,int shift,int,ring dstRing)3278 poly p_CopyEmbed(poly p, ring srcRing, int shift, int /*par_shift*/, ring dstRing)
3279 {
3280 if (dstRing == srcRing)
3281 {
3282 return(p_Copy(p,dstRing));
3283 }
3284 nMapFunc nMap=n_SetMap(srcRing->cf, dstRing->cf);
3285 poly q;
3286 // if ( nMap == nCopy)
3287 // {
3288 // q = prCopyR(p,srcRing);
3289 // }
3290 // else
3291 {
3292 int *perm = (int *)omAlloc0((rVar(srcRing)+1)*sizeof(int));
3293 int *par_perm = (int *)omAlloc0((rPar(srcRing)+1)*sizeof(int));
3294 // int *par_perm = (int *)omAlloc0((rPar(srcRing)+1)*sizeof(int));
3295 int i;
3296 // if (srcRing->P > 0)
3297 // {
3298 // for (i=0; i<srcRing->P; i++)
3299 // par_perm[i]=-i;
3300 // }
3301 if ((shift<0) || (shift > rVar(srcRing))) // ???
3302 {
3303 WerrorS("bad shifts in p_CopyEmbed");
3304 return(0);
3305 }
3306 for (i=1; i<= srcRing->N; i++)
3307 {
3308 perm[i] = shift+i;
3309 }
3310 q = p_PermPoly(p,perm,srcRing, dstRing, nMap,par_perm, rPar(srcRing));
3311 }
3312 return(q);
3313 }
3314
rIsLikeOpposite(ring rBase,ring rCandidate)3315 BOOLEAN rIsLikeOpposite(ring rBase, ring rCandidate)
3316 {
3317 /* the same basefield */
3318 int diagnose = TRUE;
3319 nMapFunc nMap = n_SetMap(rCandidate->cf, rBase->cf); // reverse?
3320
3321 ////// if (nMap != nCopy) diagnose = FALSE;
3322 if (nMap == NULL) diagnose = FALSE;
3323
3324
3325 /* same number of variables */
3326 if (rBase->N != rCandidate->N) diagnose = FALSE;
3327 /* nc and comm ring */
3328 if ( rIsPluralRing(rBase) != rIsPluralRing(rCandidate) ) diagnose = FALSE;
3329 /* both are qrings */
3330 /* NO CHECK, since it is used in building opposite qring */
3331 /* if ( ((rBase->qideal != NULL) && (rCandidate->qideal == NULL)) */
3332 /* || ((rBase->qideal == NULL) && (rCandidate->qideal != NULL)) ) */
3333 /* diagnose = FALSE; */
3334 /* TODO: varnames are e->E etc */
3335 return diagnose;
3336 }
3337
3338
3339
3340
3341 /// opposes a vector p from Rop to currRing (dst!)
pOppose(ring Rop,poly p,const ring dst)3342 poly pOppose(ring Rop, poly p, const ring dst)
3343 {
3344 /* the simplest case:*/
3345 if ( Rop == dst ) return(p_Copy(p, dst));
3346 /* check Rop == rOpposite(currRing) */
3347
3348
3349 if ( !rIsLikeOpposite(dst, Rop) )
3350 {
3351 WarnS("an opposite ring should be used");
3352 return NULL;
3353 }
3354
3355 nMapFunc nMap = n_SetMap(Rop->cf, dst->cf); // reverse?
3356
3357 /* nMapFunc nMap = nSetMap(Rop);*/
3358 /* since we know that basefields coinside! */
3359
3360 // coinside???
3361
3362 int *perm=(int *)omAlloc0((Rop->N+1)*sizeof(int));
3363 if (!p_IsConstant(p, Rop))
3364 {
3365 /* we know perm exactly */
3366 int i;
3367 for(i=1; i<=Rop->N; i++)
3368 {
3369 perm[i] = Rop->N+1-i;
3370 }
3371 }
3372 poly res = p_PermPoly(p, perm, Rop, dst, nMap);
3373 omFreeSize((ADDRESS)perm,(Rop->N+1)*sizeof(int));
3374
3375 p_Test(res, dst);
3376
3377 return res;
3378 }
3379
3380 /// opposes a module I from Rop to currRing(dst)
idOppose(ring Rop,ideal I,const ring dst)3381 ideal idOppose(ring Rop, ideal I, const ring dst)
3382 {
3383 /* the simplest case:*/
3384 if ( Rop == dst ) return id_Copy(I, dst);
3385
3386 /* check Rop == rOpposite(currRing) */
3387 if (!rIsLikeOpposite(dst, Rop))
3388 {
3389 WarnS("an opposite ring should be used");
3390 return NULL;
3391 }
3392 int i;
3393 ideal idOp = idInit(I->ncols, I->rank);
3394 for (i=0; i< (I->ncols)*(I->nrows); i++)
3395 {
3396 idOp->m[i] = pOppose(Rop,I->m[i], dst);
3397 }
3398 id_Test(idOp, dst);
3399 return idOp;
3400 }
3401
3402
nc_SetupQuotient(ring rGR,const ring rG,bool bCopy)3403 bool nc_SetupQuotient(ring rGR, const ring rG, bool bCopy)
3404 {
3405 if( rGR->qideal == NULL )
3406 return false; // no quotient = no work! done!? What about factors of SCA?
3407
3408 bool ret = true;
3409 // currently only super-commutative extension deals with factors.
3410
3411 if( ncExtensions(SCAMASK) )
3412 {
3413 bool sca_ret = sca_SetupQuotient(rGR, rG, bCopy);
3414
3415 if(sca_ret) // yes it was dealt with!
3416 ret = false;
3417 }
3418
3419 if( bCopy )
3420 {
3421 assume(rIsPluralRing(rGR) == rIsPluralRing(rG));
3422 assume((rGR->qideal==NULL) == (rG->qideal==NULL));
3423 assume(rIsSCA(rGR) == rIsSCA(rG));
3424 assume(ncRingType(rGR) == ncRingType(rG));
3425 }
3426
3427 return ret;
3428 }
3429
3430
3431
3432 // int Commutative_Context(ring r, leftv expression)
3433 // /* returns 1 if expression consists */
3434 // /* of commutative elements */
3435 // {
3436 // /* crucial: poly -> ideal, module, matrix */
3437 // }
3438
3439 // int Comm_Context_Poly(ring r, poly p)
3440 // {
3441 // poly COMM=r->GetNC()->COMM;
3442 // poly pp=pOne();
3443 // memset(pp->exp,0,r->ExpL_Size*sizeof(long));
3444 // while (p!=NULL)
3445 // {
3446 // for (i=0;i<=r->ExpL_Size;i++)
3447 // {
3448 // if ((p->exp[i]) && (pp->exp[i])) return(FALSE);
3449 // /* nonzero exponent of non-comm variable */
3450 // }
3451 // pIter(p);
3452 // }
3453 // return(TRUE);
3454 // }
3455
3456 #endif
3457
3458
3459
3460