source: git/kernel/gring.cc @ e5a4ba

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