source: git/kernel/gring.cc @ 5a9e7b

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