source: git/kernel/gring.cc @ 9f73706

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