source: git/kernel/gring.cc @ 84d05f8

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