source: git/libpolys/polys/nc/gring.cc @ e5d267

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