source: git/kernel/gring.cc @ 3c8a31

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