source: git/libpolys/polys/monomials/ring.cc @ 01c1d0

spielwiese
Last change on this file since 01c1d0 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: 142.4 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id$ */
5
6/*
7* ABSTRACT - the interpreter related ring operations
8*/
9
10/* includes */
11#include <math.h>
12
13// #include <polys/options.h>
14#include <omalloc/omalloc.h>
15#include <misc/options.h>
16#include <polys/monomials/ring.h>
17#include <polys/monomials/p_polys.h>
18#include <polys/simpleideals.h>
19#include <coeffs/numbers.h>
20// #include <???/febase.h>
21// #include <???/intvec.h>
22#include <polys/ext_fields/longalg.h>
23#include <polys/ext_fields/longtrans.h>
24#include <coeffs/ffields.h>
25#include <polys/monomials/ring.h>
26#include <polys/monomials/maps.h>
27#include <polys/prCopy.h>
28// #include "../Singular/ipshell.h"
29#include <polys/templates/p_Procs.h>
30
31#include <polys/matpol.h>
32
33#ifdef HAVE_PLURAL
34#include <polys/nc/nc.h>
35#include <polys/nc/sca.h>
36#endif
37// #include <???/maps.h>
38// #include <???/matpol.h>
39#ifdef HAVE_FACTORY
40#define SI_DONT_HAVE_GLOBAL_VARS
41#  include <factory/factory.h>
42#endif
43
44#define BITS_PER_LONG 8*SIZEOF_LONG
45
46omBin sip_sring_bin = omGetSpecBin(sizeof(ip_sring));
47
48static const char * const ringorder_name[] =
49{
50  " ?", ///< ringorder_no = 0,
51  "a", ///< ringorder_a,
52  "A", ///< ringorder_a64,
53  "c", ///< ringorder_c,
54  "C", ///< ringorder_C,
55  "M", ///< ringorder_M,
56  "S", ///< ringorder_S,
57  "s", ///< ringorder_s,
58  "lp", ///< ringorder_lp,
59  "dp", ///< ringorder_dp,
60  "rp", ///< ringorder_rp,
61  "Dp", ///< ringorder_Dp,
62  "wp", ///< ringorder_wp,
63  "Wp", ///< ringorder_Wp,
64  "ls", ///< ringorder_ls,
65  "ds", ///< ringorder_ds,
66  "Ds", ///< ringorder_Ds,
67  "ws", ///< ringorder_ws,
68  "Ws", ///< ringorder_Ws,
69  "L", ///< ringorder_L,
70  "aa", ///< ringorder_aa
71  "rs", ///< ringorder_rs,
72  "IS", ///<  ringorder_IS
73  " _" ///< ringorder_unspec
74};
75
76#define  rInternalChar(r) ((r)->cf->ch)
77
78const char * rSimpleOrdStr(int ord)
79{
80  return ringorder_name[ord];
81}
82
83/// unconditionally deletes fields in r
84void rDelete(ring r);
85/// set r->VarL_Size, r->VarL_Offset, r->VarL_LowIndex
86static void rSetVarL(ring r);
87/// get r->divmask depending on bits per exponent
88static unsigned long rGetDivMask(int bits);
89/// right-adjust r->VarOffset
90static void rRightAdjustVarOffset(ring r);
91static void rOptimizeLDeg(ring r);
92
93/*0 implementation*/
94//BOOLEAN rField_is_R(ring r)
95//{
96//  if (r->cf->ch== -1)
97//  {
98//    if (r->float_len==(short)0) return TRUE;
99//  }
100//  return FALSE;
101//}
102
103/*
104
105/// internally changes the gloabl ring and resets the relevant
106/// global variables:
107/// SHOULD BE DEPRECATED NOW...?
108void rChangeCurrRing(ring r)
109{
110 // if ((currRing!=NULL) && (currRing->cf->minpoly!=NULL))
111 // {
112 //   omCheckAddr(currRing->cf->minpoly);
113 // }
114  //------------ set global ring vars --------------------------------
115  //currRing = r;
116  //currQuotient=NULL;
117  if (r != NULL)
118  {
119    rTest(r);
120    //------------ set global ring vars --------------------------------
121    //currQuotient=r->qideal;
122
123    //------------ global variables related to coefficients ------------
124    nSetChar(r->cf);
125
126    //------------ global variables related to polys -------------------
127    p_SetGlobals(r);
128    //------------ global variables related to factory -----------------
129#ifdef HAVE_FACTORY
130    //int c=ABS(nGetChar());
131    //if (c==1) c=0;
132    //setCharacteristic( c );
133#endif
134  }
135}
136*/
137
138ring rDefault(coeffs cf, int N, char **n,int ord_size, int *ord, int *block0, int *block1)
139{
140  assume( cf != NULL);
141  ring r=(ring) omAlloc0Bin(sip_sring_bin);
142  r->cf->ch    = n_GetChar(cf);
143  r->N     = N;
144  r->cf = cf;
145  /*rPar(r)  = 0; Alloc0 */
146  /*names*/
147  r->names = (char **) omAlloc0(N * sizeof(char *));
148  int i;
149  for(i=0;i<N;i++)
150  {
151    r->names[i]  = omStrDup(n[i]);
152  }
153  /*weights: entries for 2 blocks: NULL*/
154  r->wvhdl = (int **)omAlloc0((ord_size+1) * sizeof(int *));
155  r->order = ord;
156  r->block0 = block0;
157  r->block1 = block1;
158  /*polynomial ring*/
159  r->OrdSgn    = 1;
160
161  /* complete ring intializations */
162  rComplete(r);
163  return r;
164}
165ring rDefault(int ch, int N, char **n,int ord_size, int *ord, int *block0, int *block1)
166{
167  coeffs cf;
168  if (ch==0) cf=nInitChar(n_Q,NULL);
169  else       cf=nInitChar(n_Zp,(void*)(long)ch);
170  assume( cf != NULL);
171  return rDefault(cf,N,n,ord_size,ord,block0,block1);
172}
173ring   rDefault(coeffs cf, int N, char **n)
174{ 
175  assume( cf != NULL);
176  /*order: lp,0*/
177  int *order = (int *) omAlloc(2* sizeof(int));
178  int *block0 = (int *)omAlloc0(2 * sizeof(int));
179  int *block1 = (int *)omAlloc0(2 * sizeof(int));
180  /* ringorder dp for the first block: var 1..N */
181  order[0]  = ringorder_lp;
182  block0[0] = 1;
183  block1[0] = N;
184  /* the last block: everything is 0 */
185  order[1]  = 0;
186
187  return rDefault(cf,N,n,2,order,block0,block1);
188}
189
190ring rDefault(int ch, int N, char **n)
191{
192  coeffs cf;
193  if (ch==0) cf=nInitChar(n_Q,NULL);
194  else       cf=nInitChar(n_Zp,(void*)(long)ch);
195  assume( cf != NULL);
196  return rDefault(cf,N,n);
197}
198
199///////////////////////////////////////////////////////////////////////////
200//
201// rInit: define a new ring from sleftv's
202//
203//-> ipshell.cc
204
205/////////////////////////////
206// Auxillary functions
207//
208
209// check intvec, describing the ordering
210BOOLEAN rCheckIV(intvec *iv)
211{
212  if ((iv->length()!=2)&&(iv->length()!=3))
213  {
214    WerrorS("weights only for orderings wp,ws,Wp,Ws,a,M");
215    return TRUE;
216  }
217  return FALSE;
218}
219
220int rTypeOfMatrixOrder(intvec * order)
221{
222  int i=0,j,typ=1;
223  int sz = (int)sqrt((double)(order->length()-2));
224  if ((sz*sz)!=(order->length()-2))
225  {
226    WerrorS("Matrix order is not a square matrix");
227    typ=0;
228  }
229  while ((i<sz) && (typ==1))
230  {
231    j=0;
232    while ((j<sz) && ((*order)[j*sz+i+2]==0)) j++;
233    if (j>=sz)
234    {
235      typ = 0;
236      WerrorS("Matrix order not complete");
237    }
238    else if ((*order)[j*sz+i+2]<0)
239      typ = -1;
240    else
241      i++;
242  }
243  return typ;
244}
245
246/*2
247 * set a new ring from the data:
248 s: name, chr: ch, varnames: rv, ordering: ord, typ: typ
249 */
250
251int r_IsRingVar(const char *n, ring r)
252{
253  if ((r!=NULL) && (r->names!=NULL))
254  {
255    for (int i=0; i<r->N; i++)
256    {
257      if (r->names[i]==NULL) return -1;
258      if (strcmp(n,r->names[i]) == 0) return (int)i;
259    }
260  }
261  return -1;
262}
263
264
265void rWrite(ring r)
266{
267  if ((r==NULL)||(r->order==NULL))
268    return; /*to avoid printing after errors....*/
269
270  int nblocks=rBlocks(r);
271
272  // omCheckAddrSize(r,sizeof(ip_sring));
273  omCheckAddrSize(r->order,nblocks*sizeof(int));
274  omCheckAddrSize(r->block0,nblocks*sizeof(int));
275  omCheckAddrSize(r->block1,nblocks*sizeof(int));
276  omCheckAddrSize(r->wvhdl,nblocks*sizeof(int *));
277  omCheckAddrSize(r->names,r->N*sizeof(char *));
278
279  nblocks--;
280
281  n_CoeffWrite(r->cf);
282#if 0
283  if (rField_is_GF(r))
284  {
285    Print("//   # ground field : %d\n",r->cf->ch);
286    Print("//   primitive element : %s\n", r->parameter[0]);
287    StringSetS("//   minpoly        : ");
288    nfShowMipo(r->cf);PrintS(StringAppendS("\n"));
289  }
290#ifdef HAVE_RINGS
291  else if (rField_is_Ring(r))
292  {
293    PrintS("//   coeff. ring is : ");
294    if (rField_is_Ring_Z(r)) PrintS("Integers\n");
295    long l = (long)mpz_sizeinbase(r->cf->modBase, 10) + 2;
296    char* s = (char*) omAlloc(l);
297    mpz_get_str(s,10,r->cf->modBase);
298    if (rField_is_Ring_ModN(r)) Print("Z/%s\n", s);
299    if (rField_is_Ring_2toM(r)) Print("Z/2^%lu\n", r->cf->modExponent);
300    if (rField_is_Ring_PtoM(r)) Print("Z/%s^%lu\n", s, r->cf->modExponent);
301    omFreeSize((ADDRESS)s, l);
302  }
303#endif
304  else
305  {
306    PrintS("//   characteristic : ");
307    if ( rField_is_R(r) )             PrintS("0 (real)\n");  /* R */
308    else if ( rField_is_long_R(r) )
309      Print("0 (real:%d digits, additional %d digits)\n",
310             r->float_len,r->float_len2);  /* long R */
311    else if ( rField_is_long_C(r) )
312      Print("0 (complex:%d digits, additional %d digits)\n",
313             r->float_len, r->float_len2);  /* long C */
314    else
315      Print ("%d\n",rChar(r)); /* Fp(a) */
316    if (r->cf->parameter!=NULL)
317    {
318      Print ("//   %d parameter    : ",rPar(r));
319      char **sp=r->cf->parameter;
320      int nop=0;
321      while (nop<rPar(r))
322      {
323        PrintS(*sp);
324        PrintS(" ");
325        sp++; nop++;
326      }
327      PrintS("\n//   minpoly        : ");
328      if ( rField_is_long_C(r) )
329      {
330        // i^2+1:
331        Print("(%s^2+1)\n",r->cf->parameter[0]);
332      }
333      else if (r->cf->minpoly==NULL)
334      {
335        PrintS("0\n");
336      }
337      else
338      {
339        StringSetS(""); n_Write(r->cf->minpoly,r->cf); PrintS(StringAppendS("\n"));
340      }
341      //if (r->minideal!=NULL)
342      //{
343      //  iiWriteMatrix((matrix)r->minideal,"//   minpolys",1,0);
344      //  PrintLn();
345      //}
346    }
347  }
348#endif 
349  Print("//   number of vars : %d",r->N);
350
351  //for (nblocks=0; r->order[nblocks]; nblocks++);
352  nblocks=rBlocks(r)-1;
353
354  for (int l=0, nlen=0 ; l<nblocks; l++)
355  {
356    int i;
357    Print("\n//        block %3d : ",l+1);
358
359    Print("ordering %s", rSimpleOrdStr(r->order[l]));
360
361
362    if (r->order[l] == ringorder_s)
363    {
364      assume( l == 0 );
365#ifndef NDEBUG
366      Print("  syzcomp at %d",r->typ[l].data.syz.limit);
367#endif
368      continue;
369    }
370    else if (r->order[l] == ringorder_IS)
371    {
372      assume( r->block0[l] == r->block1[l] );
373      const int s = r->block0[l];
374      assume( (-2 < s) && (s < 2) );
375      Print("(%d)", s); // 0 => prefix! +/-1 => suffix!
376      continue;
377    }
378    else if (
379    (  (r->order[l] >= ringorder_lp)
380    ||(r->order[l] == ringorder_M)
381    ||(r->order[l] == ringorder_a)
382    ||(r->order[l] == ringorder_a64)
383    ||(r->order[l] == ringorder_aa) ) && (r->order[l] < ringorder_IS) )
384    {
385      PrintS("\n//                  : names   ");
386      for (i = r->block0[l]-1; i<r->block1[l]; i++)
387      {
388        nlen = strlen(r->names[i]);
389        Print(" %s",r->names[i]);
390      }
391    }
392
393    if (r->wvhdl[l]!=NULL)
394    {
395      for (int j= 0;
396           j<(r->block1[l]-r->block0[l]+1)*(r->block1[l]-r->block0[l]+1);
397           j+=i)
398      {
399        PrintS("\n//                  : weights ");
400        for (i = 0; i<=r->block1[l]-r->block0[l]; i++)
401        {
402          if (r->order[l] == ringorder_a64)
403          {
404            int64 *w=(int64 *)r->wvhdl[l];
405            #if SIZEOF_LONG == 4
406            Print("%*lld " ,nlen,w[i+j]);
407            #else
408            Print(" %*ld"  ,nlen,w[i+j]);
409            #endif
410          }
411          else
412            Print(" %*d" ,nlen,r->wvhdl[l][i+j]);
413        }
414        if (r->order[l]!=ringorder_M) break;
415      }
416    }
417  }
418#ifdef HAVE_PLURAL
419  if(rIsPluralRing(r))
420  {
421    PrintS("\n//   noncommutative relations:");
422    poly pl=NULL;
423    int nl;
424    int i,j;
425    for (i = 1; i<r->N; i++)
426    {
427      for (j = i+1; j<=r->N; j++)
428      {
429        nl = n_IsOne(p_GetCoeff(MATELEM(r->GetNC()->C,i,j),r->GetNC()->basering), r->cf);
430        if ( (MATELEM(r->GetNC()->D,i,j)!=NULL) || (!nl) )
431        {
432          Print("\n//    %s%s=",r->names[j-1],r->names[i-1]);
433          pl = MATELEM(r->GetNC()->MT[UPMATELEM(i,j,r->N)],1,1);
434          p_Write0(pl, r, r);
435        }
436      }
437    }
438#if 0  /*Singularg should not differ from Singular except in error case*/
439    Print("\n//   noncommutative type:%d", (int)ncRingType(r));
440    Print("\n//      is skew constant:%d",r->GetNC()->IsSkewConstant);
441    if( rIsSCA(r) )
442    {
443      Print("\n//   alternating variables: [%d, %d]", scaFirstAltVar(r), scaLastAltVar(r));
444      const ideal Q = SCAQuotient(r); // resides within r!
445      PrintS("\n//   quotient of sca by ideal");
446
447      if (Q!=NULL)
448      {
449//        if (r==currRing)
450//        {
451//          PrintLn();
452//          iiWriteMatrix((matrix)Q,"scaQ",1);
453//        }
454//        else
455            PrintS(" ...");
456      }
457      else
458        PrintS(" (NULL)");
459    }
460#endif
461  }
462#endif
463  //if (r->qideal!=NULL)
464  //{
465  //  PrintS("\n// quotient ring from ideal\n");
466  //  iiWriteMatrix((matrix)r->qideal,"_",1);
467  //}
468}
469
470void rDelete(ring r)
471{
472  int i, j;
473
474  if (r == NULL) return;
475
476#ifdef HAVE_PLURAL
477  if (rIsPluralRing(r))
478    nc_rKill(r);
479#endif
480
481  nKillChar(r->cf); r->cf = NULL;
482  rUnComplete(r);
483  // delete order stuff
484  if (r->order != NULL)
485  {
486    i=rBlocks(r);
487    assume(r->block0 != NULL && r->block1 != NULL && r->wvhdl != NULL);
488    // delete order
489    omFreeSize((ADDRESS)r->order,i*sizeof(int));
490    omFreeSize((ADDRESS)r->block0,i*sizeof(int));
491    omFreeSize((ADDRESS)r->block1,i*sizeof(int));
492    // delete weights
493    for (j=0; j<i; j++)
494    {
495      if (r->wvhdl[j]!=NULL)
496        omFree(r->wvhdl[j]);
497    }
498    omFreeSize((ADDRESS)r->wvhdl,i*sizeof(int *));
499  }
500  else
501  {
502    assume(r->block0 == NULL && r->block1 == NULL && r->wvhdl == NULL);
503  }
504
505  // delete varnames
506  if(r->names!=NULL)
507  {
508    for (i=0; i<r->N; i++)
509    {
510      if (r->names[i] != NULL) omFree((ADDRESS)r->names[i]);
511    }
512    omFreeSize((ADDRESS)r->names,r->N*sizeof(char *));
513  }
514
515  // delete parameter
516  if (r->cf->parameter!=NULL)
517  {
518    char **s=r->cf->parameter;
519    j = 0;
520    while (j < rPar(r))
521    {
522      if (*s != NULL) omFree((ADDRESS)*s);
523      s++;
524      j++;
525    }
526    omFreeSize((ADDRESS)r->cf->parameter,rPar(r)*sizeof(char *));
527  }
528  omFreeBin(r, sip_sring_bin);
529}
530
531int rOrderName(char * ordername)
532{
533  int order=ringorder_unspec;
534  while (order!= 0)
535  {
536    if (strcmp(ordername,rSimpleOrdStr(order))==0)
537      break;
538    order--;
539  }
540  if (order==0) Werror("wrong ring order `%s`",ordername);
541  omFree((ADDRESS)ordername);
542  return order;
543}
544
545char * rOrdStr(ring r)
546{
547  if ((r==NULL)||(r->order==NULL)) return omStrDup("");
548  int nblocks,l,i;
549
550  for (nblocks=0; r->order[nblocks]; nblocks++);
551  nblocks--;
552
553  StringSetS("");
554  for (l=0; ; l++)
555  {
556    StringAppendS((char *)rSimpleOrdStr(r->order[l]));
557    if (
558           (r->order[l] != ringorder_c)
559        && (r->order[l] != ringorder_C)
560        && (r->order[l] != ringorder_s)
561        && (r->order[l] != ringorder_S)
562        && (r->order[l] != ringorder_IS)
563       )
564    {
565      if (r->wvhdl[l]!=NULL)
566      {
567        StringAppendS("(");
568        for (int j= 0;
569             j<(r->block1[l]-r->block0[l]+1)*(r->block1[l]-r->block0[l]+1);
570             j+=i+1)
571        {
572          char c=',';
573          if(r->order[l]==ringorder_a64)
574          {
575            int64 * w=(int64 *)r->wvhdl[l];
576            for (i = 0; i<r->block1[l]-r->block0[l]; i++)
577            {
578              StringAppend("%lld," ,w[i]);
579            }
580            StringAppend("%lld)" ,w[i]);
581            break;
582          }
583          else
584          {
585            for (i = 0; i<r->block1[l]-r->block0[l]; i++)
586            {
587              StringAppend("%d," ,r->wvhdl[l][i+j]);
588            }
589          }
590          if (r->order[l]!=ringorder_M)
591          {
592            StringAppend("%d)" ,r->wvhdl[l][i+j]);
593            break;
594          }
595          if (j+i+1==(r->block1[l]-r->block0[l]+1)*(r->block1[l]-r->block0[l]+1))
596            c=')';
597          StringAppend("%d%c" ,r->wvhdl[l][i+j],c);
598        }
599      }
600      else
601        StringAppend("(%d)",r->block1[l]-r->block0[l]+1);
602    }
603    else if (r->order[l] == ringorder_IS)
604    {
605      assume( r->block0[l] == r->block1[l] );
606      const int s = r->block0[l];
607      assume( (-2 < s) && (s < 2) );
608
609      StringAppend("(%d)", s);
610    }
611
612    if (l==nblocks) return omStrDup(StringAppendS(""));
613    StringAppendS(",");
614  }
615}
616
617char * rVarStr(ring r)
618{
619  if ((r==NULL)||(r->names==NULL)) return omStrDup("");
620  int i;
621  int l=2;
622  char *s;
623
624  for (i=0; i<r->N; i++)
625  {
626    l+=strlen(r->names[i])+1;
627  }
628  s=(char *)omAlloc((long)l);
629  s[0]='\0';
630  for (i=0; i<r->N-1; i++)
631  {
632    strcat(s,r->names[i]);
633    strcat(s,",");
634  }
635  strcat(s,r->names[i]);
636  return s;
637}
638
639char * rCharStr(ring r)
640{
641  char *s;
642  int i;
643
644#ifdef HAVE_RINGS
645  if (rField_is_Ring_Z(r))
646  {
647    s=omStrDup("integer");                   // Z
648    return s;
649  }
650  if(rField_is_Ring_2toM(r))
651  {
652    char* s = (char*) omAlloc(7+10+2);
653    sprintf(s,"integer,%lu",r->cf->modExponent);
654    return s;
655  }
656  if(rField_is_Ring_ModN(r))
657  {
658    long l = (long)mpz_sizeinbase(r->cf->modBase, 10) + 2+7;
659    char* s = (char*) omAlloc(l);
660    gmp_sprintf(s,"integer,%Zd",r->cf->modBase);
661    return s;
662  }
663  if(rField_is_Ring_PtoM(r))
664  {
665    long l = (long)mpz_sizeinbase(r->cf->modBase, 10) + 2+7+10;
666    char* s = (char*) omAlloc(l);
667    gmp_sprintf(s,"integer,%Zd^%lu",r->cf->modBase,r->cf->modExponent);
668    return s;
669  }
670#endif
671  if (r->cf->parameter==NULL)
672  {
673    i=r->cf->ch;
674    if(i==-1)
675      s=omStrDup("real");                    /* R */
676    else
677    {
678      s=(char *)omAlloc(MAX_INT_LEN+1);
679      sprintf(s,"%d",i);                   /* Q, Z/p */
680    }
681    return s;
682  }
683  if (rField_is_long_C(r))
684  {
685    s=(char *)omAlloc(21+strlen(r->cf->parameter[0]));
686    sprintf(s,"complex,%d,%s",r->float_len,r->cf->parameter[0]);   /* C */
687    return s;
688  }
689  int l=0;
690  for(i=0; i<rPar(r);i++)
691  {
692    l+=(strlen(r->cf->parameter[i])+1);
693  }
694  s=(char *)omAlloc((long)(l+MAX_INT_LEN+1));
695  s[0]='\0';
696  if (r->cf->ch<0)       sprintf(s,"%d",-r->cf->ch); /* Fp(a) */
697  else if (r->cf->ch==1) sprintf(s,"0");         /* Q(a)  */
698  else
699  {
700    sprintf(s,"%d,%s",r->cf->ch,r->cf->parameter[0]); /* GF(q)  */
701    return s;
702  }
703  char tt[2];
704  tt[0]=',';
705  tt[1]='\0';
706  for(i=0; i<rPar(r);i++)
707  {
708    strcat(s,tt);
709    strcat(s,r->cf->parameter[i]);
710  }
711  return s;
712}
713
714char * rParStr(ring r)
715{
716  if ((r==NULL)||(r->cf->parameter==NULL)) return omStrDup("");
717
718  int i;
719  int l=2;
720
721  for (i=0; i<rPar(r); i++)
722  {
723    l+=strlen(r->cf->parameter[i])+1;
724  }
725  char *s=(char *)omAlloc((long)l);
726  s[0]='\0';
727  for (i=0; i<rPar(r)-1; i++)
728  {
729    strcat(s,r->cf->parameter[i]);
730    strcat(s,",");
731  }
732  strcat(s,r->cf->parameter[i]);
733  return s;
734}
735
736char * rString(ring r)
737{
738  char *ch=rCharStr(r);
739  char *var=rVarStr(r);
740  char *ord=rOrdStr(r);
741  char *res=(char *)omAlloc(strlen(ch)+strlen(var)+strlen(ord)+9);
742  sprintf(res,"(%s),(%s),(%s)",ch,var,ord);
743  omFree((ADDRESS)ch);
744  omFree((ADDRESS)var);
745  omFree((ADDRESS)ord);
746  return res;
747}
748
749int  rIsExtension(const ring r)
750{
751  return (r->cf->parameter!=NULL); /* R, Q, Fp: FALSE */
752}
753
754static int binaryPower (const int a, const int b)
755{
756  /* computes a^b according to the binary representation of b,
757     i.e., a^7 = a^4 * a^2 * a^1. This saves some multiplications. */
758  int result = 1;
759  int factor = a;
760  int bb = b;
761  while (bb != 0)
762  {
763    if (bb % 2 != 0) result = result * factor;
764    bb = bb / 2;
765    factor = factor * factor;
766  }
767  return result;
768}
769
770int rChar(ring r)
771{
772#ifdef HAVE_RINGS
773  if (rField_is_Ring_2toM(r))
774    return binaryPower(2, (int)(unsigned long)r->cf->modExponent);
775  if (rField_is_Ring_ModN(r))
776    return (int)mpz_get_ui(r->cf->modBase);
777  if (rField_is_Ring_PtoM(r))
778    return binaryPower((int)mpz_get_ui(r->cf->modBase),
779                       (int)(unsigned long)r->cf->modExponent);
780 
781#endif
782  if (rField_is_numeric(r))
783    return 0;
784  if (!rIsExtension(r)) /* Q, Fp */
785    return r->cf->ch;
786  if (rField_is_Zp_a(r))  /* Fp(a)  */
787    return -r->cf->ch;
788  if (rField_is_Q_a(r))   /* Q(a)  */
789    return 0;
790  /*else*/               /* GF(p,n) */
791  {
792    if ((r->cf->ch & 1)==0) return 2;
793    int i=3;
794    while ((r->cf->ch % i)!=0) i+=2;
795    return i;
796  }
797}
798
799typedef char *             char_ptr;
800omBin char_ptr_bin = omGetSpecBin(sizeof(char_ptr)); // deallocation?
801
802
803// creates a commutative nc extension; "converts" comm.ring to a Plural ring
804#ifdef HAVE_PLURAL
805ring nc_rCreateNCcomm_rCopy(ring r)
806{
807  r = rCopy(r);
808  if (rIsPluralRing(r)) 
809    return r;
810
811  matrix C = mpNew(r->N,r->N); // ring-independent!?!
812  matrix D = mpNew(r->N,r->N);
813
814  for(int i=1; i<r->N; i++)
815    for(int j=i+1; j<=r->N; j++)
816      MATELEM(C,i,j) = p_One( r);
817
818  if (nc_CallPlural(C, D, NULL, NULL, r, false, true, false, r/*??currRing??*/, TRUE)) // TODO: what about quotient ideal?
819    WarnS("Error initializing multiplication!"); // No reaction!???
820
821  return r;
822}
823#endif
824
825
826/*2
827 *returns -1 for not compatible, (sum is undefined)
828 *         1 for compatible (and sum)
829 */
830/* vartest: test for variable/paramter names
831* dp_dp: for comm. rings: use block order dp + dp/ds/wp
832*/
833int rSumInternal(ring r1, ring r2, ring &sum, BOOLEAN vartest, BOOLEAN dp_dp)
834{
835
836  ip_sring tmpR;
837  memset(&tmpR,0,sizeof(tmpR));
838  /* check coeff. field =====================================================*/
839
840  TODO(Hans, switch to the use of ->cf);
841  assume( tmpR.cf !=NULL ); // ???!
842 
843  if ((rFieldType(r1)==rFieldType(r2))
844  && (rInternalChar(r1)==rInternalChar(r2)))
845  {
846    if (rField_is_Q(r1)||rField_is_Zp(r1)||rField_is_GF(r1)) /*Q, Z/p, GF(p,n)*/
847    {
848      if (r1->cf->parameter!=NULL)
849      {
850        if (!vartest || (strcmp(r1->cf->parameter[0],r2->cf->parameter[0])==0)) /* 1 par */
851        {
852          //tmpR.parameter=(char **)omAllocBin(char *);
853          //tmpR.parameter[0]=omStrDup(r1->cf->parameter[0]);
854          //tmpR.P=1;
855        }
856        else
857        {
858          WerrorS("GF(p,n)+GF(p,n)");
859          return -1;
860        }
861      }
862      tmpR.cf->ch=r1->cf->ch;
863    }
864    else if (rField_is_Extension(r1)) /* Q(a),Z/p(a) */
865    {
866      if (r1->cf->minpoly!=NULL)
867      {
868        if (r2->cf->minpoly!=NULL)
869        {
870          // HANNES: TODO: delete nSetChar
871         
872          // TODO: eliminate rChangeCurrRing
873//        rChangeCurrRing(r1);
874           
875          if ((!vartest || (strcmp(r1->cf->parameter[0],r2->cf->parameter[0])==0)) /* 1 par */
876              && n_Equal(r1->cf->minpoly,r2->cf->minpoly, r1->cf))
877          {
878            //tmpR.parameter=(char **)omAllocBin(char *);
879            //tmpR.parameter[0]=omStrDup(r1->cf->parameter[0]);
880            //tmpR.minpoly=n_Copy(r1->cf->minpoly, r1->cf);
881            //tmpR.P=1;
882            // HANNES: TODO: delete nSetChar
883          }
884          else
885          {
886            // HANNES: TODO: delete nSetChar
887            WerrorS("different minpolys");
888            return -1;
889          }
890        }
891        else
892        {
893          if ((!vartest || (strcmp(r1->cf->parameter[0],r2->cf->parameter[0])==0)) /* 1 par */
894              && (rPar(r2)==1))
895          {
896            //tmpR.parameter=(char **)omAllocBin(char *);
897            //tmpR.parameter[0]=omStrDup(r1->cf->parameter[0]);
898            //tmpR.P=1;
899            //tmpR.minpoly=n_Copy(r1->cf->minpoly, r1->cf);
900          }
901          else
902          {
903            WerrorS("different parameters and minpoly!=0");
904            return -1;
905          }
906        }
907      }
908      else /* r1->cf->minpoly==NULL */
909      {
910        if (r2->cf->minpoly!=NULL)
911        {
912          if ((!vartest || (strcmp(r1->cf->parameter[0],r2->cf->parameter[0])==0)) /* 1 par */
913              && (rPar(r1)==1))
914          {
915            //tmpR.parameter=(char **)omAllocBin(char *);
916            //tmpR.parameter[0]=omStrDup(r1->cf->parameter[0]);
917            //tmpR.P=1;
918            //tmpR.minpoly=n_Copy(r2->cf->minpoly, r2->cf);
919          }
920          else
921          {
922            WerrorS("different parameters and minpoly!=0");
923            return -1;
924          }
925        }
926        else
927        {
928          int len=rPar(r1)+rPar(r2);
929          tmpR.cf->parameter=(char **)omAlloc0(len*sizeof(char *));
930          int i;
931          for (i=0;i<rPar(r1);i++)
932          {
933            tmpR.cf->parameter[i]=omStrDup(r1->cf->parameter[i]);
934          }
935          int j,l;
936          for(j=0;j<rPar(r2);j++)
937          {
938            if (vartest)
939            {
940              for(l=0;l<i;l++)
941              {
942                if(strcmp(tmpR.cf->parameter[l],r2->cf->parameter[j])==0)
943                  break;
944              }
945            }
946            else
947              l=i;
948            if (l==i)
949            {
950              tmpR.cf->parameter[i]=omStrDup(r2->cf->parameter[j]);
951              i++;
952            }
953          }
954          if (i!=len)
955          {
956            tmpR.cf->parameter=(char**)omReallocSize(tmpR.cf->parameter,
957                                                 len*sizeof(char *),
958                                                 i*sizeof(char *));
959          }
960          tmpR.cf->P=i;
961        }
962      }
963    }
964    #ifdef HAVE_RINGS
965    else if (rField_is_Ring(r1)||rField_is_Ring(r2))
966    {
967      if (r1->ringtype != r2->ringtype)
968      {
969        Werror("rSumInternal not yet implemented for %s",
970               "different coefficient rings");
971        return -1;
972      }
973      else
974      {
975        tmpR.ch        = rInternalChar(r1);
976        tmpR.ringtype  = r1->ringtype;
977        if (r1->ringflaga != NULL)
978        {
979          omBin tmpBin = omGetSpecBin(sizeof(mpz_t));
980          tmpR.ringflaga = (int_number)omAllocBin(tmpBin);
981          mpz_init_set(tmpR.ringflaga, (int_number)r1->ringflaga);
982        }
983        tmpR.ringflagb = r1->ringflagb;
984        tmpR.nr2mModul = r1->nr2mModul;
985        if (r1->nrnModul != NULL)
986        {
987          omBin tmpBin = omGetSpecBin(sizeof(mpz_t));
988          tmpR.nrnModul = (int_number)omAllocBin(tmpBin);
989          mpz_init_set(tmpR.nrnModul, (int_number)r1->nrnModul);
990        }
991      }
992    }
993    #endif
994  }
995  else /* r1->cf->ch!=r2->cf->ch */
996  {
997    if (r1->cf->ch<-1) /* Z/p(a) */
998    {
999      if ((r2->cf->ch==0) /* Q */
1000          || (r2->cf->ch==-r1->cf->ch)) /* Z/p */
1001      {
1002        tmpR.cf->ch=r1->cf->ch;
1003        tmpR.cf->P=rPar(r1);
1004        tmpR.cf->parameter=(char **)omAlloc(rPar(r1)*sizeof(char *));
1005        int i;
1006        for (i=0;i<rPar(r1);i++)
1007        {
1008          tmpR.cf->parameter[i]=omStrDup(r1->cf->parameter[i]);
1009        }
1010        if (r1->cf->minpoly!=NULL)
1011        {
1012          tmpR.cf->minpoly=n_Copy(r1->cf->minpoly, r1->cf);
1013        }
1014      }
1015      else  /* R, Q(a),Z/q,Z/p(a),GF(p,n) */
1016      {
1017        WerrorS("Z/p(a)+(R,Q(a),Z/q(a),GF(q,n))");
1018        return -1;
1019      }
1020    }
1021    else if (r1->cf->ch==-1) /* R */
1022    {
1023      WerrorS("R+..");
1024      return -1;
1025    }
1026    else if (r1->cf->ch==0) /* Q */
1027    {
1028      if ((r2->cf->ch<-1)||(r2->cf->ch==1)) /* Z/p(a),Q(a) */
1029      {
1030        tmpR.cf->ch=r2->cf->ch;
1031        tmpR.cf->P=rPar(r2);
1032        tmpR.cf->parameter=(char **)omAlloc(rPar(r2)*sizeof(char *));
1033        int i;
1034        for (i=0;i<rPar(r2);i++)
1035        {
1036          tmpR.cf->parameter[i]=omStrDup(r2->cf->parameter[i]);
1037        }
1038        if (r2->cf->minpoly!=NULL)
1039        {
1040          tmpR.cf->minpoly=n_Copy(r2->cf->minpoly, r2->cf);
1041        }
1042      }
1043      else if (r2->cf->ch>1) /* Z/p,GF(p,n) */
1044      {
1045        tmpR.cf->ch=r2->cf->ch;
1046        if (r2->cf->parameter!=NULL)
1047        {
1048          tmpR.cf->parameter=(char_ptr *)omAllocBin(char_ptr_bin);
1049          tmpR.cf->P=1;
1050          tmpR.cf->parameter[0]=omStrDup(r2->cf->parameter[0]);
1051        }
1052      }
1053      else
1054      {
1055        WerrorS("Q+R");
1056        return -1; /* R */
1057      }
1058    }
1059    else if (r1->cf->ch==1) /* Q(a) */
1060    {
1061      if (r2->cf->ch==0) /* Q */
1062      {
1063        tmpR.cf->ch=r1->cf->ch;
1064        tmpR.cf->P=rPar(r1);
1065        tmpR.cf->parameter=(char **)omAlloc(rPar(r1)*sizeof(char *));
1066        int i;
1067        for(i=0;i<rPar(r1);i++)
1068        {
1069          tmpR.cf->parameter[i]=omStrDup(r1->cf->parameter[i]);
1070        }
1071        if (r1->cf->minpoly!=NULL)
1072        {
1073          tmpR.cf->minpoly=n_Copy(r1->cf->minpoly, r1->cf);
1074        }
1075      }
1076      else  /* R, Z/p,GF(p,n) */
1077      {
1078        WerrorS("Q(a)+(R,Z/p,GF(p,n))");
1079        return -1;
1080      }
1081    }
1082    else /* r1->cf->ch >=2 , Z/p */
1083    {
1084      if (r2->cf->ch==0) /* Q */
1085      {
1086        tmpR.cf->ch=r1->cf->ch;
1087      }
1088      else if (r2->cf->ch==-r1->cf->ch) /* Z/p(a) */
1089      {
1090        tmpR.cf->ch=r2->cf->ch;
1091        tmpR.cf->P=rPar(r2);
1092        tmpR.cf->parameter=(char **)omAlloc(rPar(r2)*sizeof(char *));
1093        int i;
1094        for(i=0;i<rPar(r2);i++)
1095        {
1096          tmpR.cf->parameter[i]=omStrDup(r2->cf->parameter[i]);
1097        }
1098        if (r2->cf->minpoly!=NULL)
1099        {
1100          tmpR.cf->minpoly=n_Copy(r2->cf->minpoly, r2->cf);
1101        }
1102      }
1103      else
1104      {
1105        WerrorS("Z/p+(GF(q,n),Z/q(a),R,Q(a))");
1106        return -1; /* GF(p,n),Z/q(a),R,Q(a) */
1107      }
1108    }
1109  }
1110  /* variable names ========================================================*/
1111  int i,j,k;
1112  int l=r1->N+r2->N;
1113  char **names=(char **)omAlloc0(l*sizeof(char *));
1114  k=0;
1115
1116  // collect all varnames from r1, except those which are parameters
1117  // of r2, or those which are the empty string
1118  for (i=0;i<r1->N;i++)
1119  {
1120    BOOLEAN b=TRUE;
1121
1122    if (*(r1->names[i]) == '\0')
1123      b = FALSE;
1124    else if ((r2->cf->parameter!=NULL) && (strlen(r1->names[i])==1))
1125    {
1126      if (vartest)
1127      {
1128        for(j=0;j<rPar(r2);j++)
1129        {
1130          if (strcmp(r1->names[i],r2->cf->parameter[j])==0)
1131          {
1132            b=FALSE;
1133            break;
1134          }
1135        }
1136      }
1137    }
1138
1139    if (b)
1140    {
1141      //Print("name : %d: %s\n",k,r1->names[i]);
1142      names[k]=omStrDup(r1->names[i]);
1143      k++;
1144    }
1145    //else
1146    //  Print("no name (par1) %s\n",r1->names[i]);
1147  }
1148  // Add variables from r2, except those which are parameters of r1
1149  // those which are empty strings, and those which equal a var of r1
1150  for(i=0;i<r2->N;i++)
1151  {
1152    BOOLEAN b=TRUE;
1153
1154    if (*(r2->names[i]) == '\0')
1155      b = FALSE;
1156    else if ((r1->cf->parameter!=NULL) && (strlen(r2->names[i])==1))
1157    {
1158      if (vartest)
1159      {
1160        for(j=0;j<rPar(r1);j++)
1161        {
1162          if (strcmp(r2->names[i],r1->cf->parameter[j])==0)
1163          {
1164            b=FALSE;
1165            break;
1166          }
1167        }
1168      }
1169    }
1170
1171    if (b)
1172    {
1173      if (vartest)
1174      {
1175        for(j=0;j<r1->N;j++)
1176        {
1177          if (strcmp(r1->names[j],r2->names[i])==0)
1178          {
1179            b=FALSE;
1180            break;
1181          }
1182        }
1183      }
1184      if (b)
1185      {
1186        //Print("name : %d : %s\n",k,r2->names[i]);
1187        names[k]=omStrDup(r2->names[i]);
1188        k++;
1189      }
1190      //else
1191      //  Print("no name (var): %s\n",r2->names[i]);
1192    }
1193    //else
1194    //  Print("no name (par): %s\n",r2->names[i]);
1195  }
1196  // check whether we found any vars at all
1197  if (k == 0)
1198  {
1199    names[k]=omStrDup("");
1200    k=1;
1201  }
1202  tmpR.N=k;
1203  tmpR.names=names;
1204  /* ordering *======================================================== */
1205  tmpR.OrdSgn=1;
1206  if (dp_dp
1207#ifdef HAVE_PLURAL
1208      && !rIsPluralRing(r1) && !rIsPluralRing(r2)
1209#endif
1210     )
1211  {
1212    tmpR.order=(int*)omAlloc(4*sizeof(int));
1213    tmpR.block0=(int*)omAlloc0(4*sizeof(int));
1214    tmpR.block1=(int*)omAlloc0(4*sizeof(int));
1215    tmpR.wvhdl=(int**)omAlloc0(4*sizeof(int *));
1216    tmpR.order[0]=ringorder_dp;
1217    tmpR.block0[0]=1;
1218    tmpR.block1[0]=rVar(r1);
1219    if (r2->OrdSgn==1)
1220    {
1221      if ((r2->block0[0]==1)
1222      && (r2->block1[0]==rVar(r2))
1223      && ((r2->order[0]==ringorder_wp)
1224        || (r2->order[0]==ringorder_Wp)
1225        || (r2->order[0]==ringorder_Dp))
1226     )
1227     {
1228       tmpR.order[1]=r2->order[0];
1229       if (r2->wvhdl[0]!=NULL)
1230         tmpR.wvhdl[1]=(int *)omMemDup(r2->wvhdl[0]);
1231     }
1232     else
1233        tmpR.order[1]=ringorder_dp;
1234    }
1235    else
1236    {
1237      tmpR.order[1]=ringorder_ds;
1238      tmpR.OrdSgn=-1;
1239    }
1240    tmpR.block0[1]=rVar(r1)+1;
1241    tmpR.block1[1]=rVar(r1)+rVar(r2);
1242    tmpR.order[2]=ringorder_C;
1243    tmpR.order[3]=0;
1244  }
1245  else
1246  {
1247    if ((r1->order[0]==ringorder_unspec)
1248        && (r2->order[0]==ringorder_unspec))
1249    {
1250      tmpR.order=(int*)omAlloc(3*sizeof(int));
1251      tmpR.block0=(int*)omAlloc(3*sizeof(int));
1252      tmpR.block1=(int*)omAlloc(3*sizeof(int));
1253      tmpR.wvhdl=(int**)omAlloc0(3*sizeof(int *));
1254      tmpR.order[0]=ringorder_unspec;
1255      tmpR.order[1]=ringorder_C;
1256      tmpR.order[2]=0;
1257      tmpR.block0[0]=1;
1258      tmpR.block1[0]=tmpR.N;
1259    }
1260    else if (l==k) /* r3=r1+r2 */
1261    {
1262      int b;
1263      ring rb;
1264      if (r1->order[0]==ringorder_unspec)
1265      {
1266        /* extend order of r2 to r3 */
1267        b=rBlocks(r2);
1268        rb=r2;
1269        tmpR.OrdSgn=r2->OrdSgn;
1270      }
1271      else if (r2->order[0]==ringorder_unspec)
1272      {
1273        /* extend order of r1 to r3 */
1274        b=rBlocks(r1);
1275        rb=r1;
1276        tmpR.OrdSgn=r1->OrdSgn;
1277      }
1278      else
1279      {
1280        b=rBlocks(r1)+rBlocks(r2)-2; /* for only one order C, only one 0 */
1281        rb=NULL;
1282      }
1283      tmpR.order=(int*)omAlloc0(b*sizeof(int));
1284      tmpR.block0=(int*)omAlloc0(b*sizeof(int));
1285      tmpR.block1=(int*)omAlloc0(b*sizeof(int));
1286      tmpR.wvhdl=(int**)omAlloc0(b*sizeof(int *));
1287      /* weights not implemented yet ...*/
1288      if (rb!=NULL)
1289      {
1290        for (i=0;i<b;i++)
1291        {
1292          tmpR.order[i]=rb->order[i];
1293          tmpR.block0[i]=rb->block0[i];
1294          tmpR.block1[i]=rb->block1[i];
1295          if (rb->wvhdl[i]!=NULL)
1296            WarnS("rSum: weights not implemented");
1297        }
1298        tmpR.block0[0]=1;
1299      }
1300      else /* ring sum for complete rings */
1301      {
1302        for (i=0;r1->order[i]!=0;i++)
1303        {
1304          tmpR.order[i]=r1->order[i];
1305          tmpR.block0[i]=r1->block0[i];
1306          tmpR.block1[i]=r1->block1[i];
1307          if (r1->wvhdl[i]!=NULL)
1308            tmpR.wvhdl[i] = (int*) omMemDup(r1->wvhdl[i]);
1309        }
1310        j=i;
1311        i--;
1312        if ((r1->order[i]==ringorder_c)
1313            ||(r1->order[i]==ringorder_C))
1314        {
1315          j--;
1316          tmpR.order[b-2]=r1->order[i];
1317        }
1318        for (i=0;r2->order[i]!=0;i++)
1319        {
1320          if ((r2->order[i]!=ringorder_c)
1321              &&(r2->order[i]!=ringorder_C))
1322          {
1323            tmpR.order[j]=r2->order[i];
1324            tmpR.block0[j]=r2->block0[i]+rVar(r1);
1325            tmpR.block1[j]=r2->block1[i]+rVar(r1);
1326            if (r2->wvhdl[i]!=NULL)
1327            {
1328              tmpR.wvhdl[j] = (int*) omMemDup(r2->wvhdl[i]);
1329            }
1330            j++;
1331          }
1332        }
1333        if((r1->OrdSgn==-1)||(r2->OrdSgn==-1))
1334          tmpR.OrdSgn=-1;
1335      }
1336    }
1337    else if ((k==rVar(r1)) && (k==rVar(r2))) /* r1 and r2 are "quite"
1338                                                the same ring */
1339      /* copy r1, because we have the variables from r1 */
1340    {
1341      int b=rBlocks(r1);
1342
1343      tmpR.order=(int*)omAlloc0(b*sizeof(int));
1344      tmpR.block0=(int*)omAlloc0(b*sizeof(int));
1345      tmpR.block1=(int*)omAlloc0(b*sizeof(int));
1346      tmpR.wvhdl=(int**)omAlloc0(b*sizeof(int *));
1347      /* weights not implemented yet ...*/
1348      for (i=0;i<b;i++)
1349      {
1350        tmpR.order[i]=r1->order[i];
1351        tmpR.block0[i]=r1->block0[i];
1352        tmpR.block1[i]=r1->block1[i];
1353        if (r1->wvhdl[i]!=NULL)
1354        {
1355          tmpR.wvhdl[i] = (int*) omMemDup(r1->wvhdl[i]);
1356        }
1357      }
1358      tmpR.OrdSgn=r1->OrdSgn;
1359    }
1360    else
1361    {
1362      for(i=0;i<k;i++) omFree((ADDRESS)tmpR.names[i]);
1363      omFreeSize((ADDRESS)names,tmpR.N*sizeof(char *));
1364      Werror("difficulties with variables: %d,%d -> %d",rVar(r1),rVar(r2),k);
1365      return -1;
1366    }
1367  }
1368  sum=(ring)omAllocBin(sip_sring_bin);
1369  memcpy(sum,&tmpR,sizeof(ip_sring));
1370  rComplete(sum);
1371
1372//#ifdef RDEBUG
1373//  rDebugPrint(sum);
1374//#endif
1375
1376
1377
1378#ifdef HAVE_PLURAL
1379  if(1)
1380  {
1381//    ring old_ring = currRing;
1382
1383    BOOLEAN R1_is_nc = rIsPluralRing(r1);
1384    BOOLEAN R2_is_nc = rIsPluralRing(r2);
1385
1386    if ( (R1_is_nc) || (R2_is_nc))
1387    {
1388      ring R1 = nc_rCreateNCcomm_rCopy(r1);
1389      assume( rIsPluralRing(R1) );
1390
1391#if 0
1392#ifdef RDEBUG
1393      rWrite(R1);
1394      rDebugPrint(R1);
1395#endif
1396#endif
1397      ring R2 = nc_rCreateNCcomm_rCopy(r2);
1398#if 0
1399#ifdef RDEBUG
1400      rWrite(R2);
1401      rDebugPrint(R2);
1402#endif
1403#endif
1404
1405//      rChangeCurrRing(sum); // ?
1406
1407      // Projections from R_i into Sum:
1408      /* multiplication matrices business: */
1409      /* find permutations of vars and pars */
1410      int *perm1 = (int *)omAlloc0((rVar(R1)+1)*sizeof(int));
1411      int *par_perm1 = NULL;
1412      if (rPar(R1)!=0) par_perm1=(int *)omAlloc0((rPar(R1)+1)*sizeof(int));
1413
1414      int *perm2 = (int *)omAlloc0((rVar(R2)+1)*sizeof(int));
1415      int *par_perm2 = NULL;
1416      if (rPar(R2)!=0) par_perm2=(int *)omAlloc0((rPar(R2)+1)*sizeof(int));
1417
1418      maFindPerm(R1->names,  rVar(R1),  R1->cf->parameter,  rPar(R1),
1419                 sum->names, rVar(sum), sum->cf->parameter, rPar(sum),
1420                 perm1, par_perm1, sum->cf->ch);
1421
1422      maFindPerm(R2->names,  rVar(R2),  R2->cf->parameter,  rPar(R2),
1423                 sum->names, rVar(sum), sum->cf->parameter, rPar(sum),
1424                 perm2, par_perm2, sum->cf->ch);
1425
1426
1427      matrix C1 = R1->GetNC()->C, C2 = R2->GetNC()->C;
1428      matrix D1 = R1->GetNC()->D, D2 = R2->GetNC()->D;
1429
1430      // !!!! BUG? C1 and C2 might live in different baserings!!!
1431
1432      int l = rVar(R1) + rVar(R2);
1433
1434      matrix C  = mpNew(l,l);
1435      matrix D  = mpNew(l,l);
1436
1437      for (i = 1; i <= rVar(R1); i++)
1438        for (j= rVar(R1)+1; j <= l; j++)
1439          MATELEM(C,i,j) = p_One(sum); // in 'sum'
1440
1441      id_Test((ideal)C, sum);
1442
1443      nMapFunc nMap1 = n_SetMap(R1->cf,sum->cf); /* can change something global: not usable
1444                                                    after the next nSetMap call :( */
1445      // Create blocked C and D matrices:
1446      for (i=1; i<= rVar(R1); i++)
1447        for (j=i+1; j<=rVar(R1); j++)
1448        {
1449          assume(MATELEM(C1,i,j) != NULL);
1450          MATELEM(C,i,j) = p_PermPoly(MATELEM(C1,i,j), perm1, R1, sum, nMap1, par_perm1, rPar(R1)); // need ADD + CMP ops.
1451
1452          if (MATELEM(D1,i,j) != NULL)
1453            MATELEM(D,i,j) = p_PermPoly(MATELEM(D1,i,j), perm1, R1, sum, nMap1, par_perm1, rPar(R1));
1454        }
1455
1456      id_Test((ideal)C, sum);
1457      id_Test((ideal)D, sum);
1458
1459
1460      nMapFunc nMap2 = n_SetMap(R2->cf,sum->cf); /* can change something global: not usable
1461                                                    after the next nSetMap call :( */
1462      for (i=1; i<= rVar(R2); i++)
1463        for (j=i+1; j<=rVar(R2); j++)
1464        {
1465          assume(MATELEM(C2,i,j) != NULL);
1466          MATELEM(C,rVar(R1)+i,rVar(R1)+j) = p_PermPoly(MATELEM(C2,i,j),perm2,R2,sum, nMap2,par_perm2,rPar(R2));
1467
1468          if (MATELEM(D2,i,j) != NULL)
1469            MATELEM(D,rVar(R1)+i,rVar(R1)+j) = p_PermPoly(MATELEM(D2,i,j),perm2,R2,sum, nMap2,par_perm2,rPar(R2));
1470        }
1471
1472      id_Test((ideal)C, sum);
1473      id_Test((ideal)D, sum);
1474
1475      // Now sum is non-commutative with blocked structure constants!
1476      if (nc_CallPlural(C, D, NULL, NULL, sum, false, false, true, sum))
1477        WarnS("Error initializing non-commutative multiplication!");
1478
1479      /* delete R1, R2*/
1480
1481#if 0
1482#ifdef RDEBUG
1483      rWrite(sum);
1484      rDebugPrint(sum);
1485
1486      Print("\nRefs: R1: %d, R2: %d\n", R1->GetNC()->ref, R2->GetNC()->ref);
1487
1488#endif
1489#endif
1490
1491
1492      rDelete(R1);
1493      rDelete(R2);
1494
1495      /* delete perm arrays */
1496      if (perm1!=NULL) omFree((ADDRESS)perm1);
1497      if (perm2!=NULL) omFree((ADDRESS)perm2);
1498      if (par_perm1!=NULL) omFree((ADDRESS)par_perm1);
1499      if (par_perm2!=NULL) omFree((ADDRESS)par_perm2);
1500
1501//      rChangeCurrRing(old_ring);
1502    }
1503
1504  }
1505#endif
1506
1507  ideal Q=NULL;
1508  ideal Q1=NULL, Q2=NULL;
1509  if (r1->qideal!=NULL)
1510  {
1511//    rChangeCurrRing(sum);
1512//     if (r2->qideal!=NULL)
1513//     {
1514//       WerrorS("todo: qring+qring");
1515//       return -1;
1516//     }
1517//     else
1518//     {}
1519    /* these were defined in the Plural Part above... */
1520    int *perm1 = (int *)omAlloc0((rVar(r1)+1)*sizeof(int));
1521    int *par_perm1 = NULL;
1522    if (rPar(r1)!=0) par_perm1=(int *)omAlloc0((rPar(r1)+1)*sizeof(int));
1523    maFindPerm(r1->names,  rVar(r1),  r1->cf->parameter,  rPar(r1),
1524               sum->names, rVar(sum), sum->cf->parameter, rPar(sum),
1525               perm1, par_perm1, sum->cf->ch);
1526    nMapFunc nMap1 = n_SetMap(r1->cf,sum->cf);
1527    Q1 = idInit(IDELEMS(r1->qideal),1);
1528
1529    for (int for_i=0;for_i<IDELEMS(r1->qideal);for_i++)
1530      Q1->m[for_i] = p_PermPoly(
1531                                r1->qideal->m[for_i], perm1,
1532                                r1, sum,
1533                                nMap1, 
1534                                par_perm1, rPar(r1));
1535   
1536    omFree((ADDRESS)perm1);
1537  }
1538
1539  if (r2->qideal!=NULL)
1540  {
1541    //if (currRing!=sum)
1542    //  rChangeCurrRing(sum);
1543    int *perm2 = (int *)omAlloc0((rVar(r2)+1)*sizeof(int));
1544    int *par_perm2 = NULL;
1545    if (rPar(r2)!=0) par_perm2=(int *)omAlloc0((rPar(r2)+1)*sizeof(int));
1546    maFindPerm(r2->names,  rVar(r2),  r2->cf->parameter,  rPar(r2),
1547               sum->names, rVar(sum), sum->cf->parameter, rPar(sum),
1548               perm2, par_perm2, sum->cf->ch);
1549    nMapFunc nMap2 = n_SetMap(r2->cf,sum->cf);
1550    Q2 = idInit(IDELEMS(r2->qideal),1);
1551     
1552    for (int for_i=0;for_i<IDELEMS(r2->qideal);for_i++)
1553      Q2->m[for_i] = p_PermPoly( 
1554                  r2->qideal->m[for_i], perm2, 
1555                  r2, sum, 
1556                  nMap2, 
1557                  par_perm2, rPar(r2));
1558   
1559    omFree((ADDRESS)perm2);
1560  }
1561  if ( (Q1!=NULL) || ( Q2!=NULL))
1562  {
1563    Q = id_SimpleAdd(Q1,Q2,sum);
1564  }
1565  sum->qideal = Q;
1566
1567#ifdef HAVE_PLURAL
1568  if( rIsPluralRing(sum) )
1569    nc_SetupQuotient( sum );
1570#endif
1571  return 1;
1572}
1573
1574/*2
1575 *returns -1 for not compatible, (sum is undefined)
1576 *         0 for equal, (and sum)
1577 *         1 for compatible (and sum)
1578 */
1579int rSum(ring r1, ring r2, ring &sum)
1580{
1581  if (r1==r2)
1582  {
1583    sum=r1;
1584    r1->ref++;
1585    return 0;
1586  }
1587  return rSumInternal(r1,r2,sum,TRUE,FALSE);
1588}
1589
1590/*2
1591 * create a copy of the ring r, which must be equivalent to currRing
1592 * used for qring definition,..
1593 * (i.e.: normal rings: same nCopy as currRing;
1594 *        qring:        same nCopy, same idCopy as currRing)
1595 * DOES NOT CALL rComplete
1596 */
1597ring rCopy0(const ring r, BOOLEAN copy_qideal, BOOLEAN copy_ordering)
1598{
1599  TODO(Hans, "Please consider the change to '->cf'...");
1600 
1601  if (r == NULL) return NULL;
1602  int i,j;
1603  ring res=(ring)omAllocBin(sip_sring_bin);
1604  memset(res,0,sizeof(ip_sring));
1605  //memcpy(res,r,sizeof(ip_sring));
1606  //memset: res->idroot=NULL; /* local objects */
1607  //ideal      minideal;
1608  res->options=r->options; /* ring dependent options */
1609
1610  //memset: res->ordsgn=NULL;
1611  //memset: res->typ=NULL;
1612  //memset: res->VarOffset=NULL;
1613  //memset: res->firstwv=NULL;
1614
1615  //struct omBin   PolyBin; /* Bin from where monoms are allocated */
1616  //memset: res->PolyBin=NULL; // rComplete
1617  res->cf->ch=r->cf->ch;     /* characteristic */
1618#ifdef HAVE_RINGS
1619  res->cf->ringtype = r->cf->ringtype;  // cring = 0 => coefficient field, cring = 1 => coeffs from Z/2^m
1620  if (r->cf->modBase!=NULL)
1621  {
1622    res->cf->modBase = (int_number) omAlloc(sizeof(mpz_t));
1623    mpz_init_set(res->cf->modBase,r->cf->modBase);
1624  }
1625  res->cf->modExponent=r->cf->modExponent;
1626  if (r->cf->modNumber!=NULL)
1627  {
1628    res->cf->modNumber = (int_number) omAlloc(sizeof(mpz_t));
1629    mpz_init_set(res->cf->modNumber,r->cf->modNumber);
1630  }
1631#endif
1632  //memset: res->ref=0; /* reference counter to the ring */
1633
1634  res->float_len=r->float_len; /* additional char-flags */
1635  res->float_len2=r->float_len2; /* additional char-flags */
1636
1637  res->N=rVar(r);      /* number of vars */
1638  res->cf->P=rPar(r);      /* number of pars */
1639  res->OrdSgn=r->OrdSgn; /* 1 for polynomial rings, -1 otherwise */
1640
1641  res->firstBlockEnds=r->firstBlockEnds;
1642#ifdef HAVE_PLURAL
1643  res->real_var_start=r->real_var_start;
1644  res->real_var_end=r->real_var_end;
1645#endif
1646
1647#ifdef HAVE_SHIFTBBA
1648  res->isLPring=r->isLPring; /* 0 for non-letterplace rings, otherwise the number of LP blocks, at least 1, known also as lV */
1649#endif
1650
1651  res->VectorOut=r->VectorOut;
1652  res->ShortOut=r->ShortOut;
1653  res->CanShortOut=r->CanShortOut;
1654  res->LexOrder=r->LexOrder; // TRUE if the monomial ordering has polynomial and power series blocks
1655  res->MixedOrder=r->MixedOrder; // ?? 1 for lex ordering (except ls), -1 otherwise
1656  res->ComponentOrder=r->ComponentOrder;
1657
1658  //memset: res->ExpL_Size=0;
1659  //memset: res->CmpL_Size=0;
1660  //memset: res->VarL_Size=0;
1661  //memset: res->pCompIndex=0;
1662  //memset: res->pOrdIndex=0;
1663  //memset: res->OrdSize=0;
1664  //memset: res->VarL_LowIndex=0;
1665  //memset: res->MinExpPerLong=0;
1666  //memset: res->NegWeightL_Size=0;
1667  //memset: res->NegWeightL_Offset=NULL;
1668  //memset: res->VarL_Offset=NULL;
1669
1670  // the following are set by rComplete unless predefined
1671  // therefore, we copy these values: maybe they are non-standard
1672  /* mask for getting single exponents */
1673  res->bitmask=r->bitmask;
1674  res->divmask=r->divmask;
1675  res->BitsPerExp = r->BitsPerExp;
1676  res->ExpPerLong =  r->ExpPerLong;
1677
1678  //memset: res->p_Procs=NULL;
1679  //memset: res->pFDeg=NULL;
1680  //memset: res->pLDeg=NULL;
1681  //memset: res->pFDegOrig=NULL;
1682  //memset: res->pLDegOrig=NULL;
1683  //memset: res->p_Setm=NULL;
1684  //memset: res->cf=NULL;
1685  res->options=r->options;
1686#ifdef HAVE_RINGS
1687  res->cf->ringtype=r->cf->ringtype;
1688#endif
1689  //
1690  if (r->algring!=NULL)
1691    r->algring->ref++;
1692  res->algring=r->algring;
1693  //memset: res->minideal=NULL;
1694  if (r->cf->parameter!=NULL)
1695  {
1696    if (r->cf->minpoly!=NULL) res->cf->minpoly=n_Copy(r->cf->minpoly,r->cf);
1697    int l=rPar(r);
1698    res->cf->parameter=(char **)omAlloc(l*sizeof(char *));
1699    int i;
1700    for(i=0;i<rPar(r);i++)
1701    {
1702      res->cf->parameter[i]=omStrDup(r->cf->parameter[i]);
1703    }
1704    if (r->minideal!=NULL)
1705    {
1706      res->minideal=id_Copy(r->minideal,r->algring);
1707    }
1708  }
1709  if (copy_ordering == TRUE)
1710  {
1711    i=rBlocks(r);
1712    res->wvhdl   = (int **)omAlloc(i * sizeof(int *));
1713    res->order   = (int *) omAlloc(i * sizeof(int));
1714    res->block0  = (int *) omAlloc(i * sizeof(int));
1715    res->block1  = (int *) omAlloc(i * sizeof(int));
1716    for (j=0; j<i; j++)
1717    {
1718      if (r->wvhdl[j]!=NULL)
1719      {
1720        res->wvhdl[j] = (int*) omMemDup(r->wvhdl[j]);
1721      }
1722      else
1723        res->wvhdl[j]=NULL;
1724    }
1725    memcpy(res->order,r->order,i * sizeof(int));
1726    memcpy(res->block0,r->block0,i * sizeof(int));
1727    memcpy(res->block1,r->block1,i * sizeof(int));
1728  }
1729  //memset: else
1730  //memset: {
1731  //memset:   res->wvhdl = NULL;
1732  //memset:   res->order = NULL;
1733  //memset:   res->block0 = NULL;
1734  //memset:   res->block1 = NULL;
1735  //memset: }
1736
1737  res->names   = (char **)omAlloc0(rVar(r) * sizeof(char *));
1738  for (i=0; i<rVar(res); i++)
1739  {
1740    res->names[i] = omStrDup(r->names[i]);
1741  }
1742  if (r->qideal!=NULL)
1743  {
1744    if (copy_qideal)
1745    {
1746      #ifndef NDEBUG
1747      if (!copy_ordering)
1748        WerrorS("internal error: rCopy0(Q,TRUE,FALSE)");
1749      else
1750      #endif
1751      {
1752      #ifndef NDEBUG
1753        WarnS("internal bad stuff: rCopy0(Q,TRUE,TRUE)");
1754      #endif
1755        rComplete(res);
1756        res->qideal= idrCopyR_NoSort(r->qideal, r, res);
1757        rUnComplete(res);
1758      }
1759    }
1760    //memset: else res->qideal = NULL;
1761  }
1762  //memset: else res->qideal = NULL;
1763  //memset: res->GetNC() = NULL; // copy is purely commutative!!!
1764  return res;
1765}
1766
1767/*2
1768 * create a copy of the ring r, which must be equivalent to currRing
1769 * used for qring definition,..
1770 * (i.e.: normal rings: same nCopy as currRing;
1771 *        qring:        same nCopy, same idCopy as currRing)
1772 */
1773ring rCopy(ring r)
1774{
1775  if (r == NULL) return NULL;
1776  ring res=rCopy0(r,FALSE,TRUE);
1777  rComplete(res, 1); // res is purely commutative so far
1778  if (r->qideal!=NULL) res->qideal=idrCopyR_NoSort(r->qideal, r, res);
1779
1780#ifdef HAVE_PLURAL
1781  if (rIsPluralRing(r))
1782    if( nc_rCopy(res, r, true) );
1783#endif
1784
1785  return res;
1786}
1787
1788// returns TRUE, if r1 equals r2 FALSE, otherwise Equality is
1789// determined componentwise, if qr == 1, then qrideal equality is
1790// tested, as well
1791BOOLEAN rEqual(ring r1, ring r2, BOOLEAN qr)
1792{
1793  int i, j;
1794
1795  if (r1 == r2) return TRUE;
1796
1797  if (r1 == NULL || r2 == NULL) return FALSE;
1798
1799  if ((r1->cf->type != r2->cf->type)
1800  || (r1->float_len != r2->float_len)
1801  || (r1->float_len2 != r2->float_len2)
1802  || (rVar(r1) != rVar(r2))
1803  || (r1->OrdSgn != r2->OrdSgn)
1804  || (rPar(r1) != rPar(r2)))
1805    return FALSE;
1806
1807  for (i=0; i<rVar(r1); i++)
1808  {
1809    if (r1->names[i] != NULL && r2->names[i] != NULL)
1810    {
1811      if (strcmp(r1->names[i], r2->names[i])) return FALSE;
1812    }
1813    else if ((r1->names[i] != NULL) ^ (r2->names[i] != NULL))
1814    {
1815      return FALSE;
1816    }
1817  }
1818
1819  i=0;
1820  while (r1->order[i] != 0)
1821  {
1822    if (r2->order[i] == 0) return FALSE;
1823    if ((r1->order[i] != r2->order[i])
1824    || (r1->block0[i] != r2->block0[i])
1825    || (r1->block1[i] != r2->block1[i]))
1826      return FALSE;
1827    if (r1->wvhdl[i] != NULL)
1828    {
1829      if (r2->wvhdl[i] == NULL)
1830        return FALSE;
1831      for (j=0; j<r1->block1[i]-r1->block0[i]+1; j++)
1832        if (r2->wvhdl[i][j] != r1->wvhdl[i][j])
1833          return FALSE;
1834    }
1835    else if (r2->wvhdl[i] != NULL) return FALSE;
1836    i++;
1837  }
1838  if (r2->order[i] != 0) return FALSE;
1839
1840  for (i=0; i<rPar(r1);i++)
1841  {
1842      if (strcmp(r1->cf->parameter[i], r2->cf->parameter[i])!=0)
1843        return FALSE;
1844  }
1845
1846  if (r1->cf->minpoly != NULL)
1847  {
1848    if (r2->cf->minpoly == NULL) return FALSE;
1849    if (! n_Equal(r1->cf->minpoly, r2->cf->minpoly, r1->cf)) return FALSE;
1850  }
1851  else if (r2->cf->minpoly != NULL) return FALSE;
1852
1853  if (qr)
1854  {
1855    if (r1->qideal != NULL)
1856    {
1857      ideal id1 = r1->qideal, id2 = r2->qideal;
1858      int i, n;
1859      poly *m1, *m2;
1860
1861      if (id2 == NULL) return FALSE;
1862      if ((n = IDELEMS(id1)) != IDELEMS(id2)) return FALSE;
1863
1864      {
1865        m1 = id1->m;
1866        m2 = id2->m;
1867        for (i=0; i<n; i++)
1868          if (! p_EqualPolys(m1[i],m2[i],r1)) return FALSE;
1869      }
1870    }
1871    else if (r2->qideal != NULL) return FALSE;
1872  }
1873
1874  return TRUE;
1875}
1876
1877// returns TRUE, if r1 and r2 represents the monomials in the same way
1878// FALSE, otherwise
1879// this is an analogue to rEqual but not so strict
1880BOOLEAN rSamePolyRep(ring r1, ring r2)
1881{
1882  int i, j;
1883
1884  if (r1 == r2) return TRUE;
1885
1886  if (r1 == NULL || r2 == NULL) return FALSE;
1887
1888  if ((r1->cf->type != r2->cf->type)
1889  || (r1->float_len != r2->float_len)
1890  || (r1->float_len2 != r2->float_len2)
1891  || (rVar(r1) != rVar(r2))
1892  || (r1->OrdSgn != r2->OrdSgn)
1893  || (rPar(r1) != rPar(r2)))
1894    return FALSE;
1895
1896  if (rVar(r1)!=rVar(r2)) return FALSE;
1897  if (rPar(r1)!=rPar(r2)) return FALSE;
1898
1899  i=0;
1900  while (r1->order[i] != 0)
1901  {
1902    if (r2->order[i] == 0) return FALSE;
1903    if ((r1->order[i] != r2->order[i])
1904    || (r1->block0[i] != r2->block0[i])
1905    || (r1->block1[i] != r2->block1[i]))
1906      return FALSE;
1907    if (r1->wvhdl[i] != NULL)
1908    {
1909      if (r2->wvhdl[i] == NULL)
1910        return FALSE;
1911      for (j=0; j<r1->block1[i]-r1->block0[i]+1; j++)
1912        if (r2->wvhdl[i][j] != r1->wvhdl[i][j])
1913          return FALSE;
1914    }
1915    else if (r2->wvhdl[i] != NULL) return FALSE;
1916    i++;
1917  }
1918  if (r2->order[i] != 0) return FALSE;
1919
1920  // we do not check minpoly
1921  // we do not check qideal
1922
1923  return TRUE;
1924}
1925
1926rOrderType_t rGetOrderType(ring r)
1927{
1928  // check for simple ordering
1929  if (rHasSimpleOrder(r))
1930  {
1931    if ((r->order[1] == ringorder_c)
1932    || (r->order[1] == ringorder_C))
1933    {
1934      switch(r->order[0])
1935      {
1936          case ringorder_dp:
1937          case ringorder_wp:
1938          case ringorder_ds:
1939          case ringorder_ws:
1940          case ringorder_ls:
1941          case ringorder_unspec:
1942            if (r->order[1] == ringorder_C
1943            ||  r->order[0] == ringorder_unspec)
1944              return rOrderType_ExpComp;
1945            return rOrderType_Exp;
1946
1947          default:
1948            assume(r->order[0] == ringorder_lp ||
1949                   r->order[0] == ringorder_rs ||
1950                   r->order[0] == ringorder_Dp ||
1951                   r->order[0] == ringorder_Wp ||
1952                   r->order[0] == ringorder_Ds ||
1953                   r->order[0] == ringorder_Ws);
1954
1955            if (r->order[1] == ringorder_c) return rOrderType_ExpComp;
1956            return rOrderType_Exp;
1957      }
1958    }
1959    else
1960    {
1961      assume((r->order[0]==ringorder_c)||(r->order[0]==ringorder_C));
1962      return rOrderType_CompExp;
1963    }
1964  }
1965  else
1966    return rOrderType_General;
1967}
1968
1969BOOLEAN rHasSimpleOrder(const ring r)
1970{
1971  if (r->order[0] == ringorder_unspec) return TRUE;
1972  int blocks = rBlocks(r) - 1;
1973  assume(blocks >= 1);
1974  if (blocks == 1) return TRUE;
1975
1976  int s = 0; 
1977  while( (s < blocks) && (r->order[s] == ringorder_IS) && (r->order[blocks-1] == ringorder_IS) )
1978  {
1979    s++;
1980    blocks--;
1981  }
1982
1983  if ((blocks - s) > 2)  return FALSE;
1984
1985  assume( blocks == s + 2 );
1986 
1987  if (
1988     (r->order[s] != ringorder_c)
1989  && (r->order[s] != ringorder_C)
1990  && (r->order[s+1] != ringorder_c)
1991  && (r->order[s+1] != ringorder_C)
1992     )
1993    return FALSE;
1994  if ((r->order[s+1] == ringorder_M)
1995  || (r->order[s] == ringorder_M))
1996    return FALSE;
1997  return TRUE;
1998}
1999
2000// returns TRUE, if simple lp or ls ordering
2001BOOLEAN rHasSimpleLexOrder(const ring r)
2002{
2003  return rHasSimpleOrder(r) &&
2004    (r->order[0] == ringorder_ls ||
2005     r->order[0] == ringorder_lp ||
2006     r->order[1] == ringorder_ls ||
2007     r->order[1] == ringorder_lp);
2008}
2009
2010BOOLEAN rOrder_is_DegOrdering(const rRingOrder_t order)
2011{
2012  switch(order)
2013  {
2014      case ringorder_dp:
2015      case ringorder_Dp:
2016      case ringorder_ds:
2017      case ringorder_Ds:
2018      case ringorder_Ws:
2019      case ringorder_Wp:
2020      case ringorder_ws:
2021      case ringorder_wp:
2022        return TRUE;
2023
2024      default:
2025        return FALSE;
2026  }
2027}
2028
2029BOOLEAN rOrder_is_WeightedOrdering(rRingOrder_t order)
2030{
2031  switch(order)
2032  {
2033      case ringorder_Ws:
2034      case ringorder_Wp:
2035      case ringorder_ws:
2036      case ringorder_wp:
2037        return TRUE;
2038
2039      default:
2040        return FALSE;
2041  }
2042}
2043
2044BOOLEAN rHasSimpleOrderAA(ring r)
2045{
2046  if (r->order[0] == ringorder_unspec) return TRUE;
2047  int blocks = rBlocks(r) - 1;
2048  assume(blocks >= 1);
2049  if (blocks == 1) return TRUE;
2050
2051  int s = 0; 
2052  while( (s < blocks) && (r->order[s] == ringorder_IS) && (r->order[blocks-1] == ringorder_IS) )
2053  {
2054    s++;
2055    blocks--;
2056  }
2057
2058  if ((blocks - s) > 3)  return FALSE;
2059 
2060//  if ((blocks > 3) || (blocks < 2)) return FALSE;
2061  if ((blocks - s) == 3)
2062  {
2063    return (((r->order[s] == ringorder_aa) && (r->order[s+1] != ringorder_M) &&
2064             ((r->order[s+2] == ringorder_c) || (r->order[s+2] == ringorder_C))) ||
2065            (((r->order[s] == ringorder_c) || (r->order[s] == ringorder_C)) &&
2066             (r->order[s+1] == ringorder_aa) && (r->order[s+2] != ringorder_M)));
2067  }
2068  else
2069  {
2070    return ((r->order[s] == ringorder_aa) && (r->order[s+1] != ringorder_M));
2071  }
2072}
2073
2074// return TRUE if p_SetComp requires p_Setm
2075BOOLEAN rOrd_SetCompRequiresSetm(ring r)
2076{
2077  if (r->typ != NULL)
2078  {
2079    int pos;
2080    for (pos=0;pos<r->OrdSize;pos++)
2081    {
2082      sro_ord* o=&(r->typ[pos]);
2083      if ((o->ord_typ == ro_syzcomp) || (o->ord_typ == ro_syz) || (o->ord_typ == ro_is) || (o->ord_typ == ro_isTemp)) return TRUE;
2084    }
2085  }
2086  return FALSE;
2087}
2088
2089// return TRUE if p->exp[r->pOrdIndex] holds total degree of p */
2090BOOLEAN rOrd_is_Totaldegree_Ordering(ring r)
2091{
2092  // Hmm.... what about Syz orderings?
2093  return (rVar(r) > 1 &&
2094          ((rHasSimpleOrder(r) &&
2095           (rOrder_is_DegOrdering((rRingOrder_t)r->order[0]) ||
2096            rOrder_is_DegOrdering(( rRingOrder_t)r->order[1]))) ||
2097           (rHasSimpleOrderAA(r) &&
2098            (rOrder_is_DegOrdering((rRingOrder_t)r->order[1]) ||
2099             rOrder_is_DegOrdering((rRingOrder_t)r->order[2])))));
2100}
2101
2102// return TRUE if p->exp[r->pOrdIndex] holds a weighted degree of p */
2103BOOLEAN rOrd_is_WeightedDegree_Ordering(ring r )
2104{
2105  // Hmm.... what about Syz orderings?
2106  return ((rVar(r) > 1) &&
2107          rHasSimpleOrder(r) &&
2108          (rOrder_is_WeightedOrdering((rRingOrder_t)r->order[0]) ||
2109           rOrder_is_WeightedOrdering(( rRingOrder_t)r->order[1])));
2110}
2111
2112BOOLEAN rIsPolyVar(int v, ring r)
2113{
2114  int  i=0;
2115  while(r->order[i]!=0)
2116  {
2117    if((r->block0[i]<=v)
2118    && (r->block1[i]>=v))
2119    {
2120      switch(r->order[i])
2121      {
2122        case ringorder_a:
2123          return (r->wvhdl[i][v-r->block0[i]]>0);
2124        case ringorder_M:
2125          return 2; /*don't know*/
2126        case ringorder_a64: /* assume: all weight are non-negative!*/
2127        case ringorder_lp:
2128        case ringorder_rs:
2129        case ringorder_dp:
2130        case ringorder_Dp:
2131        case ringorder_wp:
2132        case ringorder_Wp:
2133          return TRUE;
2134        case ringorder_ls:
2135        case ringorder_ds:
2136        case ringorder_Ds:
2137        case ringorder_ws:
2138        case ringorder_Ws:
2139          return FALSE;
2140        default:
2141          break;
2142      }
2143    }
2144    i++;
2145  }
2146  return 3; /* could not find var v*/
2147}
2148
2149#ifdef RDEBUG
2150// This should eventually become a full-fledge ring check, like pTest
2151BOOLEAN rDBTest(ring r, const char* fn, const int l)
2152{
2153  int i,j;
2154
2155  if (r == NULL)
2156  {
2157    dReportError("Null ring in %s:%d", fn, l);
2158    return FALSE;
2159  }
2160
2161
2162  if (r->N == 0) return TRUE;
2163
2164//  omCheckAddrSize(r,sizeof(ip_sring));
2165#if OM_CHECK > 0
2166  i=rBlocks(r);
2167  omCheckAddrSize(r->order,i*sizeof(int));
2168  omCheckAddrSize(r->block0,i*sizeof(int));
2169  omCheckAddrSize(r->block1,i*sizeof(int));
2170  if (r->wvhdl!=NULL)
2171  {
2172    omCheckAddrSize(r->wvhdl,i*sizeof(int *));
2173    for (j=0;j<i; j++)
2174    {
2175      if (r->wvhdl[j] != NULL) omCheckAddr(r->wvhdl[j]);
2176    }
2177  }
2178#endif
2179  if (r->VarOffset == NULL)
2180  {
2181    dReportError("Null ring VarOffset -- no rComplete (?) in n %s:%d", fn, l);
2182    return FALSE;
2183  }
2184  omCheckAddrSize(r->VarOffset,(r->N+1)*sizeof(int));
2185
2186  if ((r->OrdSize==0)!=(r->typ==NULL))
2187  {
2188    dReportError("mismatch OrdSize and typ-pointer in %s:%d");
2189    return FALSE;
2190  }
2191  omcheckAddrSize(r->typ,r->OrdSize*sizeof(*(r->typ)));
2192  omCheckAddrSize(r->VarOffset,(r->N+1)*sizeof(*(r->VarOffset)));
2193  // test assumptions:
2194  for(i=0;i<=r->N;i++) // for all variables (i = 0..N)
2195  {
2196    if(r->typ!=NULL)
2197    {
2198      for(j=0;j<r->OrdSize;j++) // for all ordering blocks (j =0..OrdSize-1)
2199      {
2200        if(r->typ[j].ord_typ == ro_isTemp)
2201        {
2202          const int p = r->typ[j].data.isTemp.suffixpos;
2203
2204          if(p <= j)
2205            dReportError("ordrec prefix %d is unmatched",j);
2206
2207          assume( p < r->OrdSize );
2208
2209          if(r->typ[p].ord_typ != ro_is)
2210            dReportError("ordrec prefix %d is unmatched (suffix: %d is wrong!!!)",j, p);
2211
2212          // Skip all intermediate blocks for undone variables:
2213          if(r->typ[j].data.isTemp.pVarOffset[i] != -1) // Check i^th variable
2214          {
2215            j = p - 1; // SKIP ALL INTERNAL BLOCKS...???
2216            continue; // To make for check OrdSize bound...
2217          }
2218        }
2219        else if (r->typ[j].ord_typ == ro_is)
2220        {
2221          // Skip all intermediate blocks for undone variables:
2222          if(r->typ[j].data.is.pVarOffset[i] != -1)
2223          {
2224            // ???
2225          }
2226
2227        }
2228        else
2229        {
2230          if (r->typ[j].ord_typ==ro_cp)
2231          {
2232            if(((short)r->VarOffset[i]) == r->typ[j].data.cp.place)
2233              dReportError("ordrec %d conflicts with var %d",j,i);
2234          }
2235          else
2236            if ((r->typ[j].ord_typ!=ro_syzcomp)
2237            && (r->VarOffset[i] == r->typ[j].data.dp.place))
2238              dReportError("ordrec %d conflicts with var %d",j,i);
2239        }
2240      }
2241    }
2242    int tmp;
2243      tmp=r->VarOffset[i] & 0xffffff;
2244      #if SIZEOF_LONG == 8
2245        if ((r->VarOffset[i] >> 24) >63)
2246      #else
2247        if ((r->VarOffset[i] >> 24) >31)
2248      #endif
2249          dReportError("bit_start out of range:%d",r->VarOffset[i] >> 24);
2250      if (i > 0 && ((tmp<0) ||(tmp>r->ExpL_Size-1)))
2251      {
2252        dReportError("varoffset out of range for var %d: %d",i,tmp);
2253      }
2254  }
2255  if(r->typ!=NULL)
2256  {
2257    for(j=0;j<r->OrdSize;j++)
2258    {
2259      if ((r->typ[j].ord_typ==ro_dp)
2260      || (r->typ[j].ord_typ==ro_wp)
2261      || (r->typ[j].ord_typ==ro_wp_neg))
2262      {
2263        if (r->typ[j].data.dp.start > r->typ[j].data.dp.end)
2264          dReportError("in ordrec %d: start(%d) > end(%d)",j,
2265            r->typ[j].data.dp.start, r->typ[j].data.dp.end);
2266        if ((r->typ[j].data.dp.start < 1)
2267        || (r->typ[j].data.dp.end > r->N))
2268          dReportError("in ordrec %d: start(%d)<1 or end(%d)>vars(%d)",j,
2269            r->typ[j].data.dp.start, r->typ[j].data.dp.end,r->N);
2270      }
2271    }
2272  }
2273  if (r->cf->minpoly!=NULL)
2274  {
2275    omCheckAddr(r->cf->minpoly);
2276  }
2277  //assume(r->cf!=NULL);
2278
2279  return TRUE;
2280}
2281#endif
2282
2283static void rO_Align(int &place, int &bitplace)
2284{
2285  // increment place to the next aligned one
2286  // (count as Exponent_t,align as longs)
2287  if (bitplace!=BITS_PER_LONG)
2288  {
2289    place++;
2290    bitplace=BITS_PER_LONG;
2291  }
2292}
2293
2294static void rO_TDegree(int &place, int &bitplace, int start, int end,
2295    long *o, sro_ord &ord_struct)
2296{
2297  // degree (aligned) of variables v_start..v_end, ordsgn 1
2298  rO_Align(place,bitplace);
2299  ord_struct.ord_typ=ro_dp;
2300  ord_struct.data.dp.start=start;
2301  ord_struct.data.dp.end=end;
2302  ord_struct.data.dp.place=place;
2303  o[place]=1;
2304  place++;
2305  rO_Align(place,bitplace);
2306}
2307
2308static void rO_TDegree_neg(int &place, int &bitplace, int start, int end,
2309    long *o, sro_ord &ord_struct)
2310{
2311  // degree (aligned) of variables v_start..v_end, ordsgn -1
2312  rO_Align(place,bitplace);
2313  ord_struct.ord_typ=ro_dp;
2314  ord_struct.data.dp.start=start;
2315  ord_struct.data.dp.end=end;
2316  ord_struct.data.dp.place=place;
2317  o[place]=-1;
2318  place++;
2319  rO_Align(place,bitplace);
2320}
2321
2322static void rO_WDegree(int &place, int &bitplace, int start, int end,
2323    long *o, sro_ord &ord_struct, int *weights)
2324{
2325  // weighted degree (aligned) of variables v_start..v_end, ordsgn 1
2326  while((start<end) && (weights[0]==0)) { start++; weights++; }
2327  while((start<end) && (weights[end-start]==0)) { end--; }
2328  int i;
2329  int pure_tdeg=1;
2330  for(i=start;i<=end;i++)
2331  {
2332    if(weights[i-start]!=1)
2333    {
2334      pure_tdeg=0;
2335      break;
2336    }
2337  }
2338  if (pure_tdeg)
2339  {
2340    rO_TDegree(place,bitplace,start,end,o,ord_struct);
2341    return;
2342  }
2343  rO_Align(place,bitplace);
2344  ord_struct.ord_typ=ro_wp;
2345  ord_struct.data.wp.start=start;
2346  ord_struct.data.wp.end=end;
2347  ord_struct.data.wp.place=place;
2348  ord_struct.data.wp.weights=weights;
2349  o[place]=1;
2350  place++;
2351  rO_Align(place,bitplace);
2352  for(i=start;i<=end;i++)
2353  {
2354    if(weights[i-start]<0)
2355    {
2356      ord_struct.ord_typ=ro_wp_neg;
2357      break;
2358    }
2359  }
2360}
2361
2362static void rO_WDegree64(int &place, int &bitplace, int start, int end,
2363    long *o, sro_ord &ord_struct, int64 *weights)
2364{
2365  // weighted degree (aligned) of variables v_start..v_end, ordsgn 1,
2366  // reserved 2 places
2367  rO_Align(place,bitplace);
2368  ord_struct.ord_typ=ro_wp64;
2369  ord_struct.data.wp64.start=start;
2370  ord_struct.data.wp64.end=end;
2371  ord_struct.data.wp64.place=place;
2372  ord_struct.data.wp64.weights64=weights;
2373  o[place]=1;
2374  place++;
2375  o[place]=1;
2376  place++;
2377  rO_Align(place,bitplace);
2378}
2379
2380static void rO_WDegree_neg(int &place, int &bitplace, int start, int end,
2381    long *o, sro_ord &ord_struct, int *weights)
2382{
2383  // weighted degree (aligned) of variables v_start..v_end, ordsgn -1
2384  while((start<end) && (weights[0]==0)) { start++; weights++; }
2385  while((start<end) && (weights[end-start]==0)) { end--; }
2386  rO_Align(place,bitplace);
2387  ord_struct.ord_typ=ro_wp;
2388  ord_struct.data.wp.start=start;
2389  ord_struct.data.wp.end=end;
2390  ord_struct.data.wp.place=place;
2391  ord_struct.data.wp.weights=weights;
2392  o[place]=-1;
2393  place++;
2394  rO_Align(place,bitplace);
2395  int i;
2396  for(i=start;i<=end;i++)
2397  {
2398    if(weights[i-start]<0)
2399    {
2400      ord_struct.ord_typ=ro_wp_neg;
2401      break;
2402    }
2403  }
2404}
2405
2406static void rO_LexVars(int &place, int &bitplace, int start, int end,
2407  int &prev_ord, long *o,int *v, int bits, int opt_var)
2408{
2409  // a block of variables v_start..v_end with lex order, ordsgn 1
2410  int k;
2411  int incr=1;
2412  if(prev_ord==-1) rO_Align(place,bitplace);
2413
2414  if (start>end)
2415  {
2416    incr=-1;
2417  }
2418  for(k=start;;k+=incr)
2419  {
2420    bitplace-=bits;
2421    if (bitplace < 0) { bitplace=BITS_PER_LONG-bits; place++; }
2422    o[place]=1;
2423    v[k]= place | (bitplace << 24);
2424    if (k==end) break;
2425  }
2426  prev_ord=1;
2427  if (opt_var!= -1)
2428  {
2429    assume((opt_var == end+1) ||(opt_var == end-1));
2430    if((opt_var != end+1) &&(opt_var != end-1)) WarnS("hier-2");
2431    int save_bitplace=bitplace;
2432    bitplace-=bits;
2433    if (bitplace < 0)
2434    {
2435      bitplace=save_bitplace;
2436      return;
2437    }
2438    // there is enough space for the optional var
2439    v[opt_var]=place | (bitplace << 24);
2440  }
2441}
2442
2443static void rO_LexVars_neg(int &place, int &bitplace, int start, int end,
2444  int &prev_ord, long *o,int *v, int bits, int opt_var)
2445{
2446  // a block of variables v_start..v_end with lex order, ordsgn -1
2447  int k;
2448  int incr=1;
2449  if(prev_ord==1) rO_Align(place,bitplace);
2450
2451  if (start>end)
2452  {
2453    incr=-1;
2454  }
2455  for(k=start;;k+=incr)
2456  {
2457    bitplace-=bits;
2458    if (bitplace < 0) { bitplace=BITS_PER_LONG-bits; place++; }
2459    o[place]=-1;
2460    v[k]=place | (bitplace << 24);
2461    if (k==end) break;
2462  }
2463  prev_ord=-1;
2464//  #if 0
2465  if (opt_var!= -1)
2466  {
2467    assume((opt_var == end+1) ||(opt_var == end-1));
2468    if((opt_var != end+1) &&(opt_var != end-1)) WarnS("hier-1");
2469    int save_bitplace=bitplace;
2470    bitplace-=bits;
2471    if (bitplace < 0)
2472    {
2473      bitplace=save_bitplace;
2474      return;
2475    }
2476    // there is enough space for the optional var
2477    v[opt_var]=place | (bitplace << 24);
2478  }
2479//  #endif
2480}
2481
2482static void rO_Syzcomp(int &place, int &bitplace, int &prev_ord,
2483    long *o, sro_ord &ord_struct)
2484{
2485  // ordering is derived from component number
2486  rO_Align(place,bitplace);
2487  ord_struct.ord_typ=ro_syzcomp;
2488  ord_struct.data.syzcomp.place=place;
2489  ord_struct.data.syzcomp.Components=NULL;
2490  ord_struct.data.syzcomp.ShiftedComponents=NULL;
2491  o[place]=1;
2492  prev_ord=1;
2493  place++;
2494  rO_Align(place,bitplace);
2495}
2496
2497static void rO_Syz(int &place, int &bitplace, int &prev_ord,
2498    long *o, sro_ord &ord_struct)
2499{
2500  // ordering is derived from component number
2501  // let's reserve one Exponent_t for it
2502  if ((prev_ord== 1) || (bitplace!=BITS_PER_LONG))
2503    rO_Align(place,bitplace);
2504  ord_struct.ord_typ=ro_syz;
2505  ord_struct.data.syz.place=place;
2506  ord_struct.data.syz.limit=0;
2507  ord_struct.data.syz.syz_index = NULL;
2508  ord_struct.data.syz.curr_index = 1;
2509  o[place]= -1;
2510  prev_ord=-1;
2511  place++;
2512}
2513
2514#ifndef NDEBUG
2515# define MYTEST 0
2516#else /* ifndef NDEBUG */
2517# define MYTEST 0
2518#endif /* ifndef NDEBUG */
2519
2520static void rO_ISPrefix(int &place, int &bitplace, int &prev_ord,
2521    long *o, int N, int *v, sro_ord &ord_struct)
2522{
2523  if ((prev_ord== 1) || (bitplace!=BITS_PER_LONG))
2524    rO_Align(place,bitplace);
2525  // since we add something afterwards - it's better to start with anew!?
2526
2527  ord_struct.ord_typ = ro_isTemp;
2528  ord_struct.data.isTemp.start = place;
2529  ord_struct.data.isTemp.pVarOffset = (int *)omMemDup(v);
2530  ord_struct.data.isTemp.suffixpos = -1;
2531
2532  // We will act as rO_Syz on our own!!!
2533  // Here we allocate an exponent as a level placeholder
2534  o[place]= -1;
2535  prev_ord=-1;
2536  place++;
2537}
2538static void rO_ISSuffix(int &place, int &bitplace, int &prev_ord, long *o,
2539  int N, int *v, sro_ord *tmp_typ, int &typ_i, int sgn)
2540{
2541
2542  // Let's find previous prefix:
2543  int typ_j = typ_i - 1;
2544  while(typ_j >= 0)
2545  {
2546    if( tmp_typ[typ_j].ord_typ == ro_isTemp)
2547      break;
2548    typ_j --;
2549  }
2550
2551  assume( typ_j >= 0 );
2552
2553  if( typ_j < 0 ) // Found NO prefix!!! :(
2554    return;
2555
2556  assume( tmp_typ[typ_j].ord_typ == ro_isTemp );
2557
2558  // Get saved state:
2559  const int start = tmp_typ[typ_j].data.isTemp.start;
2560  int *pVarOffset = tmp_typ[typ_j].data.isTemp.pVarOffset;
2561
2562/*
2563  // shift up all blocks
2564  while(typ_j < (typ_i-1))
2565  {
2566    tmp_typ[typ_j] = tmp_typ[typ_j+1];
2567    typ_j++;
2568  }
2569  typ_j = typ_i - 1; // No increment for typ_i
2570*/
2571  tmp_typ[typ_j].data.isTemp.suffixpos = typ_i;
2572
2573  // Let's keep that dummy for now...
2574  typ_j = typ_i; // the typ to change!
2575  typ_i++; // Just for now...
2576
2577
2578  for( int i = 0; i <= N; i++ ) // Note [0] == component !!! No Skip?
2579  {
2580    // Was i-th variable allocated inbetween?
2581    if( v[i] != pVarOffset[i] )
2582    {
2583      pVarOffset[i] = v[i]; // Save for later...
2584      v[i] = -1; // Undo!
2585      assume( pVarOffset[i] != -1 );
2586    }
2587    else
2588      pVarOffset[i] = -1; // No change here...
2589  }
2590
2591  if( pVarOffset[0] != -1 )
2592    pVarOffset[0] &= 0x0fff;
2593
2594  sro_ord &ord_struct = tmp_typ[typ_j];
2595
2596
2597  ord_struct.ord_typ = ro_is;
2598  ord_struct.data.is.start = start;
2599  ord_struct.data.is.end   = place;
2600  ord_struct.data.is.pVarOffset = pVarOffset;
2601
2602
2603  // What about component???
2604//   if( v[0] != -1 ) // There is a component already...???
2605//     if( o[ v[0] & 0x0fff ] == sgn )
2606//     {
2607//       pVarOffset[0] = -1; // NEVER USED Afterwards...
2608//       return;
2609//     }
2610
2611
2612  // Moreover: we need to allocate the module component (v[0]) here!
2613  if( v[0] == -1) // It's possible that there was module component v0 at the begining (before prefix)!
2614  {
2615    // Start with a whole long exponent
2616    if( bitplace != BITS_PER_LONG )
2617      rO_Align(place, bitplace);
2618
2619    assume( bitplace == BITS_PER_LONG );
2620    bitplace -= BITS_PER_LONG;
2621    assume(bitplace == 0);
2622    v[0] = place | (bitplace << 24); // Never mind whether pVarOffset[0] > 0!!!
2623    o[place] = sgn; // Singnum for component ordering
2624    prev_ord = sgn;
2625  }
2626}
2627
2628
2629static unsigned long rGetExpSize(unsigned long bitmask, int & bits)
2630{
2631  if (bitmask == 0)
2632  {
2633    bits=16; bitmask=0xffff;
2634  }
2635  else if (bitmask <= 1L)
2636  {
2637    bits=1; bitmask = 1L;
2638  }
2639  else if (bitmask <= 3L)
2640  {
2641    bits=2; bitmask = 3L;
2642  }
2643  else if (bitmask <= 7L)
2644  {
2645    bits=3; bitmask=7L;
2646  }
2647  else if (bitmask <= 0xfL)
2648  {
2649    bits=4; bitmask=0xfL;
2650  }
2651  else if (bitmask <= 0x1fL)
2652  {
2653    bits=5; bitmask=0x1fL;
2654  }
2655  else if (bitmask <= 0x3fL)
2656  {
2657    bits=6; bitmask=0x3fL;
2658  }
2659#if SIZEOF_LONG == 8
2660  else if (bitmask <= 0x7fL)
2661  {
2662    bits=7; bitmask=0x7fL; /* 64 bit longs only */
2663  }
2664#endif
2665  else if (bitmask <= 0xffL)
2666  {
2667    bits=8; bitmask=0xffL;
2668  }
2669#if SIZEOF_LONG == 8
2670  else if (bitmask <= 0x1ffL)
2671  {
2672    bits=9; bitmask=0x1ffL; /* 64 bit longs only */
2673  }
2674#endif
2675  else if (bitmask <= 0x3ffL)
2676  {
2677    bits=10; bitmask=0x3ffL;
2678  }
2679#if SIZEOF_LONG == 8
2680  else if (bitmask <= 0xfffL)
2681  {
2682    bits=12; bitmask=0xfff; /* 64 bit longs only */
2683  }
2684#endif
2685  else if (bitmask <= 0xffffL)
2686  {
2687    bits=16; bitmask=0xffffL;
2688  }
2689#if SIZEOF_LONG == 8
2690  else if (bitmask <= 0xfffffL)
2691  {
2692    bits=20; bitmask=0xfffffL; /* 64 bit longs only */
2693  }
2694  else if (bitmask <= 0xffffffffL)
2695  {
2696    bits=32; bitmask=0xffffffffL;
2697  }
2698  else if (bitmask <= 0x7fffffffffffffffL)
2699  {
2700    bits=63; bitmask=0x7fffffffffffffffL; /* for overflow tests*/
2701  }
2702  else
2703  {
2704    bits=63; bitmask=0x7fffffffffffffffL; /* for overflow tests*/
2705  }
2706#else
2707  else if (bitmask <= 0x7fffffff)
2708  {
2709    bits=31; bitmask=0x7fffffff; /* for overflow tests*/
2710  }
2711  else
2712  {
2713    bits=31; bitmask=0x7fffffffL; /* for overflow tests*/
2714  }
2715#endif
2716  return bitmask;
2717}
2718
2719/*2
2720* optimize rGetExpSize for a block of N variables, exp <=bitmask
2721*/
2722static unsigned long rGetExpSize(unsigned long bitmask, int & bits, int N)
2723{
2724  bitmask =rGetExpSize(bitmask, bits);
2725  int vars_per_long=BIT_SIZEOF_LONG/bits;
2726  int bits1;
2727  loop
2728  {
2729    if (bits == BIT_SIZEOF_LONG-1)
2730    {
2731      bits =  BIT_SIZEOF_LONG - 1;
2732      return LONG_MAX;
2733    }
2734    unsigned long bitmask1 =rGetExpSize(bitmask+1, bits1);
2735    int vars_per_long1=BIT_SIZEOF_LONG/bits1;
2736    if ((((N+vars_per_long-1)/vars_per_long) ==
2737         ((N+vars_per_long1-1)/vars_per_long1)))
2738    {
2739      vars_per_long=vars_per_long1;
2740      bits=bits1;
2741      bitmask=bitmask1;
2742    }
2743    else
2744    {
2745      return bitmask; /* and bits */
2746    }
2747  }
2748}
2749
2750
2751/*2
2752 * create a copy of the ring r, which must be equivalent to currRing
2753 * used for std computations
2754 * may share data structures with currRing
2755 * DOES CALL rComplete
2756 */
2757ring rModifyRing(ring r, BOOLEAN omit_degree,
2758                         BOOLEAN omit_comp,
2759                         unsigned long exp_limit)
2760{
2761  assume (r != NULL );
2762  assume (exp_limit > 1);
2763  BOOLEAN need_other_ring;
2764  BOOLEAN omitted_degree = FALSE;
2765
2766  int iNeedInducedOrderingSetup = 0; ///< How many induced ordering block do we have?
2767  int bits;
2768
2769  exp_limit=rGetExpSize(exp_limit, bits, r->N);
2770  need_other_ring = (exp_limit != r->bitmask);
2771
2772  int nblocks=rBlocks(r);
2773  int *order=(int*)omAlloc0((nblocks+1)*sizeof(int));
2774  int *block0=(int*)omAlloc0((nblocks+1)*sizeof(int));
2775  int *block1=(int*)omAlloc0((nblocks+1)*sizeof(int));
2776  int **wvhdl=(int**)omAlloc0((nblocks+1)*sizeof(int *));
2777
2778  int i=0;
2779  int j=0; /*  i index in r, j index in res */
2780
2781  for( int r_ord=r->order[i]; (r_ord != 0) && (i < nblocks); j++, r_ord=r->order[++i])
2782  {
2783    BOOLEAN copy_block_index=TRUE;
2784
2785    if (r->block0[i]==r->block1[i])
2786    {
2787      switch(r_ord)
2788      {
2789        case ringorder_wp:
2790        case ringorder_dp:
2791        case ringorder_Wp:
2792        case ringorder_Dp:
2793          r_ord=ringorder_lp;
2794          break;
2795        case ringorder_Ws:
2796        case ringorder_Ds:
2797        case ringorder_ws:
2798        case ringorder_ds:
2799          r_ord=ringorder_ls;
2800          break;
2801        default:
2802          break;
2803      }
2804    }
2805    switch(r_ord)
2806    {
2807      case ringorder_S:
2808      {
2809#ifndef NDEBUG
2810        Warn("Error: unhandled ordering in rModifyRing: ringorder_S = [%d]", r_ord);
2811#endif
2812        order[j]=r_ord; /*r->order[i];*/
2813        break;
2814      }
2815      case ringorder_C:
2816      case ringorder_c:
2817        if (!omit_comp)
2818        {
2819          order[j]=r_ord; /*r->order[i]*/;
2820        }
2821        else
2822        {
2823          j--;
2824          need_other_ring=TRUE;
2825          omit_comp=FALSE;
2826          copy_block_index=FALSE;
2827        }
2828        break;
2829      case ringorder_wp:
2830      case ringorder_dp:
2831      case ringorder_ws:
2832      case ringorder_ds:
2833        if(!omit_degree)
2834        {
2835          order[j]=r_ord; /*r->order[i]*/;
2836        }
2837        else
2838        {
2839          order[j]=ringorder_rs;
2840          need_other_ring=TRUE;
2841          omit_degree=FALSE;
2842          omitted_degree = TRUE;
2843        }
2844        break;
2845      case ringorder_Wp:
2846      case ringorder_Dp:
2847      case ringorder_Ws:
2848      case ringorder_Ds:
2849        if(!omit_degree)
2850        {
2851          order[j]=r_ord; /*r->order[i];*/
2852        }
2853        else
2854        {
2855          order[j]=ringorder_lp;
2856          need_other_ring=TRUE;
2857          omit_degree=FALSE;
2858          omitted_degree = TRUE;
2859        }
2860        break;
2861      case ringorder_IS:
2862      {
2863        if (omit_comp)
2864        {
2865#ifndef NDEBUG
2866          Warn("Error: WRONG USAGE of rModifyRing: cannot omit component due to the ordering block [%d]: %d (ringorder_IS)", i, r_ord);
2867#endif
2868          omit_comp = FALSE;
2869        }
2870        order[j]=r_ord; /*r->order[i];*/
2871        iNeedInducedOrderingSetup++;
2872        break;
2873      }
2874      case ringorder_s:
2875      {
2876        assume((i == 0) && (j == 0));
2877        if (omit_comp)
2878        {
2879#ifndef NDEBUG
2880          Warn("WRONG USAGE? of rModifyRing: omitting component due to the ordering block [%d]: %d (ringorder_s)", i, r_ord);
2881#endif
2882          omit_comp = FALSE;
2883        }
2884        order[j]=r_ord; /*r->order[i];*/
2885        break;
2886      }
2887      default:
2888        order[j]=r_ord; /*r->order[i];*/
2889        break;
2890    }
2891    if (copy_block_index)
2892    {
2893      block0[j]=r->block0[i];
2894      block1[j]=r->block1[i];
2895      wvhdl[j]=r->wvhdl[i];
2896    }
2897
2898    // order[j]=ringorder_no; //  done by omAlloc0
2899  }
2900  if(!need_other_ring)
2901  {
2902    omFreeSize(order,(nblocks+1)*sizeof(int));
2903    omFreeSize(block0,(nblocks+1)*sizeof(int));
2904    omFreeSize(block1,(nblocks+1)*sizeof(int));
2905    omFreeSize(wvhdl,(nblocks+1)*sizeof(int *));
2906    return r;
2907  }
2908  ring res=(ring)omAlloc0Bin(sip_sring_bin);
2909  *res = *r;
2910
2911#ifdef HAVE_PLURAL
2912  res->GetNC() = NULL;
2913#endif
2914
2915  // res->qideal, res->idroot ???
2916  res->wvhdl=wvhdl;
2917  res->order=order;
2918  res->block0=block0;
2919  res->block1=block1;
2920  res->bitmask=exp_limit;
2921  int tmpref=r->cf->ref;
2922  rComplete(res, 1);
2923  r->cf->ref=tmpref;
2924
2925  // adjust res->pFDeg: if it was changed globally, then
2926  // it must also be changed for new ring
2927  if (r->pFDegOrig != res->pFDegOrig &&
2928           rOrd_is_WeightedDegree_Ordering(r))
2929  {
2930    // still might need adjustment for weighted orderings
2931    // and omit_degree
2932    res->firstwv = r->firstwv;
2933    res->firstBlockEnds = r->firstBlockEnds;
2934    res->pFDeg = res->pFDegOrig = p_WFirstTotalDegree;
2935  }
2936  if (omitted_degree)
2937    res->pLDeg = res->pLDegOrig = r->pLDegOrig;
2938
2939  rOptimizeLDeg(res);
2940
2941  // set syzcomp
2942  if (res->typ != NULL)
2943  {
2944    if( res->typ[0].ord_typ == ro_syz) // "s" Always on [0] place!
2945    {
2946      res->typ[0] = r->typ[0]; // Copy struct!? + setup the same limit!
2947
2948      if (r->typ[0].data.syz.limit > 0)
2949      {
2950        res->typ[0].data.syz.syz_index
2951          = (int*) omAlloc((r->typ[0].data.syz.limit +1)*sizeof(int));
2952        memcpy(res->typ[0].data.syz.syz_index, r->typ[0].data.syz.syz_index,
2953              (r->typ[0].data.syz.limit +1)*sizeof(int));
2954      }
2955    }
2956
2957    if( iNeedInducedOrderingSetup > 0 )
2958    {
2959      for(j = 0, i = 0; (i < nblocks) && (iNeedInducedOrderingSetup > 0); i++)
2960        if( res->typ[i].ord_typ == ro_is ) // Search for suffixes!
2961        {
2962          ideal F = idrHeadR(r->typ[i].data.is.F, r, res); // Copy F from r into res!
2963          assume(
2964            rSetISReference( res, 
2965              F,  // WILL BE COPIED!
2966              r->typ[i].data.is.limit,
2967              j++,
2968              r->typ[i].data.is.componentWeights // WILL BE COPIED
2969              )
2970            );
2971          id_Delete(&F, res);
2972          iNeedInducedOrderingSetup--;
2973        }
2974    } // Process all induced Ordering blocks! ...
2975  }
2976  // the special case: homog (omit_degree) and 1 block rs: that is global:
2977  // it comes from dp
2978  res->OrdSgn=r->OrdSgn;
2979
2980
2981#ifdef HAVE_PLURAL
2982  if (rIsPluralRing(r))
2983  {
2984    if ( nc_rComplete(r, res, false) ) // no qideal!
2985    {
2986#ifndef NDEBUG
2987      WarnS("error in nc_rComplete");
2988#endif
2989      // cleanup?
2990
2991//      rDelete(res);
2992//      return r;
2993
2994      // just go on..
2995    }
2996
2997    if( rIsSCA(r) )
2998    {
2999      if( !sca_Force(res, scaFirstAltVar(r), scaLastAltVar(r)) )
3000      WarnS("error in sca_Force!");
3001    }
3002  }
3003#endif
3004
3005  return res;
3006}
3007
3008// construct Wp,C ring
3009ring rModifyRing_Wp(ring r, int* weights)
3010{
3011  ring res=(ring)omAlloc0Bin(sip_sring_bin);
3012  *res = *r;
3013#ifdef HAVE_PLURAL
3014  res->GetNC() = NULL;
3015#endif
3016
3017  /*weights: entries for 3 blocks: NULL*/
3018  res->wvhdl = (int **)omAlloc0(3 * sizeof(int *));
3019  /*order: Wp,C,0*/
3020  res->order = (int *) omAlloc(3 * sizeof(int *));
3021  res->block0 = (int *)omAlloc0(3 * sizeof(int *));
3022  res->block1 = (int *)omAlloc0(3 * sizeof(int *));
3023  /* ringorder Wp for the first block: var 1..r->N */
3024  res->order[0]  = ringorder_Wp;
3025  res->block0[0] = 1;
3026  res->block1[0] = r->N;
3027  res->wvhdl[0] = weights;
3028  /* ringorder C for the second block: no vars */
3029  res->order[1]  = ringorder_C;
3030  /* the last block: everything is 0 */
3031  res->order[2]  = 0;
3032  /*polynomial ring*/
3033  res->OrdSgn    = 1;
3034
3035  int tmpref=r->cf->ref;
3036  rComplete(res, 1);
3037  r->cf->ref=tmpref;
3038#ifdef HAVE_PLURAL
3039  if (rIsPluralRing(r))
3040  {
3041    if ( nc_rComplete(r, res, false) ) // no qideal!
3042    {
3043#ifndef NDEBUG
3044      WarnS("error in nc_rComplete");
3045#endif
3046      // cleanup?
3047
3048//      rDelete(res);
3049//      return r;
3050
3051      // just go on..
3052    }
3053  }
3054#endif
3055  return res;
3056}
3057
3058// construct lp, C ring with r->N variables, r->names vars....
3059ring rModifyRing_Simple(ring r, BOOLEAN ommit_degree, BOOLEAN ommit_comp, unsigned long exp_limit, BOOLEAN &simple)
3060{
3061  simple=TRUE;
3062  if (!rHasSimpleOrder(r))
3063  {
3064    simple=FALSE; // sorting needed
3065    assume (r != NULL );
3066    assume (exp_limit > 1);
3067    int bits;
3068
3069    exp_limit=rGetExpSize(exp_limit, bits, r->N);
3070
3071    int nblocks=1+(ommit_comp!=0);
3072    int *order=(int*)omAlloc0((nblocks+1)*sizeof(int));
3073    int *block0=(int*)omAlloc0((nblocks+1)*sizeof(int));
3074    int *block1=(int*)omAlloc0((nblocks+1)*sizeof(int));
3075    int **wvhdl=(int**)omAlloc0((nblocks+1)*sizeof(int *));
3076
3077    order[0]=ringorder_lp;
3078    block0[0]=1;
3079    block1[0]=r->N;
3080    if (!ommit_comp)
3081    {
3082      order[1]=ringorder_C;
3083    }
3084    ring res=(ring)omAlloc0Bin(sip_sring_bin);
3085    *res = *r;
3086#ifdef HAVE_PLURAL
3087    res->GetNC() = NULL;
3088#endif
3089    // res->qideal, res->idroot ???
3090    res->wvhdl=wvhdl;
3091    res->order=order;
3092    res->block0=block0;
3093    res->block1=block1;
3094    res->bitmask=exp_limit;
3095    int tmpref=r->cf->ref;
3096    rComplete(res, 1);
3097    r->cf->ref=tmpref;
3098
3099#ifdef HAVE_PLURAL
3100    if (rIsPluralRing(r))
3101    {
3102      if ( nc_rComplete(r, res, false) ) // no qideal!
3103      {
3104#ifndef NDEBUG
3105        WarnS("error in nc_rComplete");
3106#endif
3107        // cleanup?
3108
3109//      rDelete(res);
3110//      return r;
3111
3112      // just go on..
3113      }
3114    }
3115#endif
3116
3117    rOptimizeLDeg(res);
3118
3119    return res;
3120  }
3121  return rModifyRing(r, ommit_degree, ommit_comp, exp_limit);
3122}
3123
3124void rKillModifiedRing_Simple(ring r)
3125{
3126  rKillModifiedRing(r);
3127}
3128
3129
3130void rKillModifiedRing(ring r)
3131{
3132  rUnComplete(r);
3133  omFree(r->order);
3134  omFree(r->block0);
3135  omFree(r->block1);
3136  omFree(r->wvhdl);
3137  omFreeBin(r,sip_sring_bin);
3138}
3139
3140void rKillModified_Wp_Ring(ring r)
3141{
3142  rUnComplete(r);
3143  omFree(r->order);
3144  omFree(r->block0);
3145  omFree(r->block1);
3146  omFree(r->wvhdl[0]);
3147  omFree(r->wvhdl);
3148  omFreeBin(r,sip_sring_bin);
3149}
3150
3151static void rSetOutParams(ring r)
3152{
3153  r->VectorOut = (r->order[0] == ringorder_c);
3154  r->ShortOut = TRUE;
3155  {
3156    int i;
3157    if (r->cf->parameter!=NULL)
3158    {
3159      for (i=0;i<rPar(r);i++)
3160      {
3161        if(strlen(r->cf->parameter[i])>1)
3162        {
3163          r->ShortOut=FALSE;
3164          break;
3165        }
3166      }
3167    }
3168    if (r->ShortOut)
3169    {
3170      // Hmm... sometimes (e.g., from maGetPreimage) new variables
3171      // are introduced, but their names are never set
3172      // hence, we do the following awkward trick
3173      int N = omSizeWOfAddr(r->names);
3174      if (r->N < N) N = r->N;
3175
3176      for (i=(N-1);i>=0;i--)
3177      {
3178        if(r->names[i] != NULL && strlen(r->names[i])>1)
3179        {
3180          r->ShortOut=FALSE;
3181          break;
3182        }
3183      }
3184    }
3185  }
3186  r->CanShortOut = r->ShortOut;
3187}
3188
3189/*2
3190* sets r->MixedOrder and r->ComponentOrder for orderings with more than one block
3191* block of variables (ip is the block number, o_r the number of the ordering)
3192* o is the position of the orderingering in r
3193*/
3194static void rHighSet(ring r, int o_r, int o)
3195{
3196  switch(o_r)
3197  {
3198    case ringorder_lp:
3199    case ringorder_dp:
3200    case ringorder_Dp:
3201    case ringorder_wp:
3202    case ringorder_Wp:
3203    case ringorder_rp:
3204    case ringorder_a:
3205    case ringorder_aa:
3206    case ringorder_a64:
3207      if (r->OrdSgn==-1) r->MixedOrder=TRUE;
3208      break;
3209    case ringorder_ls:
3210    case ringorder_rs:
3211    case ringorder_ds:
3212    case ringorder_Ds:
3213    case ringorder_s:
3214      break;
3215    case ringorder_ws:
3216    case ringorder_Ws:
3217      if (r->wvhdl[o]!=NULL)
3218      {
3219        int i;
3220        for(i=r->block1[o]-r->block0[o];i>=0;i--)
3221          if (r->wvhdl[o][i]<0) { r->MixedOrder=TRUE; break; }
3222      }
3223      break;
3224    case ringorder_c:
3225      r->ComponentOrder=1;
3226      break;
3227    case ringorder_C:
3228    case ringorder_S:
3229      r->ComponentOrder=-1;
3230      break;
3231    case ringorder_M:
3232      r->LexOrder=TRUE;
3233      break;
3234    case ringorder_IS:
3235    { // TODO: What is r->ComponentOrder???
3236      r->MixedOrder=TRUE;
3237/*
3238      if( r->block0[o] != 0 ) // Suffix has the comonent
3239        r->ComponentOrder = r->block0[o];
3240      else // Prefix has level...
3241        r->ComponentOrder=-1;
3242*/
3243      break;
3244    }
3245
3246    default:
3247      dReportError("wrong internal ordering:%d at %s, l:%d\n",o_r,__FILE__,__LINE__);
3248  }
3249}
3250
3251static void rSetFirstWv(ring r, int i, int* order, int* block1, int** wvhdl)
3252{
3253  // cheat for ringorder_aa
3254  if (order[i] == ringorder_aa)
3255    i++;
3256  if(block1[i]!=r->N) r->LexOrder=TRUE;
3257  r->firstBlockEnds=block1[i];
3258  r->firstwv = wvhdl[i];
3259  if ((order[i]== ringorder_ws)
3260  || (order[i]==ringorder_Ws)
3261  || (order[i]== ringorder_wp)
3262  || (order[i]==ringorder_Wp)
3263  || (order[i]== ringorder_a)
3264   /*|| (order[i]==ringorder_A)*/)
3265  {
3266    int j;
3267    for(j=block1[i]-r->block0[i];j>=0;j--)
3268    {
3269      if (r->firstwv[j]<0) r->MixedOrder=TRUE;
3270      if (r->firstwv[j]==0) r->LexOrder=TRUE;
3271    }
3272  }
3273  else if (order[i]==ringorder_a64)
3274  {
3275    int j;
3276    int64 *w=rGetWeightVec(r);
3277    for(j=block1[i]-r->block0[i];j>=0;j--)
3278    {
3279      if (w[j]==0) r->LexOrder=TRUE;
3280    }
3281  }
3282}
3283
3284static void rOptimizeLDeg(ring r)
3285{
3286  if (r->pFDeg == p_Deg)
3287  {
3288    if (r->pLDeg == pLDeg1)
3289      r->pLDeg = pLDeg1_Deg;
3290    if (r->pLDeg == pLDeg1c)
3291      r->pLDeg = pLDeg1c_Deg;
3292  }
3293  else if (r->pFDeg == p_Totaldegree)
3294  {
3295    if (r->pLDeg == pLDeg1)
3296      r->pLDeg = pLDeg1_Totaldegree;
3297    if (r->pLDeg == pLDeg1c)
3298      r->pLDeg = pLDeg1c_Totaldegree;
3299  }
3300  else if (r->pFDeg == p_WFirstTotalDegree)
3301  {
3302    if (r->pLDeg == pLDeg1)
3303      r->pLDeg = pLDeg1_WFirstTotalDegree;
3304    if (r->pLDeg == pLDeg1c)
3305      r->pLDeg = pLDeg1c_WFirstTotalDegree;
3306  }
3307}
3308
3309// set pFDeg, pLDeg, MixOrder, ComponentOrder, etc
3310static void rSetDegStuff(ring r)
3311{
3312  int* order = r->order;
3313  int* block0 = r->block0;
3314  int* block1 = r->block1;
3315  int** wvhdl = r->wvhdl;
3316
3317  if (order[0]==ringorder_S ||order[0]==ringorder_s || order[0]==ringorder_IS)
3318  {
3319    order++;
3320    block0++;
3321    block1++;
3322    wvhdl++;
3323  }
3324  r->LexOrder = FALSE;
3325  r->MixedOrder = FALSE;
3326  r->ComponentOrder = 1;
3327  r->pFDeg = p_Totaldegree;
3328  r->pLDeg = (r->OrdSgn == 1 ? pLDegb : pLDeg0);
3329
3330  /*======== ordering type is (_,c) =========================*/
3331  if ((order[0]==ringorder_unspec) || (order[1] == 0)
3332      ||(
3333    ((order[1]==ringorder_c)||(order[1]==ringorder_C)
3334     ||(order[1]==ringorder_S) 
3335     ||(order[1]==ringorder_s))
3336    && (order[0]!=ringorder_M)
3337    && (order[2]==0))
3338    )
3339  {
3340    if ((order[0]!=ringorder_unspec)
3341    && ((order[1]==ringorder_C)||(order[1]==ringorder_S)||
3342        (order[1]==ringorder_s)))
3343      r->ComponentOrder=-1;
3344    if (r->OrdSgn == -1) r->pLDeg = pLDeg0c;
3345    if ((order[0] == ringorder_lp)
3346    || (order[0] == ringorder_ls)
3347    || (order[0] == ringorder_rp)
3348    || (order[0] == ringorder_rs))
3349    {
3350      r->LexOrder=TRUE;
3351      r->pLDeg = pLDeg1c;
3352      r->pFDeg = p_Totaldegree;
3353    }
3354    if ((order[0] == ringorder_a)
3355    || (order[0] == ringorder_wp)
3356    || (order[0] == ringorder_Wp)
3357    || (order[0] == ringorder_ws)
3358    || (order[0] == ringorder_Ws))
3359      r->pFDeg = p_WFirstTotalDegree;
3360    r->firstBlockEnds=block1[0];
3361    r->firstwv = wvhdl[0];
3362  }
3363  /*======== ordering type is (c,_) =========================*/
3364  else if (((order[0]==ringorder_c)
3365            ||(order[0]==ringorder_C)
3366            ||(order[0]==ringorder_S)
3367            ||(order[0]==ringorder_s))
3368  && (order[1]!=ringorder_M)
3369  &&  (order[2]==0))
3370  {
3371    if ((order[0]==ringorder_C)||(order[0]==ringorder_S)||
3372        order[0]==ringorder_s)
3373      r->ComponentOrder=-1;
3374    if ((order[1] == ringorder_lp)
3375    || (order[1] == ringorder_ls)
3376    || (order[1] == ringorder_rp)
3377    || order[1] == ringorder_rs)
3378    {
3379      r->LexOrder=TRUE;
3380      r->pLDeg = pLDeg1c;
3381      r->pFDeg = p_Totaldegree;
3382    }
3383    r->firstBlockEnds=block1[1];
3384    if (wvhdl!=NULL) r->firstwv = wvhdl[1];
3385    if ((order[1] == ringorder_a)
3386    || (order[1] == ringorder_wp)
3387    || (order[1] == ringorder_Wp)
3388    || (order[1] == ringorder_ws)
3389    || (order[1] == ringorder_Ws))
3390      r->pFDeg = p_WFirstTotalDegree;
3391  }
3392  /*------- more than one block ----------------------*/
3393  else
3394  {
3395    if ((r->VectorOut)||(order[0]==ringorder_C)||(order[0]==ringorder_S)||(order[0]==ringorder_s))
3396    {
3397      rSetFirstWv(r, 1, order, block1, wvhdl);
3398    }
3399    else
3400      rSetFirstWv(r, 0, order, block1, wvhdl);
3401
3402    /*the number of orderings:*/
3403    int i = 0;
3404    while (order[++i] != 0);
3405    do
3406    {
3407      i--;
3408      rHighSet(r, order[i],i);
3409    }
3410    while (i != 0);
3411
3412    if ((order[0]!=ringorder_c)
3413        && (order[0]!=ringorder_C)
3414        && (order[0]!=ringorder_S)
3415        && (order[0]!=ringorder_s))
3416    {
3417      r->pLDeg = pLDeg1c;
3418    }
3419    else
3420    {
3421      r->pLDeg = pLDeg1;
3422    }
3423    r->pFDeg = p_WTotaldegree; // may be improved: p_Totaldegree for lp/dp/ls/.. blocks
3424  }
3425 
3426  if (rOrd_is_Totaldegree_Ordering(r) || rOrd_is_WeightedDegree_Ordering(r))
3427    r->pFDeg = p_Deg;
3428
3429  r->pFDegOrig = r->pFDeg;
3430  r->pLDegOrig = r->pLDeg;
3431  rOptimizeLDeg(r);
3432}
3433
3434/*2
3435* set NegWeightL_Size, NegWeightL_Offset
3436*/
3437static void rSetNegWeight(ring r)
3438{
3439  int i,l;
3440  if (r->typ!=NULL)
3441  {
3442    l=0;
3443    for(i=0;i<r->OrdSize;i++)
3444    {
3445      if(r->typ[i].ord_typ==ro_wp_neg) l++;
3446    }
3447    if (l>0)
3448    {
3449      r->NegWeightL_Size=l;
3450      r->NegWeightL_Offset=(int *) omAlloc(l*sizeof(int));
3451      l=0;
3452      for(i=0;i<r->OrdSize;i++)
3453      {
3454        if(r->typ[i].ord_typ==ro_wp_neg)
3455        {
3456          r->NegWeightL_Offset[l]=r->typ[i].data.wp.place;
3457          l++;
3458        }
3459      }
3460      return;
3461    }
3462  }
3463  r->NegWeightL_Size = 0;
3464  r->NegWeightL_Offset = NULL;
3465}
3466
3467static void rSetOption(ring r)
3468{
3469  // set redthrough
3470  if (!TEST_OPT_OLDSTD && r->OrdSgn == 1 && ! r->LexOrder)
3471    r->options |= Sy_bit(OPT_REDTHROUGH);
3472  else
3473    r->options &= ~Sy_bit(OPT_REDTHROUGH);
3474
3475  // set intStrategy
3476#ifdef HAVE_RINGS
3477  if (
3478         rField_is_Extension(r)
3479      || rField_is_Q(r)
3480      || rField_is_Ring(r))
3481#else
3482  if (rField_is_Extension(r) || rField_is_Q(r))
3483#endif
3484    r->options |= Sy_bit(OPT_INTSTRATEGY);
3485  else
3486    r->options &= ~Sy_bit(OPT_INTSTRATEGY);
3487
3488  // set redTail
3489  if (r->LexOrder || r->OrdSgn == -1 || rField_is_Extension(r))
3490    r->options &= ~Sy_bit(OPT_REDTAIL);
3491  else
3492    r->options |= Sy_bit(OPT_REDTAIL);
3493}
3494
3495static void rCheckOrdSgn(ring r,int i/*current block*/);
3496
3497BOOLEAN rComplete(ring r, int force)
3498{
3499  if (r->VarOffset!=NULL && force == 0) return FALSE;
3500  rSetOutParams(r);
3501  int n=rBlocks(r)-1;
3502  int i;
3503  int bits;
3504  r->bitmask=rGetExpSize(r->bitmask,bits,r->N);
3505  r->BitsPerExp = bits;
3506  r->ExpPerLong = BIT_SIZEOF_LONG / bits;
3507  r->divmask=rGetDivMask(bits);
3508
3509  // will be used for ordsgn:
3510  long *tmp_ordsgn=(long *)omAlloc0(3*(n+r->N)*sizeof(long));
3511  // will be used for VarOffset:
3512  int *v=(int *)omAlloc((r->N+1)*sizeof(int));
3513  for(i=r->N; i>=0 ; i--)
3514  {
3515    v[i]=-1;
3516  }
3517  sro_ord *tmp_typ=(sro_ord *)omAlloc0(3*(n+r->N)*sizeof(sro_ord));
3518  int typ_i=0;
3519  int prev_ordsgn=0;
3520
3521  // fill in v, tmp_typ, tmp_ordsgn, determine typ_i (== ordSize)
3522  int j=0;
3523  int j_bits=BITS_PER_LONG;
3524
3525  BOOLEAN need_to_add_comp=FALSE; // Only for ringorder_s and ringorder_S!
3526
3527  for(i=0;i<n;i++)
3528  {
3529    tmp_typ[typ_i].order_index=i;
3530    switch (r->order[i])
3531    {
3532      case ringorder_a:
3533      case ringorder_aa:
3534        rO_WDegree(j,j_bits,r->block0[i],r->block1[i],tmp_ordsgn,tmp_typ[typ_i],
3535                   r->wvhdl[i]);
3536        typ_i++;
3537        break;
3538
3539      case ringorder_a64:
3540        rO_WDegree64(j,j_bits,r->block0[i],r->block1[i],tmp_ordsgn,
3541                     tmp_typ[typ_i], (int64 *)(r->wvhdl[i]));
3542        typ_i++;
3543        break;
3544
3545      case ringorder_c:
3546        rO_Align(j, j_bits);
3547        rO_LexVars_neg(j, j_bits, 0,0, prev_ordsgn,tmp_ordsgn,v,BITS_PER_LONG, -1);
3548        break;
3549
3550      case ringorder_C:
3551        rO_Align(j, j_bits);
3552        rO_LexVars(j, j_bits, 0,0, prev_ordsgn,tmp_ordsgn,v,BITS_PER_LONG, -1);
3553        break;
3554
3555      case ringorder_M:
3556        {
3557          int k,l;
3558          k=r->block1[i]-r->block0[i]+1; // number of vars
3559          for(l=0;l<k;l++)
3560          {
3561            rO_WDegree(j,j_bits,r->block0[i],r->block1[i],tmp_ordsgn,
3562                       tmp_typ[typ_i],
3563                       r->wvhdl[i]+(r->block1[i]-r->block0[i]+1)*l);
3564            typ_i++;
3565          }
3566          break;
3567        }
3568
3569      case ringorder_lp:
3570        rO_LexVars(j, j_bits, r->block0[i],r->block1[i], prev_ordsgn,
3571                   tmp_ordsgn,v,bits, -1);
3572        break;
3573
3574      case ringorder_ls:
3575        rO_LexVars_neg(j, j_bits, r->block0[i],r->block1[i], prev_ordsgn,
3576                       tmp_ordsgn,v, bits, -1);
3577        rCheckOrdSgn(r,i);
3578        break;
3579
3580      case ringorder_rs:
3581        rO_LexVars_neg(j, j_bits, r->block1[i],r->block0[i], prev_ordsgn,
3582                       tmp_ordsgn,v, bits, -1);
3583        rCheckOrdSgn(r,i);
3584        break;
3585
3586      case ringorder_rp:
3587        rO_LexVars(j, j_bits, r->block1[i],r->block0[i], prev_ordsgn,
3588                       tmp_ordsgn,v, bits, -1);
3589        break;
3590
3591      case ringorder_dp:
3592        if (r->block0[i]==r->block1[i])
3593        {
3594          rO_LexVars(j, j_bits, r->block0[i],r->block0[i], prev_ordsgn,
3595                     tmp_ordsgn,v, bits, -1);
3596        }
3597        else
3598        {
3599          rO_TDegree(j,j_bits,r->block0[i],r->block1[i],tmp_ordsgn,
3600                     tmp_typ[typ_i]);
3601          typ_i++;
3602          rO_LexVars_neg(j, j_bits, r->block1[i],r->block0[i]+1,
3603                         prev_ordsgn,tmp_ordsgn,v,bits, r->block0[i]);
3604        }
3605        break;
3606
3607      case ringorder_Dp:
3608        if (r->block0[i]==r->block1[i])
3609        {
3610          rO_LexVars(j, j_bits, r->block0[i],r->block0[i], prev_ordsgn,
3611                     tmp_ordsgn,v, bits, -1);
3612        }
3613        else
3614        {
3615          rO_TDegree(j,j_bits,r->block0[i],r->block1[i],tmp_ordsgn,
3616                     tmp_typ[typ_i]);
3617          typ_i++;
3618          rO_LexVars(j, j_bits, r->block0[i],r->block1[i]-1, prev_ordsgn,
3619                     tmp_ordsgn,v, bits, r->block1[i]);
3620        }
3621        break;
3622
3623      case ringorder_ds:
3624        if (r->block0[i]==r->block1[i])
3625        {
3626          rO_LexVars_neg(j, j_bits,r->block0[i],r->block1[i],prev_ordsgn,
3627                         tmp_ordsgn,v,bits, -1);
3628        }
3629        else
3630        {
3631          rO_TDegree_neg(j,j_bits,r->block0[i],r->block1[i],tmp_ordsgn,
3632                         tmp_typ[typ_i]);
3633          typ_i++;
3634          rO_LexVars_neg(j, j_bits, r->block1[i],r->block0[i]+1,
3635                         prev_ordsgn,tmp_ordsgn,v,bits, r->block0[i]);
3636        }
3637        rCheckOrdSgn(r,i);
3638        break;
3639
3640      case ringorder_Ds:
3641        if (r->block0[i]==r->block1[i])
3642        {
3643          rO_LexVars_neg(j, j_bits, r->block0[i],r->block0[i],prev_ordsgn,
3644                         tmp_ordsgn,v, bits, -1);
3645        }
3646        else
3647        {
3648          rO_TDegree_neg(j,j_bits,r->block0[i],r->block1[i],tmp_ordsgn,
3649                         tmp_typ[typ_i]);
3650          typ_i++;
3651          rO_LexVars(j, j_bits, r->block0[i],r->block1[i]-1, prev_ordsgn,
3652                     tmp_ordsgn,v, bits, r->block1[i]);
3653        }
3654        rCheckOrdSgn(r,i);
3655        break;
3656
3657      case ringorder_wp:
3658        rO_WDegree(j,j_bits,r->block0[i],r->block1[i],tmp_ordsgn,
3659                   tmp_typ[typ_i], r->wvhdl[i]);
3660        typ_i++;
3661        { // check for weights <=0
3662          int jj;
3663          BOOLEAN have_bad_weights=FALSE;
3664          for(jj=r->block1[i]-r->block0[i];jj>=0; jj--)
3665          {
3666            if (r->wvhdl[i][jj]<=0) have_bad_weights=TRUE;
3667          }
3668          if (have_bad_weights)
3669          {
3670             rO_TDegree(j,j_bits,r->block0[i],r->block1[i],tmp_ordsgn,
3671                                     tmp_typ[typ_i]);
3672             typ_i++;
3673             rCheckOrdSgn(r,i);
3674          }
3675        }
3676        if (r->block1[i]!=r->block0[i])
3677        {
3678          rO_LexVars_neg(j, j_bits,r->block1[i],r->block0[i]+1, prev_ordsgn,
3679                         tmp_ordsgn, v,bits, r->block0[i]);
3680        }
3681        break;
3682
3683      case ringorder_Wp:
3684        rO_WDegree(j,j_bits,r->block0[i],r->block1[i],tmp_ordsgn,
3685                   tmp_typ[typ_i], r->wvhdl[i]);
3686        typ_i++;
3687        { // check for weights <=0
3688          int jj;
3689          BOOLEAN have_bad_weights=FALSE;
3690          for(jj=r->block1[i]-r->block0[i];jj>=0; jj--)
3691          {
3692            if (r->wvhdl[i][jj]<=0) have_bad_weights=TRUE;
3693          }
3694          if (have_bad_weights)
3695          {
3696             rO_TDegree(j,j_bits,r->block0[i],r->block1[i],tmp_ordsgn,
3697                                     tmp_typ[typ_i]);
3698             typ_i++;
3699             rCheckOrdSgn(r,i);
3700          }
3701        }
3702        if (r->block1[i]!=r->block0[i])
3703        {
3704          rO_LexVars(j, j_bits,r->block0[i],r->block1[i]-1, prev_ordsgn,
3705                     tmp_ordsgn,v, bits, r->block1[i]);
3706        }
3707        break;
3708
3709      case ringorder_ws:
3710        rO_WDegree_neg(j,j_bits,r->block0[i],r->block1[i],tmp_ordsgn,
3711                       tmp_typ[typ_i], r->wvhdl[i]);
3712        typ_i++;
3713        if (r->block1[i]!=r->block0[i])
3714        {
3715          rO_LexVars_neg(j, j_bits,r->block1[i],r->block0[i]+1, prev_ordsgn,
3716                         tmp_ordsgn, v,bits, r->block0[i]);
3717        }
3718        rCheckOrdSgn(r,i);
3719        break;
3720
3721      case ringorder_Ws:
3722        rO_WDegree_neg(j,j_bits,r->block0[i],r->block1[i],tmp_ordsgn,
3723                       tmp_typ[typ_i], r->wvhdl[i]);
3724        typ_i++;
3725        if (r->block1[i]!=r->block0[i])
3726        {
3727          rO_LexVars(j, j_bits,r->block0[i],r->block1[i]-1, prev_ordsgn,
3728                     tmp_ordsgn,v, bits, r->block1[i]);
3729        }
3730        rCheckOrdSgn(r,i);
3731        break;
3732
3733      case ringorder_S:
3734        assume(typ_i == 1); // For LaScala3 only: on the 2nd place ([1])!
3735        // TODO: for K[x]: it is 0...?!
3736        rO_Syzcomp(j, j_bits,prev_ordsgn, tmp_ordsgn,tmp_typ[typ_i]);
3737        need_to_add_comp=TRUE;
3738        typ_i++;
3739        break;
3740
3741      case ringorder_s:
3742        assume(typ_i == 0 && j == 0);
3743        rO_Syz(j, j_bits, prev_ordsgn, tmp_ordsgn, tmp_typ[typ_i]); // set syz-limit?
3744        need_to_add_comp=TRUE;
3745        typ_i++;
3746        break;
3747
3748      case ringorder_IS:
3749      {
3750
3751        assume( r->block0[i] == r->block1[i] );
3752        const int s = r->block0[i];
3753        assume( -2 < s && s < 2);
3754
3755        if(s == 0) // Prefix IS
3756          rO_ISPrefix(j, j_bits, prev_ordsgn, tmp_ordsgn, r->N, v, tmp_typ[typ_i++]); // What about prev_ordsgn?
3757        else // s = +1 or -1 // Note: typ_i might be incrimented here inside!
3758        {
3759          rO_ISSuffix(j, j_bits, prev_ordsgn, tmp_ordsgn, r->N, v, tmp_typ, typ_i, s); // Suffix.
3760          need_to_add_comp=FALSE;
3761        }
3762
3763        break;
3764      }
3765      case ringorder_unspec:
3766      case ringorder_no:
3767      default:
3768        dReportError("undef. ringorder used\n");
3769        break;
3770    }
3771  }
3772
3773  int j0=j; // save j
3774  int j_bits0=j_bits; // save jbits
3775  rO_Align(j,j_bits);
3776  r->CmpL_Size = j;
3777
3778  j_bits=j_bits0; j=j0;
3779
3780  // fill in some empty slots with variables not already covered
3781  // v0 is special, is therefore normally already covered
3782  // now we do have rings without comp...
3783  if((need_to_add_comp) && (v[0]== -1))
3784  {
3785    if (prev_ordsgn==1)
3786    {
3787      rO_Align(j, j_bits);
3788      rO_LexVars(j, j_bits, 0,0, prev_ordsgn,tmp_ordsgn,v,BITS_PER_LONG, -1);
3789    }
3790    else
3791    {
3792      rO_Align(j, j_bits);
3793      rO_LexVars_neg(j, j_bits, 0,0, prev_ordsgn,tmp_ordsgn,v,BITS_PER_LONG, -1);
3794    }
3795  }
3796  // the variables
3797  for(i=1 ; i<=r->N ; i++)
3798  {
3799    if(v[i]==(-1))
3800    {
3801      if (prev_ordsgn==1)
3802      {
3803        rO_LexVars(j, j_bits, i,i, prev_ordsgn,tmp_ordsgn,v,bits, -1);
3804      }
3805      else
3806      {
3807        rO_LexVars_neg(j,j_bits,i,i, prev_ordsgn,tmp_ordsgn,v,bits, -1);
3808      }
3809    }
3810  }
3811
3812  rO_Align(j,j_bits);
3813  // ----------------------------
3814  // finished with constructing the monomial, computing sizes:
3815
3816  r->ExpL_Size=j;
3817  r->PolyBin = omGetSpecBin(POLYSIZE + (r->ExpL_Size)*sizeof(long));
3818  assume(r->PolyBin != NULL);
3819
3820  // ----------------------------
3821  // indices and ordsgn vector for comparison
3822  //
3823  // r->pCompHighIndex already set
3824  r->ordsgn=(long *)omAlloc0(r->ExpL_Size*sizeof(long));
3825
3826  for(j=0;j<r->CmpL_Size;j++)
3827  {
3828    r->ordsgn[j] = tmp_ordsgn[j];
3829  }
3830
3831  omFreeSize((ADDRESS)tmp_ordsgn,(3*(n+r->N)*sizeof(long)));
3832
3833  // ----------------------------
3834  // description of orderings for setm:
3835  //
3836  r->OrdSize=typ_i;
3837  if (typ_i==0) r->typ=NULL;
3838  else
3839  {
3840    r->typ=(sro_ord*)omAlloc(typ_i*sizeof(sro_ord));
3841    memcpy(r->typ,tmp_typ,typ_i*sizeof(sro_ord));
3842  }
3843  omFreeSize((ADDRESS)tmp_typ,(3*(n+r->N)*sizeof(sro_ord)));
3844
3845  // ----------------------------
3846  // indices for (first copy of ) variable entries in exp.e vector (VarOffset):
3847  r->VarOffset=v;
3848
3849  // ----------------------------
3850  // other indicies
3851  r->pCompIndex=(r->VarOffset[0] & 0xffff); //r->VarOffset[0];
3852  i=0; // position
3853  j=0; // index in r->typ
3854  if (i==r->pCompIndex) i++; // IS???
3855  while ((j < r->OrdSize)
3856         && ((r->typ[j].ord_typ==ro_syzcomp) ||
3857             (r->typ[j].ord_typ==ro_syz) || (r->typ[j].ord_typ==ro_isTemp) || (r->typ[j].ord_typ==ro_is) ||
3858             (r->order[r->typ[j].order_index] == ringorder_aa)))
3859  {
3860    i++; j++;
3861  }
3862  // No use of j anymore!!!????
3863
3864  if (i==r->pCompIndex) i++;
3865  r->pOrdIndex=i; // How came it is "i" here???!!!! exp[r->pOrdIndex] is order of a poly... This may be wrong!!! IS
3866
3867  // ----------------------------
3868  rSetDegStuff(r);
3869  rSetOption(r);
3870  // ----------------------------
3871  // r->p_Setm
3872  r->p_Setm = p_GetSetmProc(r);
3873
3874  // ----------------------------
3875  // set VarL_*
3876  rSetVarL(r);
3877
3878  //  ----------------------------
3879  // right-adjust VarOffset
3880  rRightAdjustVarOffset(r);
3881
3882  // ----------------------------
3883  // set NegWeightL*
3884  rSetNegWeight(r);
3885
3886  // ----------------------------
3887  // p_Procs: call AFTER NegWeightL
3888  r->p_Procs = (p_Procs_s*)omAlloc(sizeof(p_Procs_s));
3889  p_ProcsSet(r, r->p_Procs);
3890  return FALSE;
3891}
3892
3893static void rCheckOrdSgn(ring r,int i/*current block*/)
3894{ // set r->OrdSgn
3895  if ( r->OrdSgn==1)
3896  {
3897    int oo=-1;
3898    int jj;
3899    for(jj=i-1;jj>=0;jj--)
3900    {
3901      if(((r->order[jj]==ringorder_a)
3902        ||(r->order[jj]==ringorder_aa)
3903        ||(r->order[jj]==ringorder_a64))
3904      &&(r->block0[jj]<=r->block0[i])
3905      &&(r->block1[jj]>=r->block1[i]))
3906      { oo=1; break;}
3907    }
3908    r->OrdSgn=oo;
3909  }
3910}
3911
3912
3913void rUnComplete(ring r)
3914{
3915  if (r == NULL) return;
3916  if (r->VarOffset != NULL)
3917  {
3918    if (r->OrdSize!=0 && r->typ != NULL)
3919    {
3920      for(int i = 0; i < r->OrdSize; i++)
3921        if( r->typ[i].ord_typ == ro_is) // Search for suffixes! (prefix have the same VarOffset)
3922        {
3923          id_Delete(&r->typ[i].data.is.F, r);
3924          r->typ[i].data.is.F = NULL; // ?
3925
3926          if( r->typ[i].data.is.componentWeights != NULL )
3927          {
3928            delete r->typ[i].data.is.componentWeights;
3929            r->typ[i].data.is.componentWeights = NULL; // ?
3930          }
3931
3932          if( r->typ[i].data.is.pVarOffset != NULL )
3933          {
3934            omFreeSize((ADDRESS)r->typ[i].data.is.pVarOffset, (r->N +1)*sizeof(int));
3935            r->typ[i].data.is.pVarOffset = NULL; // ?
3936          }
3937        }
3938        else if (r->typ[i].ord_typ == ro_syz)
3939        {
3940          if(r->typ[i].data.syz.limit > 0)
3941            omFreeSize(r->typ[i].data.syz.syz_index, ((r->typ[i].data.syz.limit) +1)*sizeof(int));
3942          r->typ[i].data.syz.syz_index = NULL;
3943        }
3944        else if (r->typ[i].ord_typ == ro_syzcomp)
3945        {
3946          assume( r->typ[i].data.syzcomp.ShiftedComponents == NULL );
3947          assume( r->typ[i].data.syzcomp.Components        == NULL );
3948//          WarnS( "rUnComplete : ord_typ == ro_syzcomp was unhandled!!! Possibly memory leak!!!"  );
3949#ifndef NDEBUG
3950//          assume(0);
3951#endif
3952        }
3953
3954      omFreeSize((ADDRESS)r->typ,r->OrdSize*sizeof(sro_ord)); r->typ = NULL;
3955    }
3956
3957    if (r->order != NULL)
3958    {
3959      // delete r->order!!!???
3960    }
3961
3962    if (r->PolyBin != NULL)
3963      omUnGetSpecBin(&(r->PolyBin));
3964
3965    omFreeSize((ADDRESS)r->VarOffset, (r->N +1)*sizeof(int));
3966
3967    if (r->ordsgn != NULL && r->CmpL_Size != 0)
3968      omFreeSize((ADDRESS)r->ordsgn,r->ExpL_Size*sizeof(long));
3969    if (r->p_Procs != NULL)
3970      omFreeSize(r->p_Procs, sizeof(p_Procs_s));
3971    omfreeSize(r->VarL_Offset, r->VarL_Size*sizeof(int));
3972  }
3973  if (r->NegWeightL_Offset!=NULL)
3974  {
3975    omFreeSize(r->NegWeightL_Offset, r->NegWeightL_Size*sizeof(int));
3976    r->NegWeightL_Offset=NULL;
3977  }
3978}
3979
3980// set r->VarL_Size, r->VarL_Offset, r->VarL_LowIndex
3981static void rSetVarL(ring r)
3982{
3983  int  min = INT_MAX, min_j = -1;
3984  int* VarL_Number = (int*) omAlloc0(r->ExpL_Size*sizeof(int));
3985
3986  int i,j;
3987
3988  // count how often a var long is occupied by an exponent
3989  for (i=1; i<=r->N; i++)
3990  {
3991    VarL_Number[r->VarOffset[i] & 0xffffff]++;
3992  }
3993
3994  // determine how many and min
3995  for (i=0, j=0; i<r->ExpL_Size; i++)
3996  {
3997    if (VarL_Number[i] != 0)
3998    {
3999      if (min > VarL_Number[i])
4000      {
4001        min = VarL_Number[i];
4002        min_j = j;
4003      }
4004      j++;
4005    }
4006  }
4007
4008  r->VarL_Size = j; // number of long with exp. entries in
4009                    //  in p->exp
4010  r->VarL_Offset = (int*) omAlloc(r->VarL_Size*sizeof(int));
4011  r->VarL_LowIndex = 0;
4012
4013  // set VarL_Offset
4014  for (i=0, j=0; i<r->ExpL_Size; i++)
4015  {
4016    if (VarL_Number[i] != 0)
4017    {
4018      r->VarL_Offset[j] = i;
4019      if (j > 0 && r->VarL_Offset[j-1] != r->VarL_Offset[j] - 1)
4020        r->VarL_LowIndex = -1;
4021      j++;
4022    }
4023  }
4024  if (r->VarL_LowIndex >= 0)
4025    r->VarL_LowIndex = r->VarL_Offset[0];
4026
4027  r->MinExpPerLong = min;
4028  if (min_j != 0)
4029  {
4030    j = r->VarL_Offset[min_j];
4031    r->VarL_Offset[min_j] = r->VarL_Offset[0];
4032    r->VarL_Offset[0] = j;
4033  }
4034  omFree(VarL_Number);
4035}
4036
4037static void rRightAdjustVarOffset(ring r)
4038{
4039  int* shifts = (int*) omAlloc(r->ExpL_Size*sizeof(int));
4040  int i;
4041  // initialize shifts
4042  for (i=0;i<r->ExpL_Size;i++)
4043    shifts[i] = BIT_SIZEOF_LONG;
4044
4045  // find minimal bit shift in each long exp entry
4046  for (i=1;i<=r->N;i++)
4047  {
4048    if (shifts[r->VarOffset[i] & 0xffffff] > r->VarOffset[i] >> 24)
4049      shifts[r->VarOffset[i] & 0xffffff] = r->VarOffset[i] >> 24;
4050  }
4051  // reset r->VarOffset: set the minimal shift to 0
4052  for (i=1;i<=r->N;i++)
4053  {
4054    if (shifts[r->VarOffset[i] & 0xffffff] != 0)
4055      r->VarOffset[i]
4056        = (r->VarOffset[i] & 0xffffff) |
4057        (((r->VarOffset[i] >> 24) - shifts[r->VarOffset[i] & 0xffffff]) << 24);
4058  }
4059  omFree(shifts);
4060}
4061
4062// get r->divmask depending on bits per exponent
4063static unsigned long rGetDivMask(int bits)
4064{
4065  unsigned long divmask = 1;
4066  int i = bits;
4067
4068  while (i < BIT_SIZEOF_LONG)
4069  {
4070    divmask |= (((unsigned long) 1) << (unsigned long) i);
4071    i += bits;
4072  }
4073  return divmask;
4074}
4075
4076#ifdef RDEBUG
4077void rDebugPrint(ring r)
4078{
4079  if (r==NULL)
4080  {
4081    PrintS("NULL ?\n");
4082    return;
4083  }
4084  // corresponds to ro_typ from ring.h:
4085  const char *TYP[]={"ro_dp","ro_wp","ro_wp64","ro_wp_neg","ro_cp",
4086                     "ro_syzcomp", "ro_syz", "ro_isTemp", "ro_is", "ro_none"};
4087  int i,j;
4088
4089  Print("ExpL_Size:%d ",r->ExpL_Size);
4090  Print("CmpL_Size:%d ",r->CmpL_Size);
4091  Print("VarL_Size:%d\n",r->VarL_Size);
4092  Print("bitmask=0x%lx (expbound=%ld) \n",r->bitmask, r->bitmask);
4093  Print("BitsPerExp=%d ExpPerLong=%d MinExpPerLong=%d at L[%d]\n", r->BitsPerExp, r->ExpPerLong, r->MinExpPerLong, r->VarL_Offset[0]);
4094  PrintS("varoffset:\n");
4095  if (r->VarOffset==NULL) PrintS(" NULL\n");
4096  else
4097    for(j=0;j<=r->N;j++)
4098      Print("  v%d at e-pos %d, bit %d\n",
4099            j,r->VarOffset[j] & 0xffffff, r->VarOffset[j] >>24);
4100  Print("divmask=%lx\n", r->divmask);
4101  PrintS("ordsgn:\n");
4102  for(j=0;j<r->CmpL_Size;j++)
4103    Print("  ordsgn %ld at pos %d\n",r->ordsgn[j],j);
4104  Print("OrdSgn:%d\n",r->OrdSgn);
4105  PrintS("ordrec:\n");
4106  for(j=0;j<r->OrdSize;j++)
4107  {
4108    Print("  typ %s", TYP[r->typ[j].ord_typ]);
4109
4110
4111    if (r->typ[j].ord_typ==ro_syz)
4112    {
4113      const short place = r->typ[j].data.syz.place;
4114      const int limit = r->typ[j].data.syz.limit;
4115      const int curr_index = r->typ[j].data.syz.curr_index;
4116      const int* syz_index = r->typ[j].data.syz.syz_index;
4117
4118      Print("  limit %d (place: %d, curr_index: %d), syz_index: ", limit, place, curr_index);
4119
4120      if( syz_index == NULL )
4121        PrintS("(NULL)");
4122      else
4123      {
4124        Print("{");
4125        for( i=0; i <= limit; i++ )
4126          Print("%d ", syz_index[i]);
4127        Print("}");
4128      }
4129
4130    }
4131    else if (r->typ[j].ord_typ==ro_isTemp)
4132    {
4133      Print("  start (level) %d, suffixpos: %d, VO: ",r->typ[j].data.isTemp.start, r->typ[j].data.isTemp.suffixpos);
4134
4135    }
4136    else if (r->typ[j].ord_typ==ro_is)
4137    {
4138      Print("  start %d, end: %d: ",r->typ[j].data.is.start, r->typ[j].data.is.end);
4139
4140//      for( int k = 0; k <= r->N; k++) if (r->typ[j].data.is.pVarOffset[k] != -1) Print("[%2d]: %04x; ", k, r->typ[j].data.is.pVarOffset[k]);
4141
4142      Print("  limit %d\n",r->typ[j].data.is.limit);
4143      #ifndef NDEBUG
4144      //PrintS("  F: ");idShow(r->typ[j].data.is.F, r, r, 1);
4145      #endif
4146
4147      PrintS("weights: ");
4148
4149      if( r->typ[j].data.is.componentWeights == NULL )
4150        PrintS("NULL == [0,...0]\n");
4151      else
4152      {
4153        (r->typ[j].data.is.componentWeights)->show(); PrintLn();
4154      }
4155    }
4156    else
4157    {
4158      Print("  place %d",r->typ[j].data.dp.place);
4159
4160      if (r->typ[j].ord_typ!=ro_syzcomp  && r->typ[j].ord_typ!=ro_syz)
4161      {
4162        Print("  start %d",r->typ[j].data.dp.start);
4163        Print("  end %d",r->typ[j].data.dp.end);
4164        if ((r->typ[j].ord_typ==ro_wp)
4165        || (r->typ[j].ord_typ==ro_wp_neg))
4166        {
4167          PrintS(" w:");
4168          for(int l=r->typ[j].data.wp.start;l<=r->typ[j].data.wp.end;l++)
4169            Print(" %d",r->typ[j].data.wp.weights[l-r->typ[j].data.wp.start]);
4170        }
4171        else if (r->typ[j].ord_typ==ro_wp64)
4172        {
4173          PrintS(" w64:");
4174          int l;
4175          for(l=r->typ[j].data.wp64.start;l<=r->typ[j].data.wp64.end;l++)
4176            Print(" %ld",(long)(((int64*)r->typ[j].data.wp64.weights64)+l-r->typ[j].data.wp64.start));
4177          }
4178        }
4179    }
4180    PrintLn();
4181  }
4182  Print("pOrdIndex:%d pCompIndex:%d\n", r->pOrdIndex, r->pCompIndex);
4183  Print("OrdSize:%d\n",r->OrdSize);
4184  PrintS("--------------------\n");
4185  for(j=0;j<r->ExpL_Size;j++)
4186  {
4187    Print("L[%d]: ",j);
4188    if (j< r->CmpL_Size)
4189      Print("ordsgn %ld ", r->ordsgn[j]);
4190    else
4191      PrintS("no comp ");
4192    i=1;
4193    for(;i<=r->N;i++)
4194    {
4195      if( (r->VarOffset[i] & 0xffffff) == j )
4196      {  Print("v%d at e[%d], bit %d; ", i,r->VarOffset[i] & 0xffffff,
4197                                         r->VarOffset[i] >>24 ); }
4198    }
4199    if( r->pCompIndex==j ) PrintS("v0; ");
4200    for(i=0;i<r->OrdSize;i++)
4201    {
4202      if (r->typ[i].data.dp.place == j)
4203      {
4204        Print("ordrec:%s (start:%d, end:%d) ",TYP[r->typ[i].ord_typ],
4205          r->typ[i].data.dp.start, r->typ[i].data.dp.end);
4206      }
4207    }
4208
4209    if (j==r->pOrdIndex)
4210      PrintS("pOrdIndex\n");
4211    else
4212      PrintLn();
4213  }
4214  Print("LexOrder:%d, MixedOrder:%d\n",r->LexOrder, r->MixedOrder);
4215
4216  // p_Procs stuff
4217  p_Procs_s proc_names;
4218  const char* field;
4219  const char* length;
4220  const char* ord;
4221  p_Debug_GetProcNames(r, &proc_names); // changes p_Procs!!!
4222  p_Debug_GetSpecNames(r, field, length, ord);
4223
4224  Print("p_Spec  : %s, %s, %s\n", field, length, ord);
4225  PrintS("p_Procs :\n");
4226  for (i=0; i<(int) (sizeof(p_Procs_s)/sizeof(void*)); i++)
4227  {
4228    Print(" %s,\n", ((char**) &proc_names)[i]);
4229  }
4230
4231  {
4232#define pFDeg_CASE(A) if(r->pFDeg == A) PrintS( "" #A "" )
4233    Print("\npFDeg   : ");
4234   
4235    pFDeg_CASE(p_Totaldegree); else
4236      pFDeg_CASE(p_WFirstTotalDegree); else
4237      pFDeg_CASE(p_WTotaldegree); else
4238      pFDeg_CASE(p_Deg); else
4239      Print("(%p)", r->pFDeg); // default case
4240   
4241    PrintS("\n");
4242#undef pFDeg_CASE
4243  }
4244   
4245}
4246
4247void p_DebugPrint(poly p, const ring r)
4248{
4249  int i,j;
4250  p_Write(p,r);
4251  j=2;
4252  while(p!=NULL)
4253  {
4254    Print("\nexp[0..%d]\n",r->ExpL_Size-1);
4255    for(i=0;i<r->ExpL_Size;i++)
4256      Print("%ld ",p->exp[i]);
4257    PrintLn();
4258    Print("v0:%ld ",p_GetComp(p, r));
4259    for(i=1;i<=r->N;i++) Print(" v%d:%ld",i,p_GetExp(p,i, r));
4260    PrintLn();
4261    pIter(p);
4262    j--;
4263    if (j==0) { PrintS("...\n"); break; }
4264  }
4265}
4266
4267#endif // RDEBUG
4268
4269/// debug-print monomial poly/vector p, assuming that it lives in the ring R
4270static inline void m_DebugPrint(const poly p, const ring R)
4271{
4272  Print("\nexp[0..%d]\n", R->ExpL_Size - 1);
4273  for(int i = 0; i < R->ExpL_Size; i++)
4274    Print("%09lx ", p->exp[i]);
4275  PrintLn();
4276  Print("v0:%9ld ", p_GetComp(p, R));
4277  for(int i = 1; i <= R->N; i++) Print(" v%d:%5ld",i, p_GetExp(p, i, R));
4278  PrintLn();
4279}
4280
4281
4282#ifndef NDEBUG
4283/// debug-print at most nTerms (2 by default) terms from poly/vector p,
4284/// assuming that lt(p) lives in lmRing and tail(p) lives in tailRing.
4285void p_DebugPrint(const poly p, const ring lmRing, const ring tailRing, const int nTerms)
4286{
4287  assume( nTerms >= 0 );
4288  if( p != NULL )
4289  {
4290    assume( p != NULL );
4291
4292    p_Write(p, lmRing, tailRing);
4293
4294    if( (p != NULL) && (nTerms > 0) )
4295    {
4296      assume( p != NULL );
4297      assume( nTerms > 0 );
4298
4299      // debug pring leading term
4300      m_DebugPrint(p, lmRing);
4301
4302      poly q = pNext(p); // q = tail(p)
4303
4304      // debug pring tail (at most nTerms-1 terms from it)
4305      for(int j = nTerms - 1; (q !=NULL) && (j > 0); pIter(q), --j)
4306        m_DebugPrint(q, tailRing);
4307
4308      if (q != NULL)
4309        PrintS("...\n");
4310    }
4311  }
4312  else
4313    PrintS("0\n");
4314}
4315#endif
4316
4317
4318//    F = system("ISUpdateComponents", F, V, MIN );
4319//    // replace gen(i) -> gen(MIN + V[i-MIN]) for all i > MIN in all terms from F!
4320void pISUpdateComponents(ideal F, const intvec *const V, const int MIN, const ring r )
4321{
4322  assume( V != NULL );
4323  assume( MIN >= 0 );
4324
4325  if( F == NULL )
4326    return;
4327
4328  for( int j = (F->ncols*F->nrows) - 1; j >= 0; j-- )
4329  {
4330#ifdef PDEBUG
4331    Print("F[%d]:", j);
4332    p_DebugPrint(F->m[j], r, r, 0);
4333#endif
4334
4335    for( poly p = F->m[j]; p != NULL; pIter(p) )
4336    {
4337      int c = p_GetComp(p, r);
4338
4339      if( c > MIN )
4340      {
4341#ifdef PDEBUG
4342        Print("gen[%d] -> gen(%d)\n", c, MIN + (*V)[ c - MIN - 1 ]);
4343#endif
4344
4345        p_SetComp( p, MIN + (*V)[ c - MIN - 1 ], r );
4346      }
4347    }
4348#ifdef PDEBUG
4349    Print("new F[%d]:", j);
4350    p_Test(F->m[j], r);
4351    p_DebugPrint(F->m[j], r, r, 0);
4352#endif
4353  }
4354
4355}
4356
4357
4358
4359
4360/*2
4361* asssume that rComplete was called with r
4362* assume that the first block ist ringorder_S
4363* change the block to reflect the sequence given by appending v
4364*/
4365
4366#ifdef PDEBUG
4367void rDBChangeSComps(int* currComponents,
4368                     long* currShiftedComponents,
4369                     int length,
4370                     ring r)
4371{
4372  assume(r->typ[1].ord_typ == ro_syzcomp);
4373     
4374  r->typ[1].data.syzcomp.length = length;
4375  rNChangeSComps( currComponents, currShiftedComponents, r);
4376}
4377void rDBGetSComps(int** currComponents,
4378                 long** currShiftedComponents,
4379                 int *length,
4380                 ring r)
4381{
4382  assume(r->typ[1].ord_typ == ro_syzcomp);
4383 
4384  *length = r->typ[1].data.syzcomp.length;
4385  rNGetSComps( currComponents, currShiftedComponents, r);
4386}
4387#endif
4388
4389void rNChangeSComps(int* currComponents, long* currShiftedComponents, ring r)
4390{
4391  assume(r->typ[1].ord_typ == ro_syzcomp);
4392
4393  r->typ[1].data.syzcomp.ShiftedComponents = currShiftedComponents;
4394  r->typ[1].data.syzcomp.Components = currComponents;
4395}
4396
4397void rNGetSComps(int** currComponents, long** currShiftedComponents, ring r)
4398{
4399  assume(r->typ[1].ord_typ == ro_syzcomp);
4400
4401  *currShiftedComponents = r->typ[1].data.syzcomp.ShiftedComponents;
4402  *currComponents =   r->typ[1].data.syzcomp.Components;
4403}
4404
4405/////////////////////////////////////////////////////////////////////////////
4406//
4407// The following routines all take as input a ring r, and return R
4408// where R has a certain property. R might be equal r in which case r
4409// had already this property
4410//
4411ring rAssure_SyzComp(const ring r, BOOLEAN complete)
4412{
4413  if ( (r->order[0] == ringorder_s) ) return r;
4414
4415  if ( (r->order[0] == ringorder_IS) )
4416  {
4417#ifndef NDEBUG
4418    WarnS("rAssure_SyzComp: input ring has an IS-ordering!");
4419#endif
4420//    return r;
4421  }
4422  ring res=rCopy0(r, FALSE, FALSE);
4423  int i=rBlocks(r);
4424  int j;
4425
4426  res->order=(int *)omAlloc((i+1)*sizeof(int));
4427  res->block0=(int *)omAlloc0((i+1)*sizeof(int));
4428  res->block1=(int *)omAlloc0((i+1)*sizeof(int));
4429  int ** wvhdl =(int **)omAlloc0((i+1)*sizeof(int**));
4430  for(j=i;j>0;j--)
4431  {
4432    res->order[j]=r->order[j-1];
4433    res->block0[j]=r->block0[j-1];
4434    res->block1[j]=r->block1[j-1];
4435    if (r->wvhdl[j-1] != NULL)
4436    {
4437      wvhdl[j] = (int*) omMemDup(r->wvhdl[j-1]);
4438    }
4439  }
4440  res->order[0]=ringorder_s;
4441
4442  res->wvhdl = wvhdl;
4443
4444  if (complete)
4445  {
4446    rComplete(res, 1);
4447
4448#ifdef HAVE_PLURAL
4449    if (rIsPluralRing(r))
4450    {
4451      if ( nc_rComplete(r, res, false) ) // no qideal!
4452      {
4453#ifndef NDEBUG
4454        WarnS("error in nc_rComplete");      // cleanup?//      rDelete(res);//      return r;      // just go on..
4455#endif
4456      }
4457    }
4458    assume(rIsPluralRing(r) == rIsPluralRing(res));
4459#endif
4460
4461
4462#ifdef HAVE_PLURAL
4463    ring old_ring = r;
4464#endif
4465
4466    if (r->qideal!=NULL)
4467    {
4468      res->qideal= idrCopyR_NoSort(r->qideal, r, res);
4469
4470      assume(id_RankFreeModule(res->qideal, res) == 0);
4471
4472#ifdef HAVE_PLURAL
4473      if( rIsPluralRing(res) )
4474        if( nc_SetupQuotient(res, r, true) )
4475        {
4476//          WarnS("error in nc_SetupQuotient"); // cleanup?      rDelete(res);       return r;  // just go on...?
4477        }
4478
4479#endif
4480      assume(id_RankFreeModule(res->qideal, res) == 0);
4481    }
4482
4483#ifdef HAVE_PLURAL
4484    assume((res->qideal==NULL) == (old_ring->qideal==NULL));
4485    assume(rIsPluralRing(res) == rIsPluralRing(old_ring));
4486    assume(rIsSCA(res) == rIsSCA(old_ring));
4487    assume(ncRingType(res) == ncRingType(old_ring));
4488#endif
4489  }
4490
4491  return res;
4492}
4493
4494ring rAssure_TDeg(ring r, int start_var, int end_var, int &pos)
4495{
4496  int i;
4497  if (r->typ!=NULL)
4498  {
4499    for(i=r->OrdSize-1;i>=0;i--)
4500    {
4501      if ((r->typ[i].ord_typ==ro_dp)
4502      && (r->typ[i].data.dp.start==start_var)
4503      && (r->typ[i].data.dp.end==end_var))
4504      {
4505        pos=r->typ[i].data.dp.place;
4506        //printf("no change, pos=%d\n",pos);
4507        return r;
4508      }
4509    }
4510  }
4511
4512#ifdef HAVE_PLURAL
4513  nc_struct* save=r->GetNC();
4514  r->GetNC()=NULL;
4515#endif
4516  ring res=rCopy(r);
4517
4518  i=rBlocks(r);
4519  int j;
4520
4521  res->ExpL_Size=r->ExpL_Size+1; // one word more in each monom
4522  res->PolyBin=omGetSpecBin(POLYSIZE + (res->ExpL_Size)*sizeof(long));
4523  omFree((ADDRESS)res->ordsgn);
4524  res->ordsgn=(long *)omAlloc0(res->ExpL_Size*sizeof(long));
4525  for(j=0;j<r->CmpL_Size;j++)
4526  {
4527    res->ordsgn[j] = r->ordsgn[j];
4528  }
4529  res->OrdSize=r->OrdSize+1;   // one block more for pSetm
4530  if (r->typ!=NULL)
4531    omFree((ADDRESS)res->typ);
4532  res->typ=(sro_ord*)omAlloc0(res->OrdSize*sizeof(sro_ord));
4533  if (r->typ!=NULL)
4534    memcpy(res->typ,r->typ,r->OrdSize*sizeof(sro_ord));
4535  // the additional block for pSetm: total degree at the last word
4536  // but not included in the compare part
4537  res->typ[res->OrdSize-1].ord_typ=ro_dp;
4538  res->typ[res->OrdSize-1].data.dp.start=start_var;
4539  res->typ[res->OrdSize-1].data.dp.end=end_var;
4540  res->typ[res->OrdSize-1].data.dp.place=res->ExpL_Size-1;
4541  pos=res->ExpL_Size-1;
4542  //if ((start_var==1) && (end_var==res->N)) res->pOrdIndex=pos;
4543  extern void p_Setm_General(poly p, ring r);
4544  res->p_Setm=p_Setm_General;
4545  // ----------------------------
4546  omFree((ADDRESS)res->p_Procs);
4547  res->p_Procs = (p_Procs_s*)omAlloc(sizeof(p_Procs_s));
4548
4549  p_ProcsSet(res, res->p_Procs);
4550  if (res->qideal!=NULL) id_Delete(&res->qideal,res);
4551#ifdef HAVE_PLURAL
4552  r->GetNC()=save;
4553  if (rIsPluralRing(r))
4554  {
4555    if ( nc_rComplete(r, res, false) ) // no qideal!
4556    {
4557#ifndef NDEBUG
4558      WarnS("error in nc_rComplete");
4559#endif
4560      // just go on..
4561    }
4562  }
4563#endif
4564  if (r->qideal!=NULL)
4565  {
4566     res->qideal=idrCopyR_NoSort(r->qideal,r, res);
4567#ifdef HAVE_PLURAL
4568     if (rIsPluralRing(res))
4569     {
4570//       nc_SetupQuotient(res, currRing);
4571       nc_SetupQuotient(res, r); // ?
4572     }
4573     assume((res->qideal==NULL) == (r->qideal==NULL));
4574#endif
4575  }
4576
4577#ifdef HAVE_PLURAL
4578  assume(rIsPluralRing(res) == rIsPluralRing(r));
4579  assume(rIsSCA(res) == rIsSCA(r));
4580  assume(ncRingType(res) == ncRingType(r));
4581#endif
4582
4583  return res;
4584}
4585
4586ring rAssure_HasComp(const ring r)
4587{
4588  int last_block;
4589  int i=0;
4590  do
4591  {
4592     if (r->order[i] == ringorder_c ||
4593        r->order[i] == ringorder_C) return r;
4594     if (r->order[i] == 0)
4595        break;
4596     i++;
4597  } while (1);
4598  //WarnS("re-creating ring with comps");
4599  last_block=i-1;
4600
4601  ring new_r = rCopy0(r, FALSE, FALSE);
4602  i+=2;
4603  new_r->wvhdl=(int **)omAlloc0(i * sizeof(int *));
4604  new_r->order   = (int *) omAlloc0(i * sizeof(int));
4605  new_r->block0   = (int *) omAlloc0(i * sizeof(int));
4606  new_r->block1   = (int *) omAlloc0(i * sizeof(int));
4607  memcpy(new_r->order,r->order,(i-1) * sizeof(int));
4608  memcpy(new_r->block0,r->block0,(i-1) * sizeof(int));
4609  memcpy(new_r->block1,r->block1,(i-1) * sizeof(int));
4610  for (int j=0; j<=last_block; j++)
4611  {
4612    if (r->wvhdl[j]!=NULL)
4613    {
4614      new_r->wvhdl[j] = (int*) omMemDup(r->wvhdl[j]);
4615    }
4616  }
4617  last_block++;
4618  new_r->order[last_block]=ringorder_C;
4619  //new_r->block0[last_block]=0;
4620  //new_r->block1[last_block]=0;
4621  //new_r->wvhdl[last_block]=NULL;
4622
4623  rComplete(new_r, 1);
4624
4625#ifdef HAVE_PLURAL
4626  if (rIsPluralRing(r))
4627  {
4628    if ( nc_rComplete(r, new_r, false) ) // no qideal!
4629    {
4630#ifndef NDEBUG
4631      WarnS("error in nc_rComplete");      // cleanup?//      rDelete(res);//      return r;      // just go on..
4632#endif
4633    }
4634  }
4635  assume(rIsPluralRing(r) == rIsPluralRing(new_r));
4636#endif
4637
4638  return new_r;
4639}
4640
4641static ring rAssure_CompLastBlock(ring r, BOOLEAN complete = TRUE)
4642{
4643  int last_block = rBlocks(r) - 2;
4644  if (r->order[last_block] != ringorder_c &&
4645      r->order[last_block] != ringorder_C)
4646  {
4647    int c_pos = 0;
4648    int i;
4649
4650    for (i=0; i< last_block; i++)
4651    {
4652      if (r->order[i] == ringorder_c || r->order[i] == ringorder_C)
4653      {
4654        c_pos = i;
4655        break;
4656      }
4657    }
4658    if (c_pos != -1)
4659    {
4660      ring new_r = rCopy0(r, FALSE, TRUE);
4661      for (i=c_pos+1; i<=last_block; i++)
4662      {
4663        new_r->order[i-1] = new_r->order[i];
4664        new_r->block0[i-1] = new_r->block0[i];
4665        new_r->block1[i-1] = new_r->block1[i];
4666        new_r->wvhdl[i-1] = new_r->wvhdl[i];
4667      }
4668      new_r->order[last_block] = r->order[c_pos];
4669      new_r->block0[last_block] = r->block0[c_pos];
4670      new_r->block1[last_block] = r->block1[c_pos];
4671      new_r->wvhdl[last_block] = r->wvhdl[c_pos];
4672      if (complete)
4673      {
4674        rComplete(new_r, 1);
4675
4676#ifdef HAVE_PLURAL
4677        if (rIsPluralRing(r))
4678        {
4679          if ( nc_rComplete(r, new_r, false) ) // no qideal!
4680          {
4681#ifndef NDEBUG
4682            WarnS("error in nc_rComplete");   // cleanup?//      rDelete(res);//      return r;      // just go on..
4683#endif
4684          }
4685        }
4686        assume(rIsPluralRing(r) == rIsPluralRing(new_r));
4687#endif
4688      }
4689      return new_r;
4690    }
4691  }
4692  return r;
4693}
4694
4695// Moves _c or _C ordering to the last place AND adds _s on the 1st place
4696ring rAssure_SyzComp_CompLastBlock(const ring r)
4697{
4698  ring new_r_1 = rAssure_CompLastBlock(r, FALSE); // due to this FALSE - no completion!
4699  ring new_r = rAssure_SyzComp(new_r_1, FALSE); // new_r_1 is used only here!!!
4700
4701  if (new_r != r)
4702  {
4703    ring old_r = r;
4704    if (new_r_1 != new_r && new_r_1 != old_r) rDelete(new_r_1);
4705    rComplete(new_r, 1);
4706#ifdef HAVE_PLURAL
4707    if (rIsPluralRing(old_r))
4708    {
4709      if ( nc_rComplete(old_r, new_r, false) ) // no qideal!
4710      {
4711#ifndef NDEBUG
4712        WarnS("error in nc_rComplete"); // cleanup?      rDelete(res);       return r;  // just go on...?
4713#endif
4714        }
4715    }
4716    assume(rIsPluralRing(new_r) == rIsPluralRing(old_r));
4717#endif
4718///?    rChangeCurrRing(new_r);
4719    if (old_r->qideal != NULL)
4720    {
4721      new_r->qideal = idrCopyR(old_r->qideal, old_r, new_r);
4722      //currQuotient = new_r->qideal;
4723
4724#ifdef HAVE_PLURAL
4725      if( rIsPluralRing(old_r) )
4726        if( nc_SetupQuotient(new_r, old_r, true) )
4727        {
4728#ifndef NDEBUG
4729          WarnS("error in nc_SetupQuotient"); // cleanup?      rDelete(res);       return r;  // just go on...?
4730#endif
4731        }
4732#endif
4733    }
4734
4735#ifdef HAVE_PLURAL
4736    assume((new_r->qideal==NULL) == (old_r->qideal==NULL));
4737    assume(rIsPluralRing(new_r) == rIsPluralRing(old_r));
4738    assume(rIsSCA(new_r) == rIsSCA(old_r));
4739    assume(ncRingType(new_r) == ncRingType(old_r));
4740#endif
4741
4742    rTest(new_r);
4743    rTest(old_r);
4744  }
4745  return new_r;
4746}
4747
4748// use this for global orderings consisting of two blocks
4749static ring rAssure_Global(rRingOrder_t b1, rRingOrder_t b2, const ring r)
4750{
4751  int r_blocks = rBlocks(r);
4752
4753  assume(b1 == ringorder_c || b1 == ringorder_C ||
4754         b2 == ringorder_c || b2 == ringorder_C ||
4755         b2 == ringorder_S);
4756  if ((r_blocks == 3) &&
4757      (r->order[0] == b1) &&
4758      (r->order[1] == b2) &&
4759      (r->order[2] == 0))
4760    return r;
4761  ring res = rCopy0(r, TRUE, FALSE);
4762  res->order = (int*)omAlloc0(3*sizeof(int));
4763  res->block0 = (int*)omAlloc0(3*sizeof(int));
4764  res->block1 = (int*)omAlloc0(3*sizeof(int));
4765  res->wvhdl = (int**)omAlloc0(3*sizeof(int*));
4766  res->order[0] = b1;
4767  res->order[1] = b2;
4768  if (b1 == ringorder_c || b1 == ringorder_C)
4769  {
4770    res->block0[1] = 1;
4771    res->block1[1] = r->N;
4772  }
4773  else
4774  {
4775    res->block0[0] = 1;
4776    res->block1[0] = r->N;
4777  }
4778  // HANNES: This sould be set in rComplete
4779  res->OrdSgn = 1;
4780  rComplete(res, 1);
4781#ifdef HAVE_PLURAL
4782  if (rIsPluralRing(r))
4783  {
4784    if ( nc_rComplete(r, res, false) ) // no qideal!
4785    {
4786#ifndef NDEBUG
4787      WarnS("error in nc_rComplete");
4788#endif
4789    }
4790  }
4791#endif
4792//  rChangeCurrRing(res);
4793  return res;
4794}
4795
4796ring rAssure_InducedSchreyerOrdering(const ring r, BOOLEAN complete = TRUE, int sgn = 1)
4797{ // TODO: ???? Add leading Syz-comp ordering here...????
4798
4799#if MYTEST
4800    Print("rAssure_InducedSchreyerOrdering(r, complete = %d, sgn = %d): r: \n", complete, sgn);
4801    rWrite(r);
4802#ifdef RDEBUG
4803    rDebugPrint(r);
4804#endif
4805    PrintLn();
4806#endif
4807  assume((sgn == 1) || (sgn == -1));
4808
4809  ring res=rCopy0(r, FALSE, FALSE); // No qideal & ordering copy.
4810
4811  int n = rBlocks(r); // Including trailing zero!
4812
4813  // Create 2 more blocks for prefix/suffix:
4814  res->order=(int *)omAlloc0((n+2)*sizeof(int)); // 0  ..  n+1
4815  res->block0=(int *)omAlloc0((n+2)*sizeof(int));
4816  res->block1=(int *)omAlloc0((n+2)*sizeof(int));
4817  int ** wvhdl =(int **)omAlloc0((n+2)*sizeof(int**));
4818
4819  // Encapsulate all existing blocks between induced Schreyer ordering markers: prefix and suffix!
4820  // Note that prefix and suffix have the same ringorder marker and only differ in block[] parameters!
4821
4822  // new 1st block
4823  int j = 0;
4824  res->order[j] = ringorder_IS; // Prefix
4825  res->block0[j] = res->block1[j] = 0;
4826  // wvhdl[j] = NULL;
4827  j++;
4828
4829  for(int i = 0; (i <= n) && (r->order[i] != 0); i++, j++) // i = [0 .. n-1] <- non-zero old blocks
4830  {
4831    res->order [j] = r->order [i];
4832    res->block0[j] = r->block0[i];
4833    res->block1[j] = r->block1[i];
4834
4835    if (r->wvhdl[i] != NULL)
4836    {
4837      wvhdl[j] = (int*) omMemDup(r->wvhdl[i]);
4838    } // else wvhdl[j] = NULL;
4839  }
4840
4841  // new last block
4842  res->order [j] = ringorder_IS; // Suffix
4843  res->block0[j] = res->block1[j] = sgn; // Sign of v[o]: 1 for C, -1 for c
4844  // wvhdl[j] = NULL;
4845  j++;
4846
4847  // res->order [j] = 0; // The End!
4848  res->wvhdl = wvhdl;
4849
4850  // j == the last zero block now!
4851  assume(j == (n+1));
4852  assume(res->order[0]==ringorder_IS);
4853  assume(res->order[j-1]==ringorder_IS);
4854  assume(res->order[j]==0);
4855
4856
4857  if (complete)
4858  {
4859    rComplete(res, 1);
4860
4861#ifdef HAVE_PLURAL
4862    if (rIsPluralRing(r))
4863    {
4864      if ( nc_rComplete(r, res, false) ) // no qideal!
4865      {
4866#ifndef NDEBUG
4867        WarnS("error in nc_rComplete");      // cleanup?//      rDelete(res);//      return r;      // just go on..
4868#endif
4869      }
4870    }
4871    assume(rIsPluralRing(r) == rIsPluralRing(res));
4872#endif
4873
4874
4875#ifdef HAVE_PLURAL
4876    ring old_ring = r;
4877#endif
4878
4879    if (r->qideal!=NULL)
4880    {
4881      res->qideal= idrCopyR_NoSort(r->qideal, r, res);
4882
4883      assume(id_RankFreeModule(res->qideal, res) == 0);
4884
4885#ifdef HAVE_PLURAL
4886      if( rIsPluralRing(res) )
4887        if( nc_SetupQuotient(res, r, true) )
4888        {
4889//          WarnS("error in nc_SetupQuotient"); // cleanup?      rDelete(res);       return r;  // just go on...?
4890        }
4891
4892#endif
4893      assume(id_RankFreeModule(res->qideal, res) == 0);
4894    }
4895
4896#ifdef HAVE_PLURAL
4897    assume((res->qideal==NULL) == (old_ring->qideal==NULL));
4898    assume(rIsPluralRing(res) == rIsPluralRing(old_ring));
4899    assume(rIsSCA(res) == rIsSCA(old_ring));
4900    assume(ncRingType(res) == ncRingType(old_ring));
4901#endif
4902  }
4903
4904  return res;
4905}
4906
4907ring rAssure_dp_S(const ring r)
4908{
4909  return rAssure_Global(ringorder_dp, ringorder_S,r);
4910}
4911
4912ring rAssure_dp_C(const ring r)
4913{
4914  return rAssure_Global(ringorder_dp, ringorder_C,r);
4915}
4916
4917ring rAssure_C_dp(const ring r)
4918{
4919  return rAssure_Global(ringorder_C, ringorder_dp,r);
4920}
4921
4922
4923
4924/// Finds p^th IS ordering, and returns its position in r->typ[]
4925/// returns -1 if something went wrong!
4926int rGetISPos(const int p, const ring r)
4927{
4928  // Put the reference set F into the ring -ordering -recor
4929#if MYTEST
4930  Print("rIsIS(p: %d)\nF:", p);
4931  PrintLn();
4932#endif
4933
4934  if (r->typ==NULL)
4935  {
4936    dReportError("'rIsIS:' Error: wrong ring! (typ == NULL)");
4937    return -1;
4938  }
4939
4940  int j = p; // Which IS record to use...
4941  for( int pos = 0; pos < r->OrdSize; pos++ )
4942    if( r->typ[pos].ord_typ == ro_is)
4943      if( j-- == 0 )
4944      {
4945        return pos;
4946      }
4947
4948  return -1;
4949}
4950
4951
4952
4953
4954
4955
4956/// Changes r by setting induced ordering parameters: limit and reference leading terms
4957/// F belong to r, we will DO a copy! (same to componentWeights)
4958/// We will use it AS IS!
4959/// returns true is everything was allright!
4960BOOLEAN rSetISReference(const ring r, const ideal F, const int i, const int p, const intvec * componentWeights)
4961{
4962  // Put the reference set F into the ring -ordering -recor
4963
4964  // TEST THAT THERE ARE DEGs!!!
4965  // assume( componentWeights == NULL  ); // ???
4966  if( componentWeights != NULL )
4967  {
4968//    assure that the ring r has degrees!!!
4969//    Add weights to degrees of F[i]
4970  }
4971
4972  if (r->typ==NULL)
4973  {
4974    dReportError("Error: WRONG USE of rSetISReference: wrong ring! (typ == NULL)");
4975    return FALSE;
4976  }
4977
4978
4979  int pos = rGetISPos(p, r);
4980
4981  if( pos == -1 )
4982  {
4983    dReportError("Error: WRONG USE of rSetISReference: specified ordering block was not found!!!" );
4984    return FALSE;
4985  }
4986
4987#if MYTEST
4988  if( i != r->typ[pos].data.is.limit )
4989    Print("Changing record on pos: %d\nOld limit: %d --->> New Limit: %d\n", pos, r->typ[pos].data.is.limit, i);
4990#endif
4991
4992  const ideal FF = id_Copy(F, r); // idrHeadR(F, r, r);
4993
4994
4995  if( r->typ[pos].data.is.F != NULL)
4996  {
4997#if MYTEST
4998    PrintS("Deleting old reference set F... \n");        // idShow(r->typ[pos].data.is.F, r);         PrintLn();
4999#endif
5000    id_Delete(&r->typ[pos].data.is.F, r);
5001    r->typ[pos].data.is.F = NULL;
5002  }
5003
5004  assume(r->typ[pos].data.is.F == NULL);
5005
5006  r->typ[pos].data.is.F = FF; // F is owened by ring now! TODO: delete at the end!
5007
5008  if(r->typ[pos].data.is.componentWeights != NULL)
5009  {
5010#if MYTEST
5011    PrintS("Deleting old componentWeights: "); r->typ[pos].data.is.componentWeights->show(); PrintLn();
5012#endif
5013    delete r->typ[pos].data.is.componentWeights;
5014    r->typ[pos].data.is.componentWeights = NULL;
5015  }
5016
5017
5018  assume(r->typ[pos].data.is.componentWeights == NULL);
5019
5020  if( componentWeights != NULL )
5021    componentWeights = ivCopy(componentWeights); // componentWeights is owened by ring now! TODO: delete at the end!
5022
5023  r->typ[pos].data.is.componentWeights = componentWeights;
5024
5025  r->typ[pos].data.is.limit = i; // First induced component
5026
5027#if MYTEST
5028  PrintS("New reference set FF : \n");        idShow(FF, r, r, 1);         PrintLn();
5029#endif
5030
5031  return TRUE;
5032}
5033
5034
5035void rSetSyzComp(int k, const ring r)
5036{
5037  if(k < 0)
5038  {
5039    dReportError("rSetSyzComp with negative limit!");
5040    return;
5041  }
5042
5043  assume( k >= 0 );
5044  if (TEST_OPT_PROT) Print("{%d}", k);
5045  if ((r->typ!=NULL) && (r->typ[0].ord_typ==ro_syz))
5046  {
5047    if( k == r->typ[0].data.syz.limit )
5048      return; // nothing to do
5049
5050    int i;
5051    if (r->typ[0].data.syz.limit == 0)
5052    {
5053      r->typ[0].data.syz.syz_index = (int*) omAlloc0((k+1)*sizeof(int));
5054      r->typ[0].data.syz.syz_index[0] = 0;
5055      r->typ[0].data.syz.curr_index = 1;
5056    }
5057    else
5058    {
5059      r->typ[0].data.syz.syz_index = (int*)
5060        omReallocSize(r->typ[0].data.syz.syz_index,
5061                (r->typ[0].data.syz.limit+1)*sizeof(int),
5062                (k+1)*sizeof(int));
5063    }
5064    for (i=r->typ[0].data.syz.limit + 1; i<= k; i++)
5065    {
5066      r->typ[0].data.syz.syz_index[i] =
5067        r->typ[0].data.syz.curr_index;
5068    }
5069    if(k < r->typ[0].data.syz.limit) // ?
5070    {
5071#ifndef NDEBUG
5072      Warn("rSetSyzComp called with smaller limit (%d) as before (%d)", k, r->typ[0].data.syz.limit);
5073#endif
5074      r->typ[0].data.syz.curr_index = 1 + r->typ[0].data.syz.syz_index[k];
5075    }
5076
5077
5078    r->typ[0].data.syz.limit = k;
5079    r->typ[0].data.syz.curr_index++;
5080  }
5081  else if(
5082            (r->typ!=NULL) &&
5083            (r->typ[0].ord_typ==ro_isTemp)
5084           )
5085  {
5086//      (r->typ[currRing->typ[0].data.isTemp.suffixpos].data.is.limit == k)
5087#ifndef NDEBUG
5088    Warn("rSetSyzComp(%d) in an IS ring! Be careful!", k);
5089#endif
5090  }
5091  else
5092  if ((r->order[0]!=ringorder_c) && (k!=0)) // ???
5093  {
5094    dReportError("syzcomp in incompatible ring");
5095  }
5096#ifdef PDEBUG
5097  extern int pDBsyzComp;
5098  pDBsyzComp=k;
5099#endif
5100}
5101
5102// return the max-comonent wchich has syzIndex i
5103int rGetMaxSyzComp(int i, const ring r)
5104{
5105  if ((r->typ!=NULL) && (r->typ[0].ord_typ==ro_syz) &&
5106      r->typ[0].data.syz.limit > 0 && i > 0)
5107  {
5108    assume(i <= r->typ[0].data.syz.limit);
5109    int j;
5110    for (j=0; j<r->typ[0].data.syz.limit; j++)
5111    {
5112      if (r->typ[0].data.syz.syz_index[j] == i  &&
5113          r->typ[0].data.syz.syz_index[j+1] != i)
5114      {
5115        assume(r->typ[0].data.syz.syz_index[j+1] == i+1);
5116        return j;
5117      }
5118    }
5119    return r->typ[0].data.syz.limit;
5120  }
5121  else
5122  {
5123    return 0;
5124  }
5125}
5126
5127BOOLEAN rRing_is_Homog(ring r)
5128{
5129  if (r == NULL) return FALSE;
5130  int i, j, nb = rBlocks(r);
5131  for (i=0; i<nb; i++)
5132  {
5133    if (r->wvhdl[i] != NULL)
5134    {
5135      int length = r->block1[i] - r->block0[i];
5136      int* wvhdl = r->wvhdl[i];
5137      if (r->order[i] == ringorder_M) length *= length;
5138      assume(omSizeOfAddr(wvhdl) >= length*sizeof(int));
5139
5140      for (j=0; j< length; j++)
5141      {
5142        if (wvhdl[j] != 0 && wvhdl[j] != 1) return FALSE;
5143      }
5144    }
5145  }
5146  return TRUE;
5147}
5148
5149BOOLEAN rRing_has_CompLastBlock(ring r)
5150{
5151  assume(r != NULL);
5152  int lb = rBlocks(r) - 2;
5153  return (r->order[lb] == ringorder_c || r->order[lb] == ringorder_C);
5154}
5155
5156n_coeffType rFieldType(ring r)
5157{
5158  return (r->cf->type);
5159  if (rField_is_Zp(r))     return n_Zp;
5160  if (rField_is_Q(r))      return n_Q;
5161  if (rField_is_R(r))      return n_R;
5162  if (rField_is_GF(r))     return n_GF;
5163  if (rField_is_long_R(r)) return n_long_R;
5164  if (rField_is_Zp_a(r))   return n_Zp_a;
5165  if (rField_is_Q_a(r))    return n_Q_a;
5166  if (rField_is_long_C(r)) return n_long_C;
5167  #ifdef HAVE_RINGS
5168   if (rField_is_Ring_Z(r)) return n_Z;
5169   if (rField_is_Ring_ModN(r)) return n_Zn;
5170   if (rField_is_Ring_PtoM(r)) return n_Zpn;
5171   if (rField_is_Ring_2toM(r)) return  n_Z2m;
5172  #endif
5173
5174  return n_unknown;
5175}
5176
5177int64 * rGetWeightVec(ring r)
5178{
5179  assume(r!=NULL);
5180  assume(r->OrdSize>0);
5181  int i=0;
5182  while((r->typ[i].ord_typ!=ro_wp64) && (r->typ[i].ord_typ>0)) i++;
5183  assume(r->typ[i].ord_typ==ro_wp64);
5184  return (int64*)(r->typ[i].data.wp64.weights64);
5185}
5186
5187void rSetWeightVec(ring r, int64 *wv)
5188{
5189  assume(r!=NULL);
5190  assume(r->OrdSize>0);
5191  assume(r->typ[0].ord_typ==ro_wp64);
5192  memcpy(r->typ[0].data.wp64.weights64,wv,r->N*sizeof(int64));
5193}
5194
5195#include <ctype.h>
5196
5197static int rRealloc1(ring r, int size, int pos)
5198{
5199  r->order=(int*)omReallocSize(r->order, size*sizeof(int), (size+1)*sizeof(int));
5200  r->block0=(int*)omReallocSize(r->block0, size*sizeof(int), (size+1)*sizeof(int));
5201  r->block1=(int*)omReallocSize(r->block1, size*sizeof(int), (size+1)*sizeof(int));
5202  r->wvhdl=(int **)omReallocSize(r->wvhdl,size*sizeof(int *), (size+1)*sizeof(int *));
5203  for(int k=size; k>pos; k--) r->wvhdl[k]=r->wvhdl[k-1];
5204  r->order[size]=0;
5205  size++;
5206  return size;
5207}
5208//static int rReallocM1(ring r, ring src, int size, int pos)
5209//{
5210//  r->order=(int*)omReallocSize(r->order, size*sizeof(int), (size-1)*sizeof(int));
5211//  r->block0=(int*)omReallocSize(r->block0, size*sizeof(int), (size-1)*sizeof(int));
5212//  r->block1=(int*)omReallocSize(r->block1, size*sizeof(int), (size-1)*sizeof(int));
5213//  r->wvhdl=(int **)omReallocSize(r->wvhdl,size*sizeof(int *), (size-1)*sizeof(int *));
5214//  for(int k=pos+1; k<size; k++) r->wvhdl[k]=r->wvhdl[k+1];
5215//  size--;
5216//  return size;
5217//}
5218static void rOppWeight(int *w, int l)
5219{
5220  int i2=(l+1)/2;
5221  for(int j=0; j<=i2; j++)
5222  {
5223    int t=w[j];
5224    w[j]=w[l-j];
5225    w[l-j]=t;
5226  }
5227}
5228
5229#define rOppVar(R,I) (rVar(R)+1-I)
5230
5231ring rOpposite(ring src)
5232  /* creates an opposite algebra of R */
5233  /* that is R^opp, where f (*^opp) g = g*f  */
5234  /* treats the case of qring */
5235{
5236  if (src == NULL) return(NULL);
5237
5238#ifdef RDEBUG
5239  rTest(src);
5240#endif
5241
5242  //rChangeCurrRing(src);
5243
5244#ifdef RDEBUG
5245  rTest(src);
5246//  rWrite(src);
5247//  rDebugPrint(src);
5248#endif
5249
5250
5251  ring r = rCopy0(src,FALSE); /* qideal will be deleted later on!!! */
5252
5253  // change vars v1..vN -> vN..v1
5254  int i;
5255  int i2 = (rVar(r)-1)/2;
5256  for(i=i2; i>=0; i--)
5257  {
5258    // index: 0..N-1
5259    //Print("ex var names: %d <-> %d\n",i,rOppVar(r,i));
5260    // exchange names
5261    char *p;
5262    p = r->names[rVar(r)-1-i];
5263    r->names[rVar(r)-1-i] = r->names[i];
5264    r->names[i] = p;
5265  }
5266//  i2=(rVar(r)+1)/2;
5267//  for(int i=i2; i>0; i--)
5268//  {
5269//    // index: 1..N
5270//    //Print("ex var places: %d <-> %d\n",i,rVar(r)+1-i);
5271//    // exchange VarOffset
5272//    int t;
5273//    t=r->VarOffset[i];
5274//    r->VarOffset[i]=r->VarOffset[rOppVar(r,i)];
5275//    r->VarOffset[rOppVar(r,i)]=t;
5276//  }
5277  // change names:
5278  for (i=rVar(r)-1; i>=0; i--)
5279  {
5280    char *p=r->names[i];
5281    if(isupper(*p)) *p = tolower(*p);
5282    else            *p = toupper(*p);
5283  }
5284  // change ordering: listing
5285  // change ordering: compare
5286//  for(i=0; i<r->OrdSize; i++)
5287//  {
5288//    int t,tt;
5289//    switch(r->typ[i].ord_typ)
5290//    {
5291//      case ro_dp:
5292//      //
5293//        t=r->typ[i].data.dp.start;
5294//        r->typ[i].data.dp.start=rOppVar(r,r->typ[i].data.dp.end);
5295//        r->typ[i].data.dp.end=rOppVar(r,t);
5296//        break;
5297//      case ro_wp:
5298//      case ro_wp_neg:
5299//      {
5300//        t=r->typ[i].data.wp.start;
5301//        r->typ[i].data.wp.start=rOppVar(r,r->typ[i].data.wp.end);
5302//        r->typ[i].data.wp.end=rOppVar(r,t);
5303//        // invert r->typ[i].data.wp.weights
5304//        rOppWeight(r->typ[i].data.wp.weights,
5305//                   r->typ[i].data.wp.end-r->typ[i].data.wp.start);
5306//        break;
5307//      }
5308//      //case ro_wp64:
5309//      case ro_syzcomp:
5310//      case ro_syz:
5311//         WerrorS("not implemented in rOpposite");
5312//         // should not happen
5313//         break;
5314//
5315//      case ro_cp:
5316//        t=r->typ[i].data.cp.start;
5317//        r->typ[i].data.cp.start=rOppVar(r,r->typ[i].data.cp.end);
5318//        r->typ[i].data.cp.end=rOppVar(r,t);
5319//        break;
5320//      case ro_none:
5321//      default:
5322//       Werror("unknown type in rOpposite(%d)",r->typ[i].ord_typ);
5323//       break;
5324//    }
5325//  }
5326  // Change order/block structures (needed for rPrint, rAdd etc.)
5327  int j=0;
5328  int l=rBlocks(src);
5329  for(i=0; src->order[i]!=0; i++)
5330  {
5331    switch (src->order[i])
5332    {
5333      case ringorder_c: /* c-> c */
5334      case ringorder_C: /* C-> C */
5335      case ringorder_no /*=0*/: /* end-of-block */
5336        r->order[j]=src->order[i];
5337        j++; break;
5338      case ringorder_lp: /* lp -> rp */
5339        r->order[j]=ringorder_rp;
5340        r->block0[j]=rOppVar(r, src->block1[i]);
5341        r->block1[j]=rOppVar(r, src->block0[i]);
5342        break;
5343      case ringorder_rp: /* rp -> lp */
5344        r->order[j]=ringorder_lp;
5345        r->block0[j]=rOppVar(r, src->block1[i]);
5346        r->block1[j]=rOppVar(r, src->block0[i]);
5347        break;
5348      case ringorder_dp: /* dp -> a(1..1),ls */
5349      {
5350        l=rRealloc1(r,l,j);
5351        r->order[j]=ringorder_a;
5352        r->block0[j]=rOppVar(r, src->block1[i]);
5353        r->block1[j]=rOppVar(r, src->block0[i]);
5354        r->wvhdl[j]=(int*)omAlloc((r->block1[j]-r->block0[j]+1)*sizeof(int));
5355        for(int k=r->block0[j]; k<=r->block1[j]; k++)
5356          r->wvhdl[j][k-r->block0[j]]=1;
5357        j++;
5358        r->order[j]=ringorder_ls;
5359        r->block0[j]=rOppVar(r, src->block1[i]);
5360        r->block1[j]=rOppVar(r, src->block0[i]);
5361        j++;
5362        break;
5363      }
5364      case ringorder_Dp: /* Dp -> a(1..1),rp */
5365      {
5366        l=rRealloc1(r,l,j);
5367        r->order[j]=ringorder_a;
5368        r->block0[j]=rOppVar(r, src->block1[i]);
5369        r->block1[j]=rOppVar(r, src->block0[i]);
5370        r->wvhdl[j]=(int*)omAlloc((r->block1[j]-r->block0[j]+1)*sizeof(int));
5371        for(int k=r->block0[j]; k<=r->block1[j]; k++)
5372          r->wvhdl[j][k-r->block0[j]]=1;
5373        j++;
5374        r->order[j]=ringorder_rp;
5375        r->block0[j]=rOppVar(r, src->block1[i]);
5376        r->block1[j]=rOppVar(r, src->block0[i]);
5377        j++;
5378        break;
5379      }
5380      case ringorder_wp: /* wp -> a(...),ls */
5381      {
5382        l=rRealloc1(r,l,j);
5383        r->order[j]=ringorder_a;
5384        r->block0[j]=rOppVar(r, src->block1[i]);
5385        r->block1[j]=rOppVar(r, src->block0[i]);
5386        r->wvhdl[j]=r->wvhdl[j+1]; r->wvhdl[j+1]=NULL;
5387        rOppWeight(r->wvhdl[j], r->block1[j]-r->block0[j]);
5388        j++;
5389        r->order[j]=ringorder_ls;
5390        r->block0[j]=rOppVar(r, src->block1[i]);
5391        r->block1[j]=rOppVar(r, src->block0[i]);
5392        j++;
5393        break;
5394      }
5395      case ringorder_Wp: /* Wp -> a(...),rp */
5396      {
5397        l=rRealloc1(r,l,j);
5398        r->order[j]=ringorder_a;
5399        r->block0[j]=rOppVar(r, src->block1[i]);
5400        r->block1[j]=rOppVar(r, src->block0[i]);
5401        r->wvhdl[j]=r->wvhdl[j+1]; r->wvhdl[j+1]=NULL;
5402        rOppWeight(r->wvhdl[j], r->block1[j]-r->block0[j]);
5403        j++;
5404        r->order[j]=ringorder_rp;
5405        r->block0[j]=rOppVar(r, src->block1[i]);
5406        r->block1[j]=rOppVar(r, src->block0[i]);
5407        j++;
5408        break;
5409      }
5410      case ringorder_M: /* M -> M */
5411      {
5412        r->order[j]=ringorder_M;
5413        r->block0[j]=rOppVar(r, src->block1[i]);
5414        r->block1[j]=rOppVar(r, src->block0[i]);
5415        int n=r->block1[j]-r->block0[j];
5416        /* M is a (n+1)x(n+1) matrix */
5417        for (int nn=0; nn<=n; nn++)
5418        {
5419          rOppWeight(&(r->wvhdl[j][nn*(n+1)]), n /*r->block1[j]-r->block0[j]*/);
5420        }
5421        j++;
5422        break;
5423      }
5424      case ringorder_a: /*  a(...),ls -> wp/dp */
5425      {
5426        r->block0[j]=rOppVar(r, src->block1[i]);
5427        r->block1[j]=rOppVar(r, src->block0[i]);
5428        rOppWeight(r->wvhdl[j], r->block1[j]-r->block0[j]);
5429        if (src->order[i+1]==ringorder_ls)
5430        {
5431          r->order[j]=ringorder_wp;
5432          i++;
5433          //l=rReallocM1(r,src,l,j);
5434        }
5435        else
5436        {
5437          r->order[j]=ringorder_a;
5438        }
5439        j++;
5440        break;
5441      }
5442      // not yet done:
5443      case ringorder_ls:
5444      case ringorder_rs:
5445      case ringorder_ds:
5446      case ringorder_Ds:
5447      case ringorder_ws:
5448      case ringorder_Ws:
5449      // should not occur:
5450      case ringorder_S:
5451      case ringorder_IS:
5452      case ringorder_s:
5453      case ringorder_aa:
5454      case ringorder_L:
5455      case ringorder_unspec:
5456        Werror("order %s not (yet) supported", rSimpleOrdStr(src->order[i]));
5457        break;
5458    }
5459  }
5460  rComplete(r);
5461
5462
5463#ifdef RDEBUG
5464  rTest(r);
5465#endif
5466
5467  //rChangeCurrRing(r);
5468
5469#ifdef RDEBUG
5470  rTest(r);
5471//  rWrite(r);
5472//  rDebugPrint(r);
5473#endif
5474
5475
5476#ifdef HAVE_PLURAL
5477  // now, we initialize a non-comm structure on r
5478  if (rIsPluralRing(src))
5479  {
5480//    assume( currRing == r);
5481
5482    int *perm       = (int *)omAlloc0((rVar(r)+1)*sizeof(int));
5483    int *par_perm   = NULL;
5484    nMapFunc nMap   = n_SetMap(src->cf,r->cf);
5485    int ni,nj;
5486    for(i=1; i<=r->N; i++)
5487    {
5488      perm[i] = rOppVar(r,i);
5489    }
5490
5491    matrix C = mpNew(rVar(r),rVar(r));
5492    matrix D = mpNew(rVar(r),rVar(r));
5493
5494    for (i=1; i< rVar(r); i++)
5495    {
5496      for (j=i+1; j<=rVar(r); j++)
5497      {
5498        ni = r->N +1 - i;
5499        nj = r->N +1 - j; /* i<j ==>   nj < ni */
5500
5501        assume(MATELEM(src->GetNC()->C,i,j) != NULL);
5502        MATELEM(C,nj,ni) = p_PermPoly(MATELEM(src->GetNC()->C,i,j),perm,src,r, nMap,par_perm,rPar(src));
5503
5504        if(MATELEM(src->GetNC()->D,i,j) != NULL)
5505          MATELEM(D,nj,ni) = p_PermPoly(MATELEM(src->GetNC()->D,i,j),perm,src,r, nMap,par_perm,rPar(src));
5506      }
5507    }
5508
5509    id_Test((ideal)C, r);
5510    id_Test((ideal)D, r);
5511
5512    if (nc_CallPlural(C, D, NULL, NULL, r, false, false, true, r)) // no qring setup!
5513      WarnS("Error initializing non-commutative multiplication!");
5514
5515#ifdef RDEBUG
5516    rTest(r);
5517//    rWrite(r);
5518//    rDebugPrint(r);
5519#endif
5520
5521    assume( r->GetNC()->IsSkewConstant == src->GetNC()->IsSkewConstant);
5522
5523    omFreeSize((ADDRESS)perm,(rVar(r)+1)*sizeof(int));
5524  }
5525#endif /* HAVE_PLURAL */
5526
5527  /* now oppose the qideal for qrings */
5528  if (src->qideal != NULL)
5529  {
5530    id_Delete(&(r->qideal), r);
5531
5532#ifdef HAVE_PLURAL
5533    r->qideal = idOppose(src, src->qideal, r); // into the currRing: r
5534#else
5535    r->qideal = id_Copy(src->qideal, r); // ?
5536#endif
5537
5538#ifdef HAVE_PLURAL
5539    if( rIsPluralRing(r) )
5540    {
5541      nc_SetupQuotient(r);
5542#ifdef RDEBUG
5543      rTest(r);
5544//      rWrite(r);
5545//      rDebugPrint(r);
5546#endif
5547    }
5548#endif
5549  }
5550#ifdef HAVE_PLURAL
5551  if( rIsPluralRing(r) )
5552    assume( ncRingType(r) == ncRingType(src) );
5553#endif
5554  rTest(r);
5555
5556  return r;
5557}
5558
5559ring rEnvelope(ring R)
5560  /* creates an enveloping algebra of R */
5561  /* that is R^e = R \tensor_K R^opp */
5562{
5563  ring Ropp = rOpposite(R);
5564  ring Renv = NULL;
5565  int stat = rSum(R, Ropp, Renv); /* takes care of qideals */
5566  if ( stat <=0 )
5567    WarnS("Error in rEnvelope at rSum");
5568  rTest(Renv);
5569  return Renv;
5570}
5571
5572#ifdef HAVE_PLURAL
5573BOOLEAN nc_rComplete(const ring src, ring dest, bool bSetupQuotient)
5574/* returns TRUE is there were errors */
5575/* dest is actualy equals src with the different ordering */
5576/* we map src->nc correctly to dest->src */
5577/* to be executed after rComplete, before rChangeCurrRing */
5578{
5579// NOTE: Originally used only by idElimination to transfer NC structure to dest
5580// ring created by dirty hack (without nc_CallPlural)
5581  rTest(src);
5582
5583  assume(!rIsPluralRing(dest)); // destination must be a newly constructed commutative ring
5584
5585  if (!rIsPluralRing(src))
5586  {
5587    return FALSE;
5588  }
5589
5590  const int N = dest->N;
5591
5592  assume(src->N == N);
5593
5594//  ring save = currRing;
5595
5596//  if (dest != save)
5597//    rChangeCurrRing(dest);
5598
5599  const ring srcBase = src;
5600
5601  assume( n_SetMap(srcBase->cf,dest->cf) == n_SetMap(dest->cf,dest->cf) ); // currRing is important here!
5602
5603  matrix C = mpNew(N,N); // ring independent
5604  matrix D = mpNew(N,N);
5605
5606  matrix C0 = src->GetNC()->C;
5607  matrix D0 = src->GetNC()->D;
5608
5609  // map C and D into dest
5610  for (int i = 1; i < N; i++)
5611  {
5612    for (int j = i + 1; j <= N; j++)
5613    {
5614      const number n = n_Copy(p_GetCoeff(MATELEM(C0,i,j), srcBase), srcBase->cf); // src, mapping for coeffs into currRing = dest!
5615      const poly   p = p_NSet(n, dest);
5616      MATELEM(C,i,j) = p;
5617      if (MATELEM(D0,i,j) != NULL)
5618        MATELEM(D,i,j) = prCopyR(MATELEM(D0,i,j), srcBase, dest); // ?
5619    }
5620  }
5621  /* One must test C and D _only_ in r->GetNC()->basering!!! not in r!!! */
5622
5623  id_Test((ideal)C, dest); 
5624  id_Test((ideal)D, dest);
5625
5626  if (nc_CallPlural(C, D, NULL, NULL, dest, bSetupQuotient, false, true, dest)) // also takes care about quotient ideal
5627  {
5628    //WarnS("Error transferring non-commutative structure");
5629    // error message should be in the interpreter interface
5630
5631    mp_Delete(&C, dest);
5632    mp_Delete(&D, dest);
5633
5634//    if (currRing != save)
5635//       rChangeCurrRing(save);
5636
5637    return TRUE;
5638  }
5639
5640//  mp_Delete(&C, dest); // used by nc_CallPlural!
5641//  mp_Delete(&D, dest);
5642
5643//  if (dest != save)
5644//    rChangeCurrRing(save);
5645
5646  assume(rIsPluralRing(dest));
5647  return FALSE;
5648}
5649#endif
5650
5651void rModify_a_to_A(ring r)
5652// to be called BEFORE rComplete:
5653// changes every Block with a(...) to A(...)
5654{
5655   int i=0;
5656   int j;
5657   while(r->order[i]!=0)
5658   {
5659      if (r->order[i]==ringorder_a)
5660      {
5661        r->order[i]=ringorder_a64;
5662        int *w=r->wvhdl[i];
5663        int64 *w64=(int64 *)omAlloc((r->block1[i]-r->block0[i]+1)*sizeof(int64));
5664        for(j=r->block1[i]-r->block0[i];j>=0;j--)
5665                w64[j]=(int64)w[j];
5666        r->wvhdl[i]=(int*)w64;
5667        omFreeSize(w,(r->block1[i]-r->block0[i]+1)*sizeof(int));
5668      }
5669      i++;
5670   }
5671}
Note: See TracBrowser for help on using the repository browser.