source: git/polys/monomials/ring.cc @ 1b816a3

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