source: git/kernel/gring.cc @ a41623

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