source: git/libpolys/polys/monomials/ring.cc @ 028192

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