source: git/kernel/gring.cc @ 0312c5

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