source: git/kernel/gring.cc @ 45d41f

spielwiese
Last change on this file since 45d41f was 45d41f, checked in by Michael Brickenstein <bricken@…>, 19 years ago
*bricken: memory in kBucketPolyRed_Z fixed, optimized git-svn-id: file:///usr/local/Singular/svn/trunk@8433 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 56.8 KB
Line 
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.27 2005-07-18 12:41:55 bricken 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 : */
31poly 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
36poly 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
41poly 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 */
47poly 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--------------------------
57poly _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
82poly 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#ifdef PDEBUG
140        const char* s;
141        if (side==1) s="nc_p_Mult_mm";
142        else s="nc_mm_Mult_p";
143        Print("%s: exponent mismatch %d and %d\n",s,expP,expM);
144#endif
145        expOut=0;
146      }
147    }
148    p_GetExpV(p,P,r);
149    cP=p_GetCoeff(p,r);
150    cOut=n_Mult(cP,cM,r);
151    if (side==1)
152    {
153      v = nc_mm_Mult_nn(P, M, r);     
154    }
155    else
156    {
157      v = nc_mm_Mult_nn(M, P, r);
158    }
159    v = p_Mult_nn(v,cOut,r);
160    p_SetCompP(v,expOut,r);
161    if (UseBuckets) sBucket_Add_p(bu_out,v,pLength(v));
162    else out = p_Add_q(out,v,r);
163    p_DeleteLm(&p,r);
164  }
165  freeT(P,rN);
166  freeT(M,rN);
167  if (UseBuckets)
168  {
169    out = NULL;
170    int len = pLength(out);
171    sBucketDestroyAdd(bu_out, &out, &len);
172  }
173#ifdef PDEBUG
174  p_Test(out,r);
175#endif
176  return(out);
177}
178
179poly nc_mm_Mult_nn(int *F0, int *G0, const ring r)
180/* destroys nothing, no coeffs and exps */
181{
182  poly out=NULL;
183  int i,j;
184  int iF,jG,iG;
185  int rN=r->N;
186  int ExpSize=(((rN+1)*sizeof(int)+sizeof(long)-1)/sizeof(long))*sizeof(long);
187
188  int *F=(int *)omAlloc0((rN+1)*sizeof(int));
189  int *G=(int *)omAlloc0((rN+1)*sizeof(int));
190
191  memcpy(F, F0,(rN+1)*sizeof(int));
192  // pExpVectorCopy(F,F0);
193  memcpy(G, G0,(rN+1)*sizeof(int));
194  //  pExpVectorCopy(G,G0);
195  F[0]=0; /* important for p_MemAdd */
196  G[0]=0;
197
198  iF=rN;
199  while ((F[iF]==0)&&(iF>=1)) iF--; /* last exp_num of F */
200  if (iF==0) /* F0 is zero vector */
201  {
202    out=pOne();
203    p_SetExpV(out,G0,r);
204    p_Setm(out,r);
205    freeT(F,rN);
206    freeT(G,rN);
207    return(out);
208  }
209  jG=1;
210  while ((G[jG]==0)&&(jG<rN)) jG++;  /* first exp_num of G */
211  iG=rN;
212  while ((G[iG]==0)&&(iG>1)) iG--;  /* last exp_num of G */
213
214  out=pOne();
215
216  if (iF<=jG)
217    /* i.e. no mixed exp_num , MERGE case */
218  {
219    p_MemAdd_LengthGeneral(F, G, ExpSize/sizeof(long));
220    p_SetExpV(out,F,r);
221    p_Setm(out,r);
222    //    omFreeSize((ADDRESS)F,ExpSize);
223    freeT(F,rN);
224    freeT(G,rN);
225    return(out);
226  }
227
228  number cff=n_Init(1,r);
229  number tmp_num=NULL;
230  int cpower=0;
231
232  if (r->nc->type==nc_skew)
233  {
234    if (r->nc->IsSkewConstant==1)
235    {
236      int tpower=0;
237      for(j=jG; j<=iG; j++)
238      {
239        if (G[j]!=0)
240        {
241          cpower = 0;
242          for(i=j+1; i<=iF; i++)
243          {
244            cpower = cpower + F[i];
245          }
246          cpower = cpower*G[j];
247          tpower = tpower + cpower;
248        }
249      }
250      cff = n_Copy(p_GetCoeff(MATELEM(r->nc->COM,1,2),r),r);
251      nPower(cff,tpower,&tmp_num);
252      n_Delete(&cff,r);
253      cff = tmp_num;
254    }
255    else /* skew commutative with nonequal coeffs */
256    {
257      number totcff=n_Init(1,r);
258      for(j=jG; j<=iG; j++)
259      {
260        if (G[j]!=0)
261        {
262          cpower = 0;
263          for(i=j+1; i<=iF; i++)
264          {
265            if (F[i]!=0)
266            {
267              cpower = F[i]*G[j];
268              cff = n_Copy(p_GetCoeff(MATELEM(r->nc->COM,j,i),r),r);
269              nPower(cff,cpower,&tmp_num);
270              cff = nMult(totcff,tmp_num);
271              nDelete(&totcff);
272              nDelete(&tmp_num);
273              totcff = n_Copy(cff,r);
274              n_Delete(&cff,r);
275            }
276          } /* end 2nd for */
277        }
278      }
279      cff=totcff;
280    }
281    p_MemAdd_LengthGeneral(F, G, ExpSize/sizeof(long));
282    p_SetExpV(out,F,r);
283    p_Setm(out,r);
284    p_SetCoeff(out,cff,r);
285    //    p_MemAdd_NegWeightAdjust(p, r); ??? do we need this?
286    freeT(F,rN);
287    freeT(G,rN);
288    return(out);
289  } /* end nc_skew */
290   
291  /* now we have to destroy out! */
292  p_Delete(&out,r); 
293  out = NULL; 
294
295  if (iG==jG)
296    /* g is univariate monomial */
297  {
298    /*    if (ri->nc->type==nc_skew) -- postpone to TU */
299    out = nc_mm_Mult_uu(F,jG,G[jG],r);
300    freeT(F,rN);
301    freeT(G,rN);
302    return(out);
303  }
304
305  number n1=n_Init(1,r);
306  int *Prv=(int *)omAlloc0((rN+1)*sizeof(int));
307  int *Nxt=(int *)omAlloc0((rN+1)*sizeof(int));
308
309  int *log=(int *)omAlloc0((rN+1)*sizeof(int));
310  int cnt=0; int cnf=0;
311
312  /* splitting F wrt jG */
313  for (i=1;i<=jG;i++)
314  {
315    Prv[i]=F[i]; Nxt[i]=0; /* mult at the very end */
316    if (F[i]!=0) cnf++;
317  }
318
319  if (cnf==0) freeT(Prv,rN);
320
321  for (i=jG+1;i<=rN;i++)
322  {
323    Nxt[i]=F[i];
324    /*    if (cnf!=0)  Prv[i]=0; */
325    if (F[i]!=0)
326    {
327      cnt++;
328    }              /* effective part for F */
329  }
330  freeT(F,rN);
331  cnt=0;
332
333  for (i=1;i<=rN;i++)
334  {
335    if (G[i]!=0)
336    {
337     cnt++;
338     log[cnt]=i;
339     }               /* lG for G */
340   }
341
342/* ---------------------- A C T I O N ------------------------ */
343  poly D=NULL;
344  poly Rout=NULL;
345  number *c=(number *)omAlloc0((rN+1)*sizeof(number));
346  c[0]=n_Init(1,r);
347
348  int *Op=Nxt;
349  int *On=G;
350  int *U=(int *)omAlloc0((rN+1)*sizeof(int));
351
352  for (i=jG;i<=rN;i++) U[i]=Nxt[i]+G[i];  /* make leadterm */
353  Nxt=NULL;
354  G=NULL;
355  cnt=1;
356  int t=0;
357  poly w=NULL;
358  poly Pn=pOne();
359  p_SetExpV(Pn,On,r);
360  p_Setm(Pn,r);
361
362  while (On[iG]!=0)
363  {
364     t=log[cnt];
365
366     w=nc_mm_Mult_uu(Op,t,On[t],r);
367     c[cnt]=n_Mult(c[cnt-1],p_GetCoeff(w,r),r);
368     D = pNext(w);  /* getting coef and rest D */
369     p_DeleteLm(&w,r);
370     w=NULL;
371
372     Op[t] += On[t];   /* update exp_vectors */
373     On[t] = 0;
374
375     if (t!=iG)    /* not the last step */
376     {
377       p_SetExpV(Pn,On,r);
378       p_Setm(Pn,r);
379#ifdef PDEBUG
380       p_Test(Pn,r);
381#endif
382
383//       if (pNext(D)==0)
384// is D a monomial? could be postponed higher
385//       {
386//       Rout=nc_mm_Mult_nn(D,Pn,r);
387//       }
388//       else
389//       {
390       Rout=nc_p_Mult_mm(D,Pn,r);
391//       }
392     }
393     else
394     {
395       Rout=D;
396       D=NULL;
397     }
398
399     if (Rout!=NULL)
400     {
401       Rout=p_Mult_nn(Rout,c[cnt-1],r); /* Rest is ready */
402       out=p_Add_q(out,Rout,r);
403       Rout=NULL;
404     }
405     cnt++;
406  }
407  freeT(On,rN);
408  freeT(Op,rN);
409  p_Delete(&Pn,r);
410  omFreeSize((ADDRESS)log,(rN+1)*sizeof(int));
411
412  /* leadterm and Prv-part */
413
414  Rout=pOne();
415  /* U is lead.monomial */
416  U[0]=0;
417  p_SetExpV(Rout,U,r);
418  p_Setm(Rout,r);  /* use again this name Rout */
419#ifdef PDEBUG
420  p_Test(Rout,r);
421#endif
422  p_SetCoeff(Rout,c[cnt-1],r);
423  out=p_Add_q(out,Rout,r);
424  freeT(U,rN);
425  freeN(c,rN+1);
426  if (cnf!=0)  /* Prv is non-zero vector */
427  {
428    Rout=pOne();
429    Prv[0]=0;
430    p_SetExpV(Rout,Prv,r);
431    p_Setm(Rout,r);
432#ifdef PDEBUG
433    p_Test(Rout,r);
434#endif
435    out=nc_mm_Mult_p(Rout,out,r); /* getting the final result */
436    freeT(Prv,rN);
437    p_Delete(&Rout,r);
438  }
439  return (out);
440}
441
442
443poly nc_mm_Mult_uu(int *F,int jG,int bG, const ring r)
444/* f=mono(F),g=(x_iG)^bG */
445{
446  poly out=NULL;
447  int i;
448  number num=NULL;
449
450  int rN=r->N;
451  int iF=r->N;
452  while ((F[iF]==0)&&(iF>0)) iF-- ;   /* last exponent_num of F */
453
454  if (iF==0)  /* F==zero vector in other words */
455  {
456   out=pOne();
457   p_SetExp(out,jG,bG,r);
458   p_Setm(out,r);
459   return(out);
460  }
461
462  int jF=1;
463  while ((F[jF]==0)&&(jF<=rN)) jF++;  /* first exp of F */
464
465  if (iF<=jG)                       /* i.e. no mixed exp_num */
466  {
467    out=pOne();
468    F[jG]=F[jG]+bG;
469    p_SetExpV(out,F,r);
470    p_Setm(out,r);
471    return(out);
472  }
473
474  if (iF==jF)              /* uni times uni */
475  {
476   out=nc_uu_Mult_ww(iF,F[iF],jG,bG,r);
477   return(out);
478  }
479
480  /* Now: F is mono with >=2 exponents, jG<iF */
481  /* check the quasi-commutative case */
482//   matrix LCOM=r->nc->COM;
483//   number rescoef=n_Init(1,r);
484//   number tmpcoef=n_Init(1,r);
485//   int tmpint;
486//   i=iF;
487//   while (i>=jG+1)
488//     /* all the non-zero exponents */
489//   {
490//     if (MATELEM(LCOM,jG,i)!=NULL)
491//     {
492//       tmpcoef=pGetCoeff(MATELEM(LCOM,jG,i));
493//       tmpint=(int)F[i];
494//       nPower(tmpcoef,F[i],&tmpcoef);
495//       rescoef=nMult(rescoef,tmpcoef);
496//       i--;
497//     }
498//     else
499//     {
500//       if (F[i]!=0) break;
501//     }
502//   }
503//   if (iF==i)
504//   /* no action took place*/
505//   {
506
507//   }
508//   else /* power the result up to bG */
509//   {
510//     nPower(rescoef,bG,&rescoef);
511//     /* + cleanup, post-processing */
512//   }
513
514  int *Prv=(int*)omAlloc0((rN+1)*sizeof(int));
515  int *Nxt=(int*)omAlloc0((rN+1)*sizeof(int));
516  int *lF=(int *)omAlloc0((rN+1)*sizeof(int));
517  int cnt=0; int cnf=0;
518  /* splitting F wrt jG */
519  for (i=1;i<=jG;i++) /* mult at the very end */
520  {
521    Prv[i]=F[i]; Nxt[i]=0;
522    if (F[i]!=0) cnf++;
523  }
524  if (cnf==0)  freeT(Prv,rN);
525  for (i=jG+1;i<=rN;i++)
526  {
527    Nxt[i]=F[i];
528    if (cnf!=0) { Prv[i]=0;}
529    if (F[i]!=0)
530    {
531      cnt++;
532      lF[cnt]=i;
533    }                 /* eff_part,lF_for_F */
534  }
535
536  if (cnt==1) /* Nxt consists of 1 nonzero el-t only */
537  {
538    int q=lF[1];
539    poly Rout=pOne();
540    out=nc_uu_Mult_ww(q,Nxt[q],jG,bG,r);
541    freeT(Nxt,rN);
542
543    if (cnf!=0)
544    {
545       Prv[0]=0;
546       p_SetExpV(Rout,Prv,r);
547       p_Setm(Rout,r);
548#ifdef PDEBUG
549       p_Test(Rout,r);
550#endif
551       freeT(Prv,rN);
552       out=nc_mm_Mult_p(Rout,out,r); /* getting the final result */
553    }
554
555    omFreeSize((ADDRESS)lF,(rN+1)*sizeof(int));
556    p_Delete(&Rout,r);
557    return (out);
558  }
559/* -------------------- MAIN ACTION --------------------- */
560
561  poly D=NULL;
562  poly Rout=NULL;
563  number *c=(number *)omAlloc0((cnt+2)*sizeof(number));
564  c[cnt+1]=n_Init(1,r);
565  i=cnt+2;         /* later in freeN */
566  int *Op=Nxt;
567  int *On=(int *)omAlloc0((rN+1)*sizeof(int));
568  int *U=(int *)omAlloc0((rN+1)*sizeof(int));
569
570
571  //  pExpVectorCopy(U,Nxt);
572  memcpy(U, Nxt,(rN+1)*sizeof(int));
573  U[jG] = U[jG] + bG;
574
575  /* Op=Nxt and initial On=(0); */
576  Nxt=NULL;
577
578  poly Pp;
579  poly Pn;
580  int t=0;
581  int first=lF[1];
582  int nlast=lF[cnt];
583  int kk=0;
584  /*  cnt--;   */
585  /* now lF[cnt] should be <=iF-1 */
586
587  while (Op[first]!=0)
588  {
589     t=lF[cnt];   /* cnt as it was computed */
590
591     poly w=nc_uu_Mult_ww(t,Op[t],jG,bG,r);
592     c[cnt]=n_Copy(p_GetCoeff(w,r),r);
593     D = pNext(w);  /* getting coef and rest D */
594     p_DeleteLm(&w,r);
595     w=NULL;
596
597     Op[t]= 0;
598     Pp=pOne();
599     p_SetExpV(Pp,Op,r);
600     p_Setm(Pp,r);
601
602     if (t<nlast)
603     {
604       kk=lF[cnt+1];
605       On[kk]=F[kk];
606
607       Pn=pOne();
608       p_SetExpV(Pn,On,r);
609       p_Setm(Pn,r);
610
611       if (t!=first)   /* typical expr */
612       {
613         w=nc_p_Mult_mm(D,Pn,r);
614         Rout=nc_mm_Mult_p(Pp,w,r);
615         w=NULL;
616       }
617       else                   /* last step */
618       {
619         On[t]=0;
620         p_SetExpV(Pn,On,r);
621         p_Setm(Pn,r);
622         Rout=nc_p_Mult_mm(D,Pn,r);
623       }
624#ifdef PDEBUG
625       p_Test(Pp,r);
626#endif
627       p_Delete(&Pn,r);
628     }
629     else                     /* first step */
630     {
631       Rout=nc_mm_Mult_p(Pp,D,r);
632     }
633#ifdef PDEBUG
634     p_Test(Pp,r);
635#endif
636     p_Delete(&Pp,r);
637     num=n_Mult(c[cnt+1],c[cnt],r);
638     n_Delete(&c[cnt],r);
639     c[cnt]=num;
640     Rout=p_Mult_nn(Rout,c[cnt+1],r); /* Rest is ready */
641     out=p_Add_q(out,Rout,r);
642     Pp=NULL;
643     cnt--;
644  }
645  /* only to feel safe:*/
646  Pn=Pp=NULL;
647  freeT(On,rN);
648  freeT(Op,rN);
649
650/* leadterm and Prv-part with coef 1 */
651/*  U[0]=exp; */
652/*  U[jG]=U[jG]+bG;  */
653/* make leadterm */
654/* ??????????? we have done it already :-0 */
655  Rout=pOne();
656  p_SetExpV(Rout,U,r);
657  p_Setm(Rout,r);  /* use again this name */
658  p_SetCoeff(Rout,c[cnt+1],r);  /* last computed coef */
659  out=p_Add_q(out,Rout,r);
660  Rout=NULL;
661  freeT(U,rN);
662  freeN(c,i);
663  omFreeSize((ADDRESS)lF,(rN+1)*sizeof(int));
664
665  if (cnf!=0)
666  {
667    Rout=pOne();
668    p_SetExpV(Rout,Prv,r);
669    p_Setm(Rout,r);
670    freeT(Prv,rN);
671    out=nc_mm_Mult_p(Rout,out,r); /* getting the final result */
672    p_Delete(&Rout,r);
673  }
674  return (out);
675}
676
677poly nc_uu_Mult_ww_vert (int i, int a, int j, int b, const ring r)
678{
679  int k,m;
680  int rN=r->N;
681  matrix cMT=r->nc->MT[UPMATELEM(j,i,rN)];         /* cMT=current MT */
682
683  poly x=pOne();p_SetExp(x,j,1,r);p_Setm(x,r);
684/* var(j); */
685  poly y=pOne();p_SetExp(y,i,1,r);p_Setm(y,r);
686/*var(i);  for convenience */
687#ifdef PDEBUG
688  p_Test(x,r);
689  p_Test(y,r);
690#endif
691  poly t=NULL;
692/* ------------ Main Cycles ----------------------------*/
693
694  for (k=2;k<=a;k++)
695  {
696     t = nc_p_CopyGet(MATELEM(cMT,k,1),r);
697
698     if (t==NULL)   /* not computed yet */
699     {
700       t = nc_p_CopyGet(MATELEM(cMT,k-1,1),r);
701       //        t=p_Copy(MATELEM(cMT,k-1,1),r);
702       t = nc_mm_Mult_p(y,t,r);
703       MATELEM(cMT,k,1) = nc_p_CopyPut(t,r);
704       //        omCheckAddr(cMT->m);
705       p_Delete(&t,r);
706     }
707     t=NULL;
708  }
709
710  for (m=2;m<=b;m++)
711  {
712    t = nc_p_CopyGet(MATELEM(cMT,a,m),r);
713    //     t=MATELEM(cMT,a,m);
714    if (t==NULL)   //not computed yet
715    {
716      t = nc_p_CopyGet(MATELEM(cMT,a,m-1),r);
717      //      t=p_Copy(MATELEM(cMT,a,m-1),r);
718      t = nc_p_Mult_mm(t,x,r);
719      MATELEM(cMT,a,m) = nc_p_CopyPut(t,r);
720      //      MATELEM(cMT,a,m) = t;
721      //        omCheckAddr(cMT->m);
722      p_Delete(&t,r);
723    }
724    t=NULL;
725  }
726  p_Delete(&x,r);
727  p_Delete(&y,r);
728  //  t=MATELEM(cMT,a,b);
729  t= nc_p_CopyGet(MATELEM(cMT,a,b),r);
730  //  return(p_Copy(t,r));
731  /* since the last computed element was cMT[a,b] */
732  return(t);
733}
734
735poly nc_uu_Mult_ww (int i, int a, int j, int b, const ring r)
736  /* (x_i)^a times (x_j)^b */
737  /* x_i = y,  x_j = x ! */
738{
739  /* Check zero exceptions, (q-)commutativity and is there something to do? */
740  assume(a!=0);
741  assume(b!=0);
742  poly out=pOne();
743  if (i<=j)
744  {
745    p_SetExp(out,i,a,r);
746    p_AddExp(out,j,b,r);
747    p_Setm(out,r);
748    return(out);
749  }/* zero exeptions and usual case */
750  /*  if ((a==0)||(b==0)||(i<=j)) return(out); */
751
752  if (MATELEM(r->nc->COM,j,i)!=NULL)
753    /* commutative or quasicommutative case */
754  {
755    p_SetExp(out,i,a,r);
756    p_AddExp(out,j,b,r);
757    p_Setm(out,r);
758    if (r->cf->nIsOne(p_GetCoeff(MATELEM(r->nc->COM,j,i),r))) /* commutative case */
759    {
760      return(out);
761    }
762    else
763    {
764      number tmp_number=p_GetCoeff(MATELEM(r->nc->COM,j,i),r); /* quasicommutative case */
765      nPower(tmp_number,a*b,&tmp_number);
766      p_SetCoeff(out,tmp_number,r);
767      return(out);
768    }
769  }/* end_of commutative or quasicommutative case */
770  p_Delete(&out,r);
771
772  /* we are here if  i>j and variables do not commute or quasicommute */
773  /* in fact, now a>=1 and b>=1; and j<i */
774  /* now check whether the polynomial is already computed */
775  int rN=r->N;
776  int vik = UPMATELEM(j,i,rN);
777  int cMTsize=r->nc->MTsize[vik];
778  int newcMTsize=0;
779  newcMTsize=si_max(a,b);
780
781  if (newcMTsize<=cMTsize)
782  {
783    out =  nc_p_CopyGet(MATELEM(r->nc->MT[vik],a,b),r);
784    if (out !=NULL) return (out);
785  }
786  int k,m;
787  if (newcMTsize > cMTsize)
788  {
789    int inM=(((newcMTsize+6)/7)*7);
790    assume (inM>=newcMTsize);
791    newcMTsize = inM;
792    //    matrix tmp = (matrix)omAlloc0(inM*inM*sizeof(poly));
793    matrix tmp = mpNew(newcMTsize,newcMTsize);
794
795    for (k=1;k<=cMTsize;k++)
796    {
797      for (m=1;m<=cMTsize;m++)
798      {
799        out = MATELEM(r->nc->MT[UPMATELEM(j,i,rN)],k,m);
800        if ( out != NULL )
801        {
802          MATELEM(tmp,k,m) = out;/*MATELEM(r->nc->MT[UPMATELEM(j,i,rN)],k,m)*/
803          //           omCheckAddr(tmp->m);
804          MATELEM(r->nc->MT[UPMATELEM(j,i,rN)],k,m)=NULL;
805          //           omCheckAddr(r->nc->MT[UPMATELEM(j,i,rN)]->m);
806        }
807      }
808    }
809    id_Delete((ideal *)&(r->nc->MT[UPMATELEM(j,i,rN)]),r);
810    r->nc->MT[UPMATELEM(j,i,rN)] = tmp;
811    tmp=NULL;
812    r->nc->MTsize[UPMATELEM(j,i,rN)] = newcMTsize;
813  }
814  /* The update of multiplication matrix is finished */
815    pDelete(&out);
816    out = nc_uu_Mult_ww_vert(i, a, j, b, r);
817    //    out = nc_uu_Mult_ww_horvert(i, a, j, b, r);
818    return(out);
819}
820
821poly nc_uu_Mult_ww_horvert (int i, int a, int j, int b, const ring r)
822
823{
824  int k,m;
825  int rN=r->N;
826  matrix cMT=r->nc->MT[UPMATELEM(j,i,rN)];         /* cMT=current MT */
827
828  poly x=pOne();p_SetExp(x,j,1,r);p_Setm(x,r);/* var(j); */
829  poly y=pOne();p_SetExp(y,i,1,r);p_Setm(y,r); /*var(i);  for convenience */
830#ifdef PDEBUG
831  p_Test(x,r);
832  p_Test(y,r);
833#endif
834
835  poly t=NULL;
836
837  int toXY;
838  int toYX;
839
840  if (a==1) /* y*x^b, b>=2 */
841  {
842    toXY=b-1;
843    while ( (MATELEM(cMT,1,toXY)==NULL) && (toXY>=2)) toXY--;
844    for (m=toXY+1;m<=b;m++)
845    {
846      t=MATELEM(cMT,1,m);
847      if (t==NULL)   /* remove after debug */
848      {
849        t = p_Copy(MATELEM(cMT,1,m-1),r);
850        t = nc_p_Mult_mm(t,x,r);
851        MATELEM(cMT,1,m) = t;
852        /*        omCheckAddr(cMT->m); */
853      }
854      else
855      {
856        /* Error, should never get there */
857        WarnS("Error: a=1; MATELEM!=0");
858      }
859      t=NULL;
860    }
861    return(p_Copy(MATELEM(cMT,1,b),r));
862  }
863
864  if (b==1) /* y^a*x, a>=2 */
865  {
866    toYX=a-1;
867    while ( (MATELEM(cMT,toYX,1)==NULL) && (toYX>=2)) toYX--;
868    for (m=toYX+1;m<=a;m++)
869    {
870      t=MATELEM(cMT,m,1);
871      if (t==NULL)   /* remove after debug */
872      {
873        t = p_Copy(MATELEM(cMT,m-1,1),r);
874        t = nc_mm_Mult_p(y,t,r);
875        MATELEM(cMT,m,1) = t;
876        /*        omCheckAddr(cMT->m); */
877      }
878      else
879      {
880        /* Error, should never get there */
881        WarnS("Error: b=1, MATELEM!=0");
882      }
883      t=NULL;
884    }
885    return(p_Copy(MATELEM(cMT,a,1),r));
886  }
887
888/* ------------ Main Cycles ----------------------------*/
889  /*            a>1, b>1              */
890
891  int dXY=0; int dYX=0;
892  /* dXY = distance for computing x-mult, then y-mult */
893  /* dYX = distance for computing y-mult, then x-mult */
894  int toX=a-1; int toY=b-1; /* toX = to axe X, toY = to axe Y */
895  toXY=b-1; toYX=a-1;
896  /* if toX==0, toXY = dist. to computed y * x^toXY */
897  /* if toY==0, toYX = dist. to computed y^toYX * x */
898  while ( (MATELEM(cMT,toX,b)==NULL) && (toX>=1)) toX--;
899  if (toX==0) /* the whole column is not computed yet */
900  {
901    while ( (MATELEM(cMT,1,toXY)==NULL) && (toXY>=1)) toXY--;
902    /* toXY >=1 */
903    dXY=b-1-toXY;
904  }
905  dXY=dXY+a-toX; /* the distance to nearest computed y^toX x^b */
906
907  while ( (MATELEM(cMT,a,toY)==NULL) && (toY>=1)) toY--;
908  if (toY==0) /* the whole row is not computed yet */
909  {
910    while ( (MATELEM(cMT,toYX,1)==NULL) && (toYX>=1)) toYX--;
911    /* toYX >=1 */
912    dYX=a-1-toYX;
913  }
914  dYX=dYX+b-toY; /* the distance to nearest computed y^a x^toY */
915
916  if (dYX>=dXY)
917  {
918    /* first x, then y */
919    if (toX==0) /* start with the row*/
920    {
921      for (m=toXY+1;m<=b;m++)
922      {
923        t=MATELEM(cMT,1,m);
924        if (t==NULL)   /* remove after debug */
925        {
926          t = p_Copy(MATELEM(cMT,1,m-1),r);
927          t = nc_p_Mult_mm(t,x,r);
928          MATELEM(cMT,1,m) = t;
929          /*        omCheckAddr(cMT->m); */
930        }
931        else
932        {
933          /* Error, should never get there */
934          WarnS("dYX>=dXY,toXY; MATELEM==0");
935        }
936        t=NULL;
937      }
938      toX=1; /* y*x^b is computed */
939    }
940    /* Now toX>=1 */
941    for (k=toX+1;k<=a;k++)
942    {
943      t=MATELEM(cMT,k,b);
944      if (t==NULL)   /* remove after debug */
945      {
946        t = p_Copy(MATELEM(cMT,k-1,b),r);
947        t = nc_mm_Mult_p(y,t,r);
948        MATELEM(cMT,k,b) = t;
949        /*        omCheckAddr(cMT->m); */
950      }
951      else
952      {
953        /* Error, should never get there */
954        WarnS("dYX>=dXY,toX; MATELEM==0");
955      }
956      t=NULL;
957    }
958  } /* endif (dYX>=dXY) */
959
960
961  if (dYX<dXY)
962  {
963    /* first y, then x */
964    if (toY==0) /* start with the column*/
965    {
966      for (m=toYX+1;m<=a;m++)
967      {
968        t=MATELEM(cMT,m,1);
969        if (t==NULL)   /* remove after debug */
970        {
971          t = p_Copy(MATELEM(cMT,m-1,1),r);
972          t = nc_mm_Mult_p(y,t,r);
973          MATELEM(cMT,m,1) = t;
974          /*        omCheckAddr(cMT->m); */
975        }
976        else
977        {
978          /* Error, should never get there */
979          WarnS("dYX<dXY,toYX; MATELEM==0");
980        }
981        t=NULL;
982      }
983      toY=1; /* y^a*x is computed */
984    }
985    /* Now toY>=1 */
986    for (k=toY+1;k<=b;k++)
987    {
988      t=MATELEM(cMT,a,k);
989      if (t==NULL)   /* remove after debug */
990      {
991        t = p_Copy(MATELEM(cMT,a,k-1),r);
992        t = nc_p_Mult_mm(t,x,r);
993        MATELEM(cMT,a,k) = t;
994        /*        omCheckAddr(cMT->m); */
995      }
996      else
997      {
998        /* Error, should never get there */
999        WarnS("dYX<dXY,toY; MATELEM==0");
1000      }
1001      t=NULL;
1002    }
1003  } /* endif (dYX<dXY) */
1004
1005  p_Delete(&x,r);
1006  p_Delete(&y,r);
1007  t=p_Copy(MATELEM(cMT,a,b),r);
1008  return(t);  /* since the last computed element was cMT[a,b] */
1009}
1010
1011
1012/* ----------------------------- Syzygies ---------------------- */
1013
1014/*2
1015* reduction of p2 with p1
1016* do not destroy p1, but p2
1017* p1 divides p2 -> for use in NF algorithm
1018*/
1019
1020poly nc_ReduceSpoly(poly p1, poly p2,poly spNoether, const ring r)
1021{
1022  if (p_GetComp(p1,r)!=p_GetComp(p2,r)
1023  && (p_GetComp(p1,r)!=0)
1024  && (p_GetComp(p2,r)!=0))
1025  {
1026#ifdef PDEBUG
1027    Print("nc_ReduceSpoly: different components");
1028#endif
1029    return(NULL);
1030  }
1031  poly m = pOne();
1032  p_ExpVectorDiff(m,p2,p1,r);
1033  //p_Setm(m,r);
1034#ifdef PDEBUG
1035  p_Test(m,r);
1036#endif
1037  /* pSetComp(m,r)=0? */
1038  poly   N  = nc_mm_Mult_p(m, p_Head(p1,r), r);
1039  number C  = n_Copy( p_GetCoeff(N,  r), r);
1040  number cF = n_Copy( p_GetCoeff(p2, r),r);
1041  /* GCD stuff */
1042  number cG = nGcd(C, cF, r);
1043  if ( !nEqual(cG, n_Init(1,r) ) )
1044  {
1045    cF = nDiv(cF, cG);
1046    C  = nDiv(C,  cG);
1047  }
1048  p2 = p_Mult_nn(p2, C, r);
1049  poly out = nc_mm_Mult_p(m, p_Copy(pNext(p1),r), r);
1050  N = p_Add_q(N, out, r);
1051  p_Test(p2,r);
1052  p_Test(N,r);
1053  number MinusOne = n_Init(-1,r);
1054  if (!n_Equal(cF,MinusOne,r))
1055  {
1056    cF = n_Neg(cF,r);
1057    N  = p_Mult_nn(N, cF, r);
1058    p_Test(N,r);
1059  }
1060  out = p_Add_q(p2,N,r);
1061  p_Test(out,r);
1062  if ( out!=NULL ) pContent(out);
1063  p_Delete(&m,r);
1064  n_Delete(&cF,r);
1065  n_Delete(&C,r);
1066  n_Delete(&MinusOne,r);
1067  return(out);
1068}
1069
1070
1071/*3
1072* reduction of p2 with p1
1073* do not destroy p1 and p2
1074* p1 divides p2 -> for use in NF algorithm
1075*/
1076poly nc_ReduceSpolyNew(poly p1, poly p2,poly spNoether, const ring r)
1077{
1078  return(nc_ReduceSpoly(p1,p_Copy(p2,r),spNoether,r));
1079}
1080
1081/*4
1082* creates the S-polynomial of p1 and p2
1083* do not destroy p1 and p2
1084*/
1085poly nc_CreateSpoly(poly p1, poly p2,poly spNoether, const ring r)
1086{
1087  if ((p_GetComp(p1,r)!=p_GetComp(p2,r))
1088  && (p_GetComp(p1,r)!=0)
1089  && (p_GetComp(p2,r)!=0))
1090  {
1091#ifdef PDEBUG
1092    Print("nc_CreateSpoly : different components!");
1093#endif
1094    return(NULL);
1095  }
1096  if ((r->nc->type==nc_lie) && pHasNotCF(p1,p2)) /* prod crit */
1097  {
1098    return(nc_p_Bracket_qq(pCopy(p2),p1));
1099  }
1100  poly pL=pOne();
1101  poly m1=pOne();
1102  poly m2=pOne();
1103  pLcm(p1,p2,pL);
1104  p_Setm(pL,r);
1105#ifdef PDEBUG
1106  p_Test(pL,r);
1107#endif
1108  p_ExpVectorDiff(m1,pL,p1,r);
1109  //p_SetComp(m1,0,r);
1110  //p_Setm(m1,r);
1111#ifdef PDEBUG
1112  p_Test(m1,r);
1113#endif
1114  p_ExpVectorDiff(m2,pL,p2,r);
1115  //p_SetComp(m2,0,r);
1116  //p_Setm(m2,r);
1117#ifdef PDEBUG
1118  p_Test(m2,r);
1119#endif
1120  p_Delete(&pL,r);
1121  /* zero exponents ! */
1122  poly M1    = nc_mm_Mult_p(m1,p_Head(p1,r),r);
1123  number C1  = n_Copy(p_GetCoeff(M1,r),r);
1124  poly M2    = nc_mm_Mult_p(m2,p_Head(p2,r),r);
1125  number C2  = n_Copy(p_GetCoeff(M2,r),r);
1126  /* GCD stuff */
1127  number C = nGcd(C1,C2,r);
1128  if (!nEqual(C,n_Init(1,r)))
1129  {
1130    C1=nDiv(C1,C);
1131    C2=nDiv(C2,C);
1132  }
1133  M1=p_Mult_nn(M1,C2,r);
1134  p_SetCoeff(m1,C2,r);
1135  number MinusOne=n_Init(-1,r);
1136  if (n_Equal(C1,MinusOne,r))
1137  {
1138    M2=p_Add_q(M1,M2,r);
1139  }
1140  else
1141  {
1142    C1=n_Neg(C1,r);
1143    M2=p_Mult_nn(M2,C1,r);
1144    M2=p_Add_q(M1,M2,r);
1145    p_SetCoeff(m2,C1,r);
1146  }
1147  /* M1 is killed, M2=res = C2 M1 - C1 M2 */
1148  poly tmp=p_Copy(p1,r);
1149  tmp=p_LmDeleteAndNext(tmp,r);
1150  M1=nc_mm_Mult_p(m1,tmp,r);
1151  tmp=p_Copy(p2,r);
1152  tmp=p_LmDeleteAndNext(tmp,r);
1153  M2=p_Add_q(M2,M1,r);
1154  M1=nc_mm_Mult_p(m2,tmp,r);
1155  M2=p_Add_q(M2,M1,r);
1156  p_Delete(&m1,r);
1157  p_Delete(&m2,r);
1158  //  n_Delete(&C1,r);
1159  //  n_Delete(&C2,r);
1160  n_Delete(&MinusOne,r);
1161#ifdef PDEBUG
1162  p_Test(M2,r);
1163#endif
1164  if (M2!=NULL) pContent(M2);
1165  return(M2);
1166}
1167
1168/*5
1169* reduction of tail(q) with p1
1170* lead(p1) divides lead(pNext(q2)) and pNext(q2) is reduced
1171* do not destroy p1, but tail(q)
1172*/
1173void nc_ReduceSpolyTail(poly p1, poly q, poly q2, poly spNoether, const ring r)
1174{
1175  poly a1=p_Head(p1,r);
1176  poly Q=pNext(q2);
1177  number cQ=p_GetCoeff(Q,r);
1178  poly m=pOne();
1179  p_ExpVectorDiff(m,Q,p1,r);
1180  //  p_SetComp(m,0,r);
1181  //p_Setm(m,r);
1182#ifdef PDEBUG
1183  p_Test(m,r);
1184#endif
1185  /* pSetComp(m,r)=0? */
1186  poly M=nc_mm_Mult_p(m,p_Copy(p1,r),r);
1187  number C=p_GetCoeff(M,r);
1188  M=p_Add_q(M,nc_mm_Mult_p(m,p_LmDeleteAndNext(p_Copy(p1,r),r),r),r);
1189  q=p_Mult_nn(q,C,r);
1190  number MinusOne=n_Init(-1,r);
1191  if (!n_Equal(cQ,MinusOne,r))
1192  {
1193    cQ=nNeg(cQ);
1194    M=p_Mult_nn(M,cQ,r);
1195  }
1196  Q=p_Add_q(Q,M,r);
1197  pNext(q2)=Q;
1198
1199  p_Delete(&m,r);
1200  n_Delete(&C,r);
1201  n_Delete(&cQ,r);
1202  n_Delete(&MinusOne,r);
1203  /*  return(q); */
1204}
1205
1206/*6
1207* creates the commutative lcm(lm(p1),lm(p2))
1208* do not destroy p1 and p2
1209*/
1210poly nc_CreateShortSpoly(poly p1, poly p2, const ring r)
1211{
1212  if (p_GetComp(p1,r)!=p_GetComp(p2,r))
1213  {
1214#ifdef PDEBUG
1215    Print("spShort:exponent mismatch!");
1216#endif
1217    return(NULL);
1218  }
1219  poly m=pOne();
1220  pLcm(p1,p2,m);
1221  p_Setm(m,r);
1222#ifdef PDEBUG
1223  p_Test(m,r);
1224#endif
1225  return(m);
1226}
1227
1228void nc_kBucketPolyRed(kBucket_pt b, poly p, number *c)
1229{
1230  // b will not be multiplied by any constant in this impl.
1231  // ==> *c=1
1232  *c=nInit(1);
1233  poly m=pOne();
1234  pExpVectorDiff(m,kBucketGetLm(b),p);
1235  //pSetm(m);
1236#ifdef PDEBUG
1237  pTest(m);
1238#endif
1239  poly pp=nc_mm_Mult_p(m,pCopy(p),currRing);
1240  pDelete(&m);
1241  number n=nCopy(pGetCoeff(pp));
1242  number MinusOne=nInit(-1);
1243  number nn;
1244  if (!nEqual(n,MinusOne))
1245  {
1246    nn=nNeg(nInvers(n));
1247  }
1248  else nn=nInit(1);
1249  nDelete(&n);
1250  n=nMult(nn,pGetCoeff(kBucketGetLm(b)));
1251  nDelete(&nn);
1252  pp=p_Mult_nn(pp,n,currRing);
1253  nDelete(&n);
1254  nDelete(&MinusOne);
1255  int l=pLength(pp);
1256  kBucket_Add_q(b,pp,&l);
1257}
1258
1259void nc_kBucketPolyRed_Z(kBucket_pt b, poly p, number *c)
1260{
1261  // b is multiplied by a constant in this impl.
1262  poly m=pOne();
1263  pExpVectorDiff(m,kBucketGetLm(b),p);
1264  //pSetm(m);
1265#ifdef PDEBUG
1266  pTest(m);
1267#endif
1268  if(p_IsConstant(m,currRing)){
1269    pDelete(&m);
1270    *c = kBucketPolyRed(b,p,pLength(p),NULL);
1271    return;
1272  }
1273  poly pp = nc_mm_Mult_p(m,pCopy(p),currRing);
1274  pDelete(&m);
1275  *c = kBucketPolyRed(b,pp,pLength(pp),NULL);
1276  pDelete(&pp);
1277}
1278
1279void nc_PolyPolyRed(poly &b, poly p, number *c)
1280  // reduces b with p, do not delete both
1281{
1282  // b will not by multiplied by any constant in this impl.
1283  // ==> *c=1
1284  *c=nInit(1);
1285  poly m=pOne();
1286  pExpVectorDiff(m,pHead(b),p);
1287  //pSetm(m);
1288#ifdef PDEBUG
1289  pTest(m);
1290#endif
1291  poly pp=nc_mm_Mult_p(m,pCopy(p),currRing);
1292  pDelete(&m);
1293  number n=nCopy(pGetCoeff(pp));
1294  number MinusOne=nInit(-1);
1295  number nn;
1296  if (!nEqual(n,MinusOne))
1297  {
1298    nn=nNeg(nInvers(n));
1299  }
1300  else nn=nInit(1);
1301  nDelete(&n);
1302  n=nMult(nn,pGetCoeff(b));
1303  nDelete(&nn);
1304  pp=p_Mult_nn(pp,n,currRing);
1305  nDelete(&n);
1306  nDelete(&MinusOne);
1307  b=p_Add_q(b,pp,currRing);
1308}
1309
1310poly nc_p_Bracket_qq(poly p, poly q)
1311  /* returns [p,q], destroys p */
1312{
1313  if (!rIsPluralRing(currRing)) return(NULL);
1314  if (pComparePolys(p,q)) return(NULL);
1315  /* Components !? */
1316  poly Q=NULL;
1317  number coef=NULL;
1318  poly res=NULL;
1319  poly pres=NULL;
1320  int UseBuckets=1;
1321  if ((pLength(p)< MIN_LENGTH_BUCKET/2) && (pLength(q)< MIN_LENGTH_BUCKET/2) || TEST_OPT_NOT_BUCKETS) UseBuckets=0;
1322  sBucket_pt bu_out;
1323  if (UseBuckets) bu_out=sBucketCreate(currRing);
1324  while (p!=NULL)
1325  {
1326    Q=q;
1327    while(Q!=NULL)
1328    {
1329      pres=nc_mm_Bracket_nn(p,Q); /* since no coeffs are taken into account there */
1330      if (pres!=NULL)
1331      {
1332        coef = nMult(pGetCoeff(p),pGetCoeff(Q));
1333        pres = p_Mult_nn(pres,coef,currRing);
1334        if (UseBuckets) sBucket_Add_p(bu_out,pres,pLength(pres));
1335        else res=p_Add_q(res,pres,currRing);
1336        nDelete(&coef);
1337      }
1338      pIter(Q);
1339    }
1340    p=pLmDeleteAndNext(p);
1341  }
1342  if (UseBuckets)
1343  {
1344    res = NULL;
1345    int len = pLength(res);
1346    sBucketDestroyAdd(bu_out, &res, &len);
1347  }
1348  return(res);
1349}
1350
1351poly nc_mm_Bracket_nn(poly m1, poly m2)
1352  /*returns [m1,m2] for two monoms, destroys nothing */
1353  /* without coeffs */
1354{
1355  if (pLmIsConstant(m1) || pLmIsConstant(m1)) return(NULL);
1356  if (pLmCmp(m1,m2)==0) return(NULL);
1357  int rN=currRing->N;
1358  int *M1=(int *)omAlloc0((rN+1)*sizeof(int));
1359  int *M2=(int *)omAlloc0((rN+1)*sizeof(int));
1360  int *PREFIX=(int *)omAlloc0((rN+1)*sizeof(int));
1361  int *SUFFIX=(int *)omAlloc0((rN+1)*sizeof(int));
1362  pGetExpV(m1,M1);
1363  pGetExpV(m2,M2);
1364  poly res=NULL;
1365  poly ares=NULL;
1366  poly bres=NULL;
1367  poly prefix=NULL;
1368  poly suffix=NULL;
1369  int nMin,nMax;
1370  number nTmp=NULL;
1371  int i,j,k;
1372  for (i=1;i<=rN;i++)
1373  {
1374    if (M2[i]!=0)
1375    {
1376      ares=NULL;
1377      for (j=1;j<=rN;j++)
1378      {
1379        if (M1[j]!=0)
1380        {
1381          bres=NULL;
1382          /* compute [ x_j^M1[j],x_i^M2[i] ] */
1383          if (i<j) {nMax=j;  nMin=i;} else {nMax=i;  nMin=j;}
1384          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)*/
1385          { bres=NULL; }
1386          else
1387          {
1388            if (i<j) { bres=nc_uu_Mult_ww(j,M1[j],i,M2[i],currRing); }
1389            else bres=nc_uu_Mult_ww(i,M2[i],j,M1[j],currRing);
1390            if (nIsOne(pGetCoeff(bres)))
1391            {
1392              bres=pLmDeleteAndNext(bres);
1393            }
1394            else
1395            {
1396              nTmp=nSub(pGetCoeff(bres),nInit(1));
1397              pSetCoeff(bres,nTmp); /* only lc ! */
1398            }
1399#ifdef PDEBUG
1400            pTest(bres);
1401#endif
1402            if (i>j)  bres=p_Neg(bres, currRing);
1403          }
1404          if (bres!=NULL)
1405          {
1406            /* now mult (prefix, bres, suffix) */
1407            memcpy(SUFFIX, M1,(rN+1)*sizeof(int));
1408            memcpy(PREFIX, M1,(rN+1)*sizeof(int));
1409            for (k=1;k<=j;k++) SUFFIX[k]=0;
1410            for (k=j;k<=rN;k++) PREFIX[k]=0;
1411            SUFFIX[0]=0;
1412            PREFIX[0]=0;
1413            prefix=pOne();
1414            suffix=pOne();
1415            pSetExpV(prefix,PREFIX);
1416            pSetm(prefix);
1417            pSetExpV(suffix,SUFFIX);
1418            pSetm(suffix);
1419            if (!pLmIsConstant(prefix)) bres = nc_mm_Mult_p(prefix, bres,currRing);
1420            if (!pLmIsConstant(suffix)) bres = nc_p_Mult_mm(bres, suffix,currRing);
1421            ares=p_Add_q(ares, bres,currRing);
1422            /* What to give free? */
1423            /* Do we have to free PREFIX/SUFFIX? it seems so */
1424            pDelete(&prefix);
1425            pDelete(&suffix);
1426          }
1427        }
1428      }
1429      if (ares!=NULL)
1430      {
1431        /* now mult (prefix, bres, suffix) */
1432        memcpy(SUFFIX, M2,(rN+1)*sizeof(int));
1433        memcpy(PREFIX, M2,(rN+1)*sizeof(int));
1434        for (k=1;k<=i;k++) SUFFIX[k]=0;
1435        for (k=i;k<=rN;k++) PREFIX[k]=0;
1436        SUFFIX[0]=0;
1437        PREFIX[0]=0;
1438        prefix=pOne();
1439        suffix=pOne();
1440        pSetExpV(prefix,PREFIX);
1441        pSetm(prefix);
1442        pSetExpV(suffix,SUFFIX);
1443        pSetm(suffix);
1444        bres=ares;
1445        if (!pLmIsConstant(prefix)) bres = nc_mm_Mult_p(prefix, bres,currRing);
1446        if (!pLmIsConstant(suffix)) bres = nc_p_Mult_mm(bres, suffix,currRing);
1447        res=p_Add_q(res, bres,currRing);
1448        pDelete(&prefix);
1449        pDelete(&suffix);
1450      }
1451    }
1452  }
1453  freeT(M1, rN);
1454  freeT(M2, rN);
1455  freeT(PREFIX, rN);
1456  freeT(SUFFIX, rN);
1457  pTest(res);
1458  return(res);
1459}
1460
1461ideal twostd(ideal I)
1462{
1463  int i;
1464  int j;
1465  int s;
1466  int flag;
1467  poly p=NULL;
1468  poly q=NULL;
1469  ideal J=kStd(I, currQuotient,testHomog,NULL,NULL,0,0,NULL);
1470  idSkipZeroes(J);
1471  ideal K=NULL;
1472  poly varj=NULL;
1473  ideal Q=NULL;
1474  ideal id_tmp=NULL;
1475  int rN=currRing->N;
1476  int iSize=0;
1477  loop
1478  {
1479    flag=0;
1480    K=NULL;
1481    s=idElem(J);
1482    for (i=0;i<=s-1;i++)
1483    {
1484      p=J->m[i];
1485      for (j=1;j<=rN;j++)
1486      {
1487        varj = pOne();
1488        pSetExp(varj,j,1);
1489        pSetm(varj);
1490        q = nc_p_Mult_mm(pCopy(p),varj,currRing);
1491        pDelete(&varj);
1492        q = nc_ReduceSpoly(p,q,NULL,currRing);
1493        q = kNF(J,currQuotient,q,0,0);
1494        if (q!=NULL)
1495        {
1496          if (pIsConstant(q))
1497          {
1498            Q=idInit(1,1);
1499            Q->m[0]=pOne();
1500            idDelete(&J);
1501            pDelete(&q);
1502            if (K!=NULL) idDelete(&K);
1503            return(Q);
1504          }
1505          flag=1;
1506          Q=idInit(1,1);
1507          Q->m[0]=q;
1508          id_tmp=idSimpleAdd(K,Q);
1509          idDelete(&K);
1510          K=id_tmp;
1511          idDelete(&Q);
1512        }
1513      }
1514    }
1515    if (flag==0)
1516      /* i.e. all elements are two-sided */
1517    {
1518      idDelete(&K);
1519      return(J);
1520    }
1521    /* now we update GrBasis J with K */
1522    //    iSize=IDELEMS(J);
1523    iSize=idElem(J);
1524    id_tmp=idSimpleAdd(J,K);
1525    idDelete(&K);
1526    idDelete(&J); 
1527    BITSET save_test=test;
1528    test|=Sy_bit(OPT_SB_1);
1529    J=kStd(id_tmp, currQuotient, testHomog,NULL,NULL,0,iSize);
1530    test=save_test;
1531    idSkipZeroes(J);
1532  }
1533}
1534
1535matrix nc_PrintMat(int a, int b, ring r, int metric)
1536  /* returns matrix with the info on noncomm multiplication */
1537{
1538
1539  if ( (a==b) || !rIsPluralRing(r) ) return(NULL);
1540  int i;
1541  int j;
1542  if (a>b) {j=b; i=a;}
1543  else {j=a; i=b;}
1544  /* i<j */
1545  int rN=r->N;
1546  int size=r->nc->MTsize[UPMATELEM(i,j,rN)];
1547  matrix M = r->nc->MT[UPMATELEM(i,j,rN)];
1548  /*  return(M); */
1549  int sizeofres;
1550  if (metric==0)
1551  {
1552    sizeofres=sizeof(int);
1553  }
1554  if (metric==1)
1555  {
1556    sizeofres=sizeof(number);
1557  }
1558  matrix res=mpNew(size,size);
1559  int s;
1560  int t;
1561  int length;
1562  long totdeg;
1563  poly p;
1564  for(s=1;s<=size;s++)
1565  {
1566    for(t=1;t<=size;t++)
1567    {
1568      p=MATELEM(M,s,t);
1569      if (p==NULL)
1570      {
1571        MATELEM(res,s,t)=0;
1572      }
1573      else
1574      {
1575        length = pLength(p);
1576        if (metric==0) /* length */
1577        {
1578          MATELEM(res,s,t)= p_ISet(length,r);
1579        }
1580        else if (metric==1) /* sum of deg divided by the length */
1581        {
1582          totdeg=0;
1583          while (p!=NULL)
1584          {
1585            totdeg=totdeg+pDeg(p,r);
1586            pIter(p);
1587          }
1588          number ntd = nInit(totdeg);
1589          number nln = nInit(length);
1590          number nres=nDiv(ntd,nln);
1591          nDelete(&ntd);
1592          nDelete(&nln);
1593          MATELEM(res,s,t)=p_NSet(nres,r);
1594        }
1595      }
1596    }
1597  }
1598  return(res);
1599}
1600
1601void ncKill(ring r)
1602  /* kills the nc extension of ring r */
1603{
1604  int i,j;
1605  int rN=r->N;
1606  if ( rN > 1 )
1607  {
1608    for(i=1;i<rN;i++)
1609    {
1610      for(j=i+1;j<=rN;j++)
1611      {
1612        id_Delete((ideal *)&(r->nc->MT[UPMATELEM(i,j,rN)]),r->nc->basering);
1613      }
1614    }
1615    omFreeSize((ADDRESS)r->nc->MT,rN*(rN-1)/2*sizeof(matrix));
1616    omFreeSize((ADDRESS)r->nc->MTsize,rN*(rN-1)/2*sizeof(int));
1617    id_Delete((ideal *)&(r->nc->COM),r->nc->basering);
1618  }
1619  id_Delete((ideal *)&(r->nc->C),r->nc->basering);
1620  id_Delete((ideal *)&(r->nc->D),r->nc->basering);
1621  r->nc->basering->ref--;
1622  if (r->nc->basering<=0) 
1623  {
1624    rKill(r->nc->basering);
1625  }
1626  omFreeSize((ADDRESS)r->nc,sizeof(nc_struct));
1627  r->nc=NULL;
1628}
1629
1630void ncCleanUp(ring r)
1631{
1632  /* small CleanUp of r->nc */
1633  omFreeSize((ADDRESS)r->nc,sizeof(nc_struct));
1634  r->nc = NULL;
1635}
1636
1637poly nc_p_CopyGet(poly a, const ring r)
1638/* for use in getting the mult. matrix elements*/
1639/* ring r must be a currRing! */
1640/* for consistency, copies a poly from the comm. r->nc->basering */
1641/* to its image in NC ring */
1642{
1643  if (r != currRing)
1644  {
1645#ifdef PDEBUF
1646    Werror("nc_p_CopyGet call not in currRing");
1647#endif
1648    return(NULL);
1649  }
1650  if (!rIsPluralRing(r)) return(p_Copy(a,r));
1651  if (r==r->nc->basering) return(p_Copy(a,r));
1652  else
1653  {
1654    return(prCopyR_NoSort(a,r->nc->basering,r));
1655  }
1656}
1657
1658poly nc_p_CopyPut(poly a, const ring r)
1659/* for use in defining the mult. matrix elements*/
1660/* ring r must be a currRing! */
1661/* for consistency, puts a polynomial from the NC ring */
1662/* to its presentation in the comm. r->nc->basering */
1663{
1664  if (r != currRing)
1665  {
1666#ifdef PDEBUF
1667    Werror("nc_p_CopyGet call not in currRing");
1668#endif
1669    return(NULL);
1670  }
1671  if (!rIsPluralRing(r)) return(p_Copy(a,r));
1672  if (r==r->nc->basering) return(p_Copy(a,r));
1673  else
1674  {
1675    return(prCopyR_NoSort(a,r,r->nc->basering));
1676  }
1677}
1678
1679BOOLEAN nc_CheckSubalgebra(poly PolyVar, ring r)
1680  /* returns TRUE if there were errors */
1681  /* checks whether product of vars from PolyVar defines */
1682  /* an admissible subalgebra of r */
1683  /* r is indeed currRing and assumed to be noncomm. */
1684{
1685  ring save = currRing;
1686  int WeChangeRing = 0;
1687  if (currRing != r)
1688  {
1689    rChangeCurrRing(r);
1690    WeChangeRing = 1;
1691  }
1692  int rN=r->N;
1693  int *ExpVar=(int*)omAlloc0((rN+1)*sizeof(int));
1694  int *ExpTmp=(int*)omAlloc0((rN+1)*sizeof(int));
1695  p_GetExpV(PolyVar, ExpVar, r);
1696  int i; int j; int k;
1697  poly test=NULL;
1698  int OK=1;
1699  for (i=1; i<rN; i++)
1700  {
1701    if (ExpVar[i]==0) /* i.e. not in PolyVar */
1702    { 
1703      for (j=i+1; j<=rN; j++)
1704      {
1705        if (ExpVar[j]==0)
1706        {
1707          test = nc_p_CopyGet(MATELEM(r->nc->D,i,j),r);
1708          while (test!=NULL)
1709          {
1710            p_GetExpV(test, ExpTmp, r);
1711            OK=1;
1712            for (k=1;k<=rN;k++)
1713            {
1714              if (ExpTmp[k]!=0)
1715              {
1716                if (ExpVar[k]!=0) OK=0;
1717              }
1718            }
1719            if (!OK) return(TRUE);
1720            pIter(test);
1721          }
1722        }
1723      }
1724    }
1725  }
1726  p_Delete(&test,r);
1727  freeT(ExpVar,rN);
1728  freeT(ExpTmp,rN);
1729  if ( WeChangeRing )
1730    rChangeCurrRing(save);
1731  return(FALSE);
1732}
1733
1734BOOLEAN nc_CheckOrdCondition(matrix D, ring r)
1735/* returns TRUE if there were errors */
1736/* checks whether the current ordering */
1737/* is admissible for r and D == r->nc->D */
1738/* to be executed in a currRing */
1739{
1740  /* analyze D: an upper triangular matrix of polys */ 
1741  /* check the ordering condition for D */
1742  ring save = currRing;
1743  int WeChangeRing = 0;
1744  if (r != currRing)
1745  {
1746    rChangeCurrRing(r);
1747    WeChangeRing = 1;
1748  }
1749  poly p,q;
1750  int i,j;
1751  int report = 0;
1752  for(i=1; i<r->N; i++)
1753  {
1754    for(j=i+1; j<=r->N; j++)
1755    { 
1756      p = nc_p_CopyGet(MATELEM(D,i,j),r);
1757      if ( p != NULL)
1758      {
1759        q = p_ISet(1,r); // replaces pOne();
1760        p_SetExp(q,i,1,r);
1761        p_SetExp(q,j,1,r);
1762        p_Setm(q,r);
1763        if (p_LmCmp(q,p,r) != 1) /* i.e. lm(p)==xy < lm(q)==D_ij  */
1764        {
1765          Print("Bad ordering at %d,%d\n",i,j);
1766#ifdef PDEBUG
1767          p_Write(p,r);
1768          p_Write(q,r);
1769#endif
1770          report = 1;
1771        }
1772        p_Delete(&q,r);
1773        p_Delete(&p,r);
1774        p = NULL;
1775      }
1776    }
1777  }
1778  if ( WeChangeRing )
1779    rChangeCurrRing(save);
1780  return(report);
1781}
1782
1783
1784
1785BOOLEAN nc_CallPlural(matrix CCC, matrix DDD, poly CCN, poly DDN, ring r)
1786  /* returns TRUE if there were errors */
1787  /* analyze inputs, check them for consistency */
1788  /* detects nc_type, DO NOT initialize multiplication but call for it at the end*/
1789  /* checks the ordering condition and evtl. NDC */
1790{
1791  matrix CC = NULL; 
1792  matrix DD = NULL;
1793  poly CN = NULL;
1794  poly DN = NULL;
1795  matrix C;
1796  matrix D;
1797  number nN,pN,qN;
1798  int tmpIsSkewConstant;
1799  int i,j;
1800  if (r->nc != NULL)
1801  {
1802    WarnS("redefining algebra structure");
1803    if (r->nc->ref>1) /* in use by somebody else */
1804    {
1805      r->nc->ref--;
1806    }
1807    else  /* kill the previous nc data */
1808    {
1809      ncKill(r); 
1810    }
1811  }
1812  ring save = currRing;
1813  int WeChangeRing = 0;
1814  if (currRing!=r)
1815  {
1816    rChangeCurrRing(r);
1817    WeChangeRing = 1;
1818  }
1819  r->nc = (nc_struct *)omAlloc0(sizeof(nc_struct));
1820  r->nc->ref = 1;
1821  r->nc->basering = r;
1822  r->ref++;
1823  r->nc->type = nc_undef;
1824
1825  /* initialition of the matrix C */
1826  /* check the correctness of arguments */
1827
1828  if ((CCC != NULL) && ( (MATCOLS(CCC)==1) || MATROWS(CCC)==1 ) )
1829  {
1830    CN = MATELEM(CCC,1,1);
1831  }
1832  else 
1833  {
1834    if ((CCC != NULL) && ( (MATCOLS(CCC)!=r->N) || (MATROWS(CCC)!=r->N) ))
1835    {
1836      Werror("Square %d x %d  matrix expected",r->N,r->N);
1837      ncCleanUp(r);
1838      if (WeChangeRing)
1839        rChangeCurrRing(save);
1840      return TRUE;
1841    }
1842  }
1843  if (( CCC != NULL) && (CC == NULL)) CC = mpCopy(CCC);
1844  if (( CCN != NULL) && (CN == NULL)) CN = CCN;
1845
1846  /* initialition of the matrix D */
1847  /* check the correctness of arguments */
1848
1849  if ((DDD != NULL) && ( (MATCOLS(DDD)==1) || MATROWS(DDD)==1 ) )
1850  {
1851    DN = MATELEM(DDD,1,1);
1852  }
1853  else 
1854  {
1855    if ((DDD != NULL) && ( (MATCOLS(DDD)!=r->N) || (MATROWS(DDD)!=r->N) ))
1856    {
1857      Werror("Square %d x %d  matrix expected",r->N,r->N);
1858      ncCleanUp(r);
1859      if (WeChangeRing)
1860        rChangeCurrRing(save);
1861      return TRUE;
1862    }
1863  }
1864  if (( DDD != NULL) && (DD == NULL)) DD = mpCopy(DDD);
1865  if (( DDN != NULL) && (DN == NULL)) DN = DDN;
1866
1867  /* further checks */
1868
1869  if (CN != NULL)       /* create matrix C = CN * Id */
1870  {
1871    nN = p_GetCoeff(CN,r);
1872    if (n_IsZero(nN,r))
1873    {
1874      Werror("Incorrect input : zero coefficients are not allowed");
1875      ncCleanUp(r);
1876      if (WeChangeRing)
1877        rChangeCurrRing(save);
1878      return TRUE;
1879    }
1880    if (nIsOne(nN)) 
1881    {
1882      r->nc->type = nc_lie; 
1883    }
1884    else 
1885    {
1886      r->nc->type = nc_general;
1887    }
1888    r->nc->IsSkewConstant = 1;
1889    C = mpNew(r->N,r->N);
1890    for(i=1; i<r->N; i++)
1891    {
1892      for(j=i+1; j<=r->N; j++)
1893      {
1894        MATELEM(C,i,j) = nc_p_CopyPut(CN,r);
1895      }
1896    }
1897  }
1898  if ( (CN == NULL) && (CC != NULL) ) /* copy matrix C */
1899  {
1900    C = mpCopy(CC);
1901    /* analyze C */
1902    if ( MATELEM(C,1,2) == NULL ) 
1903      pN = NULL; /* check the consistency later */
1904    else 
1905      pN = p_GetCoeff(MATELEM(C,1,2),r);
1906    tmpIsSkewConstant = 1;
1907    for(i=1; i<r->N; i++)
1908    {
1909      for(j=i+1; j<=r->N; j++)
1910      { 
1911        if (MATELEM(C,i,j) == NULL)
1912          qN = NULL;
1913        else
1914          qN = p_GetCoeff(MATELEM(C,i,j),r);
1915        if ( qN == NULL )   /* check the consistency: Cij!=0 */
1916        // find also illegal pN
1917        {
1918          Werror("Incorrect input : matrix of coefficients contains zeros in the upper triangle");
1919          ncCleanUp(r);
1920          if (WeChangeRing)
1921            rChangeCurrRing(save);
1922          return TRUE;
1923        }
1924        if (!nEqual(pN,qN)) tmpIsSkewConstant = 0;
1925      }
1926    }
1927    r->nc->IsSkewConstant=tmpIsSkewConstant;
1928    if ( (tmpIsSkewConstant) && (nIsOne(pN)) ) 
1929    {
1930      r->nc->type = nc_lie;
1931    }
1932    else 
1933    {
1934      r->nc->type = nc_general;
1935    }
1936  }
1937
1938  /* initialition of the matrix D */
1939  if ( DD == NULL ) 
1940    /* we treat DN only (it could also be NULL) */
1941  {
1942    D = mpNew(r->N,r->N);
1943    if (DN  == NULL)
1944    {
1945      if ( (r->nc->type == nc_lie) || (r->nc->type == nc_undef) ) 
1946      {
1947        r->nc->type = nc_comm; /* it was nc_skew earlier */
1948      }
1949      else /* nc_general, nc_skew */
1950      {
1951        r->nc->type = nc_skew;
1952      }
1953    }
1954    else /* DN  != NULL */
1955    { 
1956      for(i=1; i<r->N; i++)
1957      {
1958        for(j=i+1; j<=r->N; j++)
1959        {
1960          MATELEM(D,i,j) = nc_p_CopyPut(DN,r);
1961        }
1962      }
1963    }
1964  }
1965  else /* DD != NULL */
1966  { 
1967    D = mpCopy(DD); 
1968  }
1969  /* analyze D */ 
1970  /* check the ordering condition for D (both matrix and poly cases) */
1971
1972  if ( nc_CheckOrdCondition(D, r) )
1973  {
1974    ncCleanUp(r);
1975    if (WeChangeRing)
1976      rChangeCurrRing(save);
1977    Werror("Matrix of polynomials violates the ordering condition");
1978    return TRUE;
1979  }
1980  r->nc->C = C;
1981  r->nc->D = D;
1982  if (WeChangeRing)
1983    rChangeCurrRing(save);
1984  return nc_InitMultiplication(r);
1985}
1986
1987BOOLEAN nc_InitMultiplication(ring r)
1988{
1989  /* returns TRUE if there were errors */
1990  /* initialize the multiplication: */
1991  /*  r->nc->MTsize, r->nc->MT, r->nc->COM, */
1992  /* and r->nc->IsSkewConstant for the skew case */
1993  if (rVar(r)==1)
1994  {
1995    r->nc->type=nc_comm;
1996    r->nc->IsSkewConstant=1;
1997    return FALSE;
1998  }
1999  ring save = currRing;
2000  int WeChangeRing = 0;
2001  if (currRing!=r)
2002  {
2003    rChangeCurrRing(r);
2004    WeChangeRing = 1;
2005  }
2006  int i,j;
2007  r->nc->MT = (matrix *)omAlloc0((r->N*(r->N-1))/2*sizeof(matrix));
2008  r->nc->MTsize = (int *)omAlloc0((r->N*(r->N-1))/2*sizeof(int));
2009  idTest(((ideal)r->nc->C));
2010  matrix COM = mpCopy(r->nc->C);
2011  poly p,q;
2012  short DefMTsize=7;
2013  int IsNonComm=0;
2014  int tmpIsSkewConstant;
2015 
2016  for(i=1; i<r->N; i++)
2017  {
2018    for(j=i+1; j<=r->N; j++)
2019    {
2020      if ( MATELEM(r->nc->D,i,j) == NULL ) /* quasicommutative case */
2021      {
2022        /* 1x1 mult.matrix */
2023        r->nc->MTsize[UPMATELEM(i,j,r->N)] = 1;
2024        r->nc->MT[UPMATELEM(i,j,r->N)] = mpNew(1,1);
2025      }
2026      else /* pure noncommutative case */
2027      {
2028        /* TODO check the special multiplication properties */
2029        IsNonComm = 1;
2030        p_Delete(&(MATELEM(COM,i,j)),r);
2031        //MATELEM(COM,i,j) = NULL; // done by p_Delete
2032        r->nc->MTsize[UPMATELEM(i,j,r->N)] = DefMTsize; /* default sizes */
2033        r->nc->MT[UPMATELEM(i,j,r->N)] = mpNew(DefMTsize, DefMTsize);
2034      }
2035      /* set MT[i,j,1,1] to c_i_j*x_i*x_j + D_i_j */
2036      p = p_ISet(1,r); /* instead of     p = pOne(); */
2037      p_SetCoeff(p,n_Copy(pGetCoeff(MATELEM(r->nc->C,i,j)),r),r);
2038      p_SetExp(p,i,1,r);
2039      p_SetExp(p,j,1,r);
2040      p_Setm(p,r);
2041      p_Test(MATELEM(r->nc->D,i,j),r);
2042      q =  nc_p_CopyGet(MATELEM(r->nc->D,i,j),r);
2043      p = p_Add_q(p,q,r);
2044      MATELEM(r->nc->MT[UPMATELEM(i,j,r->N)],1,1) = nc_p_CopyPut(p,r);
2045      p_Delete(&p,r);
2046      // p = NULL;// done by p_Delete
2047    }
2048  }
2049  if (r->nc->type==nc_undef)
2050  {
2051    if (IsNonComm==1)
2052    {
2053      //      assume(pN!=NULL);
2054      //      if ((tmpIsSkewConstant==1) && (nIsOne(pGetCoeff(pN)))) r->nc->type=nc_lie;
2055      //      else r->nc->type=nc_general;
2056    }
2057    if (IsNonComm==0) 
2058    {
2059      r->nc->type=nc_skew; /* TODO: check whether it is commutative */
2060      r->nc->IsSkewConstant=tmpIsSkewConstant;
2061    }
2062  }
2063  r->nc->COM=COM;
2064  if (WeChangeRing)
2065  {
2066    rChangeCurrRing(save);
2067  }
2068  return FALSE;
2069}
2070
2071/* substitute the n-th variable by e in p
2072* destroy p
2073* e is not a constant
2074*/
2075poly nc_pSubst(poly p, int n, poly e)
2076{
2077  int rN=currRing->N;
2078  int *PRE = (int *)omAlloc0((rN+1)*sizeof(int));
2079  int *SUF = (int *)omAlloc0((rN+1)*sizeof(int));
2080  int i,j,pow;
2081  number C;
2082  poly suf,pre;
2083  poly res = NULL;
2084  poly out = NULL;
2085  while ( p!= NULL )
2086  {
2087    C =  pGetCoeff(p);
2088    pGetExpV(p, PRE); /* faster splitting? */
2089    pow = PRE[n]; PRE[n]=0;
2090    res = NULL;
2091    if (pow!=0)
2092    {
2093      for (i=n+1; i<=rN; i++)
2094      {
2095        SUF[i] = PRE[i];
2096        PRE[i] = 0;
2097      }
2098      res =  pPower(pCopy(e),pow);
2099      /* multiply with prefix */
2100      pre = pOne();
2101      pSetExpV(pre,PRE);
2102      pSetm(pre);
2103      res = nc_mm_Mult_p(pre,res,currRing);
2104      /* multiply with suffix */
2105      suf = pOne();
2106      pSetExpV(suf,SUF);
2107      pSetm(suf);
2108      res = nc_p_Mult_mm(res,suf,currRing);
2109      res = p_Mult_nn(res,C,currRing);
2110      pSetComp(res,PRE[0]);
2111    }
2112    else /* pow==0 */
2113    {
2114      res = pHead(p);
2115    }
2116    p   = pLmDeleteAndNext(p);
2117    out = pAdd(out,res);
2118  }
2119  freeT(PRE,rN);
2120  freeT(SUF,rN);
2121  return(out);
2122}
2123
2124static ideal idPrepareStd(ideal T, ideal s,  int k)
2125{
2126  /* T is a left SB, without zeros, s is a list with zeros */
2127#ifdef PDEBUG
2128  if (IDELEMS(s)!=IDELEMS(T))
2129  {
2130    Print("ideals of diff. size!!!");
2131  }
2132#endif
2133  ideal t = idCopy(T);
2134  int j,rs=idRankFreeModule(s),rt=idRankFreeModule(t);
2135  poly p,q;
2136
2137  ideal res = idInit(2*idElem(t),1+idElem(t));
2138  if (rs == 0)
2139  {
2140    for (j=0; j<IDELEMS(t); j++)
2141    {
2142      if (s->m[j]!=NULL) pSetCompP(s->m[j],1);
2143      if (t->m[j]!=NULL) pSetCompP(t->m[j],1);
2144    }
2145    k = si_max(k,1);
2146  }
2147  for (j=0; j<IDELEMS(t); j++)
2148  {
2149    if (s->m[j]!=NULL)
2150    {
2151      p = s->m[j];
2152      q = pOne();
2153      pSetComp(q,k+1+j);
2154      pSetmComp(q);
2155#if 0     
2156      while (pNext(p)) pIter(p);
2157      pNext(p) = q;
2158#else
2159      p = pAdd(p,q);
2160      s->m[j] = p;
2161#ifdef PDEBUG
2162    pTest(p);
2163#endif
2164#endif
2165    }
2166  }
2167  res = idSimpleAdd(t,s);
2168  idDelete(&t);
2169  res->rank = 1+idElem(T);
2170  return(res);
2171}
2172
2173ideal Approx_Step(ideal L)
2174{
2175  int N=currRing->N;
2176  int i,j; // k=syzcomp
2177  int flag, flagcnt, syzcnt=0;
2178  int syzcomp = 0;
2179  int k=1; /* for ideals not modules */
2180  ideal I = kStd(L, currQuotient,testHomog,NULL,NULL,0,0,NULL);
2181  idSkipZeroes(I);
2182  ideal s_I;
2183  int idI = idElem(I);
2184  ideal trickyQuotient,s_trickyQuotient;
2185  if (currQuotient !=NULL)
2186  {
2187    trickyQuotient = idSimpleAdd(currQuotient,I);
2188  }
2189  else
2190    trickyQuotient = I;
2191  idSkipZeroes(trickyQuotient);
2192  poly *var = (poly *)omAlloc0((N+1)*sizeof(poly));
2193  //  poly *W = (poly *)omAlloc0((2*N+1)*sizeof(poly));
2194  resolvente S = (resolvente)omAlloc0((N+1)*sizeof(ideal));
2195  ideal SI, res;
2196  matrix MI;
2197  poly x=pOne();
2198  var[0]=x;
2199  ideal   h2, h3, s_h2, s_h3;
2200  poly    p,q,qq;
2201  /* init vars */
2202  for (i=1; i<=N; i++ )
2203  {
2204    x = pOne();
2205    pSetExp(x,i,1);
2206    pSetm(x);
2207    var[i]=pCopy(x);
2208  }
2209  /* init NF's */
2210  for (i=1; i<=N; i++ )
2211  {
2212    h2 = idInit(idI,1);
2213    flag = 0;
2214    for (j=0; j< idI; j++ )
2215    {
2216      q = nc_p_Mult_mm(pCopy(I->m[j]),var[i],currRing);
2217      q = kNF(I,currQuotient,q,0,0);
2218      if (q!=0)
2219      {
2220        h2->m[j]=pCopy(q);
2221        //      pShift(&(h2->m[flag]),1);
2222        flag++;
2223        pDelete(&q);
2224      }
2225      else
2226        h2->m[j]=0;
2227    }
2228    /* W[1..idElems(I)] */
2229    if (flag >0)
2230    {
2231      /* compute syzygies with values in I*/
2232      //      idSkipZeroes(h2);
2233      //      h2 = idSimpleAdd(h2,I);
2234      //      h2->rank=flag+idI+1;
2235      idTest(h2);
2236      idShow(h2);
2237      ring orig_ring=currRing;
2238      ring syz_ring=rCurrRingAssure_SyzComp();
2239      syzcomp = 1;
2240      rSetSyzComp(syzcomp);
2241      if (orig_ring != syz_ring)
2242      {
2243        s_h2=idrCopyR_NoSort(h2,orig_ring);
2244        //      s_trickyQuotient=idrCopyR_NoSort(trickyQuotient,orig_ring);
2245        //      rDebugPrint(syz_ring);
2246        s_I=idrCopyR_NoSort(I,orig_ring);
2247      }
2248      else
2249      {
2250        s_h2 = h2;
2251        s_I  = I;
2252        //      s_trickyQuotient=trickyQuotient;
2253      }
2254      idTest(s_h2);
2255      //      idTest(s_trickyQuotient);
2256      Print(".proceeding with the variable %d\n",i);
2257      s_h3 = idPrepareStd(s_I, s_h2, 1);
2258      BITSET save_test=test;
2259      test|=Sy_bit(OPT_SB_1);
2260      idTest(s_h3);
2261      idDelete(&s_h2);
2262      s_h2=idCopy(s_h3);
2263      idDelete(&s_h3);
2264      Print("...computing Syz");
2265      s_h3 = kStd(s_h2, currQuotient,(tHomog)FALSE,NULL,NULL,syzcomp,idI);
2266      test=save_test;
2267      idShow(s_h3);
2268      if (orig_ring != syz_ring)
2269      {
2270        idDelete(&s_h2);
2271        for (j=0; j<IDELEMS(s_h3); j++)
2272        {
2273          if (s_h3->m[j] != NULL)
2274          {
2275            if (p_MinComp(s_h3->m[j],syz_ring) > syzcomp) /* i.e. it is a syzygy */
2276              pShift(&s_h3->m[j], -syzcomp);
2277            else
2278              pDelete(&s_h3->m[j]);
2279          }
2280        }
2281        idSkipZeroes(s_h3);
2282        s_h3->rank -= syzcomp;
2283        rChangeCurrRing(orig_ring);
2284        //      s_h3 = idrMoveR_NoSort(s_h3, syz_ring);
2285        s_h3 = idrMoveR_NoSort(s_h3, syz_ring);
2286        rKill(syz_ring);
2287      }
2288      idTest(s_h3);
2289      S[syzcnt]=kStd(s_h3,currQuotient,(tHomog)FALSE,NULL,NULL);
2290      syzcnt++;
2291      idDelete(&s_h3);
2292    } /* end if flag >0 */
2293    else 
2294    {
2295      flagcnt++;
2296    }
2297  }
2298  if (flagcnt == N) 
2299  {
2300    Print("the input is a two--sided ideal");
2301    return(I);
2302  }
2303  if (syzcnt >0)
2304  {
2305    Print("..computing Intersect of %d modules\n",syzcnt);
2306    if (syzcnt == 1)
2307      SI = S[0];
2308    else
2309      SI = idMultSect(S, syzcnt);
2310    idShow(SI);
2311    MI = idModule2Matrix(SI);
2312    res= idInit(MATCOLS(MI),1);
2313    for (i=1; i<= MATCOLS(MI); i++)
2314    {   
2315      p = NULL;
2316      for (j=0; j< idElem(I); j++)
2317      { 
2318        q = pCopy(MATELEM(MI,j+1,i));
2319        if (q!=NULL)
2320        {
2321          q = pMult(q,pCopy(I->m[j]));
2322          p = pAdd(p,q);
2323        }
2324      }
2325      res->m[i-1]=p;
2326    }
2327    Print("final std");
2328    res = kStd(res, currQuotient,testHomog,NULL,NULL,0,0,NULL);
2329    idSkipZeroes(res);
2330    return(res);
2331  }
2332  else
2333  {
2334    Print("No syzygies");
2335    return(I);
2336  }
2337}
2338
2339
2340ring nc_rCreateNCcomm(ring r)
2341  /* creates a commutative nc extension; "converts" comm.ring to a Plural ring */
2342{
2343  if (rIsPluralRing(r)) return r;
2344  ring save = currRing;
2345  int WeChangeRing = 0;
2346  if (currRing!=r)
2347  {
2348    rChangeCurrRing(r);
2349    WeChangeRing = 1;
2350  }
2351  r->nc = (nc_struct *)omAlloc0(sizeof(nc_struct));
2352  r->nc->ref = 1;
2353  r->nc->basering = r;
2354  r->nc->type = nc_comm;
2355  r->nc->IsSkewConstant = 1;
2356  matrix C = mpNew(r->N,r->N);
2357  matrix D = mpNew(r->N,r->N);
2358  int i,j;
2359  for(i=1; i<r->N; i++)
2360  {
2361    for(j=i+1; j<=r->N; j++)
2362    {
2363      MATELEM(C,i,j) = pOne();
2364    }
2365  }
2366  r->nc->C = C;
2367  r->nc->D = D;
2368  if (nc_InitMultiplication(r))
2369  {
2370    WarnS("Error initializing multiplication!");
2371  }
2372  if (WeChangeRing)
2373  {
2374    rChangeCurrRing(save);
2375  }
2376  return r;
2377}
2378
2379poly p_CopyEmbed(poly p, ring srcRing, int shift, int par_shift)
2380  /* NOT USED ANYMORE: replaced by maFindPerm in ring.cc */
2381  /* for use with embeddings: currRing is a sum of smaller rings */
2382  /* and srcRing is one of such smaller rings */
2383  /* shift defines the position of a subring in srcRing */
2384  /* par_shift defines the position of a subfield in basefield of CurrRing */
2385{
2386  if (currRing == srcRing)
2387  {
2388    return(p_Copy(p,currRing));
2389  }
2390  nMapFunc nMap=nSetMap(srcRing);
2391  poly q;
2392  //  if ( nMap == nCopy)
2393  //  {
2394  //    q = prCopyR(p,srcRing);
2395  //  }
2396  //  else
2397  {
2398    int *perm = (int *)omAlloc0((srcRing->N+1)*sizeof(int));
2399    int *par_perm = (int *)omAlloc0((srcRing->P+1)*sizeof(int));
2400    //    int *par_perm = (int *)omAlloc0((srcRing->P+1)*sizeof(int));
2401    int i;
2402    //    if (srcRing->P > 0)
2403    //    {
2404    //      for (i=0; i<srcRing->P; i++)
2405    //  par_perm[i]=-i;
2406    //    }
2407    if ((shift<0) || (shift > currRing->N))
2408    {
2409      Werror("bad shifts in p_CopyEmbed");
2410      return(0);
2411    }
2412    for (i=1; i<= srcRing->N; i++)
2413    {
2414      perm[i] = shift+i;
2415    }
2416    q = pPermPoly(p,perm,srcRing,nMap,par_perm,srcRing->P);
2417  }
2418  return(q);
2419}
2420
2421poly pOppose(ring Rop, poly p)
2422  /* opposes a vector p from Rop to currRing */
2423{
2424  /* the simplest case:*/
2425  if (  Rop == currRing )  return(pCopy(p));
2426  /* check Rop == rOpposite(currRing) */
2427  if ( !rIsLikeOpposite(currRing, Rop) )
2428  {
2429    WarnS("an opposite ring should be used");
2430    return NULL;
2431  }
2432  /* nMapFunc nMap = nSetMap(Rop);*/
2433  /* since we know that basefields coinside! */
2434  int *perm=(int *)omAlloc0((Rop->N+1)*sizeof(int));
2435  if (!p_IsConstantPoly(p, Rop))
2436  {
2437    /* we know perm exactly */
2438    int i;
2439    for(i=1; i<=Rop->N; i++)
2440    {
2441      perm[i] = Rop->N+1-i;
2442    }
2443  }
2444  poly res = pPermPoly(p, perm, Rop, nCopy);
2445  omFreeSize((ADDRESS)perm,(Rop->N+1)*sizeof(int));
2446  return res;
2447}
2448
2449ideal idOppose(ring Rop, ideal I)
2450  /* opposes a module I from Rop to currRing */
2451{
2452  /* the simplest case:*/
2453  if ( Rop == currRing ) return idCopy(I);
2454  /* check Rop == rOpposite(currRing) */
2455  if (!rIsLikeOpposite(currRing, Rop))
2456  {
2457    WarnS("an opposite ring should be used");
2458    return NULL;
2459  }
2460  int i;
2461  ideal idOp = idInit(I->ncols, I->rank);
2462  for (i=0; i< (I->ncols)*(I->nrows); i++)
2463  { 
2464    idOp->m[i] = pOppose(Rop,I->m[i]); 
2465  }
2466  idTest(idOp);
2467  return idOp;
2468}
2469
2470BOOLEAN rIsLikeOpposite(ring rBase, ring rCandidate)
2471  /* checks whether rings rBase and rCandidate */
2472  /* could be opposite to each other */
2473  /* returns TRUE if it is so */
2474{
2475  /* the same basefield */
2476  int diagnose = TRUE;
2477  ring save = currRing;
2478  rChangeCurrRing(rBase);
2479  nMapFunc nMap = nSetMap(rCandidate);
2480  if (nMap != nCopy) diagnose = FALSE;
2481  rChangeCurrRing(save);
2482  /* same number of variables */
2483  if (rBase->N != rCandidate->N) diagnose = FALSE;
2484  /* nc and comm ring */
2485  if ( rIsPluralRing(rBase) != rIsPluralRing(rCandidate) ) diagnose = FALSE;
2486  /* both are qrings */
2487  /* NO CHECK, since it is used in building opposite qring */
2488  /*  if ( ((rBase->qideal != NULL) && (rCandidate->qideal == NULL)) */
2489  /*       || ((rBase->qideal == NULL) && (rCandidate->qideal != NULL)) ) */
2490  /*  diagnose = FALSE; */
2491  /* TODO: varnames are e->E etc */
2492  return diagnose;
2493}
2494
2495#endif
2496
2497
2498// int Commutative_Context(ring r, leftv expression)
2499//   /* returns 1 if expression consists */
2500//   /*  of commutative elements */
2501// {
2502//   /* crucial: poly -> ideal, module, matrix  */
2503// }
2504
2505// int Comm_Context_Poly(ring r, poly p)
2506// {
2507//   poly COMM=r->nc->COMM;
2508//   poly pp=pOne();
2509//   memset(pp->exp,0,r->ExpL_Size*sizeof(long));
2510//   while (p!=NULL)
2511//   {
2512//     for (i=0;i<=r->ExpL_Size;i++)
2513//     {
2514//       if ((p->exp[i]) && (pp->exp[i]))  return(FALSE);
2515//       /* nonzero exponent of non-comm variable */
2516//     }
2517//     pIter(p);
2518//   }
2519//   return(TRUE);
2520// }
Note: See TracBrowser for help on using the repository browser.