source: git/kernel/gring.cc @ 1c473f

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