Impressum fplsa4.c

/************************************************************************/
/*      File fplsa4.c   Last touch January 5, 1999                      */
/*      Calculation of finitely presented Lie superalgebras             */
/*      Version 4 of April 15, 1997.                                    */
/*      e-mail: kornyak@jinr.ru, gerdt@jinr.ru                          */
/*                      Contents                                        */
/*_0    Choice of compilation                                           */
/*_1    System constants                                                */
/*_2    Type definitions                                                */
/*_3    Constants and enumerations                                      */
/*_4    Macrodefinitions                                                */
/*_5    Global variables and arrays                                     */
/*_6    Function descriptions                                           */
/*_6_0    Main and top level functions                                  */
/*_6_1    Pairing functions                                             */
/*_6_2    Substitution (replacing) functions                            */
/*_6_3    Lie and scalar algebra functions                              */
/*_6_4    Scalar polynomial algebraic functions                         */
/*_6_5    Big number functions                                          */
/*_6_6    Copy and delete functions                                     */
/*_6_7    Technical functions                                           */
/*_6_8    Input functions                                               */
/*_6_9    Output functions                                              */
/*_6_10   Debugging functions                                           */
/************************************************************************/


/*_0    Choice of compilation============================================*/

#define RATIONAL_FIELD  /*  Working over the field R ??
                            otherwise over the ring Z */
#define ECHO_TO_SCREEN        /* Echo session to screen ?? */
#define RELATION_N_TO_SCREEN  /* Watch increasing RelationN ?? */
//#define SPP_2000  /* Big ending computer ?? */
#define SPACE_STATISTICS  /* Space statistics ?? */
//#define INTEGER_MAX_SIZE  /* Multiprecision number maximum size ?? */
//#define INTEGER_ALLOCATION_CHECK  /* Control of integer allocations ?? */
//#define POLY_ARRAY_ALLOCATION_CHECK /* Control of allocations of ?? polynomial arrays in stack */

/* GAP output ?? */
/* Avoid message file, session file,                              */
/*       compulsory suffix '.in' for input file,                  */
/*       and printing comments to the screen ?                    */
#define GAP
/**/

//#define DEBUG /* Debugging ?? */
//#define MEMORY /* Check memory balance ?? */

/* Include files */

#include <stdio.h>
#include <time.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>
#if defined(SPP_2000)
#include <alloca.h>  /* This file the genuine SPP compiler requires */
#endif


#if defined(__GNUC__) || defined(__clang__)
#define ATTRIBUTE_NORETURN  __attribute__ ((noreturn))
#else
#define ATTRIBUTE_NORETURN
#endif


#if defined(DEBUG)  /* Debug definitions ===============================*/

/* Set condition for debug output ??
*/
#define D_CONDITION   if(Debug>=5873)
/*            Examples of D_CONDITION
      `empty' if(Debug>=11951211) current!!
      if(Debug>=25283LU)*e7* *e70*
      if(Debug%10 == 0)  if(Debug>=17566)    */
/* Set condition to stop debugging ??
*/
#define D_EXIT  if(Debug > 5882) EXIT;
/*            Examples of D_EXIT
      `empty' if(Debug > 21420LU) EXIT;
       if(Debug > 26969LU) EXIT;*e7* if(Debug > 21420LU) EXIT;*e70*
       if(Debug > 30881LU) EXIT; if(Debug >= 18399) EXIT;
*/

/* Switches for debugging particular functions ?? */

#define D_CHECK_EXACTNESS_OF_DIVISION
//#define D_ADD_PAIR_TO_LIE_MONOMIAL
#define D_GENERATE_RELATIONS
//#if defined(D_GENERATE_RELATIONS)  /* Put tables ?? */
//#define D_PUT_RELATIONS
//#define D_PUT_LIE_MONOMIAL
//#define D_GET_LIE_MONOMIAL
//#define D_GET_LIE_SUM
//#define D_GET_LIE_TERM
#define D_INTEGER_CANCELLATION
#define D_INTEGER_GCD
#define D_INTEGER_PRODUCT
#define D_INTEGER_QUOTIENT
#define D_INTEGER_SUM
#define D_LIE_SUM_ADDITION
#define D_LIE_SUM_DIV_INTEGER
//#define D_LIE_SUM_DIV_SCALAR_SUM
#define D_LIE_SUM_MULT_INTEGER
//#define D_LIE_SUM_MULT_SCALAR_SUM
#define D_MAKE_RELATION_RHS
//#define D_NEW_JACOBI_RELATION
#define D_NORMALIZE_RELATION
//#define D_PAIR_MONOMIAL_MONOMIAL
//#define D_PAIR_MONOMIAL_SUM
//#define D_PAIR_SUM_MONOMIAL
//#define D_PAIR_SUM_SUM
//#define D_POLY_CONTENT
//#define D_POLY_GCD
//#define D_POLY_QUOTIENT
//#define D_POLY_PSEUDO_REMAINDER
//#define D_SCALAR_SUM_CANCELLATION
#define D_SUBSTITUTE_RELATION_IN_RELATION
#define D_IN_SET     uint debug=Debug;\
                     D_CONDITION\
                       PutDebugHeader(debug, f_name, "in"),
#define D_IN_CLOSE   Debug++; D_EXIT

#define D_OUT_OPEN   D_CONDITION PutDebugHeader(debug, f_name, "out"),

#else      /* Empty definitions for DEBUG off */

#define D_IN_SET
#define D_IN_CLOSE
#define D_OUT_OPEN

#endif

/* Check memory balance definitions */

#if defined(MEMORY)              /* `n_...' are unique */
#define IN_SET_N_LT      int n_lt  = -CurrentNLT;
#define IN_SET_N_INT     int n_int = -CurrentNINT;
#define IN_SET_N_ST      int n_st  = -CurrentNST;
#define IN_SET_N_SF      int n_sf  = -CurrentNSF;

#define OUT_SET_N_LT     n_lt  += CurrentNLT;
#define OUT_SET_N_INT    n_int += CurrentNINT;
#define OUT_SET_N_ST     n_st  += CurrentNST;
#define OUT_SET_N_SF     n_sf  += CurrentNSF;

#define ADD_LIE_SUM_NS(a) \
    AddLieSumNs(a, PLUS, &n_lt, &n_int, &n_st, &n_sf);

#define SUBTRACT_LIE_SUM_NS(a) \
    AddLieSumNs(a, MINUS, &n_lt, &n_int, &n_st, &n_sf);

#define ADD_SCALAR_SUM_NS(a) \
    AddScalarSumNs(a, PLUS, &n_int, &n_st, &n_sf);

#define SUBTRACT_SCALAR_SUM_NS(a) \
    AddScalarSumNs(a, MINUS, &n_int, &n_st, &n_sf);

#define CHECK_NS \
  if(n_lt != 0) {PutNodeBalance("\nNodeLT (Lie terms)",\
                                f_name, n_lt); EXIT;}\
  if(n_int != 0) {PutIntegerBalance(f_name, n_int); EXIT;}\
  if(n_st != 0) {PutNodeBalance("\nNodeST (scalar terms)",\
                                f_name, n_st); EXIT;}\
  if(n_sf != 0) {PutNodeBalance("\nNodeSF (scalar factors)",\
                                f_name, n_sf); EXIT;}

/* Particular functions */
/*----------------------------------------------------------*/
#define M_OUT_GET_LIE_MONOMIAL \
  OUT_SET_NS\
  SUBTRACT_LIE_SUM_NS(a)\
  CHECK_NS
/*----------------------------------------------------------*/
#define M_OUT_GET_LIE_SUM \
  OUT_SET_NS\
  SUBTRACT_LIE_SUM_NS(lsum)\
  CHECK_NS
/*----------------------------------------------------------*/
#define M_OUT_GET_LIE_TERM \
  OUT_SET_NS\
  SUBTRACT_LIE_SUM_NS(lterm)\
  CHECK_NS
/*----------------------------------------------------------*/
#define M_IN_LIE_SUM_ADDITION  \
  ADD_LIE_SUM_NS(a)\
  ADD_LIE_SUM_NS(b)
#define M_OUT_LIE_SUM_ADDITION  \
  OUT_SET_NS\
  SUBTRACT_LIE_SUM_NS(sum)\
  CHECK_NS
/*----------------------------------------------------------*/
/*   No sense to print debug information, memory check only */
#define IN_LIE_SUM_COPY  \
    IN_SET_FUNCTION_NAME("LieSumCopy...")\
    IN_SET_NS
#define OUT_LIE_SUM_COPY  \
    OUT_SET_NS\
    SUBTRACT_LIE_SUM_NS(ca)\
    CHECK_NS
/*----------------------------------------------------------*/
/*   No sense to print debug information, memory check only */
#define IN_LIE_SUM_KILL  \
    IN_SET_FUNCTION_NAME("LieSumKill...")\
    IN_SET_NS\
    ADD_LIE_SUM_NS(a)
#define OUT_LIE_SUM_KILL \
    OUT_SET_NS\
    CHECK_NS
/*----------------------------------------------------------*/
#define M_IN_LIE_SUM_DIV_INTEGER  \
  ADD_LIE_SUM_NS(b)
#define M_OUT_LIE_SUM_DIV_INTEGER  \
  OUT_SET_NS\
  SUBTRACT_LIE_SUM_NS(b)\
  CHECK_NS
/*----------------------------------------------------------*/
#define M_IN_LIE_SUM_DIV_SCALAR_SUM  \
  ADD_LIE_SUM_NS(b)\
  ADD_SCALAR_SUM_NS(den)
#define M_OUT_LIE_SUM_DIV_SCALAR_SUM  \
  OUT_SET_NS\
  SUBTRACT_LIE_SUM_NS(b)\
  CHECK_NS
/*----------------------------------------------------------*/
#define M_IN_LIE_SUM_MULT_SCALAR_SUM  \
  ADD_LIE_SUM_NS(b)\
  ADD_SCALAR_SUM_NS(num)
#define M_OUT_LIE_SUM_MULT_SCALAR_SUM  \
  OUT_SET_NS\
  SUBTRACT_LIE_SUM_NS(b)\
  CHECK_NS
/*----------------------------------------------------------*/
#define M_IN_LIE_SUM_MULT_INTEGER  \
  ADD_LIE_SUM_NS(b)
#define M_OUT_LIE_SUM_MULT_INTEGER  M_OUT_LIE_SUM_DIV_INTEGER
/*----------------------------------------------------------*/
#define M_OUT_MAKE_RELATION_RHS    M_OUT_GET_LIE_MONOMIAL
/*----------------------------------------------------------*/
#define M_OUT_NEW_JACOBI_RELATION  M_OUT_GET_LIE_MONOMIAL
/*----------------------------------------------------------*/
#define M_IN_NORMALIZE_RELATION  \
  ADD_LIE_SUM_NS(a)
#define M_OUT_NORMALIZE_RELATION   M_OUT_GET_LIE_MONOMIAL
/*----------------------------------------------------------*/
#define M_OUT_PAIR_MONOMIAL_MONOMIAL  M_OUT_GET_LIE_MONOMIAL
/*----------------------------------------------------------*/
#define M_IN_PAIR_MONOMIAL_SUM  \
  ADD_LIE_SUM_NS(a)
#define M_OUT_PAIR_MONOMIAL_SUM  \
  OUT_SET_NS\
  SUBTRACT_LIE_SUM_NS(s)\
  CHECK_NS
/*----------------------------------------------------------*/
#define M_IN_PAIR_SUM_MONOMIAL  \
  ADD_LIE_SUM_NS(a)
#define M_OUT_PAIR_SUM_MONOMIAL  \
  OUT_SET_NS\
  SUBTRACT_LIE_SUM_NS(s)\
  CHECK_NS
/*----------------------------------------------------------*/
#define M_IN_PAIR_SUM_SUM  \
  ADD_LIE_SUM_NS(a)\
  ADD_LIE_SUM_NS(b)
#define M_OUT_PAIR_SUM_SUM  \
  OUT_SET_NS\
  SUBTRACT_LIE_SUM_NS(a)\
  CHECK_NS
/*----------------------------------------------------------*/
#define M_IN_POLY_CONTENT
#define M_OUT_POLY_CONTENT  \
  OUT_SET_NS\
  SUBTRACT_SCALAR_SUM_NS(b)\
  CHECK_NS
/*----------------------------------------------------------*/
#define M_IN_POLY_GCD
#define M_OUT_POLY_GCD  \
  OUT_SET_NS\
  SUBTRACT_SCALAR_SUM_NS(b)\
  CHECK_NS
/*----------------------------------------------------------*/
#define M_IN_POLY_QUOTIENT  \
  ADD_SCALAR_SUM_NS(a)\
  ADD_SCALAR_SUM_NS(b)
#define M_OUT_POLY_QUOTIENT  \
  OUT_SET_NS\
  SUBTRACT_SCALAR_SUM_NS(c)\
  CHECK_NS
/*----------------------------------------------------------*/
#define M_IN_POLY_PSEUDO_REMAINDER  \
  ADD_SCALAR_SUM_NS(a)\
  ADD_SCALAR_SUM_NS(b)
#define M_OUT_POLY_PSEUDO_REMAINDER  \
  OUT_SET_NS\
  SUBTRACT_SCALAR_SUM_NS(a)\
  CHECK_NS
/*----------------------------------------------------------*/
/*   No sense to print debug information, memory check only */
#define IN_REDUCE_RELATIONS  \
    IN_SET_FUNCTION_NAME("ReduceRelations")\
    IN_SET_NS\
    {\
      int o;\
      for(o = 0; o < RelationN; o++)\
        ADD_LIE_SUM_NS(RELATION_LIE_SUM(o))\
    }
#define OUT_REDUCE_RELATIONS  \
    OUT_SET_NS\
    {\
      int o;\
      for(o = 0; o < RelationN; o++)\
        SUBTRACT_LIE_SUM_NS(RELATION_LIE_SUM(o))\
    }\
    CHECK_NS
/*----------------------------------------------------------*/
#define M_IN_SCALAR_SUM_CANCELLATION  \
  ADD_SCALAR_SUM_NS(*pnum)\
  ADD_SCALAR_SUM_NS(*pden)
#define M_OUT_SCALAR_SUM_CANCELLATION  \
  OUT_SET_NS\
  SUBTRACT_SCALAR_SUM_NS(*pnum)\
  SUBTRACT_SCALAR_SUM_NS(*pden)\
  CHECK_NS
/*----------------------------------------------------------*/
#define M_IN_SUBSTITUTE_RELATION_IN_RELATION  \
  ADD_LIE_SUM_NS(a)
#define M_OUT_SUBSTITUTE_RELATION_IN_RELATION  \
  OUT_SET_NS\
  SUBTRACT_LIE_SUM_NS(a)\
  CHECK_NS

#define  MM_CURRENT_N_LT    ,--CurrentNLT
#define  PP_CURRENT_N_LT    ++CurrentNLT;

#define  MM_CURRENT_N_SF    ,--CurrentNSF
#define  PP_CURRENT_N_SF    ++CurrentNSF;

#define  MM_CURRENT_N_ST    ,--CurrentNST
#define  PP_CURRENT_N_ST    ++CurrentNST;

#define  MM_CURRENT_N_INT   ,--CurrentNINT
#define  PP_CURRENT_N_INT   ;++CurrentNINT

#else       /* MEMORY off */

#define  MM_CURRENT_N_LT
#define  PP_CURRENT_N_LT
#define  MM_CURRENT_N_SF
#define  PP_CURRENT_N_SF
#define  MM_CURRENT_N_ST
#define  PP_CURRENT_N_ST
#define  MM_CURRENT_N_INT
#define  PP_CURRENT_N_INT

#define IN_SET_N_LT
#define IN_SET_N_INT
#define IN_SET_N_ST
#define IN_SET_N_SF

#define OUT_SET_N_LT
#define OUT_SET_N_INT
#define OUT_SET_N_ST
#define OUT_SET_N_SF

#define M_OUT_GET_LIE_MONOMIAL
#define M_OUT_GET_LIE_SUM
#define M_OUT_GET_LIE_TERM
#define M_IN_LIE_SUM_ADDITION
#define M_OUT_LIE_SUM_ADDITION
#define IN_LIE_SUM_COPY
#define OUT_LIE_SUM_COPY
#define IN_LIE_SUM_KILL
#define OUT_LIE_SUM_KILL
#define M_IN_LIE_SUM_DIV_INTEGER
#define M_OUT_LIE_SUM_DIV_INTEGER
#define M_IN_LIE_SUM_DIV_SCALAR_SUM
#define M_OUT_LIE_SUM_DIV_SCALAR_SUM
#define M_IN_LIE_SUM_MULT_SCALAR_SUM
#define M_OUT_LIE_SUM_MULT_SCALAR_SUM
#define M_IN_LIE_SUM_MULT_INTEGER
#define M_OUT_LIE_SUM_MULT_INTEGER
#define M_OUT_MAKE_RELATION_RHS
#define M_OUT_NEW_JACOBI_RELATION
#define M_IN_NORMALIZE_RELATION
#define M_OUT_NORMALIZE_RELATION
#define M_OUT_PAIR_MONOMIAL_MONOMIAL
#define M_IN_PAIR_MONOMIAL_SUM
#define M_OUT_PAIR_MONOMIAL_SUM
#define M_IN_PAIR_SUM_MONOMIAL
#define M_OUT_PAIR_SUM_MONOMIAL
#define M_IN_PAIR_SUM_SUM
#define M_OUT_PAIR_SUM_SUM
#define M_IN_POLY_CONTENT
#define M_OUT_POLY_CONTENT
#define M_IN_POLY_GCD
#define M_OUT_POLY_GCD
#define M_IN_POLY_QUOTIENT
#define M_OUT_POLY_QUOTIENT
#define M_IN_POLY_PSEUDO_REMAINDER
#define M_OUT_POLY_PSEUDO_REMAINDER
#define IN_REDUCE_RELATIONS
#define OUT_REDUCE_RELATIONS
#define M_IN_SCALAR_SUM_CANCELLATION
#define M_OUT_SCALAR_SUM_CANCELLATION
#define M_IN_SUBSTITUTE_RELATION_IN_RELATION
#define M_OUT_SUBSTITUTE_RELATION_IN_RELATION
#endif

#define IN_SET_NS IN_SET_N_LT IN_SET_N_INT IN_SET_N_ST IN_SET_N_SF
#define OUT_SET_NS OUT_SET_N_LT OUT_SET_N_INT OUT_SET_N_ST OUT_SET_N_SF


#if defined(DEBUG) || defined(MEMORY)             /* `f_name' is unique */
#define IN_SET_FUNCTION_NAME(fname) char * f_name = fname;
#else
#define IN_SET_FUNCTION_NAME(fname)
#endif

#if defined(D_ADD_PAIR_TO_LIE_MONOMIAL) /*------------------------------*/
                                        /*   Only debugging makes sense */
#define IN_ADD_PAIR_TO_LIE_MONOMIAL \
    IN_SET_FUNCTION_NAME("AddPairToLieMonomial")\
    D_IN_SET\
      PutDebugLieMonomial("i", i),\
      PutDebugLieMonomial("j", j),\
      PutDebugLieMonomialTable(-1);\
    D_IN_CLOSE

#define OUT_ADD_PAIR_TO_LIE_MONOMIAL_OLD \
    D_OUT_OPEN\
      PutDebugLieMonomial("Old monomial", ijp);
#define OUT_ADD_PAIR_TO_LIE_MONOMIAL_NEW  \
    D_OUT_OPEN\
      PutDebugLieMonomialTable(ijp);
#else
#define IN_ADD_PAIR_TO_LIE_MONOMIAL     D_IN_CLOSE
#define OUT_ADD_PAIR_TO_LIE_MONOMIAL_OLD
#define OUT_ADD_PAIR_TO_LIE_MONOMIAL_NEW
#endif
#if defined(D_GENERATE_RELATIONS) /*-----------------------------------*/
                                       /*   Only debugging makes sense */
#define IN_GENERATE_RELATIONS  \
    IN_SET_FUNCTION_NAME("GenerateRelations")\
    D_IN_SET\
      IN_PUT_RELATIONS\
      IN_PUT_LIE_MONOMIAL\
      PutDebugLieSum("rel", a),\
      PutDebugLieMonomial("gen", gen);\
    D_IN_CLOSE

#define OUT_GENERATE_RELATIONS \
    D_OUT_OPEN\
      PutDebugLieSum("d_rel", a);
#else
#define  IN_GENERATE_RELATIONS     D_IN_CLOSE
#define  OUT_GENERATE_RELATIONS
#endif
#if defined(D_INTEGER_CANCELLATION) /*----- (k*n)/(k*d) -> n/d --------*/
                                    /* Only debugging makes sense      */
#define IN_INTEGER_CANCELLATION \
    IN_SET_FUNCTION_NAME("IntegerCancellation")\
    D_IN_SET\
      PutDebugInteger("num", num),\
      PutDebugInteger("den", den);\
    D_IN_CLOSE
#define OUT_INTEGER_CANCELLATION \
    D_OUT_OPEN\
      PutDebugInteger("num", num),\
      PutDebugInteger("den", den);
#else
#define IN_INTEGER_CANCELLATION    D_IN_CLOSE
#define OUT_INTEGER_CANCELLATION
#endif
#if defined(D_INTEGER_GCD) /*------------------------------------------*/
                                    /* Only debugging makes sense      */
#define IN_INTEGER_GCD  \
    IN_SET_FUNCTION_NAME("IntegerGCD")\
    D_IN_SET\
      PutDebugInteger("u", u),\
      PutDebugInteger("v", v);\
    D_IN_CLOSE
#define OUT_INTEGER_GCD  \
    D_OUT_OPEN\
      PutDebugInteger("GCD", u);
#else
#define IN_INTEGER_GCD       D_IN_CLOSE
#define OUT_INTEGER_GCD
#endif
#if defined(D_INTEGER_PRODUCT)  /*------------- n*m -------------------*/
                                    /* Only debugging makes sense      */
#define IN_INTEGER_PRODUCT  \
    IN_SET_FUNCTION_NAME("IntegerProduct")\
    D_IN_SET\
      PutDebugInteger("u", u),\
      PutDebugInteger("v", v);\
    D_IN_CLOSE
#define OUT_INTEGER_PRODUCT  \
    D_OUT_OPEN\
      PutDebugInteger("u*v", w0);
#else
#define IN_INTEGER_PRODUCT    D_IN_CLOSE
#define OUT_INTEGER_PRODUCT
#endif
#if defined(D_INTEGER_QUOTIENT)  /*------------------------------------*/
                                    /* Only debugging makes sense      */
#define IN_INTEGER_QUOTIENT  \
    IN_SET_FUNCTION_NAME("IntegerQuotient")\
    D_IN_SET\
      PutDebugInteger("a", a),\
      PutDebugInteger("b", b);\
    D_IN_CLOSE
#define OUT_INTEGER_QUOTIENT  \
    D_OUT_OPEN\
      PutDebugInteger("a/b", pm);
#else
#define IN_INTEGER_QUOTIENT    D_IN_CLOSE
#define OUT_INTEGER_QUOTIENT
#endif
#if defined(D_INTEGER_SUM)       /*------------------------------------*/
                                    /* Only debugging makes sense      */
#define IN_INTEGER_SUM  \
    IN_SET_FUNCTION_NAME("IntegerSum")\
    D_IN_SET\
      PutDebugInteger("a", a),\
      PutDebugInteger("b", b);\
    D_IN_CLOSE
#define OUT_INTEGER_SUM  \
    D_OUT_OPEN\
      PutDebugInteger("c", c);
#else
#define IN_INTEGER_SUM    D_IN_CLOSE
#define OUT_INTEGER_SUM
#endif

#if defined(D_GET_LIE_MONOMIAL) /*-------------------------------------*/
#define D_IN_GET_LIE_MONOMIAL   D_IN_SET\
                                  PutDebugString("*pstr", *pstr);\
                                D_IN_CLOSE
#define D_OUT_GET_LIE_MONOMIAL  D_OUT_OPEN\
                                  PutDebugLieSum("a", a);
#else
#define D_IN_GET_LIE_MONOMIAL   D_IN_CLOSE
#define D_OUT_GET_LIE_MONOMIAL
#endif
#if defined(D_GET_LIE_SUM) /*------------------------------------------*/
#define D_IN_GET_LIE_SUM     D_IN_SET\
                               PutDebugString("*pstr", *pstr);\
                             D_IN_CLOSE
#define D_OUT_GET_LIE_SUM    D_OUT_OPEN\
                               PutDebugLieSum("lsum", lsum);
#else
#define D_IN_GET_LIE_SUM     D_IN_CLOSE
#define D_OUT_GET_LIE_SUM
#endif
#if defined(D_GET_LIE_TERM) /*-----------------------------------------*/
#define D_IN_GET_LIE_TERM  D_IN_SET\
                             PutDebugString("*pstr", *pstr);\
                           D_IN_CLOSE
#define D_OUT_GET_LIE_TERM D_OUT_OPEN\
                             PutDebugLieSum("lterm", lterm);
#else
#define D_IN_GET_LIE_TERM  D_IN_CLOSE
#define D_OUT_GET_LIE_TERM
#endif
#if defined(D_LIE_SUM_ADDITION) /*---------- a + b (Lie) --------------*/
#define D_IN_LIE_SUM_ADDITION  D_IN_SET\
                                 PutDebugLieSum("a", a),\
                                 PutDebugLieSum("b", b);\
                               D_IN_CLOSE
#define D_OUT_LIE_SUM_ADDITION D_OUT_OPEN\
                                 PutDebugLieSum("sum", sum);
#else
#define D_IN_LIE_SUM_ADDITION  D_IN_CLOSE
#define D_OUT_LIE_SUM_ADDITION
#endif
#if defined(D_LIE_SUM_DIV_INTEGER)  /*--------- uint/BIGINT ----------------*/
#define D_IN_LIE_SUM_DIV_INTEGER   D_IN_SET\
                                     PutDebugLieSum("lsum", b),\
                                     PutDebugInteger("den", den);\
                                   D_IN_CLOSE
#define D_OUT_LIE_SUM_DIV_INTEGER  D_OUT_OPEN\
                                     PutDebugLieSum("lsum", b);
#else
#define D_IN_LIE_SUM_DIV_INTEGER   D_IN_CLOSE
#define D_OUT_LIE_SUM_DIV_INTEGER
#endif
#if defined(D_LIE_SUM_DIV_SCALAR_SUM)  /*--------- uint/uint --------------*/
#define D_IN_LIE_SUM_DIV_SCALAR_SUM   D_IN_SET\
                                        PutDebugLieSum("lsum", b),\
                                        PutDebugScalarSum("den", den);\
                                      D_IN_CLOSE
#define D_OUT_LIE_SUM_DIV_SCALAR_SUM  D_OUT_OPEN\
                                        PutDebugLieSum("lsum", b);
#else
#define D_IN_LIE_SUM_DIV_SCALAR_SUM   D_IN_CLOSE
#define D_OUT_LIE_SUM_DIV_SCALAR_SUM
#endif
#if defined(D_LIE_SUM_MULT_SCALAR_SUM)  /*--------- uint*uint --------------*/
#define D_IN_LIE_SUM_MULT_SCALAR_SUM   D_IN_SET\
                                         PutDebugLieSum("lsum", b),\
                                         PutDebugScalarSum("num", num);\
                                       D_IN_CLOSE
#define D_OUT_LIE_SUM_MULT_SCALAR_SUM  D_OUT_OPEN\
                                         PutDebugLieSum("lsum", b);
#else
#define D_IN_LIE_SUM_MULT_SCALAR_SUM   D_IN_CLOSE
#define D_OUT_LIE_SUM_MULT_SCALAR_SUM
#endif
#if defined(D_LIE_SUM_MULT_INTEGER)  /*------- uint*BIGINT -----------------*/
#define D_IN_LIE_SUM_MULT_INTEGER   D_IN_SET\
                                      PutDebugLieSum("lsum", b),\
                                      PutDebugInteger("num", num);\
                                    D_IN_CLOSE
#define D_OUT_LIE_SUM_MULT_INTEGER  D_OUT_OPEN\
                                      PutDebugLieSum("lsum", b);
#else
#define D_IN_LIE_SUM_MULT_INTEGER   D_IN_CLOSE
#define D_OUT_LIE_SUM_MULT_INTEGER
#endif
#if defined(D_MAKE_RELATION_RHS)   /*----------------------------------*/
#define D_IN_MAKE_RELATION_RHS     D_IN_SET\
                                     PutDebugLieSum("rel",\
                                         RELATION_LIE_SUM(i));\
                                   D_IN_CLOSE
#define D_OUT_MAKE_RELATION_RHS    D_OUT_OPEN\
                                     PutDebugLieSum("r.h.s", a);
#else
#define D_IN_MAKE_RELATION_RHS     D_IN_CLOSE
#define D_OUT_MAKE_RELATION_RHS
#endif
#if defined(D_NEW_JACOBI_RELATION) /*----------------------------------*/
#define D_IN_NEW_JACOBI_RELATION   D_IN_SET\
                                     PutDebugLieMonomial("l", l);\
                                   D_IN_CLOSE
#define D_OUT_NEW_JACOBI_RELATION  D_OUT_OPEN\
                                     PutDebugLieSum("a", a);
#else
#define D_IN_NEW_JACOBI_RELATION   D_IN_CLOSE
#define D_OUT_NEW_JACOBI_RELATION
#endif
#if defined(D_NORMALIZE_RELATION) /*-----------------------------------*/
#define D_IN_NORMALIZE_RELATION    D_IN_SET\
                                     PutDebugLieSum("a", a);\
                                   D_IN_CLOSE
#define D_OUT_NORMALIZE_RELATION   D_OUT_OPEN\
                                     PutDebugLieSum("a", a);
#else
#define D_IN_NORMALIZE_RELATION    D_IN_CLOSE
#define D_OUT_NORMALIZE_RELATION
#endif
#if defined(D_PAIR_MONOMIAL_MONOMIAL) /*---------[mona, monb]----------*/
#define D_IN_PAIR_MONOMIAL_MONOMIAL    D_IN_SET\
                                         PutDebugLieMonomial("i", i),\
                                         PutDebugLieMonomial("j", j);\
                                       D_IN_CLOSE
#define D_OUT_PAIR_MONOMIAL_MONOMIAL   D_OUT_OPEN\
                                         PutDebugLieSum("a", a);
#else
#define D_IN_PAIR_MONOMIAL_MONOMIAL    D_IN_CLOSE
#define D_OUT_PAIR_MONOMIAL_MONOMIAL
#endif
#if defined(D_PAIR_MONOMIAL_SUM) /*-----------[mon, a]-----------------*/
#define D_IN_PAIR_MONOMIAL_SUM  D_IN_SET\
                                  PutDebugLieMonomial("mon", mon),\
                                  PutDebugLieSum("a", a);\
                                D_IN_CLOSE
#define D_OUT_PAIR_MONOMIAL_SUM D_OUT_OPEN\
                                  PutDebugLieSum("s", s);
#else
#define D_IN_PAIR_MONOMIAL_SUM  D_IN_CLOSE
#define D_OUT_PAIR_MONOMIAL_SUM
#endif
#if defined(D_PAIR_SUM_MONOMIAL) /*-----------[a, mon]-----------------*/
#define D_IN_PAIR_SUM_MONOMIAL  D_IN_SET\
                                  PutDebugLieSum("a", a),\
                                  PutDebugLieMonomial("mon", mon);\
                                D_IN_CLOSE
#define D_OUT_PAIR_SUM_MONOMIAL D_OUT_OPEN\
                                  PutDebugLieSum("s", s);
#else
#define D_IN_PAIR_SUM_MONOMIAL  D_IN_CLOSE
#define D_OUT_PAIR_SUM_MONOMIAL
#endif
#if defined(D_PAIR_SUM_SUM)    /*----------------[a, b]---------------*/
#define D_IN_PAIR_SUM_SUM  D_IN_SET\
                             PutDebugLieSum("a", a),\
                             PutDebugLieSum("b", b);\
                           D_IN_CLOSE
#define D_OUT_PAIR_SUM_SUM D_OUT_OPEN\
                             PutDebugLieSum("a", a);
#else
#define D_IN_PAIR_SUM_SUM    D_IN_CLOSE
#define D_OUT_PAIR_SUM_SUM
#endif
#if defined(D_POLY_CONTENT) /*--------------------------------------------*/
#define D_IN_POLY_CONTENT  D_IN_SET\
                             PutDebugScalarSum("a", a),\
                             PutDebugString("mp", \
                                ParameterName + mp*NameLength1);\
                           D_IN_CLOSE
#define D_OUT_POLY_CONTENT D_OUT_OPEN\
                             PutDebugScalarSum("b", b);
#else
#define D_IN_POLY_CONTENT  D_IN_CLOSE
#define D_OUT_POLY_CONTENT
#endif
#if defined(D_POLY_GCD) /*--------------------------------------------*/
#define D_IN_POLY_GCD  D_IN_SET\
                         PutDebugScalarSum("a", a),\
                         PutDebugScalarSum("b", b);\
                       D_IN_CLOSE
#define D_OUT_POLY_GCD D_OUT_OPEN\
                         PutDebugScalarSum("b", b);
#else
#define D_IN_POLY_GCD    D_IN_CLOSE
#define D_OUT_POLY_GCD
#endif
#if defined(D_POLY_QUOTIENT) /*--------------------------------------------*/
#define D_IN_POLY_QUOTIENT  D_IN_SET\
                              PutDebugScalarSum("a", a),\
                              PutDebugScalarSum("b", b);\
                            D_IN_CLOSE
#define D_OUT_POLY_QUOTIENT D_OUT_OPEN\
                              PutDebugScalarSum("c", c);
#else
#define D_IN_POLY_QUOTIENT    D_IN_CLOSE
#define D_OUT_POLY_QUOTIENT
#endif
#if defined(D_POLY_PSEUDO_REMAINDER) /*------------------------------------*/
#define D_IN_POLY_PSEUDO_REMAINDER  D_IN_SET\
                                      PutDebugScalarSum("a", a),\
                                      PutDebugScalarSum("b", b);\
                                    D_IN_CLOSE
#define D_OUT_POLY_PSEUDO_REMAINDER D_OUT_OPEN\
                                      PutDebugScalarSum("a", a);
#else
#define D_IN_POLY_PSEUDO_REMAINDER  D_IN_CLOSE
#define D_OUT_POLY_PSEUDO_REMAINDER
#endif
#if defined(D_SCALAR_SUM_CANCELLATION) /*----- (k*n)/(k*d) -> n/d ----*/
#define D_IN_SCALAR_SUM_CANCELLATION  D_IN_SET\
                                        PutDebugScalarSum("*pnum", *pnum),\
                                        PutDebugScalarSum("*pden", *pden);\
                                      D_IN_CLOSE
#define D_OUT_SCALAR_SUM_CANCELLATION D_OUT_OPEN\
                                        PutDebugScalarSum("*pnum", *pnum),\
                                        PutDebugScalarSum("*pden", *pden);
#else
#define D_IN_SCALAR_SUM_CANCELLATION    D_IN_CLOSE
#define D_OUT_SCALAR_SUM_CANCELLATION
#endif
#if defined(D_SUBSTITUTE_RELATION_IN_RELATION)  /*--------------------*/
#define D_IN_SUBSTITUTE_RELATION_IN_RELATION  D_IN_SET\
                                                PutDebugLieSum("r", r),\
                                                PutDebugLieSum("a", a);\
                                              D_IN_CLOSE
#define D_OUT_SUBSTITUTE_RELATION_IN_RELATION  D_OUT_OPEN\
                                                 PutDebugLieSum("a", a);
#else
#define D_IN_SUBSTITUTE_RELATION_IN_RELATION  D_IN_CLOSE
#define D_OUT_SUBSTITUTE_RELATION_IN_RELATION
#endif

/* Set INs and OUTs for particular functions */

/*--------------------------------------------------------------------*/
#define IN_GET_LIE_MONOMIAL \
    IN_SET_FUNCTION_NAME("GetLieMonomial")\
    IN_SET_NS\
    D_IN_GET_LIE_MONOMIAL
#define OUT_GET_LIE_MONOMIAL \
    D_OUT_GET_LIE_MONOMIAL\
    M_OUT_GET_LIE_MONOMIAL
/*--------------------------------------------------------------------*/
#define IN_GET_LIE_SUM \
    IN_SET_FUNCTION_NAME("GetLieSum")\
    IN_SET_NS\
    D_IN_GET_LIE_SUM
#define OUT_GET_LIE_SUM \
    D_OUT_GET_LIE_SUM\
    M_OUT_GET_LIE_SUM
/*--------------------------------------------------------------------*/
#define IN_GET_LIE_TERM \
    IN_SET_FUNCTION_NAME("GetLieTerm")\
    IN_SET_NS\
    D_IN_GET_LIE_TERM
#define OUT_GET_LIE_TERM \
    D_OUT_GET_LIE_TERM\
    M_OUT_GET_LIE_TERM
/*--------------------------------------------------------------------*/
#define IN_LIE_SUM_ADDITION  \
    IN_SET_FUNCTION_NAME("LieSumAddition")\
    IN_SET_NS\
    D_IN_LIE_SUM_ADDITION\
    M_IN_LIE_SUM_ADDITION
#define OUT_LIE_SUM_ADDITION  \
    D_OUT_LIE_SUM_ADDITION\
    M_OUT_LIE_SUM_ADDITION
/*--------------------------------------------------------------------*/
#define IN_LIE_SUM_DIV_INTEGER  \
    LS_B_LSUM_DIV_INT\
    IN_SET_FUNCTION_NAME("LieSumDivInteger")\
    IN_SET_NS\
    D_IN_LIE_SUM_DIV_INTEGER\
    M_IN_LIE_SUM_DIV_INTEGER
#define OUT_LIE_SUM_DIV_INTEGER  \
    D_OUT_LIE_SUM_DIV_INTEGER\
    M_OUT_LIE_SUM_DIV_INTEGER
/*--------------------------------------------------------------------*/
#define IN_LIE_SUM_DIV_SCALAR_SUM  \
    LS_B_LSUM_DIV_SS\
    IN_SET_FUNCTION_NAME("LieSumDivScalarSum")\
    IN_SET_NS\
    D_IN_LIE_SUM_DIV_SCALAR_SUM\
    M_IN_LIE_SUM_DIV_SCALAR_SUM
#define OUT_LIE_SUM_DIV_SCALAR_SUM  \
    D_OUT_LIE_SUM_DIV_SCALAR_SUM\
    M_OUT_LIE_SUM_DIV_SCALAR_SUM
/*--------------------------------------------------------------------*/
#define IN_LIE_SUM_MULT_SCALAR_SUM  \
    LS_B_LSUM_MULT_SS\
    IN_SET_FUNCTION_NAME("LieSumMultScalarSum")\
    IN_SET_NS\
    D_IN_LIE_SUM_MULT_SCALAR_SUM\
    M_IN_LIE_SUM_MULT_SCALAR_SUM
#define OUT_LIE_SUM_MULT_SCALAR_SUM  \
    D_OUT_LIE_SUM_MULT_SCALAR_SUM\
    M_OUT_LIE_SUM_MULT_SCALAR_SUM
/*--------------------------------------------------------------------*/
#define IN_LIE_SUM_MULT_INTEGER  \
    LS_B_LSUM_MULT_INT\
    IN_SET_FUNCTION_NAME("LieSumMultInteger")\
    IN_SET_NS\
    D_IN_LIE_SUM_MULT_INTEGER\
    M_IN_LIE_SUM_MULT_INTEGER
#define OUT_LIE_SUM_MULT_INTEGER  \
    D_OUT_LIE_SUM_MULT_INTEGER\
    M_OUT_LIE_SUM_MULT_INTEGER
/*--------------------------------------------------------------------*/
#define IN_MAKE_RELATION_RHS \
    IN_SET_FUNCTION_NAME("MakeRelationRHS...")\
    IN_SET_NS\
    D_IN_MAKE_RELATION_RHS
#define OUT_MAKE_RELATION_RHS \
    D_OUT_MAKE_RELATION_RHS\
    M_OUT_MAKE_RELATION_RHS
/*--------------------------------------------------------------------*/
#define IN_NEW_JACOBI_RELATION  \
    IN_SET_FUNCTION_NAME("NewJacobiRelation")\
    IN_SET_NS\
    D_IN_NEW_JACOBI_RELATION
#define OUT_NEW_JACOBI_RELATION  \
    D_OUT_NEW_JACOBI_RELATION\
    M_OUT_NEW_JACOBI_RELATION
/*--------------------------------------------------------------------*/
#define IN_NORMALIZE_RELATION  \
    IN_SET_FUNCTION_NAME("NormalizeRelation...")\
    IN_SET_NS\
    D_IN_NORMALIZE_RELATION\
    M_IN_NORMALIZE_RELATION
#define OUT_NORMALIZE_RELATION  \
    D_OUT_NORMALIZE_RELATION\
    M_OUT_NORMALIZE_RELATION
/*--------------------------------------------------------------------*/
#define IN_PAIR_MONOMIAL_MONOMIAL  \
    IN_SET_FUNCTION_NAME("PairMonomialMonomial...")\
    IN_SET_NS\
    D_IN_PAIR_MONOMIAL_MONOMIAL
#define OUT_PAIR_MONOMIAL_MONOMIAL  \
    D_OUT_PAIR_MONOMIAL_MONOMIAL\
    M_OUT_PAIR_MONOMIAL_MONOMIAL
/*--------------------------------------------------------------------*/
#define IN_PAIR_MONOMIAL_SUM  \
    IN_SET_FUNCTION_NAME("PairMonomialSum...")\
    IN_SET_NS\
    D_IN_PAIR_MONOMIAL_SUM\
    M_IN_PAIR_MONOMIAL_SUM
#define OUT_PAIR_MONOMIAL_SUM  \
    D_OUT_PAIR_MONOMIAL_SUM\
    M_OUT_PAIR_MONOMIAL_SUM
/*--------------------------------------------------------------------*/
#define IN_PAIR_SUM_MONOMIAL  \
    IN_SET_FUNCTION_NAME("PairSumMonomial...")\
    IN_SET_NS\
    D_IN_PAIR_SUM_MONOMIAL\
    M_IN_PAIR_SUM_MONOMIAL
#define OUT_PAIR_SUM_MONOMIAL  \
    D_OUT_PAIR_SUM_MONOMIAL\
    M_OUT_PAIR_SUM_MONOMIAL
/*--------------------------------------------------------------------*/
#define IN_PAIR_SUM_SUM  \
    IN_SET_FUNCTION_NAME("PairSumSum...")\
    IN_SET_NS\
    D_IN_PAIR_SUM_SUM\
    M_IN_PAIR_SUM_SUM
#define OUT_PAIR_SUM_SUM  \
    D_OUT_PAIR_SUM_SUM\
    M_OUT_PAIR_SUM_SUM
/*--------------------------------------------------------------------*/
#define IN_POLY_CONTENT  \
    IN_SET_FUNCTION_NAME("PolyContent")\
    IN_SET_NS\
    D_IN_POLY_CONTENT\
    M_IN_POLY_CONTENT
#define OUT_POLY_CONTENT  \
    D_OUT_POLY_CONTENT\
    M_OUT_POLY_CONTENT
/*--------------------------------------------------------------------*/
#define IN_POLY_GCD  \
    IN_SET_FUNCTION_NAME("PolyGCD")\
    IN_SET_NS\
    D_IN_POLY_GCD\
    M_IN_POLY_GCD
#define OUT_POLY_GCD  \
    D_OUT_POLY_GCD\
    M_OUT_POLY_GCD
/*--------------------------------------------------------------------*/
#define IN_POLY_QUOTIENT  \
    IN_SET_FUNCTION_NAME("PolyQuotient")\
    IN_SET_NS\
    D_IN_POLY_QUOTIENT\
    M_IN_POLY_QUOTIENT
#define OUT_POLY_QUOTIENT  \
    D_OUT_POLY_QUOTIENT\
    M_OUT_POLY_QUOTIENT
/*--------------------------------------------------------------------*/
#define IN_POLY_PSEUDO_REMAINDER  \
    IN_SET_FUNCTION_NAME("PolyPseudoRemainder")\
    IN_SET_NS\
    D_IN_POLY_PSEUDO_REMAINDER\
    M_IN_POLY_PSEUDO_REMAINDER
#define OUT_POLY_PSEUDO_REMAINDER  \
    D_OUT_POLY_PSEUDO_REMAINDER\
    M_OUT_POLY_PSEUDO_REMAINDER
/*--------------------------------------------------------------------*/
#define IN_SCALAR_SUM_CANCELLATION  \
    IN_SET_FUNCTION_NAME("ScalarSumCancellation")\
    IN_SET_NS\
    D_IN_SCALAR_SUM_CANCELLATION\
    M_IN_SCALAR_SUM_CANCELLATION
#define OUT_SCALAR_SUM_CANCELLATION  \
    D_OUT_SCALAR_SUM_CANCELLATION\
    M_OUT_SCALAR_SUM_CANCELLATION
/*--------------------------------------------------------------------*/
#define IN_SUBSTITUTE_RELATION_IN_RELATION  \
    IN_SET_FUNCTION_NAME("SubstituteRelationInRelation...")\
    IN_SET_NS\
    D_IN_SUBSTITUTE_RELATION_IN_RELATION\
    M_IN_SUBSTITUTE_RELATION_IN_RELATION
#define OUT_SUBSTITUTE_RELATION_IN_RELATION  \
    D_OUT_SUBSTITUTE_RELATION_IN_RELATION\
    M_OUT_SUBSTITUTE_RELATION_IN_RELATION
/* Conditional print of Relation and Monomial tables */

#if defined(D_PUT_RELATIONS)   /*-------------------------------------*/
#define IN_PUT_RELATIONS       PutDebugRelations(),
#define OUT_PUT_RELATIONS      D_CONDITION PutDebugRelations();
#else
#define IN_PUT_RELATIONS
#define OUT_PUT_RELATIONS
#endif

#if defined(D_PUT_LIE_MONOMIAL)   /*-------------------------------------*/
#define IN_PUT_LIE_MONOMIAL       PutDebugLieMonomialTable(-1),
#define OUT_PUT_LIE_MONOMIAL      D_CONDITION PutDebugLieMonomialTable(-1);
#else
#define IN_PUT_LIE_MONOMIAL
#define OUT_PUT_LIE_MONOMIAL
#endif

#if defined(D_LIE_SUM_DIV_INTEGER) || defined(MEMORY)
#define LS_B_LSUM_DIV_INT  uint b = lsum; /* For "LieSumDivInteger" */
#else
#define LS_B_LSUM_DIV_INT
#endif
#if defined(D_LIE_SUM_MULT_INTEGER) || defined(MEMORY)
#define LS_B_LSUM_MULT_INT  uint b = lsum; /* For "LieSumMultInteger" */
#else
#define LS_B_LSUM_MULT_INT
#endif
#if defined(D_LIE_SUM_DIV_SCALAR_SUM) || defined(MEMORY)
#define LS_B_LSUM_DIV_SS  uint b = lsum; /* For "LieSumDivScalarSum" */
#else
#define LS_B_LSUM_DIV_SS
#endif
#if defined(D_LIE_SUM_MULT_SCALAR_SUM) || defined(MEMORY)
#define LS_B_LSUM_MULT_SS  uint b = lsum; /* For "LieSumMultScalarSum" */
#else
#define LS_B_LSUM_MULT_SS
#endif

/* Test of char functions ?? */
//#define TEST_FUNCTION

#if defined(INTEGER_ALLOCATION_CHECK)
#define INTEGER_IN_STACK(n) ;if((n)==NULL) Error(E_A_STACK_INTEGER)
#define INTEGER_IN_HEAP(n)  ;if((n)==NULL) Error(E_A_HEAP_INTEGER) \
                            PP_CURRENT_N_INT
#else
#define INTEGER_IN_STACK(n)
#define INTEGER_IN_HEAP(n)  PP_CURRENT_N_INT
#endif

#if defined(POLY_ARRAY_ALLOCATION_CHECK)
#define POLY_ARRAY_IN_STACK(a) ;if((a)==NULL) Error(E_A_STACK_POLY_ARRAY)
#else
#define POLY_ARRAY_IN_STACK(a)
#endif

/*_1    System constants=================================================*/

#define NIL                                           (0u)
#define NOTHING                                     (~NIL)

#define INTEGER_SIGN_MASK                   ((LIMB)0x8000)
                      /* 1 << 15 or 0000 0000 0000 0001 */
#define INTEGER_N_LIMBS_MASK  ((LIMB)(~INTEGER_SIGN_MASK))
                                 /* 1111 1111 1111 1110 */
#define BITS_PER_LIMB                                 (16)
#define BASE_LIMB                              (0x10000lu)
#define MAX_LIMB                                 (0xFFFFu)
#define FIRST_GENUINE_PARAMETER  1 /* Skip i number */

/*_2    Type definitions================================================*/

typedef unsigned char byte;
typedef unsigned short LIMB; /* Limb of big integer */
typedef LIMB * BIGINT;  /* (Pointer to) big integer array */
typedef unsigned int uint;

/* Element of LieMonomial table */

typedef struct
{
  int order;       /* Order of this position element */
  int position;    /* Position of this order */
  int left;        /* (Position of) left submonomial of Lie bracket */
                   /* or ~(index of generator) if not commutator */
  int right;       /* (Position of) right submonomial of Lie bracket  */
                   /* or 1 for 1st generator and 0 for nexts */
  struct
  {
    uint parity : 1;
    int index : 23;  /* Index if basis element or */
                   /* ~(index of relation with leading ordinal) */
                   /* Number of relations < 4194303 */
    uint weight : 8;  /* Weight of Lie monomial < 256 */
  } info;
} LIE_MON;

/* Element of NodeLT pool for Lie terms */

typedef struct
{
  int monomial;   /* (Position of) Lie monomial in LieMonomial table */
  union
  {
    uint  scalar_sum;
    BIGINT integer;
  } numerator;
  union
  {
    uint  scalar_sum;
    BIGINT integer;
  } denominator;  /* Pointer to next Lie term */
  uint rptr;
} NODE_LT;

/* Element of NodeSF pool for scalar factors */

typedef struct
{
#if defined(SPP_2000)
  byte parameter; /* Parameter ordinal   */
  byte degree;    /* Degree of parameter */
#else
  byte degree;    /* Degree of parameter */
  byte parameter; /* Parameter ordinal   */
#endif
  uint   rptr;      /* Pointer to next parametric factor */
} NODE_SF;

/* Element of NodeST pool for scalar terms */

typedef struct
{
  uint  monomial;  /* Scalar monomial */
  BIGINT numerator; /* Integer coefficient */
  uint  rptr;      /* Pointer to next parametric term */
} NODE_ST;

/* Element of Relation table */

typedef struct
{
  uint    lie_sum;           /* Expression of relation (Lie sum) */
  byte  min_generator;     /* Minimal generator for differentiation */
  byte to_be_substituted; /* YES if relation must be substituted */
} REL;                     /* into higher relations */

/*_3    Constants and enumerations======================================*/

/* Enumeration constants */

enum boolean {NO = 0, YES = 1};
enum signs {PLUS = 0, MINUS = 1};
enum parity {EVEN = 0, ODD = 1};
enum orders {ORDER12 = -1, ORDER11 = 0, ORDER21 = 1};
enum scalar_types   /* Types of scalar factors */
{
  I_NUMBER  =  0
};
enum init_file_cases   /* Initiating file fplsa4.ini (fplsa416.ini) */
{
  COEFFICIENT_SUM_TABLE_SIZE =  0,
  CRUDE_TIME                 =  1,
  ECHO_INPUT_FILE            =  2,
  EVEN_BASIS_SYMBOL          =  3,
  GAP_ALGEBRA_NAME           =  4,
  GAP_BASIS_NAME             =  5,
  GAP_RELATIONS_NAME         =  6,
  GAP_OUTPUT_BASIS           =  7,
  GAP_OUTPUT_COMMUTATORS     =  8,
  GAP_OUTPUT_RELATIONS       =  9,
  GENERATOR_MAX_N            = 10,
  INPUT_DIRECTORY            = 11,
  INPUT_INTEGER_SIZE         = 12,
  INPUT_STRING_SIZE          = 13,
  LEFT_NORMED_OUTPUT         = 14,
  LIE_MONOMIAL_SIZE          = 15,
  LINE_LENGTH                = 16,
  NAME_LENGTH                = 17,
  NODE_LT_SIZE               = 18,
  NODE_SF_SIZE               = 19,
  NODE_ST_SIZE               = 20,
  ODD_BASIS_SYMBOL           = 21,
  OUT_LINE_SIZE              = 22,
  PARAMETER_MAX_N            = 23,
  PUT_BASIS_ELEMENTS         = 24,
  PUT_COMMUTATORS            = 25,
  PUT_HILBERT_SERIES         = 26,
  PUT_INITIAL_RELATIONS      = 27,
  PUT_NON_ZERO_COEFFICIENTS  = 28,
  PUT_PROGRAM_HEADING        = 29,
  PUT_REDUCED_RELATIONS      = 30,
  PUT_STATISTICS             = 31,
  RELATION_SIZE              = 32,
  N_INIT_CASES               = 33
};
enum input_file_cases   /* Items of input files */
{
  GENERATORS      = 0,
  LIMITING_WEIGHT = 1,
  PARAMETERS      = 2,
  RELATIONS       = 3,
  WEIGHTS         = 4,
  N_INPUT_CASES   = 5
};
enum messages
{
  /* Head messages */

  H_PROGRAM               =  0,
  H_ENTER_FILE            =  1,
  H_INPUT_FILE            =  2,
  H_CREATE_NEW_FILE       =  3,
  H_ENTER_GENERATORS      =  4,
  H_ENTER_WEIGHTS_IN_FILE =  5,
  H_ENTER_LIMITING_WEIGHT =  6,
  H_ENTER_PARAMETERS      =  7,
  H_ENTER_RELATIONS       =  8,
  H_SHOW_INPUT            =  9,
  H_IN_RELATIONS          = 10,
  H_NON_ZERO_COEFFICIENTS = 11,
  H_REDUCED_RELATIONS     = 12,
  H_BASIS_ELEMENTS        = 13,
  H_HILBERT_SERIES        = 14,
  H_COMMUTATORS           = 15,
  H_NO_PUT_COMMUTATORS    = 16,

  /* Error messages */

  ERROR                            = 17,
  E_WRONG_INI_CASE                 = 18,
  E_WRONG_INPUT_CASE               = 19,
  E_CANCEL_PROGRAM                 = 20,
  E_UNEXPECTED_EOF                 = 21,
  E_A_GENERATOR_NAME               = 22,
  E_A_PARAMETER_NAME               = 23,
  E_A_OUT_LINE                     = 24,
  E_A_RELATION                     = 25,
  E_A_LIE_MONOMIAL                 = 26,
  E_A_NODE_LT                      = 27,
  E_A_NODE_ST                      = 28,
  E_A_NODE_SF                      = 29,
  E_A_HEAP_INTEGER                 = 30,
  E_A_COEFF_PARA_TABLE             = 31,
  E_A_COEFF_SUM_TABLE              = 32,
  E_ALLOC                          = 33,
  E_A_STACK_INPUT_STRING           = 34,
  E_A_STACK_INTEGER                = 35,
  E_A_STACK_INTEGER_DECIMAL_STRING = 36,
  E_A_STACK_POLY_ARRAY             = 37,
  E_INPUT_STRING_SIZE              = 38,
  E_OUT_LINE_SIZE                  = 39,
  E_LIE_MONOMIAL_SIZE              = 40,
  E_RELATION_SIZE                  = 41,
  E_NODE_LT_SIZE                   = 42,
  E_NODE_SF_SIZE                   = 43,
  E_NODE_ST_SIZE                   = 44,
  E_COEFF_SUM_TABLE_SIZE           = 45,
  E_GENERATOR_MAX_N                = 46,
  E_PARAMETER_MAX_N                = 47,
  E_TOO_MUCH_INPUT_WEIGHTS         = 48,
  E_NON_NUM_INPUT_WEIGHT           = 49,
  E_NO_R_PARENTHESIS               = 50,
  E_NO_GENERAL_POWER               = 51,
  E_UNDECLARED_GENERATOR           = 52,
  E_NO_COMMUTATOR_COMMA            = 53,
  E_NO_COMMUTATOR_BRACKET          = 54,
  E_INVALID_CHARACTER              = 55,
  E_MESSAGE                        = 56
};

/* Constants for input and output */

#define LEFT_COMMENT              '<'  /* In *.ini and input files */
#define RIGHT_COMMENT             '>'
#define SUBSCRIPT_INPUT_SIGN   '_'
#define ODD_GENERATOR_INPUT_SIGN  '-'  /* At input */

#define LEVEL                   '\xFF'
#define MAIN_LEVEL              '\x2'
#define MARGIN                  (LEVEL-1)

#define CASE_STRING_SIZE      256    /* Size of string to match case */
#define GAP_NAME_SIZE          64    /* Including ending '\0' */
#define GAP_WIDTH              79    /* Width of GAP page */

/*_4    Macrodefinitions==================================================*/
#define COUNT_LEADING_ZERO_BITS_IN_LIMB(n,w) n=CountLeadingZeroBitsInLimb((LIMB)(w))
#if 0 /* ?? Exclude the macro: the above function is slightly faster */
#define COUNT_LEADING_ZERO_BITS_IN_LIMB(n, w)  (n) = (w);\
                                          if((n) >= 0x100) \
                                            if((n) >= 0x1000) \
                                              if((n) >= 0x4000) \
                                                if((n) >= 0x8000) \
                                                  (n) = 0; \
                                                else \
                                                  (n) = 1; \
                                              else \
                                                if((n) >= 0x2000) \
                                                  (n) = 2; \
                                                else \
                                                  (n) = 3; \
                                            else \
                                              if((n) >= 0x400) \
                                                if((n) >= 0x800) \
                                                  (n) = 4; \
                                                else \
                                                  (n) = 5; \
                                              else \
                                                if((n) >= 0x200) \
                                                  (n) = 6; \
                                                else \
                                                  (n) = 7; \
                                          else \
                                            if((n) >= 0x10) \
                                              if((n) >= 0x40) \
                                                if((n) >= 0x80) \
                                                  (n) = 8; \
                                                else \
                                                  (n) = 9; \
                                              else \
                                                if((n) >= 0x20) \
                                                  (n) = 10; \
                                                else \
                                                  (n) = 11; \
                                            else \
                                              if((n) >= 0x4) \
                                                if((n) >= 0x8) \
                                                  (n) = 12; \
                                                else \
                                                  (n) = 13; \
                                              else \
                                                if((n) >= 0x2) \
                                                  (n) = 14; \
                                                else \
                                                  if((n)) \
                                                    (n) = 15; \
                                                  else \
                                                    (n) = 16
#endif

#define CUT_ARRAY(arr, type, n) (arr)=(type *)realloc(arr,sizeof(type)*(n))

#define EXIT     do { TIME_OFF; PutStatistics(); exit(1); } while(0)

#define IN_LINE_MARGIN    OutLine[++PosOutLine]=MARGIN

#define INTEGER_MINUS(ia)  do {  if(INTEGER_IS_NEGATIVE(ia))\
                                   (ia)[0] &= INTEGER_N_LIMBS_MASK;\
                                 else\
                                   (ia)[0] |= INTEGER_SIGN_MASK; } while(0)
#define INTEGER_IS_NEGATIVE(ia)  (((ia)[0]&INTEGER_SIGN_MASK)!=0)
#define INTEGER_IS_POSITIVE(ia)  (((ia)[0]&INTEGER_SIGN_MASK)==0)

#define INTEGER_IS_UNIT(ia)  (((ia)[0]==1) && ((ia)[1]==1))
#define INTEGER_IS_UNIT_ABS(ia)  ((((ia)[0]&INTEGER_N_LIMBS_MASK)==1) &&\
                              ((ia)[1]==1))
#define INTEGER_IS_NOT_UNIT(ia)  (((ia)[0]!=1) || ((ia)[1]!=1))
#define INTEGER_IS_NOT_UNIT_ABS(ia)  ((((ia)[0]&INTEGER_N_LIMBS_MASK)!=1) ||\
                                  ((ia)[1]!=1))
#define INTEGER_N_LIMBS(ia)  ((ia)[0]&INTEGER_N_LIMBS_MASK)

#define INTEGER_SET_MINUS(ia)       (ia)[0] |= INTEGER_SIGN_MASK
#define INTEGER_SET_PLUS(ia)        (ia)[0] &= INTEGER_N_LIMBS_MASK
#define INTEGER_SIGN(ia)            ((ia)[0]&INTEGER_SIGN_MASK)

#define INTEGER_HEAP_NEW(n,i)       (n)=(BIGINT)malloc(sizeof(LIMB)*(i))\
                                    INTEGER_IN_HEAP(n)
#define INTEGER_HEAP_COPY(n,o,i)    (i)=INTEGER_N_LIMBS(o);\
                                    INTEGER_HEAP_NEW(n,++(i));\
                                    do{(i)--; (n)[i] = (o)[i];}while(i)
#define INTEGER_HEAP_COPY_DOUBLE_1(n1,n2,o,i) \
                                    (i)=INTEGER_N_LIMBS(o);\
                                    INTEGER_HEAP_NEW(n1,++(i)+1);\
                                    INTEGER_HEAP_NEW(n2,i);\
                                    do{(i)--; (n1)[i] = (n2)[i] = (o)[i];}\
                                    while(i)
#define INTEGER_KILL(bn)            free(bn) MM_CURRENT_N_INT


#define INTEGER_STACK_NEW(n,i)      (n)=(BIGINT)alloca(sizeof(LIMB)*(i))\
                                    INTEGER_IN_STACK(n)

#define INTEGER_STACK_COPY(n,o,i)   (i)=INTEGER_N_LIMBS(o);\
                                    INTEGER_STACK_NEW(n,++(i));\
                                    do{(i)--; (n)[i] = (o)[i];}while(i)

#define INTEGER_STACK_COPY_1(n,o,i) (i)=INTEGER_N_LIMBS(o);\
                                    INTEGER_STACK_NEW(n,++(i)+1);\
                                    do{(i)--; (n)[i] = (o)[i];}while(i)

#define LIE_MONOMIAL_I_RELATION(pos)    (~LIE_MONOMIAL_INDEX(pos))
#define LIE_MONOMIAL_INDEX_BY_ORDER(ord) \
                         LIE_MONOMIAL_INDEX(LIE_MONOMIAL_POSITION(ord))
#define LIE_MONOMIAL_IS_EVEN(pos)       (LIE_MONOMIAL_PARITY(pos)==EVEN)
#define LIE_MONOMIAL_IS_LEADING_BY_ORDER(ord) \
                    LIE_MONOMIAL_IS_LEADING(LIE_MONOMIAL_POSITION(ord))
#define LIE_MONOMIAL_IS_OCCUPIED(pos)   LIE_MONOMIAL_WEIGHT(pos)
#define LIE_MONOMIAL_IS_GENERATOR(pos)  ((pos) < GeneratorN)
#define LIE_MONOMIAL_IS_GENERATOR_BY_ORDER(ord) \
                              (LIE_MONOMIAL_POSITION(ord) < GeneratorN)
#define LIE_MONOMIAL_IS_COMMUTATOR(pos) ((pos) >= GeneratorN)
#define LIE_MONOMIAL_IS_ODD(pos)        LIE_MONOMIAL_PARITY(pos)
#define LIE_MONOMIAL_IS_SQUARE(pos)     (LIE_MONOMIAL_LEFT(pos)==\
                                         LIE_MONOMIAL_RIGHT(pos))
#define LIE_MONOMIAL_IS_NOT_SQUARE(pos)  (LIE_MONOMIAL_LEFT(pos)!=\
                                          LIE_MONOMIAL_RIGHT(pos))
#define LIE_MONOMIAL_LEFT(pos)          (LieMonomial[pos].left)
#define LIE_MONOMIAL_RIGHT(pos)         (LieMonomial[pos].right)
#define LIE_MONOMIAL_ORDER(pos)         (LieMonomial[pos].order)
#define LIE_MONOMIAL_LEFT_ORDER(pos)   \
                      LIE_MONOMIAL_ORDER(LIE_MONOMIAL_LEFT(pos))
#define LIE_MONOMIAL_RIGHT_ORDER(pos)   \
                      LIE_MONOMIAL_ORDER(LIE_MONOMIAL_RIGHT(pos))
#define LIE_MONOMIAL_POSITION(ord)      (LieMonomial[ord].position)
#define LIE_MONOMIAL_WEIGHT_BY_ORDER(ord) \
                      LIE_MONOMIAL_WEIGHT(LIE_MONOMIAL_POSITION(ord))
#define LIE_MONOMIAL_INDEX(pos)         (LieMonomial[pos].info.index)
#define LIE_MONOMIAL_IS_BASIS(pos)      (LieMonomial[pos].info.index >= 0)
#define LIE_MONOMIAL_IS_LEADING(pos)    (LieMonomial[pos].info.index < 0)
#define LIE_MONOMIAL_PARITY(pos)        (LieMonomial[pos].info.parity)
#define LIE_MONOMIAL_WEIGHT(pos)        (LieMonomial[pos].info.weight)

#define LIE_TERM_MONOMIAL(a)               (NodeLT[a].monomial)
#define LIE_TERM_MONOMIAL_ORDER(a)         LIE_MONOMIAL_ORDER(\
                                             (NodeLT[a].monomial))
#define LIE_TERM_R(a)                      (NodeLT[a].rptr)
#define LIE_TERM_NUMERATOR_INTEGER(a)      (NodeLT[a].numerator.integer)
#define LIE_TERM_MINUS_INTEGER(a) \
                        INTEGER_MINUS(NodeLT[a].numerator.integer)
#define LIE_TERM_NUMERATOR_SCALAR_SUM(a)   (NodeLT[a].numerator.scalar_sum)
#define LIE_TERM_DENOMINATOR_INTEGER(a)    (NodeLT[a].denominator.integer)
#define LIE_TERM_DENOMINATOR_SCALAR_SUM(a) (NodeLT[a].denominator.scalar_sum)

#define MAX(a,b)    (((a) > (b)) ? (a) : (b))

#define NODE_LT_KILL(a)    LIE_TERM_R(a)=NodeLTTop,NodeLTTop=(a)\
                                                     MM_CURRENT_N_LT
#define NODE_SF_KILL(a)    SCALAR_FACTOR_R(a)=NodeSFTop,NodeSFTop=(a)\
                                                        MM_CURRENT_N_SF
#define NODE_ST_KILL(a)    SCALAR_TERM_R(a)=NodeSTTop,NodeSTTop=(a)\
                                                      MM_CURRENT_N_ST

#define RELATION_LIE_SUM(i)           (Relation[i].lie_sum)
#define RELATION_MIN_GENERATOR(i)     (Relation[i].min_generator)
#define RELATION_TO_BE_SUBSTITUTED(i) (Relation[i].to_be_substituted)

#define SCALAR_SUM_IS_UNIT(a)  (SCALAR_TERM_MONOMIAL(a)==NIL&&\
                                INTEGER_IS_UNIT(SCALAR_TERM_NUMERATOR(a)))

#define SCALAR_SUM_IS_UNIT_ABS(a)  (SCALAR_TERM_MONOMIAL(a)==NIL&&\
                              INTEGER_IS_UNIT_ABS(SCALAR_TERM_NUMERATOR(a)))

#define SCALAR_SUM_IS_NOT_UNIT(a)  (SCALAR_TERM_MONOMIAL(a)!=NIL||\
                                    INTEGER_IS_NOT_UNIT(\
                                          SCALAR_TERM_NUMERATOR(a)))

#define SCALAR_TERM_MONOMIAL(a)  (NodeST[a].monomial)
#define SCALAR_TERM_NUMERATOR(a) (NodeST[a].numerator)
#define SCALAR_TERM_R(a)         (NodeST[a].rptr)

#define SCALAR_FACTOR_PARAMETER(a)    (NodeSF[a].parameter)
#define SCALAR_FACTOR_IS_I_NUMBER(a)  (NodeSF[a].parameter==I_NUMBER)
#define SCALAR_FACTOR_DEGREE(a)       (NodeSF[a].degree)
#define SCALAR_FACTOR_WORD(a)         (*(unsigned short *)(NodeSF+(a)))
#define SCALAR_FACTOR_R(a)            (NodeSF[a].rptr)

#define SCALAR_TERM_MINUS(a)  INTEGER_MINUS(NodeST[a].numerator)
#define SCALAR_TERM_MAIN_PARAMETER(a) \
                           SCALAR_FACTOR_PARAMETER(SCALAR_TERM_MONOMIAL(a))
#define SCALAR_TERM_MAIN_PARAMETER_WORD(a) \
                                SCALAR_FACTOR_WORD(SCALAR_TERM_MONOMIAL(a))
#define SCALAR_TERM_MAIN_DEGREE(a) \
                              SCALAR_FACTOR_DEGREE(SCALAR_TERM_MONOMIAL(a))
#define POLY_MAIN_PARAMETER(a) ((SCALAR_TERM_MONOMIAL(a)==NIL) ? -1 :\
                          SCALAR_FACTOR_PARAMETER(SCALAR_TERM_MONOMIAL(a)))
#define POLY_ARRAY_STACK_NEW(a,n)   (a)=(uint*)alloca(sizeof(uint)*(n))\
                                    POLY_ARRAY_IN_STACK(a)

#define TIME_OFF  TimeC += (CrudeTime ? time(NULL) : clock()) - TimeA
#define TIME_ON   TimeA = (CrudeTime ? time(NULL) : clock())

/*_5    Global variables and arrays=====================================*/

/* Files */

#if !defined(GAP)
FILE *MessageFile;
FILE *SessionFile;
#endif

/* Single variables */

int IncompletedBasis;
int IncompletedRelations;
int IsParametric;
int LieMonomialIsNew;
int SubstitutionIsDone;
uint LimitingWeight;
int GeneratorN;
int ParameterN;

/* Arrays */

LIE_MON *LieMonomial;    /* Set of Lie monomials */
int LieMonomialSize;
int LieMonomialN;
int LieMonomialFreePosition; /* Start search of free position */
#if defined(SPACE_STATISTICS)
int LieMonomialMaxN;
#endif


NODE_LT *NodeLT;         /* Pool of nodes for Lie terms */
uint NodeLTSize;
uint NodeLTTop = 1;
#if defined(SPACE_STATISTICS)
uint NodeLTTopMax;
#endif

NODE_SF *NodeSF;         /* Pool of nodes for scalar factors */
uint NodeSFSize;
uint NodeSFTop = 1;
#if defined(SPACE_STATISTICS)
uint NodeSFTopMax;
#endif

NODE_ST *NodeST;         /* Pool of nodes for scalar terms */
uint NodeSTSize;
uint NodeSTTop = 1;
#if defined(SPACE_STATISTICS)
uint NodeSTTopMax;
#endif

REL *Relation;           /* Ordered set of relations */
int RelationSize;
int RelationN;
#if defined(SPACE_STATISTICS)
int MaxNRelation;
#endif
#if defined(INTEGER_MAX_SIZE)
LIMB IntegerMaxSize;
#endif

int CoeffSumTableSize; /* Non-zero coefficient tables variables */
int CoeffSumTableN;
uint *CoeffSumTable;     /* Non-zero parametric sums */
int *CoeffParamTable;   /* Table for memorizing single */
                         /* non-zero parameters */

/* Input and output variables */

char BasisSymbolEven;
char BasisSymbolOdd;
int CurrentLevel;
char * GeneratorName;     /* Input names */
char * ParameterName;
int GeneratorMaxN;     /* Maximum number of input generators */
uint InputIntegerSize;  /* Maximum size of input integer in LIMBs */
int InputStringSize;   /* Size of string for reading input */
int LastItemEnd;
int LineLength;
int Margin;
int MaxLevel;
int MinLevel;
int NameLength1;       /* Maximum length of input name (with ending '\0') */
int NewMargin;
int ParameterMaxN = 1; /* Maximum number of input parameters (i at least) */
int BasisElementsPut;
int CommutatorsPut;
int CrudeTime;         /* Prevent time variable wrapping for large tasks */
int EchoInput;
int GAPOutputCommutators;
int GAPOutputBasis;
int GAPOutputRelations;
char GAPAlgebraName[GAP_NAME_SIZE];
char GAPBasisName[GAP_NAME_SIZE];
char GAPRelationsName[GAP_NAME_SIZE];
int HeadingPut;
int HilbertSeriesPut;
int InitialRelationsPut;
int NonZeroCoefficientsPut;
int ReducedRelationsPut;
int StatisticsPut;
char * OutLine;           /* String for preparation of output block */
int OutLineSize;
int PosOutLine;
int PreviousEnd;
int PrintEnd;

uint TimeA, TimeC;

#if defined(DEBUG)
uint Debug;
#endif
#if defined(MEMORY)
int CurrentNLT;
int CurrentNSF;
int CurrentNST;
int CurrentNINT;

#endif

/*_6    Function descriptions===========================================*/

/*_6_0          Main and top level functions============================*/

void ConstructFreeAlgebraBasis(void);       /* 1 call! */
int FindNewPositionInRelation(int lmo);
void GenerateRelations(void);               /* 1 call! */
uint NewJacobiRelation(int l);           /* 1 call! */
int ReduceRelations(int i);

/*_6_1          Pairing functions=======================================*/

int AddPairToLieMonomial(int i, int j);
uint MakeRelationRHSInteger(int i);             /* 1 call!! */
uint MakeRelationRHSParametric(int i);          /* 1 call!! */
uint PairMonomialMonomialInteger(int i, int j);
uint PairMonomialMonomialParametric(int i, int j);
uint PairMonomialSumInteger(int mon, uint a);
uint PairMonomialSumParametric(int mon, uint a);
uint PairSumMonomialInteger(uint a, int mon);
uint PairSumMonomialParametric(uint a, int mon);
uint PairSumSumInteger(uint a, uint b);
uint PairSumSumParametric(uint a, uint b);

/*_6_2          Substitution (replacing) functions======================*/

int IsMonomialInMonomial(int submon, int mon);
uint SubstituteRelationInRelationInteger(uint r, uint a);
uint SubstituteRelationInRelationParametric(uint r, uint a);
uint SubstituteRHSInMonomialInteger(int mon, int lmonr, uint r);
uint SubstituteRHSInMonomialParametric(int mon, int lmonr, uint r);

/*_6_3          Lie and scalar algebra functions========================*/

int LieLikeTermsCollectionInteger(uint a, uint b);
int LieLikeTermsCollectionParametric(uint a, uint b);
uint LieSumAddition(uint a, uint b);
void LieSumDivInteger(uint lsum, BIGINT den);
void LieSumDivScalarSum(uint lsum, uint den);
void LieSumMinusInteger(uint a);
void LieSumMinusParametric(uint a);
void LieSumMultInteger(uint lsum, BIGINT num);
#if defined(RATIONAL_FIELD)
void LieSumMultRationalInteger(int a, BIGINT num, BIGINT den);
#endif
void LieSumMultScalarSum(uint lsum, uint num);
void NormalizeRelationInteger(uint a);
void NormalizeRelationParametric(uint a);
uint ScalarMonomialMultiplication(int *pchange_sign, uint ma, uint mb);
uint ScalarSumAddition(uint a, uint b);
void ScalarSumCancellation(uint *pnum, uint *pden);
void ScalarSumMinus(uint a);
uint ScalarSumMultiplication(uint a, uint b);
void ScalarTermMultiplication(uint a, uint b);     /* 1 call! */

/*_6_4          Scalar polynomial algebraic functions===================*/

uint ContentOfScalarSum(uint cont, uint a);
void InCoeffParamTable(uint cont);
void InCoeffSumTable(uint sum);
void InCoeffTable(uint coe);
uint PolyCoeffAtMainParameter(uint *pa, int mp);
uint PolyContent(uint a, int mp);
uint PolyGCD(uint a, uint b);
uint PolyMainParameterTerm(uint *pa, int mp, int mpdeg);
int PolynomialsAreEqual(uint a, uint b);
uint PolyPseudoRemainder(uint a, uint b, int mp);  /* 1 call! */
uint PolyTermGCD(uint a, uint b);
void PolyTermQuotient(uint a, uint b);
uint PolyQuotient(uint a, uint b);

/*_6_5          Big number functions====================================*/

int BigNMinusBigN(BIGINT a, int na, BIGINT b, int nb);
LIMB BigNShiftLeft(BIGINT bign, int n, int cnt);
int BigNShiftRight(BIGINT bign, int n, int cnt);
int CountLeadingZeroBitsInLimb(LIMB w);
void IntegerCancellation(BIGINT num, BIGINT den);
BIGINT IntegerGCD(BIGINT u, BIGINT v);
void IntegerProduct(BIGINT w, BIGINT u, BIGINT v);
void IntegerQuotient(BIGINT c, BIGINT a, BIGINT b);
void IntegerSum(BIGINT c, BIGINT a, BIGINT b);

/*_6_6          Copy and delete functions===============================*/

uint LieSumCopyInteger(uint a);
uint LieSumCopyIntegerNegative(uint a);
uint LieSumCopyParametric(uint a);
void LieSumKillInteger(uint a);
void LieSumKillParametric(uint a);
uint LieTermFromMonomialInteger(int mon);
uint LieTermFromMonomialParametric(int mon);
uint ScalarSumCopy(uint a);
void ScalarSumKill(uint a);
uint ScalarTermCopy(uint a);

/*_6_7          Technical functions=====================================*/

void Error(int i_message) ATTRIBUTE_NORETURN;
void Initialization(void);
void *NewArray(uint n, uint size, int i_message);
uint NodeLTNew(void);
uint NodeSFNew(void);
uint NodeSTNew(void);
FILE *OpenFile(char * file_name, char * file_type);

/*_6_8          Input functions=========================================*/

int BinaryQuestion(int i_message);
int FindNameInTable(char * name, char * nametab, int n_nametab);
void GetGenerator(char * str);
void GetInput(int n, char * fin);
void GetInteger(BIGINT a, char **pstr);
uint GetLieMonomial(char **pstr);
uint GetLieSum(char **pstr);
uint GetLieTerm(char **pstr);
uint GetUInteger(char **pstr);
void GetParameter(char * str);
void GetRelation(char * str);
uint GetScalarSum(char **pstr);
uint GetScalarTerm(char **pstr);
void GetWeight(char * str);
int KeyBoardBytesToString(char * str);
int KeyBoardStringToFile(int i_m, char * prefix, char * str, FILE *file);
void ReadAndProcessStringsFromFile(void (*proc_func)(char * str), FILE *inf,
                                                        char sep, char end);
int ReadBooleanFromFile(FILE *file);
int ReadCaseFromFile(FILE * file, char * case_str[], int n_cases);
uint ReadDecimalFromFile(FILE *file);
short ReadStringFromFile(char * str, FILE *file);
short SkipCommentInFile(FILE *file);
void SkipName(char **pstr);
void SkipSpaces(char **pstr);
short SkipSpacesInFile(FILE *file);

/*_6_9          Output functions========================================*/

void AddSymbolToOutLine(char c, int position);
void InLineLevel(int level);
void InLineNumberInBrackets(uint n);
void InLineString(char * str);
void InLineSubscript(char * s);
void InLineSymbol(char symbol);
void InLineTableName(char * name);
char * UToString(uint n);
void PutBasis(void);
#if defined(GAP)
void PutBasisGAP(void);
#endif
void PutBlock(void);
void PutCharacter(char c);
#if defined(GAP)
void PutCharacterGAP(char c);
#endif
void PutCoefficientTable(void);
void PutCommutators(void);
#if defined(GAP)
void PutCommutatorsGAP(void);
#endif
void PutDegree(uint deg);
void PutDimensions(void);
void PutDots(void);
void PutEnd(void);
void PutFormattedU(char * format, uint i);
void PutIntegerUnsigned(BIGINT bn);
#if defined(GAP)
void PutIntegerUnsignedGAP(BIGINT bn);
#endif
void PutLieBareTerm(void (*put_lie_mon)(int a), uint a);
void PutLieBasisElement(int pos);
void PutLieMonomialLeftNormed(int pos);
void PutLieMonomialStandard(int pos);
#if defined(GAP)
void PutLieMonomialGAP(int pos);
#endif
void PutLieSum(void (*put_lie_mon)(int a), uint a);
void PutMessage(int i_message);
void PutRelations(int i);
#if defined(GAP)
void PutRelationsGAP(void);
#endif
void PutScalarBareTerm(uint a);
void PutScalarFactor(uint a);
void PutScalarSum(uint a);
void PutStart(void);
void PutStatistics(void);
void PutString(char * str);
#if defined(GAP)
void PutStringGAP(char * str);
#endif
void PutStringStandard(char * str);
void PutSymbol(char c);

/* Global function variables */

int (*LieLikeTermsCollection)(uint a, uint b) = LieLikeTermsCollectionInteger;
uint (*LieSumCopy)(uint a) = LieSumCopyInteger;
void (*LieSumKill)(uint a) = LieSumKillInteger;
void (*LieSumMinus)(uint a) = LieSumMinusInteger;
uint (*LieTermFromMonomial)(int mon) = LieTermFromMonomialInteger;
void (*NormalizeRelation)(uint a) = NormalizeRelationInteger;
uint (*PairMonomialMonomial)(int i, int j) = PairMonomialMonomialInteger;
uint (*PairMonomialSum)(int mon, uint a) = PairMonomialSumInteger;
uint (*PairSumMonomial)(uint a, int mon) = PairSumMonomialInteger;
uint (*PairSumSum)(uint a, uint b) = PairSumSumInteger;
void (*PutLieMonomial)(int pos) = PutLieMonomialStandard;
uint (*SubstituteRelationInRelation)(uint r, uint a) =
SubstituteRelationInRelationInteger;

/*_6_10   Debugging functions===========================================*/

#if defined(DEBUG)
void PutDebugHeader(uint debug, char * f_name, char * in_out);
void PutDebugInteger(char * name, BIGINT u);
void PutDebugLieMonomial(char * name, int a);
void PutDebugLieMonomialTable(int newmon);
void PutDebugLieSum(char * name, uint a);
void PutDebugLieTerm(char * name, uint a);
void PutDebugU(char * name, uint i);
#if defined(D_PUT_RELATIONS)
void PutDebugRelations(void);
#endif
void PutDebugScalarSum(char * name, uint a);
void PutDebugString(char * strname, char * str);
#endif
#if defined(MEMORY)
void AddLieSumNs(uint a, int minus_or_plus,
              int *pn_lt, int *pn_int, int *pn_st, int *pn_sf);
void AddScalarSumNs(uint a, int minus_or_plus, int *pn_int, int *pn_st, int *pn_sf);
void PutIntegerBalance(char * fname, int dn);
void PutNodeBalance(char * type, char * fname, int dn);
#endif

/*_6_0          Main and top level functions============================*/

#if !defined(TEST_FUNCTION)
/*=main======================================
*/
int main(int narg, char ** fin)
{
  Initialization();
  GetInput(narg, fin[1]);
  if(RelationN)
  {
    if(InitialRelationsPut)
      PutRelations(H_IN_RELATIONS);
    GenerateRelations();
    if(NonZeroCoefficientsPut)
      PutCoefficientTable();
  }
  if(RelationN)
  {
    if(ReducedRelationsPut)
      PutRelations(H_REDUCED_RELATIONS);
  }
  else /* Free algebra */
    ConstructFreeAlgebraBasis();
  PutBasis();
  if(HilbertSeriesPut)
    PutDimensions();
  if(CommutatorsPut)
    PutCommutators();
#if defined(DEBUG)
  PutDebugU("Top Debug", Debug);
#endif
  if(StatisticsPut)
  {
    TIME_OFF;
    PutStatistics();
  }
#if defined(GAP)
  if(!IsParametric && !IncompletedRelations && !IncompletedBasis)
  {
/*    if(GAPOutputCommutators || GAPOutputBasis || GAPOutputRelations) */
/*    {                                                                */
/*      fclose(SessionFile);                                           */
/*#if defined(SPP_2000)                                                */
/*      SessionFile = OpenFile("fplsa4.gap", "w");                     */
/*#else                                                                */
/*      SessionFile = OpenFile("fplsa4.gap", "wt");                    */
/*#endif                                                               */
/*    }                                                                */
    if(GAPOutputCommutators)
      PutCommutatorsGAP();
    if(GAPOutputBasis)
      PutBasisGAP();
    if(GAPOutputRelations)
      PutRelationsGAP();
  }
#endif
  exit(0);
}
#endif
/*=ConstructFreeAlgebraBasis==========================================
 Make all regular monomials up to LimitingWeight
*/
void ConstructFreeAlgebraBasis(void)
{
  int i = 0, j, moni, monj;
  while(YES)
  {
    moni = LIE_MONOMIAL_POSITION(i);
    if(LIE_MONOMIAL_IS_NOT_SQUARE(moni)) /* ?? Not so for non-standard */
    {
      j = 0;
      while(j < i)
      {
        monj = LIE_MONOMIAL_POSITION(j);
        if(LIE_MONOMIAL_IS_NOT_SQUARE(monj))
          if(LIE_MONOMIAL_IS_GENERATOR(moni) ||
             j >= LIE_MONOMIAL_RIGHT_ORDER(moni))
          {
            if(LIE_MONOMIAL_WEIGHT(moni) + LIE_MONOMIAL_WEIGHT(0) >
               LimitingWeight) /* Out of weight for all consequent i & j */
              goto incompleted;
            if(LIE_MONOMIAL_WEIGHT(moni) + LIE_MONOMIAL_WEIGHT(monj) >
               LimitingWeight) /* Out of weight for all consequent j */
              goto next_i;
            AddPairToLieMonomial(moni, monj);
          }
        j++;
      }
      if(LIE_MONOMIAL_IS_ODD(moni))
      {                           /* Add square */
        if(LIE_MONOMIAL_WEIGHT(moni) + LIE_MONOMIAL_WEIGHT(0) >
           LimitingWeight) /* Out of weight for all consequent i & j */
          goto incompleted;
        if(2*LIE_MONOMIAL_WEIGHT(moni) > LimitingWeight)
          goto next_i;         /* Out of weight for all consequent j */
        AddPairToLieMonomial(moni, moni);
      }
    }
    next_i:
    if(++i >= LieMonomialN)
      return;            /* One generator case */
  }
  incompleted:
  IncompletedBasis = YES;
}
/*=FindNewPositionInRelation======================================
 Find position of first relation with leading monomial order > lmo
 among 0, 1,..., RelationN - 1
*/
int FindNewPositionInRelation(int lmo)
{
  int left = 0;
  if(RelationN)                                /* Binary search */
  {                           /* `right' must be of signed type */
    int m, right = RelationN-1;
    do
    {
      m = (left + right)/2;
      if(lmo < LIE_TERM_MONOMIAL_ORDER(RELATION_LIE_SUM(m)))
        right = --m;
      else
        left = ++m;
    }while(left <= right);
  }
  return left;
}
/*=GenerateRelations=======================================================
 Generate and process new relations
*/
void GenerateRelations(void)
{
  int i, k, l, mon, mona,
      gen;    /* LieMonomial table position of differentiating generator */
  uint a, lim_weight_i;
  i = 0;
  while(i < RelationN)                           /* Differentiation loop */
    if(RELATION_MIN_GENERATOR(i) < GeneratorN)
      if(LIE_MONOMIAL_IS_BASIS(gen = RELATION_MIN_GENERATOR(i)))
      {
        /* Program assures the ban of differentiation OF leading generators
           in the process of new relation adding by setting negation of the
           first IF predicate */
        a = RELATION_LIE_SUM(i);
        mona = LIE_TERM_MONOMIAL(a);   /* Irregular triple criterion */
        if(LIE_MONOMIAL_IS_SQUARE(mona) ||
           LIE_MONOMIAL_ORDER(gen) < LIE_MONOMIAL_RIGHT_ORDER(mona))
        {
IN_GENERATE_RELATIONS  /*------------------------------------------------*/
          if(LIE_MONOMIAL_WEIGHT(gen) +                 /* Out of weight */
             LIE_MONOMIAL_WEIGHT(LIE_TERM_MONOMIAL(a)) > LimitingWeight)
          {
            RELATION_MIN_GENERATOR(i++) = GeneratorN;
            continue;     /* There might be lower weight next generators */
          }
          RELATION_MIN_GENERATOR(i)++;       /* For next differentiation */
          if((a = (*PairSumMonomial)((*LieSumCopy)(a), gen)) != NIL)
          {
            add_new_relation:
            if(RelationN == RelationSize)
              Error(E_RELATION_SIZE);
            gen = LIE_TERM_MONOMIAL(a);
            l = FindNewPositionInRelation(k = LIE_MONOMIAL_ORDER(gen));
#if defined(SPACE_STATISTICS)
            if(RelationN >= MaxNRelation)
              MaxNRelation = RelationN + 1;
#endif
            (*NormalizeRelation)(a);
OUT_GENERATE_RELATIONS /*------------------------------------------------*/
            LIE_MONOMIAL_INDEX(gen) = ~l;    /* Set position of relation */

            /* Shift positions of higher  relations in LieMonomial */

            while(++k < LieMonomialN)
              if(LIE_MONOMIAL_IS_LEADING_BY_ORDER(k))
                --LIE_MONOMIAL_INDEX_BY_ORDER(k);

            /* Make room for new relation */

            for(k = RelationN; k > l; k--)
              Relation[k] = Relation[k-1];

            if(LIE_MONOMIAL_IS_GENERATOR(gen))    /* Ban differentiating */
              RELATION_MIN_GENERATOR(l) = GeneratorN;  /* lead. generat. */
            else
            {
              RELATION_MIN_GENERATOR(l) = 0;
              if(l <= i)                       /* Shift min. diff. index */
                i = l;
            }
            RELATION_LIE_SUM(l) = a;                 /* Set new relation */
            if(LIE_MONOMIAL_WEIGHT(gen) + 1 > LimitingWeight)
              RELATION_MIN_GENERATOR(l) = GeneratorN;
            if(l == RelationN++)
              RELATION_TO_BE_SUBSTITUTED(l) = NO;
            else
            {                           /* New relation inside the table */
              RELATION_TO_BE_SUBSTITUTED(l) = YES;
              l = ReduceRelations(l);
              if(l <= i)
                i = l;
            }
     /* #if defined(RELATION_N_TO_SCREEN)
            TIME_OFF;
            printf("\n%10d", RelationN);
            TIME_ON;
    #endif*/
OUT_PUT_LIE_MONOMIAL /*--------------------------------------------------*/
OUT_PUT_RELATIONS /*-----------------------------------------------------*/
          }
        }
        else                   /* Any next generator >= right of leading */
          RELATION_MIN_GENERATOR(i++) = GeneratorN;
      }
      else                  /* Skip differentiation BY leading generator */
        RELATION_MIN_GENERATOR(i)++;
    else
      i++;                    /* Skip completely differentiated relation */
  IncompletedBasis =                       /* Limiting weight is reached */
  IncompletedRelations =
    (LIE_MONOMIAL_WEIGHT(LIE_TERM_MONOMIAL(RELATION_LIE_SUM(RelationN-1)))
     >= LimitingWeight);

/*??#if 0 vvvvvvvvvvvvvvvvv Off checking vvvvvvvvvvvvvvvvvvvvvvvvvvvvv*/
  /* Check regular pairs */

  k = 0;
  while(k < LieMonomialN)
  {
    gen = LIE_MONOMIAL_POSITION(k);
    if(LIE_MONOMIAL_IS_BASIS(gen))
    {
#if 0 /*?? Experiment with individual basis elements vvvvvvvvvvvvvvvv */
      if(LIE_MONOMIAL_IS_COMMUTATOR(gen))
      {
                             /* Old pairs */
        a = NewJacobiRelation(gen);
        if(a != NIL)
/*??*/{
/*??*/PutDebugLieSum("\n***New Jacobi Relation from Old Basis", a);
          goto add_new_relation;
/*??*/}
      }
#endif /*??  Experiment ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */
      if(LIE_MONOMIAL_WEIGHT(gen) + LIE_MONOMIAL_WEIGHT(0)
         > LimitingWeight)
      {
        IncompletedBasis = YES;
        goto loops_out;
      }
      if(LIE_MONOMIAL_IS_NOT_SQUARE(gen)) /* ?? non-standard square */
      {
        lim_weight_i = LimitingWeight - LIE_MONOMIAL_WEIGHT(gen);
        i = 0;
        do
        {
          mon = LIE_MONOMIAL_POSITION(i);
          if(LIE_MONOMIAL_IS_BASIS(mon) &&
             LIE_MONOMIAL_IS_NOT_SQUARE(mon) &&
             (i != k || LIE_MONOMIAL_IS_ODD(mon)) &&
             (LIE_MONOMIAL_IS_GENERATOR(gen) ||
              i >= LIE_MONOMIAL_RIGHT_ORDER(gen)))
          {
            if(LIE_MONOMIAL_WEIGHT(mon) > lim_weight_i)
              break; /* Stop considering next right monomials */
            mon = AddPairToLieMonomial(gen, mon);
            if(LieMonomialIsNew)
            {                              /* New pairs */
#if 0
/*??*/PutDebugLieMonomial("\n***New Lie Monomial", mon);
#endif
              a = NewJacobiRelation(mon);
              if(a != NIL)
#if 0
/*??*/{
/*??*/PutDebugLieSum("\n***New Jacobi Relation", a);
#endif
                goto add_new_relation;
#if 0
/*??*/}
#endif
            }
          }
        }while(++i <= k);
      }
    }
    ++k;
  }
  loops_out: ;
/*??#endif ^^^^^^^^^^^^^^^^^^^^ Off checking ^^^^^^^^^^^^^^^^^^^^*/
}
/*=NewJacobiRelation===================================================
  Construct independent Jacobi relation containing leading commutator L
*/
uint NewJacobiRelation(int l)
{
  uint a, b;
  int x, y, z;
IN_NEW_JACOBI_RELATION /*--------------------------------------------*/
  /* Try left pair in l = [[x,y],z] to make relation:  */
  /*                             p(x)p(y)              */
  /* [[x,y],z] - [x,[y,z]] + (-1)        [y,[x,z]] = 0 */

  x = LIE_MONOMIAL_LEFT(l);
  if(LIE_MONOMIAL_IS_COMMUTATOR(x))
  {
    z = LIE_MONOMIAL_RIGHT(l);
    y = LIE_MONOMIAL_RIGHT(x);
    x = LIE_MONOMIAL_LEFT(x);

    /* 1st of triple */

    a = (*LieTermFromMonomial)(l);

    /* 2nd of triple */

    if((b = (*PairMonomialMonomial)(z, y)) != NIL)
    {
      if(LIE_MONOMIAL_IS_ODD(y) && LIE_MONOMIAL_IS_ODD(z))
        b = (*PairSumMonomial)(b, x);
      else
        b = (*PairMonomialSum)(x, b);
      a = LieSumAddition(a, b);
    }

    /* 3d of triple */

    if(LIE_MONOMIAL_ORDER(x) > LIE_MONOMIAL_ORDER(z))
    {
      if((b = (*PairMonomialMonomial)(x, z)) != NIL)
      {
        if(LIE_MONOMIAL_IS_ODD(x) && LIE_MONOMIAL_IS_ODD(y))
        {
          b = (*PairSumMonomial)(b, y);
          if(LIE_MONOMIAL_IS_EVEN(z))
            (*LieSumMinus)(b);
        }
        else
          b = (*PairMonomialSum)(y, b);
        a = LieSumAddition(a, b);
      }
    }
    else
    {
      if((b = (*PairMonomialMonomial)(z, x)) != NIL)
      {
        if(LIE_MONOMIAL_IS_ODD(z) &&
           LIE_MONOMIAL_PARITY(x) != LIE_MONOMIAL_PARITY(y))
        {
          b = (*PairMonomialSum)(y, b);
          if(LIE_MONOMIAL_IS_EVEN(x))
            (*LieSumMinus)(b);
        }
        else
          b = (*PairSumMonomial)(b, y);
        a = LieSumAddition(a, b);
      }
    }

    if(a != NIL)
      goto out;
  }

  /* Try right pair in l = [x,[y,z]] to make relation: */
  /*                             p(x)p(y)              */
  /* [x,[y,z]] - [[x,y],z] - (-1)        [y,[x,z]] = 0 */

  y = LIE_MONOMIAL_RIGHT(l);
  if(LIE_MONOMIAL_IS_COMMUTATOR(y))
  {
    x = LIE_MONOMIAL_LEFT(l);
    z = LIE_MONOMIAL_RIGHT(y);
    y = LIE_MONOMIAL_LEFT(y);

    /* 1st of triple */

    a = (*LieTermFromMonomial)(l);

    /* 2nd of triple */

    if((b = (*PairMonomialMonomial)(x, y)) != NIL)
    {
      b = (*PairSumMonomial)(b, z);
      (*LieSumMinus)(b);
      a = LieSumAddition(a, b);
    }

    /* 3d of triple */

    if((b = (*PairMonomialMonomial)(x, z)) != NIL)
    {
      b = (*PairMonomialSum)(y, b);
      if(LIE_MONOMIAL_IS_EVEN(x) || LIE_MONOMIAL_IS_EVEN(y))
          (*LieSumMinus)(b);
      a = LieSumAddition(a, b);
    }

    goto out;
  }
  a = NIL;
  out:
OUT_NEW_JACOBI_RELATION /*-------------------------------------------*/
  return a;
}
/*=ReduceRelations=====================================================
  Reduce the system of relations starting from Ith one. For further
  differentiations returns lowest new position (or starting one int).
*/
int ReduceRelations(int i)
{
  int i_min, j, lordj, lordl, min_gen, lmoni, k, l, m;
  uint ai, aj;
  i_min = i;
  do           /* While relations with new leading monomials arise */
    if(RELATION_TO_BE_SUBSTITUTED(i))
    {
IN_REDUCE_RELATIONS  /*----------------------------------------------*/
      new_i:
      ai = RELATION_LIE_SUM(i);
      lmoni = LIE_TERM_MONOMIAL(ai);
      j = i + 1;
      while(j < RelationN)
      {
        aj = RELATION_LIE_SUM(j);
        lordj = LIE_TERM_MONOMIAL_ORDER(aj);
        aj = (*SubstituteRelationInRelation)(ai, aj);
        test_substitution_result:
        if(SubstitutionIsDone)
        {
          if(aj == NIL)                        /* Killed relation */
          {
            --RelationN;
            k = lordj;
            l = LIE_MONOMIAL_POSITION(k);
            LIE_MONOMIAL_INDEX(l) = 0;            /* To ban usage */
            while(++k < LieMonomialN)   /* Shift int's of relations */
            {
              l = LIE_MONOMIAL_POSITION(k);
              if(LIE_MONOMIAL_IS_LEADING(l))       /* LieMonomial */
                ++LIE_MONOMIAL_INDEX(l);
            }
            for(k = j; k < RelationN; k++) /* Remove gap shifting */
              Relation[k] = Relation[k+1]; /* down top relations  */
            continue;
          }
          else          /* Non-killing substitution has been done */
          {
            (*NormalizeRelation)(aj);
            lordl = LIE_TERM_MONOMIAL_ORDER(aj);
            if(lordl < lordj)        /* Change of leading ordinal */
            {
              l = 0;
              k = j - 1;         /* Binary search of new position */
              while(l <= k)
              {
                m = (l + k)/2;
                if(lordl < LIE_TERM_MONOMIAL_ORDER(RELATION_LIE_SUM(m)))
                  k = --m;
                else
                  l = ++m;
              }
              if(l &&
                 LIE_TERM_MONOMIAL_ORDER(RELATION_LIE_SUM(l-1)) == lordl)
              {             /* Substitute once more to avoid collision */
                aj = (*SubstituteRelationInRelation)
                                           (RELATION_LIE_SUM(l-1), aj);
                goto test_substitution_result;
              }
              /* Set index of dropped relation in LieMonomial */
              LIE_MONOMIAL_INDEX_BY_ORDER(lordl) = ~l;
              /* Set zero to ban using old leading monomial
                 of dropped relation */
              LIE_MONOMIAL_INDEX_BY_ORDER(lordj) = 0;
              min_gen = RELATION_MIN_GENERATOR(j);
              /* Shift upper relations and their indices */
              k = (j+1 == RelationN) ? LieMonomialN :
                  LIE_TERM_MONOMIAL_ORDER(RELATION_LIE_SUM(j+1));
              while(--k > lordl)
                if(LIE_MONOMIAL_IS_LEADING_BY_ORDER(k))
                  --LIE_MONOMIAL_INDEX_BY_ORDER(k); /* Means ++ */
              for(k = j; k > l; k--)
                Relation[k] = Relation[k-1];
              RELATION_MIN_GENERATOR(l) =
                LIE_MONOMIAL_IS_GENERATOR_BY_ORDER(lordl) ?
                GeneratorN : /* Don't differentiate leading generator */
                min_gen;   /* Avoiding redifferentiation (or set 0 ??) */
              RELATION_TO_BE_SUBSTITUTED(l) = (byte)(l < RelationN-1);
              RELATION_LIE_SUM(l) = aj;
              if(l <= i)
              {         /* New substituted relation with lesser ordinal */
                i = l;
                if(i < i_min)
                  i_min = i;
                goto new_i;
              }
            }
            else                  /* No change of leading ordinal */
              RELATION_LIE_SUM(j) = aj;
          }
        }
        j++;
      }
      /* Remove consequences of leading monomial of substituted relation
         from LieMonomial table */
      lordj = LieMonomialN-1;
      m = LIE_MONOMIAL_WEIGHT(lmoni);
      while(LIE_MONOMIAL_WEIGHT_BY_ORDER(lordj) > (uint)m)
      {
        l = LIE_MONOMIAL_POSITION(lordj);
        if(IsMonomialInMonomial(lmoni, l))
        {
          if(l < LieMonomialFreePosition)   /* Set possibly lowest */
            LieMonomialFreePosition = l;    /* free position */
          LIE_MONOMIAL_IS_OCCUPIED(l) = NO; /* Mark free position */
          LieMonomialN--;
          lordl = k = lordj;
          while(k < LieMonomialN)
          {
            l = LIE_MONOMIAL_POSITION(++lordl); /* Shift positions */
            LIE_MONOMIAL_POSITION(k++) = l;     /* of upper orders */
            --LIE_MONOMIAL_ORDER(l); /* Decrease orders of upper monomials */
          }
        }
        lordj--;
      }
      RELATION_TO_BE_SUBSTITUTED(i) = NO;
OUT_REDUCE_RELATIONS /*----------------------------------------------*/
    }
  while(++i < RelationN);
  return i_min;
}

/*_6_1          Pairing functions=======================================*/

/*=AddPairToLieMonomial=================================================
 Find position in LieMonomial table for regular pair [i,j].
 LieMonomialIsNew is global boolean signal:
 LieMonomialIsNew == YES if pair has to be added in table,
 if so, add the pair to the table,
 LieMonomialIsNew == NO if pair exists already in table.
 Return position for [i,j]
*/
int AddPairToLieMonomial(int i, int j)
{
  uint wt = LIE_MONOMIAL_WEIGHT(i) + LIE_MONOMIAL_WEIGHT(j);
  int ijo /* left */, r /* right */, m /* middle */, ijp;
IN_ADD_PAIR_TO_LIE_MONOMIAL  /*------------------------------------------*/
  ijo = 0;
  r = LieMonomialN - 1;
  do
  {
    m = (ijo + r)/2;
    ijp = LIE_MONOMIAL_POSITION(m);

    /* Compare wrt weights */

    if(wt > LIE_MONOMIAL_WEIGHT(ijp))
      goto shift_left;
    if(wt < LIE_MONOMIAL_WEIGHT(ijp))
      goto shift_right;

    /* Equal weights: compare lexicographically */

    if(LIE_MONOMIAL_ORDER(i) > LIE_MONOMIAL_LEFT_ORDER(ijp))
      goto shift_left;
    if(LIE_MONOMIAL_ORDER(i) < LIE_MONOMIAL_LEFT_ORDER(ijp))
      goto shift_right;
    if(LIE_MONOMIAL_ORDER(j) > LIE_MONOMIAL_RIGHT_ORDER(ijp))
      goto shift_left;
    if(LIE_MONOMIAL_ORDER(j) < LIE_MONOMIAL_RIGHT_ORDER(ijp))
      goto shift_right;
    LieMonomialIsNew = NO;
OUT_ADD_PAIR_TO_LIE_MONOMIAL_OLD  /*-------------------------------------*/
    return ijp;
    shift_left:
    ijo = ++m;
    continue;
    shift_right:
    r = --m;
  }while(ijo <= r);

  /* Add new monomial to table */

  LieMonomialIsNew = YES;
  if(LieMonomialN >= LieMonomialSize)
  {
    TIME_OFF;
    PutStatistics();                          /* No room for new element */
    Error(E_LIE_MONOMIAL_SIZE);
  }
  m = r = LieMonomialN++;
#if defined(SPACE_STATISTICS)
  if(LieMonomialN > LieMonomialMaxN)
    LieMonomialMaxN = LieMonomialN;
#endif
  while(m > ijo)
  {
    ijp = LIE_MONOMIAL_POSITION(--r); /* Shift positions of upper orders */
    LIE_MONOMIAL_POSITION(m--) = ijp;
    ++LIE_MONOMIAL_ORDER(ijp);      /* Increase orders of upper elements */
  }
  while(LIE_MONOMIAL_IS_OCCUPIED(LieMonomialFreePosition))
    LieMonomialFreePosition++;                   /* Search free position */
  ijp = LieMonomialFreePosition++;
  LIE_MONOMIAL_ORDER(ijp) = ijo;
  LIE_MONOMIAL_POSITION(ijo) = ijp;  /* Parity is 1bit field: + is mod 2 */
  LIE_MONOMIAL_PARITY(ijp) = LIE_MONOMIAL_PARITY(i) + LIE_MONOMIAL_PARITY(j);
  LIE_MONOMIAL_LEFT(ijp) = i;
  LIE_MONOMIAL_RIGHT(ijp) = j;
  LIE_MONOMIAL_WEIGHT(ijp) = wt;
  LIE_MONOMIAL_INDEX(ijp) = 0;              /* Set type of basis element */
OUT_ADD_PAIR_TO_LIE_MONOMIAL_NEW  /*-------------------------------------*/
  return ijp;
}
/*=MakeRelationRHSInteger============================================
 Make r.h.s. of i-th relation (Integer regime)
 Relation is in normalized form (no denominators, no content)
*/
uint MakeRelationRHSInteger(int i)
{
  uint c;
  if((c = LIE_TERM_R(RELATION_LIE_SUM(i))) != NIL)
  {
    BIGINT n, cn;
#if !defined(RATIONAL_FIELD)
    BIGINT ln = LIE_TERM_NUMERATOR_INTEGER(RELATION_LIE_SUM(i));
#endif
    uint a, ea;
IN_MAKE_RELATION_RHS  /*-------------------------------------------*/
    /* Copy r.h.s setting negations for numerators */
    a = ea = NodeLTNew();
    LIE_TERM_MONOMIAL(ea) = LIE_TERM_MONOMIAL(c);
    n = LIE_TERM_NUMERATOR_INTEGER(c);
    INTEGER_HEAP_COPY(cn, n, i);
    INTEGER_MINUS(cn);
    LIE_TERM_NUMERATOR_INTEGER(ea) = cn;
    if((n = LIE_TERM_DENOMINATOR_INTEGER(c)) != NULL)
    {
      INTEGER_HEAP_COPY(cn, n, i);
      LIE_TERM_DENOMINATOR_INTEGER(ea) = cn;
    }
    else
      LIE_TERM_DENOMINATOR_INTEGER(ea) = NULL;
    while((c = LIE_TERM_R(c)) != NIL)
    {
      LIE_TERM_R(ea) = NodeLTNew();
      ea = LIE_TERM_R(ea);
      LIE_TERM_MONOMIAL(ea) = LIE_TERM_MONOMIAL(c);
      n = LIE_TERM_NUMERATOR_INTEGER(c);
      INTEGER_HEAP_COPY(cn, n, i);
      INTEGER_MINUS(cn);
      LIE_TERM_NUMERATOR_INTEGER(ea)= cn;
      if((n = LIE_TERM_DENOMINATOR_INTEGER(c)) != NULL)
      {
        INTEGER_HEAP_COPY(cn, n, i);
        LIE_TERM_DENOMINATOR_INTEGER(ea) = cn;
      }
      else
        LIE_TERM_DENOMINATOR_INTEGER(ea) = NULL;
    }
#if !defined(RATIONAL_FIELD)
    /* Divide by leading coefficient */
    if(INTEGER_IS_NOT_UNIT(ln))
    {
      INTEGER_STACK_COPY(cn, ln, i);
      LieSumDivInteger(a, cn);
    }
#endif
OUT_MAKE_RELATION_RHS /*-------------------------------------------*/
    return a;
  }
  return NIL;
}
/*=MakeRelationRHSParametric============================================
 Make r.h.s. of i-th relation (Parametric regime)
 Relation is in normalized form (no denominators, no content)
*/
uint MakeRelationRHSParametric(int i)
{
  uint a;
IN_MAKE_RELATION_RHS  /*----------------------------------------------*/
  if((a = LieSumCopyParametric(LIE_TERM_R(RELATION_LIE_SUM(i)))) != NIL)
  {
    uint lc;

    /* Negation */

    LieSumMinusParametric(a);

    /* Divide by leading coefficient */

    lc = LIE_TERM_NUMERATOR_SCALAR_SUM(RELATION_LIE_SUM(i));
    if(SCALAR_SUM_IS_NOT_UNIT(lc))
    {
      lc = ScalarSumCopy(lc);
      LieSumDivScalarSum(a, lc);
    }
  }
OUT_MAKE_RELATION_RHS /*----------------------------------------------*/
  return a;
}
/*=PairMonomialMonomialInteger==============================
 Make regular expression from two monomials (Integer regime)
 Caller ensures ORDER(i) >= ORDER(j)
*/
uint PairMonomialMonomialInteger(int i, int j)
{
  uint a;
  int k;
IN_PAIR_MONOMIAL_MONOMIAL  /*-----------------------------*/
  if(LIE_MONOMIAL_IS_SQUARE(j))
  {                          /* [i,[j,j]] = 2 [[i,j],j] */
    j = LIE_MONOMIAL_LEFT(j);
    a = PairMonomialMonomialInteger(i, j);
    a = PairSumMonomialInteger(a, j);
    if(a != NIL)
    {
      LIMB two[2] = {1, 2};
      LieSumMultInteger(a, two);
    }
  }
  else if(LIE_MONOMIAL_IS_SQUARE(i))
  {
    LIMB two[2] = {1, 2};
    i = LIE_MONOMIAL_LEFT(i);
    if(i == j)
      a = NIL;                         /* [[i,i],i] = 0 */
    else
    {
      if(LIE_MONOMIAL_ORDER(i) < LIE_MONOMIAL_ORDER(j))
        a = PairMonomialMonomialInteger(j, i);
                           /* [[i,i],j] = - 2 [[j,i],i] */
      else
      {
        a = PairMonomialMonomialInteger(i, j);
        if(LIE_MONOMIAL_IS_EVEN(j))
          goto last_pairing; /* [[i,i],j] = 2 [[i,j],i] */
      }       /* IS_ODD(j) -> [[i,i],j] = - 2 [[i,j],i] */
      INTEGER_SET_MINUS(two);                /* 2 -> -2 */
      last_pairing:
      a = PairSumMonomialInteger(a, i);
      if(a != NIL)
        LieSumMultInteger(a, two);
    }
  }
  else if(LIE_MONOMIAL_IS_COMMUTATOR(i) &&  /* Irregular triple */
          LIE_MONOMIAL_ORDER(j) < LIE_MONOMIAL_RIGHT_ORDER(i))
  {
    uint b;              /* i > k > j  =>  [[i,k],j]] = */
    k = LIE_MONOMIAL_RIGHT(i);
    i = LIE_MONOMIAL_LEFT(i);
    a = PairMonomialMonomialInteger(i, j);
    a = PairSumMonomialInteger(a, k);    /* + [[i,j],k] */
    if(LIE_MONOMIAL_IS_ODD(j) && LIE_MONOMIAL_IS_ODD(k))
    {
      uint c = a;                          /* - [[i,j],k] */
      while(c != NIL)
      {
        LIE_TERM_MINUS_INTEGER(c);
        c = LIE_TERM_R(c);
      }
    }
    b = PairMonomialMonomialInteger(k, j);
    b = PairSumMonomialInteger(b, i);    /* + [[k,j],i] */
    if(LIE_MONOMIAL_IS_EVEN(i) ||
       LIE_MONOMIAL_PARITY(j) == LIE_MONOMIAL_PARITY(k))
    {
      uint c = b;                          /* - [[k,j],i] */
      while(c != NIL)
      {
        LIE_TERM_MINUS_INTEGER(c);
        c = LIE_TERM_R(c);
      }
    }
    a = LieSumAddition(a, b);
  }
  else if(i == j && LIE_MONOMIAL_IS_EVEN(i))
    a = NIL;                    /* [i,i] = 0 for even i */
  else
  {                                     /* Regular pair */
    k = AddPairToLieMonomial(i, j);
    if(LIE_MONOMIAL_IS_LEADING(k))
      a = MakeRelationRHSInteger(LIE_MONOMIAL_I_RELATION(k));
    else
    {
      BIGINT n;
      a = NodeLTNew();
      LIE_TERM_MONOMIAL(a) = k;
      INTEGER_HEAP_NEW(n, 2);
      n[0] = n[1] = 1;
      LIE_TERM_NUMERATOR_INTEGER(a) = n;
      LIE_TERM_DENOMINATOR_INTEGER(a) = NULL;
    }
  }
OUT_PAIR_MONOMIAL_MONOMIAL  /*----------------------------*/
  return a;
}
/*=PairMonomialMonomialParametric==============================
 Make regular expression from two monomials (Parametric regime)
 Caller ensures i >= j by Shirshov
*/
uint PairMonomialMonomialParametric(int i, int j)
{
  uint a;
  int k;
IN_PAIR_MONOMIAL_MONOMIAL  /*--------------------------------*/
  if(LIE_MONOMIAL_IS_SQUARE(j))
  {                               /* [i,[j,j]] = 2 [[i,j],j] */
    j = LIE_MONOMIAL_LEFT(j);
    a = PairMonomialMonomialParametric(i, j);
    a = PairSumMonomialParametric(a, j);
    if(a != NIL)
    {
      uint two;
      BIGINT n2;
      INTEGER_HEAP_NEW(n2, 2);
      n2[0] = 1;
      n2[1] = 2;
      two = NodeSTNew();
      SCALAR_TERM_MONOMIAL(two) = NIL;
      SCALAR_TERM_NUMERATOR(two) = n2;
      LieSumMultScalarSum(a, two);
    }
  }
  else if(LIE_MONOMIAL_IS_SQUARE(i))
  {
    BIGINT n2;
    i = LIE_MONOMIAL_LEFT(i);
    if(i == j)
      a = NIL;                              /* [[i,i],i] = 0 */
    else
    {
      INTEGER_HEAP_NEW(n2, 2);
      n2[0] = 1;
      n2[1] = 2;
      if(LIE_MONOMIAL_ORDER(i) < LIE_MONOMIAL_ORDER(j))
        a = PairMonomialMonomialParametric(j, i);
                                /* [[i,i],j] = - 2 [[j,i],i] */
      else
      {
        a = PairMonomialMonomialParametric(i, j);
        if(LIE_MONOMIAL_IS_EVEN(j))
          goto last_pairing;      /* [[i,i],j] = 2 [[i,j],i] */
      }            /* IS_ODD(j) -> [[i,i],j] = - 2 [[i,j],i] */
      INTEGER_SET_MINUS(n2);                      /* 2 -> -2 */
      last_pairing:
      a = PairSumMonomialParametric(a, i);
      if(a != NIL)
      {
        uint two = NodeSTNew();
        SCALAR_TERM_MONOMIAL(two) = NIL;
        SCALAR_TERM_NUMERATOR(two) = n2;
        LieSumMultScalarSum(a, two);
      }
      else
        INTEGER_KILL(n2);
    }
  }
  else if(LIE_MONOMIAL_IS_COMMUTATOR(i) &&  /* Irregular triple */
          LIE_MONOMIAL_ORDER(j) < LIE_MONOMIAL_RIGHT_ORDER(i))
  {
    uint b;                     /* i > k > j  =>  [[i,k],j]] = */
    k = LIE_MONOMIAL_RIGHT(i);
    i = LIE_MONOMIAL_LEFT(i);
    a = PairMonomialMonomialParametric(i, j);
    a = PairSumMonomialParametric(a, k);      /* + [[i,j],k] */
    if(LIE_MONOMIAL_IS_ODD(j) && LIE_MONOMIAL_IS_ODD(k))
      LieSumMinusParametric(a);               /* - [[i,j],k] */
    b = PairMonomialMonomialParametric(k, j);
    b = PairSumMonomialParametric(b, i);      /* + [[k,j],i] */
    if(LIE_MONOMIAL_IS_EVEN(i) ||
       LIE_MONOMIAL_PARITY(j) == LIE_MONOMIAL_PARITY(k))
      LieSumMinusParametric(b);               /* - [[k,j],i] */
    a = LieSumAddition(a, b);
  }
  else if(i == j && LIE_MONOMIAL_IS_EVEN(i))
    a = NIL;                         /* [i,i] = 0 for even i */
  else
  {                                          /* Regular pair */
    k = AddPairToLieMonomial(i, j);
    if(LIE_MONOMIAL_IS_LEADING(k))
      a = MakeRelationRHSParametric(LIE_MONOMIAL_I_RELATION(k));
    else
    {
      BIGINT n;
      uint c;
      a = NodeLTNew();
      LIE_TERM_MONOMIAL(a) = k;
      INTEGER_HEAP_NEW(n, 2);
      n[0] = n[1] = 1;
      c = NodeSTNew();
      SCALAR_TERM_MONOMIAL(c) = NIL;
      SCALAR_TERM_NUMERATOR(c) = n;
      LIE_TERM_NUMERATOR_SCALAR_SUM(a) = c;
      LIE_TERM_DENOMINATOR_SCALAR_SUM(a) = NIL;
    }
  }
OUT_PAIR_MONOMIAL_MONOMIAL  /*--------------------------------*/
  return a;
}
/*=PairMonomialSumInteger==========================================

 Make commutator of the form [mon, Lie_sum] (Integer regime)
*/
uint PairMonomialSumInteger(int mon, uint a)
{
  uint b, s;
  BIGINT nb, db;                        /* Sum has definite parity */
  int monb, change_sign = (LIE_MONOMIAL_IS_EVEN(mon) ||
                      LIE_MONOMIAL_IS_EVEN(LIE_TERM_MONOMIAL(a)));
IN_PAIR_MONOMIAL_SUM  /*----------------------------------------*/
  s = NIL;
  while(a != NIL)
  {
    b = a;
    a = LIE_TERM_R(a);
    monb = LIE_TERM_MONOMIAL(b);     /* Take full info from `b' */
    nb = LIE_TERM_NUMERATOR_INTEGER(b);
    db = LIE_TERM_DENOMINATOR_INTEGER(b);
    NODE_LT_KILL(b);
    if(mon != monb || LIE_MONOMIAL_IS_ODD(mon))
    {                                       /* [mon, monb] != 0 */
      if(LIE_MONOMIAL_ORDER(mon) < LIE_MONOMIAL_ORDER(monb))
      {                                       /* Swap monomials */
        if(change_sign)
          INTEGER_MINUS(nb);
        b = PairMonomialMonomialInteger(monb, mon);
      }
      else
        b = PairMonomialMonomialInteger(mon, monb);
      if(INTEGER_IS_NOT_UNIT(nb))
        LieSumMultInteger(b, nb);
      if(db != NULL)
        LieSumDivInteger(b, db);
      s = LieSumAddition(s, b);
    }
    INTEGER_KILL(nb);
    if(db != NULL)
      INTEGER_KILL(db);
  }
OUT_PAIR_MONOMIAL_SUM  /*---------------------------------------*/
  return s;
}
/*=PairMonomialSumParametric======================================
 Make commutator of the form [mon, Lie_sum] (Parametric regime)
*/
uint PairMonomialSumParametric(int mon, uint a)
{
  uint b, s, nb, db;                    /* Sum has definite parity */
  int monb, change_sign = (LIE_MONOMIAL_IS_EVEN(mon) ||
                      LIE_MONOMIAL_IS_EVEN(LIE_TERM_MONOMIAL(a)));
IN_PAIR_MONOMIAL_SUM   /*---------------------------------------*/
  s = NIL;
  while(a != NIL)
  {
    b = a;
    a = LIE_TERM_R(a);
    monb = LIE_TERM_MONOMIAL(b);     /* Take full info from `b' */
    nb = LIE_TERM_NUMERATOR_SCALAR_SUM(b);
    db = LIE_TERM_DENOMINATOR_SCALAR_SUM(b);
    NODE_LT_KILL(b);
    if(mon != monb || LIE_MONOMIAL_IS_ODD(mon))
    {                                       /* [mon, monb] != 0 */
      if(LIE_MONOMIAL_ORDER(mon) < LIE_MONOMIAL_ORDER(monb))
      {                                       /* Swap monomials */
        if(change_sign)
          ScalarSumMinus(nb);
        b = PairMonomialMonomialParametric(monb, mon);
      }
      else
        b = PairMonomialMonomialParametric(mon, monb);
      if(SCALAR_SUM_IS_NOT_UNIT(nb))
        LieSumMultScalarSum(b, nb);
      else
        ScalarSumKill(nb);
      if(db != NIL)
        LieSumDivScalarSum(b, db);
      s = LieSumAddition(s, b);
    }
    else
    {
      ScalarSumKill(nb);
      ScalarSumKill(db);
    }
  }
OUT_PAIR_MONOMIAL_SUM  /*---------------------------------------*/
  return s;
}
/*=PairSumMonomialInteger==========================================

 Make commutator of the form [Lie_sum, mon] (Integer regime)
*/
uint PairSumMonomialInteger(uint a, int mon)
{
  uint b, s;
  BIGINT nb, db;                        /* Sum has definite parity */
  int monb, change_sign = (LIE_MONOMIAL_IS_EVEN(mon) ||
                      LIE_MONOMIAL_IS_EVEN(LIE_TERM_MONOMIAL(a)));
IN_PAIR_SUM_MONOMIAL   /*---------------------------------------*/
  s = NIL;
  while(a != NIL)
  {
    b = a;
    a = LIE_TERM_R(a);
    monb = LIE_TERM_MONOMIAL(b);     /* Take full info from `b' */
    nb = LIE_TERM_NUMERATOR_INTEGER(b);
    db = LIE_TERM_DENOMINATOR_INTEGER(b);
    NODE_LT_KILL(b);
    if(mon != monb || LIE_MONOMIAL_IS_ODD(mon))
    {                                       /* [monb, mon] != 0 */
      if(LIE_MONOMIAL_ORDER(monb) < LIE_MONOMIAL_ORDER(mon))
      {                                       /* Swap monomials */
        if(change_sign)
          INTEGER_MINUS(nb);
        b = PairMonomialMonomialInteger(mon, monb);
      }
      else
        b = PairMonomialMonomialInteger(monb, mon);
      if(INTEGER_IS_NOT_UNIT(nb))
        LieSumMultInteger(b, nb);
      if(db != NULL)
        LieSumDivInteger(b, db);
      s = LieSumAddition(s, b);
    }
    INTEGER_KILL(nb);
    if(db != NULL)
      INTEGER_KILL(db);
  }
OUT_PAIR_SUM_MONOMIAL  /*---------------------------------------*/
  return s;
}
/*=PairSumMonomialParametric======================================
 Make commutator of the form [Lie_sum, mon] (Parametric regime)
*/
uint PairSumMonomialParametric(uint a, int mon)
{
  uint b, s, nb, db;                    /* Sum has definite parity */
  int monb, change_sign = (LIE_MONOMIAL_IS_EVEN(mon) ||
                      LIE_MONOMIAL_IS_EVEN(LIE_TERM_MONOMIAL(a)));
IN_PAIR_SUM_MONOMIAL   /*---------------------------------------*/
  s = NIL;
  while(a != NIL)
  {
    b = a;
    a = LIE_TERM_R(a);
    monb = LIE_TERM_MONOMIAL(b);     /* Take full info from `b' */
    nb = LIE_TERM_NUMERATOR_SCALAR_SUM(b);
    db = LIE_TERM_DENOMINATOR_SCALAR_SUM(b);
    NODE_LT_KILL(b);
    if(mon != monb || LIE_MONOMIAL_IS_ODD(mon))
    {                                       /* [monb, mon] != 0 */
      if(LIE_MONOMIAL_ORDER(monb) < LIE_MONOMIAL_ORDER(mon))
      {                                       /* Swap monomials */
        if(change_sign)
          ScalarSumMinus(nb);
        b = PairMonomialMonomialParametric(mon, monb);
      }
      else
        b = PairMonomialMonomialParametric(monb, mon);
      if(SCALAR_SUM_IS_NOT_UNIT(nb))
        LieSumMultScalarSum(b, nb);
      else
        ScalarSumKill(nb);
      if(db != NIL)
        LieSumDivScalarSum(b, db);
      s = LieSumAddition(s, b);
    }
  }
OUT_PAIR_SUM_MONOMIAL  /*---------------------------------------*/
  return s;
}
/*=PairSumSumInteger======================================
 Commutator of the form [Lie_sum,Lie_sum] (Integer regime)
*/
uint PairSumSumInteger(uint a, uint b)
{
IN_PAIR_SUM_SUM  /*-------------------------------------*/
  if(a == NIL)
    LieSumKillInteger(b);
  else if(b == NIL)
  {
    LieSumKillInteger(a);
    a = b;
  }
  else
  {
    BIGINT num, den;
    uint c, d, s;
    s = NIL;
    do
    {
      c = a;
      a = LIE_TERM_R(c);
      d = (a != NIL) ? LieSumCopyInteger(b) : b;
      d = PairMonomialSumInteger(LIE_TERM_MONOMIAL(c), d);
      num = LIE_TERM_NUMERATOR_INTEGER(c);
      den = LIE_TERM_DENOMINATOR_INTEGER(c);
      NODE_LT_KILL(c);
      if(INTEGER_IS_NOT_UNIT(num))
        LieSumMultInteger(d, num);
      INTEGER_KILL(num);
      if(den != NULL)
      {
        LieSumDivInteger(d, den);
        INTEGER_KILL(den);
      }
      s = LieSumAddition(s, d);
    }while(a != NIL);
    a = s;
  }
OUT_PAIR_SUM_SUM /*-------------------------------------*/
  return a;
}
/*=PairSumSumParametric======================================
 Commutator of the form [Lie_sum,Lie_sum] (Parametric regime)
*/
uint PairSumSumParametric(uint a, uint b)
{
IN_PAIR_SUM_SUM  /*----------------------------------------*/
  if(a == NIL)
    LieSumKillParametric(b);
  else if(b == NIL)
  {
    LieSumKillParametric(a);
    a = b;
  }
  else
  {
    uint num, den, c, d, s;
    s = NIL;
    do
    {
      c = a;
      a = LIE_TERM_R(c);
      d = (a != NIL) ? LieSumCopyParametric(b) : b;
      d = PairMonomialSumParametric(LIE_TERM_MONOMIAL(c), d);
      num = LIE_TERM_NUMERATOR_SCALAR_SUM(c);
      den = LIE_TERM_DENOMINATOR_SCALAR_SUM(c);
      NODE_LT_KILL(c);
      if(SCALAR_SUM_IS_NOT_UNIT(num))
        LieSumMultScalarSum(d, num);
      else
        ScalarSumKill(num);
      if(den != NIL)
        LieSumDivScalarSum(d, den);
      s = LieSumAddition(s, d);
    }while(a != NIL);
    a = s;
  }
OUT_PAIR_SUM_SUM /*----------------------------------------*/
  return a;
}

/*_6_2          Substitution (replacing) functions======================*/

/*=IsMonomialInMonomial==========================================

 Check whether `mon' contains `submon'
*/
int IsMonomialInMonomial(int submon, int mon)
{
  if(LIE_MONOMIAL_ORDER(submon) > LIE_MONOMIAL_ORDER(mon))
    return NO;
  if(submon == mon)
    return YES;
  if(LIE_MONOMIAL_IS_GENERATOR(mon))
    return NO;
  return IsMonomialInMonomial(submon, LIE_MONOMIAL_LEFT(mon)) ||
         IsMonomialInMonomial(submon, LIE_MONOMIAL_RIGHT(mon));
}
/*=SubstituteRelationInRelationInteger================================
 R is donor and unchanged, A is acceptor and changed. Integer regime
*/
uint SubstituteRelationInRelationInteger(uint r, uint a)
{
  uint b, bl, bf, rhs;
  BIGINT nb;
#if defined(RATIONAL_FIELD)
  BIGINT db;
#endif
  int lmon, monb, lord;
IN_SUBSTITUTE_RELATION_IN_RELATION /*-------------------------------*/
  lmon = LIE_TERM_MONOMIAL(r);
  lord = LIE_MONOMIAL_ORDER(lmon);
  SubstitutionIsDone = NO;
  bl = NIL;
  bf = b = a;
  do
  {
    monb = LIE_TERM_MONOMIAL(b);
    if(LIE_MONOMIAL_ORDER(monb) < lord)
    {
      if(SubstitutionIsDone)
      {
        a = LieSumAddition(a, bf);
        LieSumKillInteger(rhs);
      }
      goto out;
    }
    if(IsMonomialInMonomial(lmon, monb))
    {
      if(SubstitutionIsDone == NO)        /* First substituent term */
      {                                /* Make right hand side of R */
        rhs = LieSumCopyIntegerNegative(LIE_TERM_R(r));
#if !defined(RATIONAL_FIELD)
        if(rhs != NIL)
        {
          nb = LIE_TERM_NUMERATOR_INTEGER(r);
          if(INTEGER_IS_NOT_UNIT(nb))
          {
            BIGINT den;
            int i;                   /* Divide by leading coefficient */
            INTEGER_STACK_COPY(den, nb, i);
            LieSumDivInteger(rhs, den);
          }
        }
#endif
        SubstitutionIsDone = YES;
        a = NIL;
      }
      nb = LIE_TERM_NUMERATOR_INTEGER(b);
      if((r = SubstituteRHSInMonomialInteger(monb, lmon, rhs)) != NIL)
      {
#if defined(RATIONAL_FIELD)                 /* Field R case compiling */
        db = LIE_TERM_DENOMINATOR_INTEGER(b);
        if(INTEGER_IS_UNIT(nb))
        {
          if(db != NULL)
          {
            LieSumDivInteger(r, db);
            INTEGER_KILL(db);
          }
        }
        else
          if(db != NULL)
          {
            LieSumMultRationalInteger(r, nb, db);
            INTEGER_KILL(db);
          }
          else
            LieSumMultInteger(r, nb);
#else                                         /* Ring Z case compiling */
        if(INTEGER_IS_NOT_UNIT(nb))
          LieSumMultInteger(r, nb);
#endif
        a = LieSumAddition(a, r);
      }
      INTEGER_KILL(nb);
      if(bl != NIL)
      {                /* There is substitution-free sublist before */
        LIE_TERM_R(bl) = NIL;
        a = LieSumAddition(a, bf);
      }
      bl = b;
      bf = b = LIE_TERM_R(b);
      NODE_LT_KILL(bl);
      bl = NIL;
    }
    else
    {                                /* Skip substitution-free term */
      bl = b;
      b = LIE_TERM_R(b);
    }
  }while(b != NIL);
  if(SubstitutionIsDone)
  {                                /* Append substitution-free tail */
    if(bl != NIL)
      a = LieSumAddition(a, bf);
    LieSumKillInteger(rhs);
  }
  out:
OUT_SUBSTITUTE_RELATION_IN_RELATION /*------------------------------*/
  return a;
}
/*=SubstituteRelationInRelationParametric================================
 R is donor and unchanged, A is acceptor and changed. Parametric regime
*/
uint SubstituteRelationInRelationParametric(uint r, uint a)
{
  uint b, bl, bf, rhs, pb;
  int lmon, monb, lord;
IN_SUBSTITUTE_RELATION_IN_RELATION /*----------------------------------*/
  lmon = LIE_TERM_MONOMIAL(r);
  lord = LIE_MONOMIAL_ORDER(lmon);
  SubstitutionIsDone = NO;
  bl = NIL;
  bf = b = a;
  do
  {
    monb = LIE_TERM_MONOMIAL(b);
    if(LIE_MONOMIAL_ORDER(monb) < lord)
    {
      if(SubstitutionIsDone)
      {
        a = LieSumAddition(a, bf);
        LieSumKillParametric(rhs);
      }
      goto out;
    }
    if(IsMonomialInMonomial(lmon, monb))
    {
      if(SubstitutionIsDone == NO)           /* First substituent term */
      {                                   /* Make right hand side of R */
        if((rhs = LieSumCopyParametric(LIE_TERM_R(r))) != NIL)
        {
          pb = LIE_TERM_NUMERATOR_SCALAR_SUM(r);
          if(SCALAR_SUM_IS_NOT_UNIT(pb))   /* Divide by leading coeff. */
            LieSumDivScalarSum(rhs, ScalarSumCopy(pb));
          LieSumMinusParametric(rhs);                      /* Negation */
        }
        SubstitutionIsDone = YES;
        a = NIL;
      }
      pb = LIE_TERM_NUMERATOR_SCALAR_SUM(b);
      if((r = SubstituteRHSInMonomialParametric(monb, lmon, rhs)) != NIL)
      {
        if(SCALAR_SUM_IS_NOT_UNIT(pb))
          LieSumMultScalarSum(r, pb);
        else
          ScalarSumKill(pb);
        a = LieSumAddition(a, r);
      }
      if(bl != NIL)
      {                   /* There is substitution-free sublist before */
        LIE_TERM_R(bl) = NIL;
        a = LieSumAddition(a, bf);
      }
      bl = b;
      bf = b = LIE_TERM_R(b);
      NODE_LT_KILL(bl);
      bl = NIL;
    }
    else
    {                                   /* Skip substitution-free term */
      bl = b;
      b = LIE_TERM_R(b);
    }
  }while(b != NIL);
  if(SubstitutionIsDone)
  {                                   /* Append substitution-free tail */
    if(bl != NIL)
      a = LieSumAddition(a, bf);
    LieSumKillParametric(rhs);
  }
  out:
OUT_SUBSTITUTE_RELATION_IN_RELATION /*---------------------------------*/
  return a;
}
/*=SubstituteRHSInMonomialInteger=======================================
 Insert right hand side `r' of relation with leading monomial `lmonr'
 in monomial `mon'.
 Returned NOTHING means "no substitution", NIL means 0.
 Function saves input `r'. Integer regime.
*/
uint SubstituteRHSInMonomialInteger(int mon, int lmonr, uint r)
{
  /* Single monomial matching case */

  if(mon == lmonr)
    return LieSumCopyInteger(r);

  /* Possible substitution(s) in submonomial(s) */

  if(LIE_MONOMIAL_ORDER(mon) > LIE_MONOMIAL_ORDER(lmonr)
     && LIE_MONOMIAL_IS_COMMUTATOR(mon))
  {
    uint res;
    int monl = LIE_MONOMIAL_LEFT(mon);
    mon = LIE_MONOMIAL_RIGHT(mon);
    res = SubstituteRHSInMonomialInteger(monl, lmonr, r);
    if(res == NIL)
      return NIL;                                     /*  [0, x] -> 0 */
    {
      uint resr = SubstituteRHSInMonomialInteger(mon, lmonr, r);
      if(res == NOTHING)
      {
        if(resr == NOTHING)
          return NOTHING;    /* No substitutions in both submonomials */
        if(resr == NIL)
          return NIL;                                  /* [x, 0] -> 0 */
        return PairMonomialSumInteger(monl, resr); /* -> [monl, resr] */
      }
      if(resr == NOTHING)
        return PairSumMonomialInteger(res, mon);     /* -> [res, mon] */
      if(resr == NIL)
      {
        LieSumKillInteger(res);
        return NIL;
      }
      return PairSumSumInteger(res, resr);          /* -> [res, resr] */
    }
  }
  return NOTHING;
}
/*=SubstituteRHSInMonomialParametric======================================
 Insert right hand side `r' of relation with leading monomial `lmonr'
 in monomial `mon'.
 Returned NOTHING means "no substitution", NIL means 0.
 Function saves input `r'. Parametric regime.
*/
uint SubstituteRHSInMonomialParametric(int mon, int lmonr, uint r)
{
  /* Single monomial matching case */

  if(mon == lmonr)
    return LieSumCopyParametric(r);

  /* Possible substitution(s) in submonomial(s) */

  if(LIE_MONOMIAL_ORDER(mon) > LIE_MONOMIAL_ORDER(lmonr)
     && LIE_MONOMIAL_IS_COMMUTATOR(mon))
  {
    uint res;
    int monl = LIE_MONOMIAL_LEFT(mon);
    mon = LIE_MONOMIAL_RIGHT(mon);
    res = SubstituteRHSInMonomialParametric(monl, lmonr, r);
    if(res == NIL)
      return NIL;                                        /*  [0, x] -> 0 */
    {
      uint resr = SubstituteRHSInMonomialParametric(mon, lmonr, r);
      if(res == NOTHING)
      {
        if(resr == NOTHING)
          return NOTHING;       /* No substitutions in both submonomials */
        if(resr == NIL)
          return NIL;                                     /* [x, 0] -> 0 */
        return PairMonomialSumParametric(monl, resr); /* -> [monl, resr] */
      }
      if(resr == NOTHING)
        return PairSumMonomialParametric(res, mon);     /* -> [res, mon] */
      if(resr == NIL)
      {
        LieSumKillParametric(res);
        return NIL;
      }
      return PairSumSumParametric(res, resr);          /* -> [res, resr] */
    }
  }
  return NOTHING;
}

/*_6_3 Lie and scalar algebra functions========================*/

/*=LieLikeTermsCollectionInteger==========================================
 For Lie terms `a' and `b' sum rational integers  `na/da' and `nb/db'
 destructing input terms, set nonzero result in `a', return YES for
 nonzero result, otherwise kill `a' and return NO;
 `da', `db'= NIL means 1.
*/
int LieLikeTermsCollectionInteger(uint a, uint b)
{
  BIGINT na, da, nb, db;
  nb = LIE_TERM_NUMERATOR_INTEGER(b);
  db = LIE_TERM_DENOMINATOR_INTEGER(b);
  NODE_LT_KILL(b);
  na = LIE_TERM_NUMERATOR_INTEGER(a);
  da = LIE_TERM_DENOMINATOR_INTEGER(a);
  if(da != NULL) if(db != NULL)                 /* `da' != 1, `db' != 1 */
  {
    BIGINT g, h;
    int i;
    INTEGER_STACK_COPY(g, da, i);
    INTEGER_STACK_COPY(h, db, i);
    if((g = IntegerGCD(g, h)) != NULL)           /* g = GCD(da, db) > 1 */
    {
      BIGINT k, m, daa;
      INTEGER_STACK_COPY(k, g, i);                 /*  k = GCD(da, db)' */
      INTEGER_STACK_NEW(m, 2+INTEGER_N_LIMBS(da)-INTEGER_N_LIMBS(k));
      INTEGER_STACK_COPY_1(daa, da, i);
      INTEGER_KILL(da);
      IntegerQuotient(m, daa, k);               /*  m = da/GCD(da, db)' */
      INTEGER_STACK_NEW(h, 1+INTEGER_N_LIMBS(nb)+INTEGER_N_LIMBS(m));
      IntegerProduct(h, nb, m);             /*  h = nb*da'/GCD(da, db)' */
      INTEGER_KILL(nb);
      INTEGER_STACK_COPY_1(k, db, i);                       /*  k = db' */
      INTEGER_STACK_COPY(da, g, i);                /* da = GCD(da, db)' */
      INTEGER_STACK_NEW(nb, 2+INTEGER_N_LIMBS(k)-INTEGER_N_LIMBS(da));
      IntegerQuotient(nb, k, da);              /* nb = db'/GCD(da, db)' */
      INTEGER_STACK_NEW(k, 1+INTEGER_N_LIMBS(na)+INTEGER_N_LIMBS(nb));
      IntegerProduct(k, na, nb);            /*  k = na*db'/GCD(da, db)' */
      INTEGER_KILL(na);
      INTEGER_STACK_NEW(na, 3+MAX(INTEGER_N_LIMBS(h),INTEGER_N_LIMBS(k)));
      IntegerSum(na, h, k);                               /* na = h + k */
      if(INTEGER_N_LIMBS(na) != 0)
      {
        INTEGER_STACK_COPY(h, na, i);                      /* Numerator */
        if((h = IntegerGCD(h, g)) != NULL)
        {
          INTEGER_STACK_COPY(g, h, i);
          INTEGER_HEAP_NEW(k, 2+INTEGER_N_LIMBS(na)-INTEGER_N_LIMBS(g));
          IntegerQuotient(k, na, g);
          LIE_TERM_NUMERATOR_INTEGER(a) = k;
          INTEGER_STACK_NEW(k, 2+INTEGER_N_LIMBS(db)-INTEGER_N_LIMBS(h));
          INTEGER_STACK_COPY_1(daa, db, i);
          INTEGER_KILL(db);
          IntegerQuotient(k, daa, h);
          INTEGER_HEAP_NEW(db, 1+INTEGER_N_LIMBS(m)+INTEGER_N_LIMBS(k));
          IntegerProduct(db, m, k);
          if(INTEGER_IS_UNIT(db))
          {
            INTEGER_KILL(db);
            db = NULL;     /* Standard convention for unit denominators */
          }
          LIE_TERM_DENOMINATOR_INTEGER(a) = db;
        }
        else
        {
          INTEGER_HEAP_COPY(h, na, i);
          LIE_TERM_NUMERATOR_INTEGER(a) = h;
          INTEGER_HEAP_NEW(k, 1+INTEGER_N_LIMBS(m)+INTEGER_N_LIMBS(db));
          IntegerProduct(k, m, db);
          INTEGER_KILL(db);
          if(INTEGER_IS_UNIT(k))
          {
            INTEGER_KILL(k);
            k = NULL;      /* Standard convention for unit denominators */
          }
          LIE_TERM_DENOMINATOR_INTEGER(a) = k;
        }
        goto non_zero;
      }
      else
      {
        INTEGER_KILL(db);
        goto zero;
      }
    }
    else                                       /* Mutually prime da, db */
    {
      INTEGER_STACK_NEW(g, 1+INTEGER_N_LIMBS(nb)+INTEGER_N_LIMBS(da));
      IntegerProduct(g, nb, da);                           /* g = nb*da */
      INTEGER_KILL(nb);
      INTEGER_STACK_NEW(nb, 1+INTEGER_N_LIMBS(na)+INTEGER_N_LIMBS(db));
      IntegerProduct(nb, na, db);                         /* nb = na*db */
      INTEGER_KILL(na);
      INTEGER_HEAP_NEW(na, 2+MAX(INTEGER_N_LIMBS(nb),INTEGER_N_LIMBS(g)));
      IntegerSum(na, nb, g);                      /* na = nb*da + na*db */
      if(INTEGER_N_LIMBS(na) != 0)
      {
        LIE_TERM_NUMERATOR_INTEGER(a) = na;
        INTEGER_HEAP_NEW(nb, 1+INTEGER_N_LIMBS(da)+INTEGER_N_LIMBS(db));
        IntegerProduct(nb, da, db);                       /* nb = da*db */
        LIE_TERM_DENOMINATOR_INTEGER(a) = nb;
        INTEGER_KILL(da);
        INTEGER_KILL(db);
        goto non_zero;
      }
      else
      {
        INTEGER_KILL(na);
        INTEGER_KILL(da);
        INTEGER_KILL(db);
        goto zero;
      }
    }
  }
  else                                           /* `da' != 1, `db' = 1 */
  {
    INTEGER_STACK_NEW(db, 1+INTEGER_N_LIMBS(nb)+INTEGER_N_LIMBS(da));
    IntegerProduct(db, nb, da);
    INTEGER_KILL(nb);
    INTEGER_HEAP_NEW(nb, 2+MAX(INTEGER_N_LIMBS(db),INTEGER_N_LIMBS(na)));
    IntegerSum(nb, db, na);
    INTEGER_KILL(na);
    LIE_TERM_NUMERATOR_INTEGER(a) = nb;
    goto non_zero;
  }
  else if(db != NULL)                            /* `da' = 1, `db' != 1 */
  {
    INTEGER_STACK_NEW(da, 1+INTEGER_N_LIMBS(na)+INTEGER_N_LIMBS(db));
    IntegerProduct(da, na, db);
    INTEGER_KILL(na);
    INTEGER_HEAP_NEW(na, 2+MAX(INTEGER_N_LIMBS(da),INTEGER_N_LIMBS(nb)));
    IntegerSum(na, da, nb);
    INTEGER_KILL(nb);
    LIE_TERM_NUMERATOR_INTEGER(a) = na;
    LIE_TERM_DENOMINATOR_INTEGER(a) = db;
    goto non_zero;
  }
  else                                               /* `da' = `db' = 1 */
  {
    INTEGER_HEAP_NEW(da, 2+MAX(INTEGER_N_LIMBS(na),INTEGER_N_LIMBS(nb)));
    IntegerSum(da, na, nb);
    INTEGER_KILL(na);
    INTEGER_KILL(nb);
    if(INTEGER_N_LIMBS(da) != 0)
    {
      LIE_TERM_NUMERATOR_INTEGER(a) = da;
      goto non_zero;
    }
    else
    {
      INTEGER_KILL(da);
      goto zero;
    }
  }
  non_zero:
  return YES;
  zero:
  NODE_LT_KILL(a);                                   /* `na' + `nb' = 0 */
  return NO;
}
/*=LieLikeTermsCollectionParametric=======================================
 For Lie terms `a' and `b' sum rational functions  `na/da' and `nb/db'
 destructing input terms, set nonzero result in `a', return YES for
 nonzero result, otherwise kill `a' and return NO;
 `da', `db'= NIL means 1.
*/
int LieLikeTermsCollectionParametric(uint a, uint b)
{
  uint na, da, nb, db;
  nb = LIE_TERM_NUMERATOR_SCALAR_SUM(b);
  db = LIE_TERM_DENOMINATOR_SCALAR_SUM(b);
  NODE_LT_KILL(b);
  na = LIE_TERM_NUMERATOR_SCALAR_SUM(a);
  da = LIE_TERM_DENOMINATOR_SCALAR_SUM(a);
  if(da != NIL) if(db != NIL)                /* `da' != 1 and `db' != 1 */
  {
    uint g;
    if((g = PolyGCD(da, db)) != NIL)            /* g = GCD(da, db) != 1 */
    {
      da = PolyQuotient(da, g);                /*  da' = da/GCD(da, db) */
                                                       /*  nb' = nb*da' */
      nb = ScalarSumMultiplication(nb, ScalarSumCopy(da));
                                            /*  na' = na*db/GCD(da, db) */
      na = ScalarSumMultiplication(na, PolyQuotient(ScalarSumCopy(db), g));
      na = ScalarSumAddition(na, nb);               /* na'' = na' + nb' */
      if(na != NIL)
      {
        if((nb = PolyGCD(na, g)) != NIL)
        {                                      /* Set na''/GCD(na'', g) */
          LIE_TERM_NUMERATOR_SCALAR_SUM(a) = PolyQuotient(na, nb);
          db = PolyQuotient(db, nb);           /* db' = db/GCD(na'', g) */
          ScalarSumKill(nb);                       /* Kill GCD(na'', g) */
        }
        else
          LIE_TERM_NUMERATOR_SCALAR_SUM(a) = na;
        ScalarSumKill(g);                           /* Kill GCD(da, db) */
        da = ScalarSumMultiplication(da, db);
        if(SCALAR_SUM_IS_UNIT(da))
        {
          ScalarSumKill(da);
          da = NIL;
        }
        LIE_TERM_DENOMINATOR_SCALAR_SUM(a) = da;
        goto non_zero;
      }
      ScalarSumKill(g);             /* Nontrivial g at na == NIL branch */
    }
    else                                       /* Mutually prime da, db */
    {                                          /* na' = na*db' + nb*da' */
      na = ScalarSumAddition(ScalarSumMultiplication(na, ScalarSumCopy(db)),
                             ScalarSumMultiplication(nb, ScalarSumCopy(da)));
      if(na != NIL)
      {
        LIE_TERM_NUMERATOR_SCALAR_SUM(a) = na;
        LIE_TERM_DENOMINATOR_SCALAR_SUM(a) = ScalarSumMultiplication(da, db);
        goto non_zero;
      }
    }
    ScalarSumKill(da);                              /* na == NIL branch */
    ScalarSumKill(db);
    goto zero;
  }
  else                   /* `da' != 1 and `db' = 1 --> (na + nb*da')/da */
  {
    LIE_TERM_NUMERATOR_SCALAR_SUM(a) =
      ScalarSumAddition(na, ScalarSumMultiplication(nb, ScalarSumCopy(da)));
    goto non_zero;
  }
  else if(db != NIL)     /* `da' = 1 and `db' != 1 --> (nb + na*db')/db */
  {
    LIE_TERM_NUMERATOR_SCALAR_SUM(a) =
      ScalarSumAddition(nb, ScalarSumMultiplication(na, ScalarSumCopy(db)));
    LIE_TERM_DENOMINATOR_SCALAR_SUM(a) = db;
    goto non_zero;
  }
  else if((na = ScalarSumAddition(na, nb)) != NIL)
  {                           /* `da' = `db' = 1 --> (na + nb)/1 (!= 0) */
    LIE_TERM_NUMERATOR_SCALAR_SUM(a) = na;
    goto non_zero;
  }
  zero:
  NODE_LT_KILL(a);  /* `na' + `nb' = 0 */
  return NO;
  non_zero:
  return YES;
}
/*=LieSumAddition=============================================
 Sum of two Lie expressions
*/
uint LieSumAddition(uint a, uint b)
{
  uint sum = NIL, last, wa, wb;
IN_LIE_SUM_ADDITION /*---------------------------------------*/
  while(YES)
  {
    next_pair:
    if(b == NIL)
    {                   /* List b is ended, append rest of a */
      if(sum == NIL)
        sum = a;
      else
        LIE_TERM_R(last) = a;
      break;
    }
    if(a == NIL)
    {                   /* List a is ended, append rest of b */
      if(sum == NIL)
        sum = b;
      else
        LIE_TERM_R(last) = b;
      break;
    }

    /* Compare algebra terms */

    if(LIE_TERM_MONOMIAL_ORDER(a) > LIE_TERM_MONOMIAL_ORDER(b))
      goto order_12;
    if(LIE_TERM_MONOMIAL_ORDER(a) < LIE_TERM_MONOMIAL_ORDER(b))
      goto order_21;

    /* Reduce like algebra terms */

    wa = a;
    wb = b;
    a = LIE_TERM_R(a);
    b = LIE_TERM_R(b);

        /* Sum rational coefficients */

    if((*LieLikeTermsCollection)(wa, wb))
      goto append_term;
    else
      goto next_pair;

    order_12:
    wa = a;
    a = LIE_TERM_R(a);
    goto append_term;

    order_21:
    wa = b;
    b = LIE_TERM_R(b);

    append_term:
    if(sum == NIL)
      sum = wa;
    else
      LIE_TERM_R(last) = wa;
    last = wa;
  }
OUT_LIE_SUM_ADDITION /*--------------------------------------*/
  return sum;
}
/*=LieSumDivInteger=======================================================
 Divide Lie sum by integer (of unknown nature) on spot in Integer regime
 Integer `den' is spoiled
*/
void LieSumDivInteger(uint lsum, BIGINT den)
{
  if(lsum != NIL)
  {
    BIGINT d, da, dao;
    int i, n;
    uint a;
IN_LIE_SUM_DIV_INTEGER  /*----------------------------------------------*/
    n = INTEGER_N_LIMBS(den);
    INTEGER_STACK_NEW(d, 1+n);          /* Space for copies input `den' */
    do
    {
      a = lsum;
      lsum = LIE_TERM_R(lsum);
      if(lsum != NIL)
      {
        i = n;
        do
          d[i] = den[i];
        while(i--);
      }
      else
        d = den;
      IntegerCancellation(LIE_TERM_NUMERATOR_INTEGER(a), d);
      if(INTEGER_IS_NOT_UNIT(d))
      {
        if((dao = LIE_TERM_DENOMINATOR_INTEGER(a)) != NULL)
        {                           /* Nontrivial old denominator `dao' */
          INTEGER_HEAP_NEW(da, 1+INTEGER_N_LIMBS(d)+INTEGER_N_LIMBS(dao));
          IntegerProduct(da, d, dao);
          INTEGER_KILL(dao);
        }
        else
        {
          INTEGER_HEAP_COPY(da, d, i); /* Composite statement */
        }
        LIE_TERM_DENOMINATOR_INTEGER(a) = da;
      }
    }while(lsum != NIL);
OUT_LIE_SUM_DIV_INTEGER  /*---------------------------------------------*/
  }
}
/*=LieSumDivScalarSum======================================
 Divide Lie sum by scalar sum on spot in Parametric regime
 `den' is killed
*/
void LieSumDivScalarSum(uint lsum, uint den)
{
  if(lsum == NIL)
    ScalarSumKill(den);
  else
  {
    uint n, d, a;
IN_LIE_SUM_DIV_SCALAR_SUM /*-----------------------------*/
    do
    {
      a = lsum;
      lsum = LIE_TERM_R(lsum);
      d = (lsum != NIL) ? ScalarSumCopy(den) : den;
      n = LIE_TERM_NUMERATOR_SCALAR_SUM(a);
      ScalarSumCancellation(&n, &d);
      LIE_TERM_NUMERATOR_SCALAR_SUM(a) = n;
      if(SCALAR_SUM_IS_NOT_UNIT(d))
      {
        if((n = LIE_TERM_DENOMINATOR_SCALAR_SUM(a)) != NIL)
          d = ScalarSumMultiplication(d, n);   /* Absorb */
        LIE_TERM_DENOMINATOR_SCALAR_SUM(a) = d;
      }
      else
        ScalarSumKill(d);                   /* Kill unit */
    }while(lsum != NIL);
OUT_LIE_SUM_DIV_SCALAR_SUM /*----------------------------*/
  }
}
/*=LieSumMinusInteger====================
 Change signs in Lie sum (Integer regime)
*/
void LieSumMinusInteger(uint a)
{
  while(a != NIL)
  {
    LIE_TERM_MINUS_INTEGER(a);
    a = LIE_TERM_R(a);
  }
}

/*=LieSumMinusParametric====================
 Change signs in Lie sum (Parametric regime)
*/
void LieSumMinusParametric(uint a)
{
  uint b;
  while(a != NIL)
  {
    b = LIE_TERM_NUMERATOR_SCALAR_SUM(a);
    do
      SCALAR_TERM_MINUS(b);
    while((b = SCALAR_TERM_R(b)) != NIL);
    a = LIE_TERM_R(a);
  }
}
/*=LieSumMultInteger=======================================================
 Multiply Lie sum by integer (of unknown nature) on spot in Integer regime
 Integer `num' is spoiled
*/
void LieSumMultInteger(uint lsum, BIGINT num)
{
  if(lsum != NIL)
  {
    BIGINT nw, nao, da;
    int i, n;
    uint a;
IN_LIE_SUM_MULT_INTEGER  /*----------------------------------------------*/
    n = INTEGER_N_LIMBS(num);
    INTEGER_STACK_NEW(nw, 1+n);       /* Space for copies of input `num' */
    do
    {
      a = lsum;
      lsum = LIE_TERM_R(lsum);
      nao = LIE_TERM_NUMERATOR_INTEGER(a);        /* Old numerator `nao' */
      if((da = LIE_TERM_DENOMINATOR_INTEGER(a)) != NULL)
      {                                   /* Nontrivial denominator `da' */
        if(lsum != NIL)
        {                                            /* Copy if not last */
          i = n;
          do
            nw[i] = num[i];
          while(i--);
        }
        else
          nw = num;
        IntegerCancellation(nw, da);
        if(INTEGER_IS_UNIT(da))
        {                                 /* Trivialize unit denominator */
          INTEGER_KILL(da);
          LIE_TERM_DENOMINATOR_INTEGER(a) = NULL;
        }
        if(INTEGER_IS_NOT_UNIT(nw))
        {
          INTEGER_HEAP_NEW(da, 1+INTEGER_N_LIMBS(nw)+INTEGER_N_LIMBS(nao));
          IntegerProduct(da, nw, nao);
          goto stick_new;
        }
      }
      else
      {
        INTEGER_HEAP_NEW(da, 1+INTEGER_N_LIMBS(num)+INTEGER_N_LIMBS(nao));
        IntegerProduct(da, num, nao);
        stick_new:
        INTEGER_KILL(nao);
        LIE_TERM_NUMERATOR_INTEGER(a) = da;
      }
    }while(lsum != NIL);
OUT_LIE_SUM_MULT_INTEGER /*----------------------------------------------*/
  }
}
#if defined(RATIONAL_FIELD)
/*=LieSumMultRationalInteger=============================================
 num and den are non-NULL integers of unknown nature
*/
void LieSumMultRationalInteger(int a, BIGINT num, BIGINT den)
{
  BIGINT  numc, denc, numa, dena, w;
  int i, nn = INTEGER_N_LIMBS(num), nd = INTEGER_N_LIMBS(den);
  INTEGER_STACK_NEW(numc, 1+nn);
  INTEGER_STACK_NEW(denc, 1+nd);
  while(a != NIL)
  {
    for(i = 0; i <= nn; i++)                   /* Copy input numerator */
      numc[i] = num[i];
    if((dena = LIE_TERM_DENOMINATOR_INTEGER(a)) != NULL)
    {                                                         /* n1/d2 */
      IntegerCancellation(numc, dena);
      if(INTEGER_IS_UNIT(dena))
      {
        INTEGER_KILL(dena);
        dena = NULL;
      }
    }
    numa = LIE_TERM_NUMERATOR_INTEGER(a);
    for(i = 0; i <= nd; i++)               /* Copy input denominator */
      denc[i] = den[i];
    IntegerCancellation(numa, denc);                        /* n2/d1 */
    INTEGER_HEAP_NEW(w, 1+INTEGER_N_LIMBS(numc)+INTEGER_N_LIMBS(numa));
    IntegerProduct(w, numc, numa);                            /* n1*n2 */
    INTEGER_KILL(numa);
    LIE_TERM_NUMERATOR_INTEGER(a) = w;
    if(INTEGER_IS_UNIT(denc))
       if(dena == NULL)
         w = NULL;                                              /* 1*1 */
       else
       {                            /* Copy to avoid garbage tail in w */
         INTEGER_HEAP_COPY(w, dena, i);                        /* 1*d2 */
         INTEGER_KILL(dena);
       }
    else if(dena == NULL)
    {
      INTEGER_HEAP_COPY(w, denc, i);                           /* d1*1 */
    }
    else
    {
      INTEGER_HEAP_NEW(w, 1+INTEGER_N_LIMBS(denc)+INTEGER_N_LIMBS(dena));
      IntegerProduct(w, denc, dena);                          /* d1*d2 */
      INTEGER_KILL(dena);
    }
    LIE_TERM_DENOMINATOR_INTEGER(a) = w;
    a = LIE_TERM_R(a);
  }
}
#endif
/*=LieSumMultScalarSum===============================================
 Multiply Lie sum by scalar sum on spot in Parametric regime
 `num' is killed
*/
void LieSumMultScalarSum(uint lsum, uint num)
{
  if(lsum == NIL)
    ScalarSumKill(num);
  else
  {
    uint n, d, a;
IN_LIE_SUM_MULT_SCALAR_SUM /*--------------------------------------*/
    do
    {
      a = lsum;
      lsum = LIE_TERM_R(lsum);
      n = (lsum != NIL) ? ScalarSumCopy(num) : num;
      if((d = LIE_TERM_DENOMINATOR_SCALAR_SUM(a)) != NIL)
      {
        ScalarSumCancellation(&n, &d);
        if(SCALAR_SUM_IS_UNIT(d))
        {
          ScalarSumKill(d);                           /* Kill unit */
          d = NIL;
        }
        LIE_TERM_DENOMINATOR_SCALAR_SUM(a) = d;
      }
      LIE_TERM_NUMERATOR_SCALAR_SUM(a) =
        ScalarSumMultiplication(LIE_TERM_NUMERATOR_SCALAR_SUM(a), n);
    }while(lsum != NIL);
OUT_LIE_SUM_MULT_SCALAR_SUM /*-------------------------------------*/
  }
}
/*=NormalizeRelationInteger================================================
 Normalize sign, remove GCD of integer numerators,
 remove denominators for non-NIL relation
*/
void NormalizeRelationInteger(uint a)
{
  uint b;
  BIGINT n2, n1 = LIE_TERM_NUMERATOR_INTEGER(a);
IN_NORMALIZE_RELATION   /*-----------------------------------------------*/
  /* Normalize sign */

  if(INTEGER_IS_NEGATIVE(n1))
  {
    INTEGER_SET_PLUS(n1);
    b = LIE_TERM_R(a);
    while(b != NIL)
    {
      INTEGER_MINUS(LIE_TERM_NUMERATOR_INTEGER(b));
      b = LIE_TERM_R(b);
    }
  }
#if defined(RATIONAL_FIELD)         /* Field R case compiling */
  n2 = LIE_TERM_DENOMINATOR_INTEGER(a);
  if(INTEGER_IS_UNIT(n1))  /* Either 1 (nothing to do or 1/n2 */
  {
    if(n2 != NULL)        /* Leading coefficient in form 1/n2 */
    {
      LIE_TERM_DENOMINATOR_INTEGER(a) = NULL;
      LieSumMultInteger(LIE_TERM_R(a), n2);
      INTEGER_KILL(n2);              /* Free spoiled array n2 */
    }
  }
  else           /* Leading coefficient is either n1 or n1/n2 */
  {
    if(n2 != NULL)
    {
      LIE_TERM_DENOMINATOR_INTEGER(a) = NULL;
      LieSumMultRationalInteger(LIE_TERM_R(a), n2, n1);
      INTEGER_KILL(n2);              /* Free spoiled array n2 */
    }
    else                    /* Leading coefficient in form n1 */
      LieSumDivInteger(LIE_TERM_R(a), n1);
    n1[0] = n1[1] = 1;    /* Set big number unit in old array */
  }
#else                                /* Ring Z case compiling */
  {
    BIGINT gcd;
    int i;

    /* Remove GCD of numerators */

    if(INTEGER_IS_UNIT(n1))
      goto kill_denominators;
    INTEGER_STACK_COPY(gcd, n1, i);
    b = LIE_TERM_R(a);
    while(b != NIL)
    {
      n1 = LIE_TERM_NUMERATOR_INTEGER(b);
      if(INTEGER_IS_UNIT_ABS(n1))
        goto kill_denominators;
      INTEGER_HEAP_COPY(n2, n1, i);                /* Working heap integer */
      gcd = IntegerGCD(gcd, n2);
      INTEGER_KILL(n2);                               /* Free heap integer */
      if(gcd == NULL)
        goto kill_denominators;
      b = LIE_TERM_R(b);
    }
    LieSumDivInteger(a, gcd);

    /* Remove denominators */

    kill_denominators:
    b = a;
    do
      if(LIE_TERM_DENOMINATOR_INTEGER(b) != NULL)
      {
        BIGINT n3, n4, lcm;

        /* Make first LCM */

        n1 = LIE_TERM_DENOMINATOR_INTEGER(b);
        INTEGER_HEAP_COPY(lcm, n1, i);
        while((b = LIE_TERM_R(b)) != NIL)
          if((n1 = LIE_TERM_DENOMINATOR_INTEGER(b)) != NULL)
          {
            /* Make copy of previous LCM */

            INTEGER_HEAP_COPY(n3, lcm, i);

            /* Make 2 copies of current denominator */

            INTEGER_HEAP_COPY_DOUBLE_1(n2, n4, n1, i);

            /* GCD of LCM and current denominator (stored in `n3') */

            gcd = IntegerGCD(n3, n4);
            INTEGER_KILL(n4);

            /* Divide current denominator by GCD */

            if(gcd != NULL)
            {
              INTEGER_HEAP_NEW(n1, 2+INTEGER_N_LIMBS(n2)-INTEGER_N_LIMBS(gcd));
              IntegerQuotient(n1, n2, gcd);
              INTEGER_KILL(n2);
            }
            else
              n1 = n2;
            INTEGER_KILL(n3);

            /* New LCM */

            n2 = lcm;
            INTEGER_HEAP_NEW(lcm, 1+INTEGER_N_LIMBS(n1)+INTEGER_N_LIMBS(n2));
            IntegerProduct(lcm, n1, n2);
            INTEGER_KILL(n2);
            INTEGER_KILL(n1);
          }

        /* Kill denominators */

        LieSumMultInteger(a, lcm);
        INTEGER_KILL(lcm);
        break;
      }
    while((b = LIE_TERM_R(b)) != NIL);
  }
#endif
OUT_NORMALIZE_RELATION   /*----------------------------------------------*/
}
/*=NormalizeRelationParametric===========================================
 Normalize sign, remove GCD of polynomial numerators, remove denominators
 for non-NIL relation, set in table common factor and leading coefficient
*/
void NormalizeRelationParametric(uint a)
{
  uint b, c = LIE_TERM_NUMERATOR_SCALAR_SUM(a), d, e;
IN_NORMALIZE_RELATION   /*---------------------------------------------*/

  /* Normalize sign */

  if(INTEGER_IS_NEGATIVE(SCALAR_TERM_NUMERATOR(c)))
    LieSumMinusParametric(a);

  /* Remove GCD of numerators */

  if(SCALAR_SUM_IS_UNIT(c))
    goto kill_denominators;
  b = LIE_TERM_R(a);
  c = ScalarSumCopy(c);
  while(b != NIL)
  {
    d = LIE_TERM_NUMERATOR_SCALAR_SUM(b);
    if(SCALAR_SUM_IS_UNIT_ABS(d))
    {
      ScalarSumKill(c);
      goto kill_denominators;
    }
    e = c;
    c = PolyGCD(e, d);
    ScalarSumKill(e);                     /* Kill old GCD */
    if(c == NIL)
      goto kill_denominators;
    b = LIE_TERM_R(b);
  }
  if(NonZeroCoefficientsPut)
    InCoeffTable(ScalarSumCopy(c));
  LieSumDivScalarSum(a, c);

  /* Remove denominators */

  kill_denominators:
  b = a;
  do
    if(LIE_TERM_DENOMINATOR_SCALAR_SUM(b) != NIL)
    {

      /* Make first LCM */

      c = ScalarSumCopy(LIE_TERM_DENOMINATOR_SCALAR_SUM(b));
      while((b = LIE_TERM_R(b)) != NIL)
        if((d = LIE_TERM_DENOMINATOR_SCALAR_SUM(b)) != NIL)
        {
          /* GCD of LCM and current denominator */

          e = PolyGCD(c, d);

          /* Divide current denominator by GCD */

          d = ScalarSumCopy(d);
          if(e != NIL)
          {
            d = PolyQuotient(d, e);
            ScalarSumKill(e);
          }

          /* New LCM */

          c = ScalarSumMultiplication(c, d);
        }

      /* Kill denominators */

      LieSumMultScalarSum(a, c);
      if(NonZeroCoefficientsPut)
        InCoeffTable(ScalarSumCopy(LIE_TERM_NUMERATOR_SCALAR_SUM(a)));
      break;
    }
  while((b = LIE_TERM_R(b)) != NIL);
OUT_NORMALIZE_RELATION   /*--------------------------------------------*/
}
/*=ScalarMonomialMultiplication============================================
*/
uint ScalarMonomialMultiplication(int *pchange_sign, uint ma, uint mb)
{
  uint mc, wa, wb, last;
  *pchange_sign = NO;
  mc = NIL;
  while(YES)
  {
    next_pair:
    if(mb == NIL)
    {                             /* List mb is ended, append rest of ma */
      if(mc == NIL)
        mc = ma;
      else
        SCALAR_FACTOR_R(last) = ma;
      break;
    }
    if(ma == NIL)
    {                             /* List ma is ended, append rest of mb */
      if(mc == NIL)
        mc = mb;
      else
        SCALAR_FACTOR_R(last) = mb;
      break;
    }

    /* Compare scalar factors */

    if(SCALAR_FACTOR_PARAMETER(ma) > SCALAR_FACTOR_PARAMETER(mb))
      goto order_12;
    if(SCALAR_FACTOR_PARAMETER(ma) < SCALAR_FACTOR_PARAMETER(mb))
      goto order_21;

    /* Reduce like factors */

    wa = ma;
    wb = mb;
    ma = SCALAR_FACTOR_R(ma);
    mb = SCALAR_FACTOR_R(mb);
    if(SCALAR_FACTOR_IS_I_NUMBER(wa))       /* Imaginary unit i*i --> -1 */
    {
      NODE_SF_KILL(wa);
      NODE_SF_KILL(wb);
      if(*pchange_sign)
        *pchange_sign = NO;               /* Convey change of sign to up */
      else
        *pchange_sign = YES;
      goto next_pair;
    }
    SCALAR_FACTOR_DEGREE(wa) += SCALAR_FACTOR_DEGREE(wb); /* Sum degrees */
    NODE_SF_KILL(wb);
    goto append_term;

    order_12:
    wa = ma;
    ma = SCALAR_FACTOR_R(ma);
    goto append_term;

    order_21:
    wa = mb;
    mb = SCALAR_FACTOR_R(mb);

    append_term:
    if(mc == NIL)
      mc = wa;
    else
      SCALAR_FACTOR_R(last) = wa;
    last = wa;
  }
  return mc;
}
/*=ScalarSumAddition=====================================================
 Sum of two scalar (polynomial) expressions
*/
uint ScalarSumAddition(uint a, uint b)
{
  uint sum = NIL, last, wa, wb, ma, mb;
  BIGINT na, nb, nc;
  while(YES)
  {
    next_pair:
    if(b == NIL)
    {                 /* List b is ended, append rest of a */
      if(sum == NIL)
        sum = a;
      else
        SCALAR_TERM_R(last) = a;
      break;
    }
    if(a == NIL)
    {                 /* List a is ended, append rest of b */
      if(sum == NIL)
        sum = b;
      else
        SCALAR_TERM_R(last) = b;
      break;
    }

    /* Compare scalar terms */

    ma = SCALAR_TERM_MONOMIAL(a);
    mb = SCALAR_TERM_MONOMIAL(b);
    while(YES)
    {
      if(ma == NIL)
      {
        if(mb != NIL)
          goto order_21;   /* a-monomial is contained in b-monomial */
        break;             /* a-monomial == b-monomial  */
      }                    /* (including both are NILs) */
      if(mb == NIL)
        goto order_12;     /* a-monomial contains b-monomial */
      if(SCALAR_FACTOR_WORD(ma) > SCALAR_FACTOR_WORD(mb))
        goto order_12;
      if(SCALAR_FACTOR_WORD(ma) < SCALAR_FACTOR_WORD(mb))
        goto order_21;
      ma = SCALAR_FACTOR_R(ma);  /* Skip equal factors in monomials */
      mb = SCALAR_FACTOR_R(mb);
    }

    /* Reduce like scalar terms */

    wa = a;
    wb = b;
    a = SCALAR_TERM_R(a);
    b = SCALAR_TERM_R(b);

        /* Sum integer coefficients */

    na = SCALAR_TERM_NUMERATOR(wa);
    nb = SCALAR_TERM_NUMERATOR(wb);
    INTEGER_HEAP_NEW(nc, 2+MAX(INTEGER_N_LIMBS(na),INTEGER_N_LIMBS(nb)));
    IntegerSum(nc, na, nb);
    INTEGER_KILL(na);
    INTEGER_KILL(nb);

        /* Kill head of wb and its monomial */

    mb = SCALAR_TERM_MONOMIAL(wb);
    NODE_ST_KILL(wb);
    while(mb != NIL)
    {
      ma = mb;
      mb = SCALAR_FACTOR_R(mb);
      NODE_SF_KILL(ma);
    }

        /* Nonzero collection */

    if(INTEGER_N_LIMBS(nc) != 0)
    {
      SCALAR_TERM_NUMERATOR(wa) = nc;
      goto append_term;
    }

        /* Zero collection: kill nc and wa */

    INTEGER_KILL(nc);
    ma = SCALAR_TERM_MONOMIAL(wa);
    NODE_ST_KILL(wa);
    while(ma != NIL)
    {
      mb = ma;
      ma = SCALAR_FACTOR_R(ma);
      NODE_SF_KILL(mb);
    }
    goto next_pair;

    order_12:
    wa = a;
    a = SCALAR_TERM_R(a);
    goto append_term;

    order_21:
    wa = b;
    b = SCALAR_TERM_R(b);

    append_term:
    if(sum == NIL)
      sum = wa;
    else
      SCALAR_TERM_R(last) = wa;
    last = wa;
  }
  return sum;
}
/*=ScalarSumCancellation===================
*/
void ScalarSumCancellation(uint *pnum, uint *pden)
{
  uint g;
IN_SCALAR_SUM_CANCELLATION /*------------*/
  if((g = PolyGCD(*pnum, *pden)) != NIL)
  {
    *pnum = PolyQuotient(*pnum, g);
    *pden = PolyQuotient(*pden, g);
    ScalarSumKill(g);
  }
OUT_SCALAR_SUM_CANCELLATION /*-----------*/
}
/*=ScalarSumMinus===========================
 Change sign of scalar sum in Parametric regime
*/
void ScalarSumMinus(uint a)
{
  while(a != NIL)
  {
    SCALAR_TERM_MINUS(a);
    a = SCALAR_TERM_R(a);
  }
}
/*=ScalarSumMultiplication===========================================
 Expanded product of two positive nonzero general scalar expressions,
 caller ensures A != NIL and B != NIL.
*/
uint ScalarSumMultiplication(uint a, uint b)
{
  uint s, aw, bcur, bw, ac, bc;
  s = NIL;
  while(a != NIL)
  {
    aw = a;
    a = SCALAR_TERM_R(a);
    bcur = b;
    while(bcur != NIL)
    {
      bw = bcur;
      bcur = SCALAR_TERM_R(bcur);
      if(a == NIL)
      {
        bc = bw;
        SCALAR_TERM_R(bc) = NIL;
      }
      else
        bc = ScalarTermCopy(bw);
      if(bcur == NIL)
      {
        ac = aw;
        SCALAR_TERM_R(ac) = NIL;
      }
      else
        ac = ScalarTermCopy(aw);
      ScalarTermMultiplication(ac, bc);
      s = ScalarSumAddition(s, ac);
    }
  }
  return s;
}
/*=ScalarTermMultiplication===============================================
 Product of two scalar terms on place of A, B deleted
*/
void ScalarTermMultiplication(uint a, uint b)
{
  BIGINT na, nb, nc;
  uint ma, mb, mc, last, aa;

  /* Multiply integer coefficients */

  na = SCALAR_TERM_NUMERATOR(a);
  nb = SCALAR_TERM_NUMERATOR(b);
  INTEGER_HEAP_NEW(nc, 1+INTEGER_N_LIMBS(na)+INTEGER_N_LIMBS(nb));
  IntegerProduct(nc, na, nb);
  INTEGER_KILL(na);
  INTEGER_KILL(nb);
  SCALAR_TERM_NUMERATOR(a) = nc;

  /* Multiply monomials */

  ma = SCALAR_TERM_MONOMIAL(a);
  mb = SCALAR_TERM_MONOMIAL(b);
  NODE_ST_KILL(b);
  mc = NIL;
  while(YES)
  {
    next_pair:
    if(mb == NIL)
    {                 /* List mb is ended, append rest of ma */
      if(mc == NIL)
        mc = ma;
      else
        SCALAR_FACTOR_R(last) = ma;
      break;
    }
    if(ma == NIL)
    {                 /* List ma is ended, append rest of mb */
      if(mc == NIL)
        mc = mb;
      else
        SCALAR_FACTOR_R(last) = mb;
      break;
    }

    /* Compare scalar factors */

    if(SCALAR_FACTOR_PARAMETER(ma) > SCALAR_FACTOR_PARAMETER(mb))
      goto order_12;
    if(SCALAR_FACTOR_PARAMETER(ma) < SCALAR_FACTOR_PARAMETER(mb))
      goto order_21;

    /* Reduce like factors */

    aa = ma;
    b = mb;                    /* uint-uint mixion */
    ma = SCALAR_FACTOR_R(ma);
    mb = SCALAR_FACTOR_R(mb);
    if(SCALAR_FACTOR_IS_I_NUMBER(aa)) /* Imaginary unit i*i --> -1 */
    {
      NODE_SF_KILL(aa);
      NODE_SF_KILL(b);
      INTEGER_MINUS(nc);
      goto next_pair;
    }
    SCALAR_FACTOR_DEGREE(aa) += SCALAR_FACTOR_DEGREE(b); /* Sum degrees */
    NODE_SF_KILL(b);
    goto append_term;

    order_12:
    aa = ma;
    ma = SCALAR_FACTOR_R(ma);
    goto append_term;

    order_21:
    aa = mb;
    mb = SCALAR_FACTOR_R(mb);

    append_term:
    if(mc == NIL)
      mc = aa;
    else
      SCALAR_FACTOR_R(last) = aa;
    last = aa;
  }
  SCALAR_TERM_MONOMIAL(a) = mc;
}

/*_6_4          Scalar polynomial algebraic functions===================*/

/*=ContentOfScalarSum=========================================
 Returns relative (or full if initial CONT == NIL) single term
 content of scalar sum. NIL corresponds to 1.
 CONT is destroyed, A remains.
*/
uint ContentOfScalarSum(uint cont, uint a)
{
  if(cont == NIL)
  {
    cont = ScalarTermCopy(a);
    INTEGER_SET_PLUS(SCALAR_TERM_NUMERATOR(cont));
    if((a = SCALAR_TERM_R(a)) == NIL)
      goto out_cont;
  }
  while((cont = PolyTermGCD(cont, a)) != NIL
        && (a = SCALAR_TERM_R(a)) != NIL)
    ;
  out_cont:
  return cont;
}
/*=InCoeffParamTable===============================================
 Set in CoeffParamTable parameters of scalar term `cont' killing it
*/
void InCoeffParamTable(uint cont)
{
  uint a = SCALAR_TERM_MONOMIAL(cont);
  INTEGER_KILL(SCALAR_TERM_NUMERATOR(cont));
  NODE_ST_KILL(cont);

  if(a != NIL)
  {
    if(SCALAR_FACTOR_IS_I_NUMBER(a))
    {
      NODE_SF_KILL(a);
      return;
    }
    if(CoeffParamTable == NULL)
      CoeffParamTable = (int*)NewArray(ParameterN, sizeof(int),
                                             E_A_COEFF_PARA_TABLE);
    do
    {
      cont = a;
      a = SCALAR_FACTOR_R(a);
      CoeffParamTable[SCALAR_FACTOR_PARAMETER(cont)] = YES;
      NODE_SF_KILL(cont);
    }while(a != NIL);
  }
}
/*=InCoeffSumTable========================================================
 Insert parametric content-free SUM in table or delete if already exists
*/
void InCoeffSumTable(uint sum)
{
  if(SCALAR_FACTOR_IS_I_NUMBER(SCALAR_TERM_MONOMIAL(sum)))
    ScalarSumKill(sum);                  /* Kill complex number a*i + b */
  else
  {
    int i;
    uint gcd, quocoe, quosum;
    for(i = 0; i < CoeffSumTableN; i++)
      if(PolynomialsAreEqual(sum, CoeffSumTable[i]))
      {
        ScalarSumKill(sum);
        return;
      }
      else
      {
        gcd = PolyGCD(sum, CoeffSumTable[i]);
        if(gcd != NIL)
        {
          quocoe = PolyQuotient(CoeffSumTable[i], gcd);
          quosum = PolyQuotient(sum, gcd);
          --CoeffSumTableN; /* Remove gap in ith position */
          while(i < CoeffSumTableN)
          {
            CoeffSumTable[i] = CoeffSumTable[i+1];
            ++i;
          }
          InCoeffSumTable(gcd);
          if(SCALAR_TERM_R(quocoe) == NIL) /* Either sum or 1 certainly */
          {
            INTEGER_KILL(SCALAR_TERM_NUMERATOR(quocoe));      /* Kill 1 */
            NODE_ST_KILL(quocoe);
          }
          else
            InCoeffSumTable(quocoe);
          if(SCALAR_TERM_R(quosum) == NIL) /* Either sum or 1 certainly */
          {
            INTEGER_KILL(SCALAR_TERM_NUMERATOR(quosum));      /* Kill 1 */
            NODE_ST_KILL(quosum);
          }
          else
            InCoeffSumTable(quosum);
          return;
        }
      }
    if(CoeffSumTableN >= CoeffSumTableSize)
      Error(E_COEFF_SUM_TABLE_SIZE);
    if(CoeffSumTable == NULL)
      CoeffSumTable = (uint *)NewArray(CoeffSumTableSize, sizeof(uint),
                                               E_A_COEFF_SUM_TABLE);
    CoeffSumTable[CoeffSumTableN++] = sum;
  }
}
/*=InCoeffTable===========================================
 Set in tables components of non-NIL parametric polynomial
*/
void InCoeffTable(uint coe)
{
  if(SCALAR_TERM_R(coe) != NIL)
  {
    uint cont;
    if((cont = ContentOfScalarSum(NIL, coe)) != NIL)
    {
      coe = PolyQuotient(coe, cont);
      InCoeffParamTable(cont);
    }
    InCoeffSumTable(coe);
  }
  else
    InCoeffParamTable(coe);
}
/*=PolyCoeffAtMainParameter===========================================
 Polynomial coefficient with normalized sign at the current degree of
 the main parameter. *PA initially points to the start of list, at the
 end of work it points to tail of list. NIL corresponds to 1.
 Initial polynomial remains.
*/
uint PolyCoeffAtMainParameter(uint *pa, int mp)
{
  uint a;
  int isnegative = INTEGER_IS_NEGATIVE(SCALAR_TERM_NUMERATOR(*pa));
  if(POLY_MAIN_PARAMETER(*pa) < mp)
  {                                                    /* Free term */
    a = ScalarSumCopy(*pa);
    if(isnegative)                                /* Normalize sign */
    {
      *pa = a;
      do
        SCALAR_TERM_MINUS(*pa);
      while((*pa = SCALAR_TERM_R(*pa)) != NIL);
    }
    else
      *pa = NIL;
  }
  else
  {
    uint b;
    uint mf;
    int mppow = SCALAR_TERM_MAIN_DEGREE(*pa);
    a = b = ScalarTermCopy(*pa);
    if(isnegative)
      SCALAR_TERM_MINUS(a);
    mf = SCALAR_TERM_MONOMIAL(a);      /* Strike out main parameter */
    SCALAR_TERM_MONOMIAL(a) = SCALAR_FACTOR_R(mf);
    NODE_SF_KILL(mf);
    while((*pa = SCALAR_TERM_R(*pa)) != NIL
          && POLY_MAIN_PARAMETER(*pa) == mp
          && SCALAR_TERM_MAIN_DEGREE(*pa) == mppow)
    {
      SCALAR_TERM_R(b) = ScalarTermCopy(*pa);
      b = SCALAR_TERM_R(b);
      if(isnegative)
      {
        SCALAR_TERM_MINUS(b);
      }
      mf = SCALAR_TERM_MONOMIAL(b);    /* Strike out main parameter */
      SCALAR_TERM_MONOMIAL(b) = SCALAR_FACTOR_R(mf);
      NODE_SF_KILL(mf);
    }
  }
  if(SCALAR_SUM_IS_UNIT(a))
  {
    INTEGER_KILL(SCALAR_TERM_NUMERATOR(a));
    NODE_ST_KILL(a);
    a = NIL;
  }
  return a;
}
/*=PolyContent=============================================
 Polynomial content of polynomial w.r.t. main parameter MP.
 A remains unchanged.
*/
uint PolyContent(uint a, int mp)
{
  uint b;
IN_POLY_CONTENT /*---------------------------------------*/
  if((b = PolyCoeffAtMainParameter(&a, mp)) != NIL)
  {
    uint c, d;
    while(a != NIL)
    {
      if((c = PolyCoeffAtMainParameter(&a, mp)) == NIL)
      {
        ScalarSumKill(b);
        b = NIL;
        break;
      }
      d = b;
      if((b = PolyGCD(b, c)) == NIL)
      {
        ScalarSumKill(d);
        ScalarSumKill(c);
        break;
      }
    }
  }
OUT_POLY_CONTENT /*--------------------------------------*/
  return b;
}
/*=PolyGCD==============================================================
 Returns Greatest Common Divisor of two multivariate polynomials in
 the form GCD(PP(A), PP(B)) * GCD(CONT(A), CONT(B)).
 A, B unchanged.
 Returned NIL means trivial GCD = 1
*/
uint PolyGCD(uint a, uint b)
{
  uint c;
IN_POLY_GCD /*--------------------------------------------------------*/
  if(SCALAR_TERM_R(a) == NIL || SCALAR_TERM_R(b) == NIL)
  {               /* At least one of the polynomials is not a sum */
    if(SCALAR_TERM_R(a) != NIL)
    {                                    /* Set A to be a single term */
      c = a;
      a = b;
      b = c;
    }
    c = ScalarSumCopy(a);
    INTEGER_SET_PLUS(SCALAR_TERM_NUMERATOR(c));
    b = ContentOfScalarSum(c, b);
  }
  else                                 /* Both are polynomials really */
  {
    uint conta, contb;
    int mp, mpb;                                     /* Main parameters */
    mp = SCALAR_FACTOR_PARAMETER(SCALAR_TERM_MONOMIAL(a));
    mpb = SCALAR_FACTOR_PARAMETER(SCALAR_TERM_MONOMIAL(b));
    if(mpb > mp ||
       (mpb == mp &&
        SCALAR_TERM_MAIN_DEGREE(b) > SCALAR_TERM_MAIN_DEGREE(a)))
    {                           /* Parameters go in DECREASING order! */
      c = a;                                      /* Swap polynomials */
      a = b;
      b = c;
      mp = mpb;
    }
    a = ScalarSumCopy(a);
    b = ScalarSumCopy(b);
    contb = PolyContent(b, mp);
    if((conta = PolyContent(a, mp)) != NIL)
    {                     /* Make primitive parts and GCD of contents */
      if(contb != NIL)
      {
        c = PolyGCD(conta, contb);
        b = PolyQuotient(b, contb);                 /* Primitive part */
        ScalarSumKill(contb);
      }
      else
        c = NIL;
      a = PolyQuotient(a, conta);                   /* Primitive part */
      ScalarSumKill(conta);
    }
    else
    {
      if(contb != NIL)
      {
        b = PolyQuotient(b, contb);                 /* Primitive part */
        ScalarSumKill(contb);
      }
      c = NIL;
    }
    while((conta = PolyPseudoRemainder(a, ScalarSumCopy(b), mp)) != NIL)
    {
      if(SCALAR_TERM_MONOMIAL(conta) == NIL || /* Pure number */
         SCALAR_FACTOR_PARAMETER(SCALAR_TERM_MONOMIAL(conta)) != mp)
      {                 /* Zero degree with respect to main parameter */
        ScalarSumKill(b);
        ScalarSumKill(conta);
        b = c;          /* char is content ?? */
        goto out;
      }
      a = b;
      if((contb = ContentOfScalarSum(NIL, conta)) != NIL)
      {                                               /* Term content */
        conta = PolyQuotient(conta, contb);
            ScalarSumKill(contb);
      }
      if((contb = PolyContent(conta, mp)) != NIL)
      {                                         /* Polynomial content */
        b = PolyQuotient(conta, contb);             /* Primitive part */
        ScalarSumKill(contb);
      }
      else
            b = conta;
    }
    if(c != NIL)
      b = ScalarSumMultiplication(b, c);
    if(INTEGER_IS_NEGATIVE(SCALAR_TERM_NUMERATOR(b)))
    {                                      /* Standardize sign of GCD */
      c = b;
      do
        SCALAR_TERM_MINUS(c);
      while((c = SCALAR_TERM_R(c)) != NIL);
    }
    if(SCALAR_SUM_IS_UNIT(b))
    {                                      /* Standardize trivial GCD */
      INTEGER_KILL(SCALAR_TERM_NUMERATOR(b));
      NODE_ST_KILL(b);
      b = NIL;
    }
  }
  out:
OUT_POLY_GCD /*-------------------------------------------------------*/
  return b;
}
/*=PolyMainParameterTerm================================================
 Take sublist of terms, containing main parameter MP of given degree
 MPDEG, *PA points to the next term after end of A.
 This function is applied in succession starting from top degree.
 Initial expression *PA is destructed.
*/
uint PolyMainParameterTerm(uint *pa, int mp, int mpdeg)
{
  uint a;
  if(mpdeg)
  {
    unsigned short w; /* Word combining degree and index of parameter (MPDEG,MP) */
    w = mpdeg;
#if defined(SPP_2000)
    *((byte*)&w) = mp;
#else
    *((byte*)&w+1) = mp;
#endif
    if(SCALAR_TERM_MONOMIAL(*pa) != NIL &&
       SCALAR_TERM_MAIN_PARAMETER_WORD(*pa) == w)
    {                                         /* There is such degree */
      uint b;
      a = *pa;
      while(YES)
      {
        b = *pa;                            /* Remember previous term */
        *pa = SCALAR_TERM_R(*pa);
        if(*pa == NIL)
          break;                   /* Whole expression is homogeneous */
        if(SCALAR_TERM_MONOMIAL(*pa) == NIL ||
           SCALAR_TERM_MAIN_PARAMETER_WORD(*pa) != w)
        {                         /* End of homogeneous part is found */
          SCALAR_TERM_R(b) = NIL;               /* Set end of sublist */
              break;
        }
      }
    }
    else                                            /* No such degree */
      a = NIL;
  }
  else
  {                                                      /* Free term */
    a = *pa;
    *pa = NIL;
  }
  return a;
}
/*=PolynomialsAreEqual====================================
*/
int PolynomialsAreEqual(uint a, uint b)
{
  uint ma, mb;
  BIGINT na, nb;
  while(YES)
  {
    /* Compare monomials */

    ma = SCALAR_TERM_MONOMIAL(a);
    mb = SCALAR_TERM_MONOMIAL(b);
    while(YES)
    {
      if(ma == NIL)
        if(mb == NIL)
          break;
            else
          goto no;
      else if(mb == NIL)
            goto no;
      if(SCALAR_FACTOR_WORD(ma) != SCALAR_FACTOR_WORD(mb))
            goto no;
      ma = SCALAR_FACTOR_R(ma);
      mb = SCALAR_FACTOR_R(mb);
    }
    /* Compare numerators */

    na = SCALAR_TERM_NUMERATOR(a);
    nb = SCALAR_TERM_NUMERATOR(b);
    if(na[0] != nb[0])
      goto no;
    ma = INTEGER_N_LIMBS(na);
    do
      if(na[ma] != nb[ma])
        goto no;
    while(ma-- > 0);

    a = SCALAR_TERM_R(a);
    b = SCALAR_TERM_R(b);
    if(a == NIL)
      if(b == NIL)
        return YES;
      else
        goto no;
        else if(b == NIL)
        goto no;
  }
  no:
  return NO;
}
/*=PolyPseudoRemainder==================================================
 Returns pseudo-remainder of two polynomials. MP is main parameter.
 main_degree(A) >= main_degree(B). A, B destructed.
*/
uint PolyPseudoRemainder(uint a, uint b, int mp)
{
IN_POLY_PSEUDO_REMAINDER /*-------------------------------------------*/
  if(SCALAR_TERM_MAIN_PARAMETER(b) != mp)
  {                         /* B doesn't contain MP => return 0 (NIL) */
    ScalarSumKill(a);
    ScalarSumKill(b);
    a = NIL;
  }
  else
  {
    uint *u, *v, vn, c, w;
    int m, n, j, k;
    m = SCALAR_TERM_MAIN_DEGREE(a);
    n = SCALAR_TERM_MAIN_DEGREE(b);
    POLY_ARRAY_STACK_NEW(u, m+1);
    POLY_ARRAY_STACK_NEW(v, n);
    for(j = m; j >= 0; j--)                     /*          j         */
      u[j] = PolyMainParameterTerm(&a, mp, j);  /* u[j] = mp  u  etc. */
    vn = PolyMainParameterTerm(&b, mp, n);      /*             j      */
    for(j = n - 1; j >= 0; j--)
      v[j] = PolyMainParameterTerm(&b, mp, j);
    for(k = m - n; k >= 0; k--)
    {
      j = n + k - 1;
      while(j >= k)                             /*    n+j             */
      {                                         /*  mp   (u  = v  u ) */
        if(u[j] != NIL)                         /*         j    n  j  */
          u[j] = ScalarSumMultiplication(ScalarSumCopy(vn), u[j]);
        if(u[n+k] != NIL)
          if(v[j-k] != NIL)
          {
            a = ScalarSumMultiplication(ScalarSumCopy(v[j-k]),
                                        ScalarSumCopy(u[n+k]));
            b = a;
            do
              SCALAR_TERM_MINUS(b);
            while((b = SCALAR_TERM_R(b)) != NIL);     /*     n+j          */
              u[j] = ScalarSumAddition(u[j], a);      /* - mp   v    u    */
          }                                           /*         j-k  n+k */
        if(u[j] != NIL)
        {
          c = u[j];                  /* Drop degree of main parameter */
          do
            if((SCALAR_TERM_MAIN_DEGREE(c) -= n) == 0)
            {                            /* Strike out main parameter */
              w = SCALAR_TERM_MONOMIAL(c);
              SCALAR_TERM_MONOMIAL(c) = SCALAR_FACTOR_R(w);
              NODE_SF_KILL(w);
            }
          while((c = SCALAR_TERM_R(c)) != NIL);
        }
        j--;
      }
      if(u[n+k] != NIL)
        ScalarSumKill(u[n+k]);
      while(j >= 0)
      {
        if(u[j] != NIL)
        {
          u[j] = ScalarSumMultiplication(ScalarSumCopy(vn), u[j]);
          c = u[j];                  /* Drop degree of main parameter */
          do
            if((SCALAR_TERM_MAIN_DEGREE(c) -= n) == 0)
            {                            /* Strike out main parameter */
              w = SCALAR_TERM_MONOMIAL(c);
              SCALAR_TERM_MONOMIAL(c) = SCALAR_FACTOR_R(w);
              NODE_SF_KILL(w);
            }
          while((c = SCALAR_TERM_R(c)) != NIL);
        }
        j--;
      }
    }
    ScalarSumKill(vn);
    for(j = n - 1; j >= 0; j--)
      if(v[j] != NIL)
        ScalarSumKill(v[j]);
    j = n - 1;
    while(j >= 0 && u[j] == NIL)
      j--;                          /* Search first nonzero term u[j] */
    if(j >= 0)
    {                   /* Concatenate pseudoremainder from array u[] */
      a = u[j--];
      b = a;
      while(SCALAR_TERM_R(b) != NIL)
        b = SCALAR_TERM_R(b);
      while(j >= 0)
      {
        if(u[j] != NIL)
            {
              SCALAR_TERM_R(b) = u[j];
          while(SCALAR_TERM_R(b) != NIL)
            b = SCALAR_TERM_R(b);
        }
        j--;
      }
    }
    else                                        /* All u[j] are zeros */
      a = NIL;
  }
OUT_POLY_PSEUDO_REMAINDER /*------------------------------------------*/
  return a;
}
/*=PolyTermGCD==========================================================
 GCD of two single (non-NIL) terms, A is destroyed, B remains,
 caller ensures A is positive, returned NIL corresponds to 1
*/
uint PolyTermGCD(uint a, uint b)
{
  BIGINT na, naa, nb, nbb;
  int i;
  uint ma, maa;

  /* Do integer coefficients */

  na = SCALAR_TERM_NUMERATOR(a);
  nb = SCALAR_TERM_NUMERATOR(b);
  naa = na;                               /* Anyway it will be killed */
  INTEGER_STACK_COPY(nbb, nb, i);
  if((naa = IntegerGCD(naa, nbb)) != NULL)   /* naa = GCD(na, nb) > 1 */
  {
    INTEGER_HEAP_COPY(nb, naa, i);
    SCALAR_TERM_NUMERATOR(a) = nb;
  }
  INTEGER_KILL(na);

  /* Do parametric monomials (parameters go in DECREASING order) */

  maa = NIL;
  if((ma = SCALAR_TERM_MONOMIAL(a)) != NIL)
  {
    uint maw, mb, mal;
    mb = SCALAR_TERM_MONOMIAL(b);
    while(YES)
    {
      if(mb == NIL)
      {               /* MB is ended - kill tail of MA and break loop */
        while(ma != NIL)
        {
          maw = ma;
          ma = SCALAR_FACTOR_R(ma);
          NODE_SF_KILL(maw);
        }
        if(maa != NIL)
          SCALAR_FACTOR_R(mal) = NIL;
        break;
      }
      if(SCALAR_FACTOR_PARAMETER(ma) > SCALAR_FACTOR_PARAMETER(mb))
      {
        maw = ma;                        /* No match for MA - kill it */
        ma = SCALAR_FACTOR_R(ma);
        NODE_SF_KILL(maw);
        if(ma == NIL)
        {
          if(maa != NIL)
            SCALAR_FACTOR_R(mal) = NIL;
          break;
        }
      }
      else if(SCALAR_FACTOR_PARAMETER(ma) < SCALAR_FACTOR_PARAMETER(mb))
        mb = SCALAR_FACTOR_R(mb);       /* No match for MB - shift it */
      else
      {                                 /* Match - set minimum degree */
        if(SCALAR_FACTOR_DEGREE(mb) < SCALAR_FACTOR_DEGREE(ma))
          SCALAR_FACTOR_DEGREE(ma) = SCALAR_FACTOR_DEGREE(mb);
        if(maa == NIL)
          maa = ma;
        else                                    /* Append to last MAA */
          SCALAR_FACTOR_R(mal) = ma;
        mal = ma;
        ma = SCALAR_FACTOR_R(ma);
        if(ma == NIL)
          break;
      }
    }
  }
  SCALAR_TERM_MONOMIAL(a) = maa;          /* Set constructed monomial */
  if(naa == NULL)
  {
    if(maa == NIL)
    {                                                  /* Trivial GCD */
      NODE_ST_KILL(a);
      a = NIL;
    }
    else                         /* Make standard integer coefficient */
    {
      INTEGER_HEAP_NEW(na, 2);
      na[0] = na[1] = 1;
      SCALAR_TERM_NUMERATOR(a) = na;
    }
  }
  return a;
}
/*=PolyTermQuotient===================================================
 Exact division of term A by term B: A = char*B on place of A, B remains.
 Parameters go in decreasing order.
*/
void PolyTermQuotient(uint a, uint b)
{
  BIGINT na, nb, naa, nbb, nc;
  int i;
  uint mb;

  /* Divide integer numerator */

  na = SCALAR_TERM_NUMERATOR(a);
  nb = SCALAR_TERM_NUMERATOR(b);
  INTEGER_STACK_COPY(nbb, nb, i);
  INTEGER_HEAP_NEW(nc, 2+INTEGER_N_LIMBS(na)-INTEGER_N_LIMBS(nbb));
  INTEGER_STACK_COPY_1(naa, na, i);
  INTEGER_KILL(na);
  IntegerQuotient(nc, naa, nbb);
  SCALAR_TERM_NUMERATOR(a) = nc;

  /* Divide parametric monomial */

  if((mb = SCALAR_TERM_MONOMIAL(b)) != NIL)
  {
    uint ma, maa, mae, maw;
    ma = SCALAR_TERM_MONOMIAL(a);
    maa = NIL;
    do
    {
      while(SCALAR_FACTOR_PARAMETER(ma) > SCALAR_FACTOR_PARAMETER(mb))
      {
        if(maa == NIL)
          maa = ma;
        else
          SCALAR_FACTOR_R(mae) = ma;
        mae = ma;
        ma = SCALAR_FACTOR_R(ma);
      }
      if((SCALAR_FACTOR_DEGREE(ma) -= SCALAR_FACTOR_DEGREE(mb)) != 0)
      {
        if(maa == NIL)
          maa = ma;
        else
          SCALAR_FACTOR_R(mae) = ma;
        mae = ma;
        ma = SCALAR_FACTOR_R(ma);
      }
      else
      {
        maw = ma;
        ma = SCALAR_FACTOR_R(ma);
        NODE_SF_KILL(maw);
      }
    }while((mb = SCALAR_FACTOR_R(mb)) != NIL);
    if(maa == NIL)                      /* Append tail of numerator */
      maa = ma;
    else
      SCALAR_FACTOR_R(mae) = ma;
    SCALAR_TERM_MONOMIAL(a) = maa;
  }
}
/*=PolyQuotient====================================================
 Exact division of polynomial A by polynomial B: A = char*B, return char.
 Caller ensures A, B != NIL, B is positive.
 A is destructed, B remains unchanged.
*/
uint PolyQuotient(uint a, uint b)
{
  uint c;
IN_POLY_QUOTIENT /*----------------------------------------------*/
  if(SCALAR_TERM_R(b) == NIL)       /* Division by single term B */
  {
    c = a;
    if(SCALAR_SUM_IS_NOT_UNIT(b))             /* Nontrivial term */
      do
        PolyTermQuotient(a, b);
      while((a = SCALAR_TERM_R(a)) != NIL);
  }
  else                               /* Division by polynomial B */
  {
    uint aw, bw, cw;
    BIGINT n;
    bw = SCALAR_TERM_R(b);
    c = NIL;
    do
    {
      aw = a;
      a = SCALAR_TERM_R(a);
      SCALAR_TERM_R(aw) = NIL;
      PolyTermQuotient(aw, b);
      cw = ScalarTermCopy(aw);
      n = SCALAR_TERM_NUMERATOR(aw);
      INTEGER_MINUS(n);
      aw = ScalarSumMultiplication(aw, ScalarSumCopy(bw));
      a = ScalarSumAddition(a, aw);            /* Remainder of A */
      c = ScalarSumAddition(c, cw);                /* Quotient char */
    }while(a != NIL);
  }
OUT_POLY_QUOTIENT /*---------------------------------------------*/
  return c;
}

/*_6_5          Big number functions====================================*/

/*=BigNMinusBigN================================================
 Subtract two big numbers on the place of first one: `a' -= `b',
 caller provides `a' > `b', returns new size of `a'
*/
int BigNMinusBigN(BIGINT a, int na, BIGINT b, int nb)
{
  uint lw;
  LIMB k = 1;
  int i = 0;
  while(i < nb)         /* Common part */
  {
    lw = MAX_LIMB + (uint)k + (uint)a[i] - (uint)b[i];
    k = (lw > MAX_LIMB);
    a[i++] = (LIMB)lw;
  }
  while(i < na)
  {
    lw = MAX_LIMB + (uint)k + (uint)a[i];
    k = (lw > MAX_LIMB);
    a[i++] = (LIMB)lw;
  }
  while(--i >= 0)
    if(a[i] != 0)
      break;
  return (++i);
}
/*=BigNShiftLeft==============================================
 Add on spot 0 <= `cnt' < BITS_PER_LIMB lowest zero bits to
 `bign' of size `n', return the bits shifted out from the most
 significant LIMB digit
*/
LIMB BigNShiftLeft(BIGINT bign, int n, int cnt)
{
  if(cnt)
  {
    int cocnt = BITS_PER_LIMB - cnt;
    LIMB low_limb,
         high_limb = bign[--n],
         pushed_out = high_limb >> cocnt;
    while(--n >= 0)
    {
      low_limb = bign[n];
      bign[n+1] = (high_limb << cnt) | (low_limb >> cocnt);
      high_limb = low_limb;
    }
    bign[n+1] = high_limb << cnt;
    return pushed_out;
  }
  return 0;
}
/*=BigNShiftRight==========================================
 Remove on spot 0 <= `cnt' < BITS_PER_LIMB lowest bits from
 `bign' of size `n', return size of the result
*/
int BigNShiftRight(BIGINT bign, int n, int cnt)
{
  if(cnt)
  {
    BIGINT bigni;
    int high_limb, low_limb, i, cocnt = BITS_PER_LIMB - cnt;
    low_limb = *bign;
    bigni = bign;
    bigni++;
    for(i = n-1; i > 0; i--)
    {
      high_limb = *(bigni++);
      *(bign++) = (low_limb >> cnt) | (high_limb << cocnt);
      low_limb = high_limb;
    }
    low_limb >>= cnt;
    if(low_limb != 0)
      *bign = low_limb;
    else
      n--;
#if 0
    LIMB high_limb, low_limb;
    int i, cocnt = BITS_PER_LIMB - cnt;
    low_limb = bign[0];
    for(i = 1; i < n; i++)
    {
      high_limb = bign[i];
      bign[i-1] = (low_limb >> cnt) | (high_limb << cocnt);
      low_limb = high_limb;
    }
    low_limb >>= cnt;
    if(low_limb != 0)
      bign[i-1] = low_limb;
    else
      n--;
#endif
  }
  return n;
}
/*=CountLeadingZeroBitsInLimb=========================
 Count number of leading zero bits in LIMB word
*/
int CountLeadingZeroBitsInLimb(LIMB w)
{
  if(w >= 0x100u) /* [0,  7] */
    if(w >= 0x1000u)    /* [0, 3] */
      if(w >= 0x4000u)          /* [0, 1] */
        if(w >= 0x8000u)
          return 0;
        else
          return 1;
      else                      /* [2, 3] */
        if(w >= 0x2000u)
          return 2;
        else
          return 3;
    else                /* [4, 7] */
      if(w >= 0x400u)           /* [4, 5] */
        if(w >= 0x800u)
          return 4;
        else
          return 5;
      else                      /* [6, 7] */
        if(w >= 0x200u)
          return 6;
        else
          return 7;
  else            /* [8, 16] */
    if(w >= 0x10u)      /* [ 8, 11] */
      if(w >= 0x40u)            /* [ 8,  9] */
        if(w >= 0x80u)
          return 8;
        else
          return 9;
      else                      /* [10, 11] */
        if(w >= 0x20u)
          return 10;
        else
          return 11;
    else                /* [12, 16] */
      if(w >= 0x4u)             /* [12, 13] */
        if(w >= 0x8u)
          return 12;
        else
          return 13;
      else                      /* [14, 16] */
        if(w >= 0x2u)
          return 14;
        else                            /* [15, 16] */
          if(w)
            return 15;
          else
            return 16;
}
/*=IntegerCancellation========================================
 Results are placed in `num' and `den' arrays
*/
void IntegerCancellation(BIGINT num, BIGINT den)
{
  BIGINT n, d, g;
  int i;
IN_INTEGER_CANCELLATION  /*----------------------*/
  INTEGER_STACK_COPY_1(n, num, i);
  INTEGER_STACK_COPY_1(d, den, i);
  if((g = IntegerGCD(num, den)) != NULL)
  {
    BIGINT gg;
    INTEGER_STACK_COPY(gg, g, i);  /* `g' in `num' */

    /* Cancel `den' */

    IntegerQuotient(den, d, g);

    /* Cancel `num' */

    IntegerQuotient(num, n, gg);
  }
  else    /* Lozh' vzad */
  {
    i = INTEGER_N_LIMBS(n);
    do
      num[i] = n[i];
    while(i--);
    i = d[0];  /* Denominators and GCDs are positive always */
    do
      den[i] = d[i];
    while(i--);
  }
OUT_INTEGER_CANCELLATION /*----------------------*/
}
/*=IntegerGCD=========================================================
 Binary algorithm is used,
 Returns pointer to array of greatest common divisor
 or NULL interpreted as 1 by caller,
 the function spoils both arrays `u' and `v',
 result is placed in `u' array
*/
BIGINT IntegerGCD(BIGINT u, BIGINT v)
{
  int i, nu, nv, bcnt, w_bcnt;
  BIGINT u0;
  LIMB carry_digit;
IN_INTEGER_GCD  /*--------------------------------------------------*/
  nu = INTEGER_N_LIMBS(u);
  nv = INTEGER_N_LIMBS(v);
  u0 = ++u;       /* Skip size information limbs and memorize begin */
  ++v;
  i = 0;                             /* Shift down uint to make it odd */
  while(*u == 0)
  {                                              /* Skip zero limbs */
    ++i;
    ++u;
  }
  COUNT_LEADING_ZERO_BITS_IN_LIMB(bcnt, *u & -*u);
  bcnt = BITS_PER_LIMB - 1 - bcnt;
  nu = BigNShiftRight(u, nu - i, bcnt);
  bcnt += i * BITS_PER_LIMB;
  w_bcnt = bcnt;
  i = 0;                             /* Shift down void to make it odd */
  while(*v == 0)
  {                                              /* Skip zero limbs */
    ++i;
    ++v;
  }
  COUNT_LEADING_ZERO_BITS_IN_LIMB(bcnt, *v & -*v);
  bcnt = BITS_PER_LIMB - 1 - bcnt;
  nv = BigNShiftRight(v, nv - i, bcnt);
  bcnt += i * BITS_PER_LIMB;
  if(bcnt < w_bcnt)
    w_bcnt = bcnt;                  /* Number of common 2 factors.  */
  while(YES)
  {
    if(nu > nv)
      goto u_greater_v;
    if(nu < nv)
      goto v_greater_u;
    i = nu;
    while(--i >= 0)
    {
      if(u[i] > v[i])
        goto u_greater_v;
      if(v[i] > u[i])
        goto v_greater_u;
    }
    break;   /* If uint and void have become equal, we have found the GCD */
    u_greater_v:  /* Replace uint by (uint - void) >> cnt making uint odd again */
    nu = BigNMinusBigN(u, nu, v, nv);
    while(*u == 0)
    {
      --nu;
      ++u;
    }
    COUNT_LEADING_ZERO_BITS_IN_LIMB(bcnt, *u & -*u);
    bcnt = BITS_PER_LIMB - 1 - bcnt;
    nu = BigNShiftRight(u, nu, bcnt);
    continue;
    v_greater_u:  /* Replace void by (void - uint) >> cnt making void odd again */
    nv = BigNMinusBigN(v, nv, u, nu);
    while(*v == 0)
    {
      --nv;
      ++v;
    }
    COUNT_LEADING_ZERO_BITS_IN_LIMB(bcnt, *v & -*v);
    bcnt = BITS_PER_LIMB - 1 - bcnt;
    nv = BigNShiftRight(v, nv, bcnt);
  }
  /* GCD(U_IN, V_IN) now is uint * 2**W_BCNT.  */
  carry_digit = BigNShiftLeft(u, nu, w_bcnt % BITS_PER_LIMB);
  i = w_bcnt / BITS_PER_LIMB;
  u -= i;
  i += nu;
  if(carry_digit != 0)
  {
    if(u > u0)
    {
      u0 = u--;   /* Shift to left by 1 limb to make room for carry */
      for(nu = 0; nu < i; nu++)
        u[nu] = u0[nu];
    }
    u[i++] = carry_digit;
  }
  if(i == 1 && u[0] == 1)
    return NULL;
  --u;
  u[0] = i;                                    /* PLUS == 0 assumed */
OUT_INTEGER_GCD /*--------------------------------------------------*/
  return u;
}
/*=IntegerProduct======================================
 Traditional multiplication of two signed non-zero
 big numbers uint and void, result in W, W[] != uint[] or void[]
*/
void IntegerProduct(BIGINT w, BIGINT u, BIGINT v)
{
  LIMB carry;
  uint luw;
  BIGINT w0;
  int set_minus, i, j, n, m;
IN_INTEGER_PRODUCT  /*-------------------------------*/
  n = INTEGER_N_LIMBS(u);
  m = INTEGER_N_LIMBS(v);
  set_minus = (INTEGER_SIGN(u) != INTEGER_SIGN(v));
  u++;
  v++;
  w0 = w;
  w++;
  i = j = 0;
  do
    w[i] = 0;
  while(++i < n);
  do
  {
    i = carry = 0;
    do
    {
      luw = (uint)u[i]*(uint)v[j] + (uint)w[i+j] + (uint)carry;
      w[i+j] = (LIMB)luw;
      carry = (LIMB)(luw/BASE_LIMB);
    }while(++i < n);
    w[j+n] = carry;
  }while(++j < m);
  n += m - (carry == 0);
#if defined(INTEGER_MAX_SIZE)
  if(n > IntegerMaxSize)
    IntegerMaxSize = n;
#endif
  w0[0] = n;
  if(set_minus)
    INTEGER_SET_MINUS(w0);
OUT_INTEGER_PRODUCT /*-------------------------------*/
}
/*=IntegerQuotient=====================================================
 Exact division of big numbers:
 Quotient in char[*PM - N + 1 or *PM - N stored in *PM] = A[M] / B[N],
 char != A, char != B.
 Function is spoiling input A and B.
 Array for A should have 1 additional LIMB at the top
 for increasing A at normalizing B.
*/
void IntegerQuotient(BIGINT c, BIGINT a, BIGINT b)
{
  uint lw;
  LIMB q;
  int i, n, set_minus = (INTEGER_SIGN(a) != INTEGER_SIGN(b));
#if defined(D_CHECK_EXACTNESS_OF_DIVISION)
  int nr;
#endif
  BIGINT pm = c;
IN_INTEGER_QUOTIENT  /*----------------------------------------------*/
  *pm = INTEGER_N_LIMBS(a);
  n = INTEGER_N_LIMBS(b);
  a++;
  b++;
  c++;
  if(n == 1)
  {                                      /* Division by short number */
    i = *pm - 1;
    q = a[i] % *b;                                /* Carried residue */
    if((c[i] = a[i] / *b) == 0)
      --*pm;
    while(i)
    {
      lw = (uint)(q)*BASE_LIMB + (uint)a[--i];
      q = (LIMB)(lw % *b);
      c[i] = (LIMB)(lw / *b);
    }
#if defined(D_CHECK_EXACTNESS_OF_DIVISION)
    nr = (q != 0);
#endif
  }
  else
  {                                 /* Division by big number: n > 1 */
    BIGINT aw, bq;
    int k, j, shift, n1 = n + 1;
    LIMB carry, aj, aj1, aj2, bn1, bn2;
    INTEGER_STACK_NEW(bq, n1);                        /* For B[] * Q */
    COUNT_LEADING_ZERO_BITS_IN_LIMB(shift, b[n-1]);
    if(shift)             /* Normalize to make b[n-1] >= BASE_LIMB/2 */
    {
      a[*pm] = BigNShiftLeft(a, *pm, shift);
      BigNShiftLeft(b, n, shift);
    }
    else
      a[*pm] = 0;
    bn1 = b[n-1];
    bn2 = b[n-2];
    j = *pm;
    k = *pm - n;                     /* Top digit char[M-N] may be zero */
    *pm = k + 1;                /* Number of char[K] getting iterations */
    aw = a + k;  /* Start of current subarray of A of the length N+1 */
    do
    {
      aj = a[j];
      aj1 = a[j-1];
      aj2 = a[j-2];
      lw = (uint)aj * BASE_LIMB + (uint)aj1;
      q = (aj == bn1) ? MAX_LIMB : (LIMB)(lw/bn1);
      lw -= (uint)q*(uint)bn1;
      if(lw < BASE_LIMB && (uint)bn2*(uint)q > lw*BASE_LIMB + (uint)aj2)
      {                      /* Knuth's criterion shows Q is too big */
        q--;
        lw += (uint)bn1;
        if(lw < BASE_LIMB && (uint)bn2*(uint)q > lw*BASE_LIMB + (uint)aj2)
          q--;  /* Q was still too big */
      }
      if(q)
      {                                     /* Multiply and subtract */
        i = carry = 0;                /* Make copy of product B by Q */
        do
        {
          lw = (uint)q * (uint)b[i] + (uint)carry;
          carry = (LIMB)(lw/BASE_LIMB);
          bq[i] = (LIMB)lw;
        }while(++i < n);
        bq[i] = carry;               /* BQ[] padded by zero if needs */
        i = n;
        do
          if(aw[i] != bq[i])
          {                                      /* AW[] - BQ[] != 0 */
            if(aw[i] < bq[i])
            {                             /* AW[] - BQ[] is negative */
              q--;                          /* Additional correction */
              BigNMinusBigN(bq, n1, b, n);
            }
            break;
          }
        while(--i >= 0);
#if defined(D_CHECK_EXACTNESS_OF_DIVISION)
        nr =
#endif
        BigNMinusBigN(aw, n1, bq, n1);        /* AW[] - BQ[] on spot */
      }
      c[k--] = q;                    /* Set current LIMB of quotient */
      --aw;                            /* Shift subarray of A digits */
    }while(--j >= n);
    if(c[*pm-1] == 0)
      --*pm;           /* Real quotient size is M - N, not M - N + 1 */
#if defined(D_CHECK_EXACTNESS_OF_DIVISION)
    if(nr && shift)
      nr = BigNShiftRight(a, nr, shift);      /* Normalize remainder */
#endif
  }
  if(set_minus)
    INTEGER_SET_MINUS(pm);
OUT_INTEGER_QUOTIENT /*----------------------------------------------*/
#if defined(D_CHECK_EXACTNESS_OF_DIVISION)
  if(nr)
  {
    PutFormattedU("\n***Division violation at Debug==%u\n", Debug);
    EXIT;
  }
#endif
}
/*=IntegerSum=========================================
 Sum of two signed big numbers A and B, result in char
*/
void IntegerSum(BIGINT c, BIGINT a, BIGINT b)
{
  int set_minus, i, na, nb;
  uint lw;
  LIMB carry;
IN_INTEGER_SUM /*-----------------------------------*/
  if(INTEGER_N_LIMBS(a) < INTEGER_N_LIMBS(b))
  {
    BIGINT w = a;   /* Swap input numbers if necessary */
    a = b;
    b = w;
  }
  na = INTEGER_N_LIMBS(a);
  nb = INTEGER_N_LIMBS(b);
  if(INTEGER_SIGN(a) == INTEGER_SIGN(b))
  {                     /* The same signs: addition */
    set_minus = INTEGER_IS_NEGATIVE(a);
    i = 1;
    carry = 0;
    while(i <= nb)                   /* Common part */
    {
      lw = (uint)carry + (uint)a[i] + (uint)b[i];
      carry = (lw > MAX_LIMB);
      c[i++] = (LIMB)lw;
    }
    while(i <= na)                          /* Tail */
    {
      lw = (uint)carry + (uint)a[i];
      carry = (lw > MAX_LIMB);
      c[i++] = (LIMB)lw;
    }
    if(carry)
      c[i] = 1;
    else
      i--;
#if defined(INTEGER_MAX_SIZE)
    if(i > IntegerMaxSize)
      IntegerMaxSize = i;
#endif
  }
  else              /* Different signs: subtraction */
  {
    if(na == nb)
    {
      i = na;
      while(i > 0)
      {
        if(a[i] < b[i])
        {                           /* Swap numbers */
          BIGINT w = a;
          a = b;
          b = w;
          goto subtract;
        }
        else if(a[i] != b[i])
          goto subtract;
        i--;
      }
      c[0] = 0;       /* Zero result of subtraction */
      goto out;
    }
    subtract:
    set_minus = INTEGER_IS_NEGATIVE(a);
    i = carry = 1;
    while(i <= nb)                   /* Common part */
    {
      lw = MAX_LIMB + (uint)carry + (uint)a[i] - (uint)b[i];
      carry = (lw > MAX_LIMB);
      c[i++] = (LIMB)lw;
    }
    while(i <= na)
    {
      lw = MAX_LIMB + (uint)carry + (uint)a[i];
      carry = (lw > MAX_LIMB);
      c[i++] = (LIMB)lw;
    }
    while(--i > 0)
      if(c[i] != 0)
        break;
  }
  c[0] = i;               /* Number of LIMBs in sum */
  if(set_minus)
    INTEGER_SET_MINUS(c);
  out:
OUT_INTEGER_SUM /*----------------------------------*/
  return;
}

/*_6_6          Copy and delete functions===============================*/

/*=LieSumCopyInteger===================================
  Integer regime
*/
uint LieSumCopyInteger(uint a)
{
  if(a != NIL)
  {
    uint ca, eca;
    BIGINT n, cn;
    int i;
IN_LIE_SUM_COPY  /*----------------------------------*/
    ca = eca = NodeLTNew();
    LIE_TERM_MONOMIAL(eca) = LIE_TERM_MONOMIAL(a);
    n = LIE_TERM_NUMERATOR_INTEGER(a);
    INTEGER_HEAP_COPY(cn, n, i);
    LIE_TERM_NUMERATOR_INTEGER(eca) = cn;
    if((n = LIE_TERM_DENOMINATOR_INTEGER(a)) != NULL)
    {
      INTEGER_HEAP_COPY(cn, n, i);
      LIE_TERM_DENOMINATOR_INTEGER(eca) = cn;
    }
    else
      LIE_TERM_DENOMINATOR_INTEGER(eca) = NULL;
    while((a = LIE_TERM_R(a)) != NIL)
    {
      LIE_TERM_R(eca) = NodeLTNew();
      eca = LIE_TERM_R(eca);
      LIE_TERM_MONOMIAL(eca) = LIE_TERM_MONOMIAL(a);
      n = LIE_TERM_NUMERATOR_INTEGER(a);
      INTEGER_HEAP_COPY(cn, n, i);
      LIE_TERM_NUMERATOR_INTEGER(eca)= cn;
      if((n = LIE_TERM_DENOMINATOR_INTEGER(a)) != NULL)
      {
        INTEGER_HEAP_COPY(cn, n, i);
        LIE_TERM_DENOMINATOR_INTEGER(eca) = cn;
      }
      else
        LIE_TERM_DENOMINATOR_INTEGER(eca) = NULL;
    }
OUT_LIE_SUM_COPY /*----------------------------------*/
    return ca;
  }
  return NIL;
}
/*=LieSumCopyIntegerNegative===================================
  Copy changing sign. Integer regime
*/
uint LieSumCopyIntegerNegative(uint a)
{
  if(a != NIL)
  {
    uint ca, eca;
    BIGINT n, cn;
    int i;
    ca = eca = NodeLTNew();
    LIE_TERM_MONOMIAL(eca) = LIE_TERM_MONOMIAL(a);
    n = LIE_TERM_NUMERATOR_INTEGER(a);
    INTEGER_HEAP_COPY(cn, n, i);
    INTEGER_MINUS(cn);
    LIE_TERM_NUMERATOR_INTEGER(eca) = cn;
    if((n = LIE_TERM_DENOMINATOR_INTEGER(a)) != NULL)
    {
      INTEGER_HEAP_COPY(cn, n, i);
      LIE_TERM_DENOMINATOR_INTEGER(eca) = cn;
    }
    else
      LIE_TERM_DENOMINATOR_INTEGER(eca) = NULL;
    while((a = LIE_TERM_R(a)) != NIL)
    {
      LIE_TERM_R(eca) = NodeLTNew();
      eca = LIE_TERM_R(eca);
      LIE_TERM_MONOMIAL(eca) = LIE_TERM_MONOMIAL(a);
      n = LIE_TERM_NUMERATOR_INTEGER(a);
      INTEGER_HEAP_COPY(cn, n, i);
      INTEGER_MINUS(cn);
      LIE_TERM_NUMERATOR_INTEGER(eca)= cn;
      if((n = LIE_TERM_DENOMINATOR_INTEGER(a)) != NULL)
      {
        INTEGER_HEAP_COPY(cn, n, i);
        LIE_TERM_DENOMINATOR_INTEGER(eca) = cn;
      }
      else
        LIE_TERM_DENOMINATOR_INTEGER(eca) = NULL;
    }
    return ca;
  }
  return NIL;
}
/*=LieSumCopyParametric=====================================
  Parametric regime
*/
uint LieSumCopyParametric(uint a)
{
  if(a != NIL)
  {
    uint ca, eca;
IN_LIE_SUM_COPY  /*---------------------------------------*/
    ca = eca = NodeLTNew();
    LIE_TERM_MONOMIAL(eca) = LIE_TERM_MONOMIAL(a);
    LIE_TERM_NUMERATOR_SCALAR_SUM(eca) =

      ScalarSumCopy(LIE_TERM_NUMERATOR_SCALAR_SUM(a));
    if(LIE_TERM_DENOMINATOR_SCALAR_SUM(a) != NIL)
      LIE_TERM_DENOMINATOR_SCALAR_SUM(eca) =
        ScalarSumCopy(LIE_TERM_DENOMINATOR_SCALAR_SUM(a));
    else
      LIE_TERM_DENOMINATOR_SCALAR_SUM(eca) = NIL;
    while((a = LIE_TERM_R(a)) != NIL)
    {
      LIE_TERM_R(eca) = NodeLTNew();
      eca = LIE_TERM_R(eca);
      LIE_TERM_MONOMIAL(eca) = LIE_TERM_MONOMIAL(a);
      LIE_TERM_NUMERATOR_SCALAR_SUM(eca) =
        ScalarSumCopy(LIE_TERM_NUMERATOR_SCALAR_SUM(a));
      if(LIE_TERM_DENOMINATOR_SCALAR_SUM(a) != NIL)
        LIE_TERM_DENOMINATOR_SCALAR_SUM(eca) =
          ScalarSumCopy(LIE_TERM_DENOMINATOR_SCALAR_SUM(a));
      else
        LIE_TERM_DENOMINATOR_SCALAR_SUM(eca) = NIL;
    }
OUT_LIE_SUM_COPY /*---------------------------------------*/
    return ca;
  }
  return NIL;
}
/*=LieSumKillInteger================================
 a == NIL is admitted (Integer regime)
*/
void LieSumKillInteger(uint a)
{
  uint b;
  BIGINT d;
IN_LIE_SUM_KILL  /*-------------------------------*/
  while(a != NIL)
  {
    b = a;
    a = LIE_TERM_R(a);
    INTEGER_KILL(LIE_TERM_NUMERATOR_INTEGER(b));
    if((d = LIE_TERM_DENOMINATOR_INTEGER(b)) != NULL)
      INTEGER_KILL(d);
    NODE_LT_KILL(b);
  }
OUT_LIE_SUM_KILL  /*------------------------------*/
}
/*=LieSumKillParametric===============================
 a == NIL is admitted (Parametric regime)
*/
void LieSumKillParametric(uint a)
{
  uint b;
IN_LIE_SUM_KILL  /*---------------------------------*/
  while(a != NIL)
  {
    b = a;
    a = LIE_TERM_R(a);
    ScalarSumKill(LIE_TERM_NUMERATOR_SCALAR_SUM(b));
    ScalarSumKill(LIE_TERM_DENOMINATOR_SCALAR_SUM(b));
    NODE_LT_KILL(b);
  }
OUT_LIE_SUM_KILL  /*--------------------------------*/
}
/*=LieTermFromMonomialInteger============
*/
uint LieTermFromMonomialInteger(int mon)
{
  BIGINT num;
  uint a = NodeLTNew();

  LIE_TERM_MONOMIAL(a) = mon;

  INTEGER_HEAP_NEW(num, 2);
  num[0] = num[1] = 1;
  LIE_TERM_NUMERATOR_INTEGER(a) = num;

  LIE_TERM_DENOMINATOR_INTEGER(a) = NULL;

  return a;
}
/*=LieTermFromMonomialParametric===========
*/
uint LieTermFromMonomialParametric(int mon)
{
  BIGINT num;
  uint c = NodeSTNew(),
    a = NodeLTNew();

  LIE_TERM_MONOMIAL(a) = mon;

  SCALAR_TERM_MONOMIAL(c) = NIL;
  INTEGER_HEAP_NEW(num, 2);
  num[0] = num[1] = 1;
  SCALAR_TERM_NUMERATOR(c) = num;

  LIE_TERM_NUMERATOR_SCALAR_SUM(a) = c;

  LIE_TERM_DENOMINATOR_SCALAR_SUM(a) = NIL;

  return a;
}
/*=ScalarSumCopy=======================================
 Caller ensures a != NIL
*/
uint ScalarSumCopy(uint a)
{
  int i;
  BIGINT n, o;
  uint ca, bca, b, cb;
  bca = ca = NodeSTNew();
  while(YES)
  {
    /* Copy integer coefficient */

    o = SCALAR_TERM_NUMERATOR(a);
    INTEGER_HEAP_COPY(n, o, i);
    SCALAR_TERM_NUMERATOR(ca) = n;

    /* Copy scalar monomial */

    b = SCALAR_TERM_MONOMIAL(a);
    if(b != NIL)
    {
      SCALAR_TERM_MONOMIAL(ca) = cb = NodeSFNew();
      SCALAR_FACTOR_WORD(cb) = SCALAR_FACTOR_WORD(b);
      while((b = SCALAR_FACTOR_R(b)) != NIL)
      {
        SCALAR_FACTOR_R(cb) = NodeSFNew();
        cb = SCALAR_FACTOR_R(cb);
        SCALAR_FACTOR_WORD(cb) = SCALAR_FACTOR_WORD(b);
      }
    }
    else
      SCALAR_TERM_MONOMIAL(ca) = NIL;

    if((a = SCALAR_TERM_R(a)) == NIL)
      break;
    SCALAR_TERM_R(ca) = NodeSTNew();
    ca = SCALAR_TERM_R(ca);
  }
  return bca;
}
/*=ScalarSumKill=================================================
  Only at IsParametric == YES, a == NIL is admitted
*/
void ScalarSumKill(uint a)
{
  uint b, c;
  while(a != NIL)
  {
    b = a;
    a = SCALAR_TERM_R(a);
    INTEGER_KILL(SCALAR_TERM_NUMERATOR(b));
    c = SCALAR_TERM_MONOMIAL(b);
    NODE_ST_KILL(b);
    while(c != NIL)   /* Scalar monomial may be NIL, uint-uint mix */
    {
      b = c;
      c = SCALAR_FACTOR_R(c);
      NODE_SF_KILL(b);
    }
  }
}
/*=ScalarTermCopy====================================
 Caller ensures a != NIL
*/
uint ScalarTermCopy(uint a)
{
  int i;
  BIGINT cn, n;
  uint m, ca = NodeSTNew();

  /* Copy integer coefficient */

  n = SCALAR_TERM_NUMERATOR(a);
  INTEGER_HEAP_COPY(cn, n, i);
  SCALAR_TERM_NUMERATOR(ca) = cn;

  /* Copy monomial */

  if((m = SCALAR_TERM_MONOMIAL(a)) != NIL)
  {
    uint cm;
    SCALAR_TERM_MONOMIAL(ca) = cm = NodeSFNew();
    SCALAR_FACTOR_WORD(cm) = SCALAR_FACTOR_WORD(m);
    while((m = SCALAR_FACTOR_R(m)) != NIL)
    {
      SCALAR_FACTOR_R(cm) = NodeSFNew();
      cm = SCALAR_FACTOR_R(cm);
      SCALAR_FACTOR_WORD(cm) = SCALAR_FACTOR_WORD(m);
    }
  }
  else
    SCALAR_TERM_MONOMIAL(ca) = NIL;
  return ca;
}

/*_6_7          Technical functions=====================================*/

/*=Error!===============
*/
void Error(int i_message)
{
  PutMessage(ERROR);
  PutMessage(i_message);
  EXIT;
}
/*=Initialization==========================================================
*/
void Initialization(void)
{
  FILE *inif;
  short c;
  uint i, j;
  char * init_case[N_INIT_CASES];

  /* Set cases strings */

  init_case[COEFFICIENT_SUM_TABLE_SIZE] = "Co";
  init_case[CRUDE_TIME]                 = "Cr";
  init_case[ECHO_INPUT_FILE]            = "Ec";
  init_case[EVEN_BASIS_SYMBOL]          = "Ev";
  init_case[GAP_ALGEBRA_NAME]           = "GAP a";
  init_case[GAP_BASIS_NAME]             = "GAP b";
  init_case[GAP_RELATIONS_NAME]         = "GAP r";
  init_case[GAP_OUTPUT_BASIS]           = "GAP output b";
  init_case[GAP_OUTPUT_COMMUTATORS]     = "GAP output c";
  init_case[GAP_OUTPUT_RELATIONS]       = "GAP output r";
  init_case[GENERATOR_MAX_N]            = "Ge";
  init_case[INPUT_DIRECTORY]            = "Input d";
  init_case[INPUT_INTEGER_SIZE]         = "Input i";
  init_case[INPUT_STRING_SIZE]          = "Input s";
  init_case[LEFT_NORMED_OUTPUT]         = "Le";
  init_case[LIE_MONOMIAL_SIZE]          = "Lie";
  init_case[LINE_LENGTH]                = "Lin";
  init_case[NAME_LENGTH]                = "Na";
  init_case[NODE_LT_SIZE]               = "Node L";
  init_case[NODE_SF_SIZE]               = "Node scalar f";
  init_case[NODE_ST_SIZE]               = "Node scalar t";
  init_case[ODD_BASIS_SYMBOL]           = "Od";
  init_case[OUT_LINE_SIZE]              = "Ou";
  init_case[PARAMETER_MAX_N]            = "Pa";
  init_case[PUT_BASIS_ELEMENTS]         = "Put b";
  init_case[PUT_COMMUTATORS]            = "Put c";
  init_case[PUT_HILBERT_SERIES]         = "Put H";
  init_case[PUT_INITIAL_RELATIONS]      = "Put i";
  init_case[PUT_NON_ZERO_COEFFICIENTS]  = "Put n";
  init_case[PUT_PROGRAM_HEADING]        = "Put p";
  init_case[PUT_REDUCED_RELATIONS]      = "Put r";
  init_case[PUT_STATISTICS]             = "Put s";
  init_case[RELATION_SIZE]              = "Re";

#if !defined(GAP)
#if defined(SPP_2000)
  MessageFile = OpenFile("fplsa4.msg", "r");
  SessionFile = OpenFile("fplsa4.ses", "w");
  inif = OpenFile("fplsa4.ini", "r");
#else
  MessageFile = OpenFile("fplsa4.msg", "rt");
  SessionFile = OpenFile("fplsa4.ses", "wt");
  inif = OpenFile("fplsa4.ini", "rt");
#endif
#else
#if defined(SPP_2000)
  inif = OpenFile("fplsa4.ini", "r");
#else
  inif = OpenFile("fplsa4.ini", "rt");
#endif
#endif
  while(YES)
    switch(ReadCaseFromFile(inif, init_case, N_INIT_CASES))
    {
      case COEFFICIENT_SUM_TABLE_SIZE:
        CoeffSumTableSize = ReadDecimalFromFile(inif);
      break;
      case CRUDE_TIME:
        CrudeTime = ReadBooleanFromFile(inif);
      break;
      case ECHO_INPUT_FILE:
        EchoInput = ReadBooleanFromFile(inif);
      break;
      case EVEN_BASIS_SYMBOL:
        BasisSymbolEven = fgetc(inif);
      break;
      case GAP_ALGEBRA_NAME:
        if(ReadStringFromFile(GAPAlgebraName, inif) == EOF)
          goto out;
      break;
      case GAP_BASIS_NAME:
        if(ReadStringFromFile(GAPBasisName, inif) == EOF)
          goto out;
      break;
      case GAP_RELATIONS_NAME:
        if(ReadStringFromFile(GAPRelationsName, inif) == EOF)
          goto out;
      break;
      case GAP_OUTPUT_BASIS:
        GAPOutputBasis = ReadBooleanFromFile(inif);
      break;
      case GAP_OUTPUT_COMMUTATORS:
        GAPOutputCommutators = ReadBooleanFromFile(inif);
      break;
      case GAP_OUTPUT_RELATIONS:
        GAPOutputRelations = ReadBooleanFromFile(inif);
      break;
      case GENERATOR_MAX_N:
        GeneratorMaxN = ReadDecimalFromFile(inif);
      break;
      case INPUT_DIRECTORY:
        while((c = fgetc(inif)) != '\n')
        {
          if(c == LEFT_COMMENT)
          {
            ungetc(c, inif);
            break;
          }
          if(c == EOF)
            goto out;
          if(!isspace(c))
            OutLine[PosOutLine++] = (char)c;
        }
      break;
      case INPUT_INTEGER_SIZE:
        InputIntegerSize =  (uint)ReadDecimalFromFile(inif);
        InputIntegerSize++; /* For head */
      break;
      case INPUT_STRING_SIZE:
        InputStringSize =  (uint)ReadDecimalFromFile(inif);
      break;
      case LEFT_NORMED_OUTPUT:
        if(ReadBooleanFromFile(inif))
          PutLieMonomial = PutLieMonomialLeftNormed;
      break;
      case LINE_LENGTH:
        LineLength = (uint)ReadDecimalFromFile(inif);
      break;
      case LIE_MONOMIAL_SIZE:
        LieMonomialSize = (int)ReadDecimalFromFile(inif);
        LieMonomial = (LIE_MON*)NewArray(LieMonomialSize, sizeof(LIE_MON),
                                                         E_A_LIE_MONOMIAL);
      break;
      case NAME_LENGTH:
        NameLength1 = ReadDecimalFromFile(inif);
        NameLength1++;
      break;
      case NODE_LT_SIZE:
        NodeLTSize = (uint)ReadDecimalFromFile(inif);
        NodeLT =
          (NODE_LT*)NewArray(NodeLTSize, sizeof(NODE_LT), E_A_NODE_LT);
        i = 1;
        j = 2;
        while(j < NodeLTSize) /* Install NodeLT for Lie terms */
          LIE_TERM_R(i++) = j++;
#if defined(SPACE_STATISTICS)
        LIE_TERM_R(i) = NOTHING;
#else
        LIE_TERM_R(i) = NIL;
#endif
      break;
      case NODE_SF_SIZE:
        NodeSFSize = (uint)ReadDecimalFromFile(inif);
      break;
      case NODE_ST_SIZE:
        NodeSTSize = (uint)ReadDecimalFromFile(inif);
      break;
      case ODD_BASIS_SYMBOL:
        BasisSymbolOdd = fgetc(inif);
      break;
      case OUT_LINE_SIZE:
        OutLineSize = (uint)ReadDecimalFromFile(inif);
        OutLine = (char *)NewArray(OutLineSize, 1, E_A_OUT_LINE);
      break;
      case PARAMETER_MAX_N:
        ParameterMaxN = ReadDecimalFromFile(inif);
      break;
      case PUT_BASIS_ELEMENTS:
        BasisElementsPut = ReadBooleanFromFile(inif);
      break;
      case PUT_COMMUTATORS:
        CommutatorsPut = ReadBooleanFromFile(inif);
      break;
      case PUT_HILBERT_SERIES:
        HilbertSeriesPut = ReadBooleanFromFile(inif);
      break;
      case PUT_INITIAL_RELATIONS:
        InitialRelationsPut = ReadBooleanFromFile(inif);
      break;
      case PUT_NON_ZERO_COEFFICIENTS:
        NonZeroCoefficientsPut = ReadBooleanFromFile(inif);
      break;
      case PUT_PROGRAM_HEADING:
        HeadingPut = ReadBooleanFromFile(inif);
      break;
      case PUT_REDUCED_RELATIONS:
        ReducedRelationsPut = ReadBooleanFromFile(inif);
      break;
      case PUT_STATISTICS:
        StatisticsPut = ReadBooleanFromFile(inif);
      break;
      case RELATION_SIZE:
        RelationSize = (int)ReadDecimalFromFile(inif);
        Relation = (REL*)NewArray(RelationSize, sizeof(REL), E_A_RELATION);
      break;
      case EOF:
        goto out;
      default:
        Error(E_WRONG_INI_CASE);
    }
  out:
  fclose(inif);
  if(HeadingPut)
    PutMessage(H_PROGRAM);
  GeneratorName = (char *)NewArray(GeneratorMaxN*NameLength1, sizeof(char),
                                                 E_A_GENERATOR_NAME);
}
/*=NewArray=========================================
*/
void * NewArray(uint n, uint size, int i_message)
{
  void * new_pointer = (void *)calloc(n, size);
  if(new_pointer == NULL && n != 0)
  {
    const char * format = "%u elements of size %u\n%u bytes\n";
    PutMessage(E_ALLOC);
    PutMessage(i_message);
#if defined(ECHO_TO_SCREEN)
    printf(format, n,  size, n*size);
#endif
#if !defined(GAP)
    fprintf(SessionFile, format, n,  size, n*size);
#endif
    EXIT;
  }
  return new_pointer;
}
/*=NodeLTNew===================
 Get node from NodeLT pool.
*/
uint NodeLTNew(void)
{
  uint a = NodeLTTop;
#if !defined(SPACE_STATISTICS)
  if(a == NIL)
  {
    TIME_OFF;
    PutStatistics();
    Error(E_NODE_LT_SIZE);
  }
#endif
  NodeLTTop = LIE_TERM_R(a);
#if defined(SPACE_STATISTICS)
  if(NodeLTTopMax < NodeLTTop)
  {
    if(NodeLTTop > NodeLTSize)
    {
      TIME_OFF;
      PutStatistics();
      Error(E_NODE_LT_SIZE);
    }
    NodeLTTopMax = NodeLTTop;
  }
#endif
  LIE_TERM_R(a) = NIL;
PP_CURRENT_N_LT     /* MEMORY */
    return a;
}
/*=NodeSFNew=====================
 Get node from NodeSF pool.
*/
uint NodeSFNew(void)
{
  uint a = NodeSFTop;
#if !defined(SPACE_STATISTICS)
  if(a == NIL)
  {
    TIME_OFF;
    PutStatistics();
    Error(E_NODE_SF_SIZE);
  }
#endif
  NodeSFTop = SCALAR_FACTOR_R(a);
#if defined(SPACE_STATISTICS)
  if(NodeSFTopMax < NodeSFTop)
  {
    if(NodeSFTop > NodeSFSize)
    {
      TIME_OFF;
      PutStatistics();
      Error(E_NODE_SF_SIZE);
    }
    NodeSFTopMax = NodeSFTop;
  }
#endif
  SCALAR_FACTOR_R(a) = NIL;
PP_CURRENT_N_SF     /* MEMORY */
    return a;
}
/*=NodeSTNew===================
 Get node from NodeST pool.
*/
uint NodeSTNew(void)
{
  uint a = NodeSTTop;
#if !defined(SPACE_STATISTICS)
  if(a == NIL)
  {
    TIME_OFF;
    PutStatistics();
    Error(E_NODE_ST_SIZE);
  }
#endif
  NodeSTTop = SCALAR_TERM_R(a);
#if defined(SPACE_STATISTICS)
  if(NodeSTTopMax < NodeSTTop)
  {
    if(NodeSTTop > NodeSTSize)
    {
      TIME_OFF;
      PutStatistics();
      Error(E_NODE_ST_SIZE);
    }
    NodeSTTopMax = NodeSTTop;
  }
#endif
  SCALAR_TERM_R(a) = NIL;
PP_CURRENT_N_ST     /* MEMORY */
    return a;
}
/*=OpenFile================================
*/
FILE *OpenFile(char * file_name, char * file_type)
{
  FILE *file = fopen(file_name, file_type);
  if(file == NULL)
  {
    printf("\nNo file: %s", file_name);
    exit(1);
  }
  return file;
}

/*_6_8 Input functions=========================================*/

/*=BinaryQuestion!================
*/
int BinaryQuestion(int i_message)
{
  char c[2];
  get_symbol:
  PutMessage(i_message);
  scanf("%1s", c);
#if !defined(GAP)
  fputc(c[0], SessionFile);
#endif
  switch(c[0])
  {
    case 'y': case 'Y': case '\n':
      return YES;
    case 'n': case 'N':
      break;
    case 'c': case 'C':
      Error(E_CANCEL_PROGRAM);
    default:
      goto get_symbol;
  }
  return NO;
}
/*=FindNameInTable==============================================
 Find name from string in table ...NameIn
*/
int FindNameInTable(char * name, char * nametab, int n_nametab)
{
  char *w_nametab, *w_name;
  int j = 0;
  while(j < n_nametab)
  {
    w_nametab = nametab;
    w_name = name;
    while(YES)
    {
      if(*w_nametab == '\0')              /* Table name ended */
      {
        if(!isalnum(*w_name) && *w_name != SUBSCRIPT_INPUT_SIGN)
          goto out;                      /* String name ended */
        break;
      }                                 /* Table name goes on */
      if(!isalnum(*w_name) && *w_name != SUBSCRIPT_INPUT_SIGN)
        break;                           /* String name ended */
      if(*w_nametab != *w_name)
        break;                             /* Different names */
      w_nametab++;
      w_name++;
    }
    nametab += NameLength1;
    ++j;
  }
  out:
  return j;
}
/*=GetGenerator===========================================================
 Read single generator from description string
*/
void GetGenerator(char * str)
{
  char * name = GeneratorName + GeneratorN*NameLength1;
  if(GeneratorN == GeneratorMaxN)
    Error(E_GENERATOR_MAX_N);
  do
  {
    if(*str == ODD_GENERATOR_INPUT_SIGN)        /* EVEN == 0 is assumed */
      LIE_MONOMIAL_PARITY(GeneratorN) = ODD;
    else
      *(name++) = *str;
  }while(*(str++));
  LIE_MONOMIAL_ORDER(GeneratorN) =                 /* Standard settings */
  LIE_MONOMIAL_POSITION(GeneratorN) = GeneratorN;
  LIE_MONOMIAL_RIGHT(GeneratorN) = 1; /* To avoid SQUARE interpretation */
  LIE_MONOMIAL_WEIGHT(GeneratorN++) = 1;
}
/*=GetInput===============================================================
*/
void GetInput(int n, char * fin)
{
  char *instr, *sfname, *in_case[N_INPUT_CASES];
  FILE *inf;
  uint i, j;

  in_case[GENERATORS]      = "G";
  in_case[LIMITING_WEIGHT] = "L";
  in_case[PARAMETERS]      = "P";
  in_case[RELATIONS]       = "R";
  in_case[WEIGHTS]         = "W";

  sfname = OutLine + PosOutLine;
  if(n == 1)      /* No input file at call */
  {
    PutMessage(H_ENTER_FILE);
    if (fgets(sfname, OutLineSize - PosOutLine, stdin)) {
      // success
      sfname[strcspn(sfname, "\n")] = '\0';
    }
    else {
      // failure
      exit(1);
    }
  }
  else
    do
      *(sfname++) = *fin;
    while(*(fin++));
  sfname = OutLine;
  while(YES)
  {
    if(*sfname == '.')
      break;
    if(*sfname == '\0')
    {
      *sfname = '.';
      *++sfname = 'i';
      *++sfname = 'n';
      *++sfname = '\0';
      break;
    }
    ++sfname;
  }
  /*
  PutMessage(H_INPUT_FILE);
  PutStringStandard(OutLine);
  */
#if defined(SPP_2000)
  if((inf = fopen(OutLine, "r")) == NULL)
#else
  if((inf = fopen(OutLine, "rt")) == NULL)
#endif
  { /* New file */
    if(!BinaryQuestion(H_CREATE_NEW_FILE))
      exit(1);
    instr = (char *)alloca(InputStringSize);
    if(instr == NULL)
      Error(E_A_STACK_INPUT_STRING);
    fgetc(stdin);
#if defined(SPP_2000)
    inf = OpenFile(OutLine, "w");
#else
    inf = OpenFile(OutLine, "wt");
#endif
    KeyBoardStringToFile(H_ENTER_GENERATORS, "Generators: ", instr, inf);
    KeyBoardStringToFile(H_ENTER_WEIGHTS_IN_FILE, "Weights: ", instr, inf);
    KeyBoardStringToFile(H_ENTER_LIMITING_WEIGHT, "Limiting weight: ",
                                                               instr, inf);
    KeyBoardStringToFile(H_ENTER_PARAMETERS, "Parameters: ", instr, inf);
    if(KeyBoardStringToFile(H_ENTER_RELATIONS, "Relations:\n", instr, inf))
    {
      n = YES; /* Now non-last */
      while(n && (j = KeyBoardBytesToString(instr)) > 0)
        for(i = 0; i <= j; i++)
        {
          if(instr[i] == '.')
            n = NO;             /* Last input */
          fputc(instr[i], inf); /* Copy entered string to file */
        }
    }
    fclose(inf);
#if defined(SPP_2000)
    inf = OpenFile(OutLine, "r");
#else
    inf = OpenFile(OutLine, "rt");
#endif
  }
  if(EchoInput)
  {
    short c;
    PutMessage(H_SHOW_INPUT);
    while((c = fgetc(inf)) != EOF)
      PutCharacter((char)c);
  }
  rewind(inf);
  while(YES)
    switch(ReadCaseFromFile(inf, in_case, N_INPUT_CASES))
    {
      case GENERATORS:
        ReadAndProcessStringsFromFile(GetGenerator, inf, ' ', ';');
        CUT_ARRAY(GeneratorName, char, NameLength1*GeneratorN);
        LieMonomialFreePosition = LieMonomialN = GeneratorN;
#if defined(SPACE_STATISTICS)
        LieMonomialMaxN = LieMonomialN;
#endif
      break;
      case WEIGHTS:    /* Should come after generators */
        GeneratorMaxN = GeneratorN;
        GeneratorN = 0;
        ReadAndProcessStringsFromFile(GetWeight, inf, ' ', ';');
        GeneratorN = GeneratorMaxN;

        /* Reorder generators in accordance with new weights */

        i = 1;
        while(i < (uint)GeneratorN)
        {
          for(j = GeneratorN - 1; j >= i; j--)
            if(LIE_MONOMIAL_WEIGHT(j-1) > LIE_MONOMIAL_WEIGHT(j))
            {
              byte wt;       /* To save swapped walues */

              /* Swap generator names */

              fin = GeneratorName + j*NameLength1; /* Next name */
              instr = fin - NameLength1;           /* Previous name */
              for(n = 1; n < NameLength1; n++)
              {
                wt = *instr;
                *(instr++) = *fin;
                *(fin++) = wt;
              }

              /* Swap weights */

              wt = LIE_MONOMIAL_WEIGHT(j-1);
              LIE_MONOMIAL_WEIGHT(j-1) = LIE_MONOMIAL_WEIGHT(j);
              LIE_MONOMIAL_WEIGHT(j) = wt;

              /* Swap parities */

              wt = LIE_MONOMIAL_PARITY(j-1);
              LIE_MONOMIAL_PARITY(j-1) = LIE_MONOMIAL_PARITY(j);
              LIE_MONOMIAL_PARITY(j) = wt;
            }
          i++;
        }
      break;
      case PARAMETERS:
        IsParametric = YES;
        LieLikeTermsCollection = LieLikeTermsCollectionParametric;
        LieSumCopy = LieSumCopyParametric;
        LieSumKill = LieSumKillParametric;
        LieSumMinus = LieSumMinusParametric;
        NormalizeRelation = NormalizeRelationParametric;
        LieTermFromMonomial = LieTermFromMonomialParametric;
        PairMonomialMonomial = PairMonomialMonomialParametric;
        PairMonomialSum = PairMonomialSumParametric;
        PairSumMonomial = PairSumMonomialParametric;
        PairSumSum = PairSumSumParametric;
        SubstituteRelationInRelation =
                            SubstituteRelationInRelationParametric;
        ParameterName = (char *)NewArray(NameLength1*ParameterMaxN,
                                    sizeof(char), E_A_PARAMETER_NAME);
        ParameterName[0] = 'i';         /* Obligatory imaginary unit */
        ParameterN = 1;
        ReadAndProcessStringsFromFile(GetParameter, inf, ' ', ';');
        CUT_ARRAY(ParameterName, char, NameLength1*ParameterN);
        NodeST =
          (NODE_ST*)NewArray(NodeSTSize, sizeof(NODE_ST), E_A_NODE_ST);
        i = 1;
        j = 2;
        while(j < NodeSTSize)     /* Install NodeST for scalar terms */
          SCALAR_TERM_R(i++) = j++;
#if defined(SPACE_STATISTICS)
        SCALAR_TERM_R(i) = NOTHING;
#else
        SCALAR_TERM_R(i) = NIL;
#endif
        NodeSF =
          (NODE_SF*)NewArray(NodeSFSize, sizeof(NODE_SF), E_A_NODE_SF);
        i = 1;
        j = 2;
        while(j < NodeSFSize)   /* Install NodeSF for scalar factors */
          SCALAR_FACTOR_R(i++) = j++;
#if defined(SPACE_STATISTICS)
        SCALAR_FACTOR_R(i) = NOTHING;
#else
        SCALAR_FACTOR_R(i) = NIL;
#endif
      break;
      case RELATIONS:
        if(!IsParametric)
          NonZeroCoefficientsPut = NO;
        ReadAndProcessStringsFromFile(GetRelation, inf, ';', '.');
      break;
      case LIMITING_WEIGHT:
        LimitingWeight = ReadDecimalFromFile(inf);
        while(fgetc(inf) != '\n')                 /* Go to next line */
          ;
      break;
      case EOF:
        goto out;
      default:
        Error(E_WRONG_INPUT_CASE);
    }
  out:
  fclose(inf);
  if(LimitingWeight == 0)
  {
    PutMessage(H_ENTER_LIMITING_WEIGHT);
    scanf("%d", &LimitingWeight);
#if !defined(GAP)
    fprintf(SessionFile, "%d\n", LimitingWeight);
#endif
  }
  TIME_ON;  /* First start of time */
}
/*=GetInteger===========================================================
 Read big integer with shift in string.
 A is already allocated array of LIMBs A[].
*/
void GetInteger(BIGINT a, char **pstr)
{
  BIGINT w;
  int i;
  LIMB digit[2], ten[2];
  digit[0] = ten[0] = 1;
  ten[1] = 10;
  INTEGER_STACK_NEW(w, InputIntegerSize);
  while(**pstr == '0') /* Skip leading zeros */
    ++*pstr;
  if(isdigit(**pstr))
  {                     /* First digit: IntegerProduct does not eat 0 */
    a[0] = 1;
    a[1] = (LIMB)(**pstr - '0');
    while(isdigit(*(++*pstr)))
    {
      IntegerProduct(w, a, ten);                              /* a*10 */
      i = w[0];
      while(i >= 0) {
        a[i] = w[i];                                 /* Copy w to a */
        i--;
      }
      if((digit[1] = (LIMB)(**pstr - '0')) != 0)
        IntegerSum(a, a, digit);                      /* a*10 + digit */
    }
  }
  else        /* Caller interprets absence as  */
    a[0] = 0; /* 1 (or 0) depending on context */
}
/*=GetLieMonomial=========================================================
 Read monomial from string with transformations and substitutions
*/
uint GetLieMonomial(char **pstr)
{
  int mon;
  uint a;
IN_GET_LIE_MONOMIAL   /*------------------------------------------------*/
  if(isalpha(**pstr))
    if((mon=FindNameInTable(*pstr,GeneratorName,GeneratorN)) < GeneratorN)
    {
      BIGINT num;
      SkipName(pstr);
      SkipSpaces(pstr);
      a = NodeLTNew();
      LIE_TERM_MONOMIAL(a) = mon;
      INTEGER_HEAP_NEW(num, 2);     /* Make integer 1 */
      num[0] = num[1] = 1;
      if(IsParametric)
      {
        uint st = NodeSTNew();
        SCALAR_TERM_MONOMIAL(st) = NIL;
        SCALAR_TERM_NUMERATOR(st) = num;
        LIE_TERM_NUMERATOR_SCALAR_SUM(a) = st;
        LIE_TERM_DENOMINATOR_SCALAR_SUM(a) = NIL;
      }
      else
      {
        LIE_TERM_NUMERATOR_INTEGER(a) = num;
        LIE_TERM_DENOMINATOR_INTEGER(a) = NULL;
      }
    }
    else
      Error(E_UNDECLARED_GENERATOR);
  else if(**pstr == '[')
  {
    uint b;
    SkipSpaces(pstr);
    ++*pstr;
    a = GetLieMonomial(pstr);
    SkipSpaces(pstr);
    if(**pstr != ',')
      Error(E_NO_COMMUTATOR_COMMA);
    ++*pstr;
    SkipSpaces(pstr);
    b = GetLieMonomial(pstr);
    SkipSpaces(pstr);
    if(**pstr != ']')
      Error(E_NO_COMMUTATOR_BRACKET);
    ++*pstr;
    a = (*PairSumSum)(a, b);
  }
  else
    Error(E_INVALID_CHARACTER);
OUT_GET_LIE_MONOMIAL  /*------------------------------------------------*/
  return a;
}
/*=GetLieSum=====================================================
 Read Lie expression from string and make internal representation
*/
uint GetLieSum(char **pstr)
{
  uint lsum, term;
  int sign = PLUS;
IN_GET_LIE_SUM   /*--------------------------------------------*/
  SkipSpaces(pstr);
  if(**pstr == '-')
  {
    sign = MINUS;
    ++*pstr;
    SkipSpaces(pstr);
  }
  lsum = GetLieTerm(pstr);
  SkipSpaces(pstr);
  if(sign)
    LieSumMinus(lsum);
  while(**pstr == '+' || **pstr == '-')
  {
    sign = (**pstr == '+') ? PLUS : MINUS;
    ++*pstr;
    SkipSpaces(pstr);
    term = GetLieTerm(pstr);
    SkipSpaces(pstr);
    if(sign)
      LieSumMinus(term);
    lsum = LieSumAddition(lsum, term);
  }
OUT_GET_LIE_SUM  /*--------------------------------------------*/
  return lsum;
}
/*=GetLieTerm===========================================================
*/
uint GetLieTerm(char **pstr)
{
  uint lterm;
IN_GET_LIE_TERM   /*--------------------------------------------------*/
  if(IsParametric)
  {
    uint num = GetScalarSum(pstr),
       den = NIL;
    SkipSpaces(pstr);
    if(**pstr == '/')
    {
      ++*pstr;
      SkipSpaces(pstr);
      den = GetScalarSum(pstr);
      ScalarSumCancellation(&num, &den);
    }
    lterm = GetLieMonomial(pstr);    /* May be sum with generic coeffs */
    if(den != NIL)
      LieSumDivScalarSum(lterm, den);
    LieSumMultScalarSum(lterm, num);
  }
  else
  {
    BIGINT num, den;
    INTEGER_STACK_NEW(num, InputIntegerSize);
    INTEGER_STACK_NEW(den, InputIntegerSize);
    GetInteger(num, pstr);
    den[0] = 0;
    SkipSpaces(pstr);
    if(**pstr == '/')
    {
      ++*pstr;
      SkipSpaces(pstr);
      GetInteger(den, pstr);
      SkipSpaces(pstr);
      IntegerCancellation(num, den);
    }
    lterm = GetLieMonomial(pstr);    /* May be sum with generic coeffs */
    if(den[0] != 0 && (den[0] != 1 || den[1] != 1))
      LieSumDivInteger(lterm, den);
    if(num[0] != 0 &&
       (INTEGER_IS_NEGATIVE(num) || num[0] != 1 || num[1] != 1))
      LieSumMultInteger(lterm, num);
  }
OUT_GET_LIE_TERM  /*--------------------------------------------------*/
  return lterm;
}
/*=GetUInteger====================================
 Read with shift long unsigned integer from string
*/
uint GetUInteger(char **pstr)
{
  uint i = 0;
  while(isdigit(**pstr))
  {
    i = i*10 + **pstr - '0';
    ++*pstr;
  }
  return i;
}
/*=GetParameter=====================================================
 Read single parameter from description string
*/
void GetParameter(char * str)
{
  if(str[0] != 'i' || str[1] != '\0') /* Skip already settled `i' */
  {
    char * name = ParameterName + ParameterN*NameLength1;
    if(ParameterN == ParameterMaxN)
      Error(E_PARAMETER_MAX_N);
    do
      *(name++) = *str;
    while(*(str++));
    ++ParameterN;
  }
}
/*=GetRelation===============================================================
 Read single relation from string, reduce and set in array
*/
void GetRelation(char * str)
{
  uint a;
  if(str[0] != '\0' && (a = GetLieSum(&str)) != NIL)
  {
    int lmonpos, pos, i, l;
    if(RelationN + 1 == RelationSize)
      Error(E_RELATION_SIZE);
    (*NormalizeRelation)(a);
    lmonpos = LIE_TERM_MONOMIAL(a);
    i = LIE_MONOMIAL_ORDER(lmonpos);
    l = FindNewPositionInRelation(i);
    LIE_MONOMIAL_INDEX(lmonpos) = ~l;  /* Set int of relation in LieMonomial */
    while(++i < LieMonomialN)     /* Shift int's of relations in LieMonomial */
      if(LIE_MONOMIAL_IS_LEADING(pos = LIE_MONOMIAL_POSITION(i)))
        --LIE_MONOMIAL_INDEX(pos);
    for(i = RelationN; i > l; i--) /* Make room moving Relation structures */
      Relation[i] = Relation[i-1];
    RELATION_MIN_GENERATOR(l) = LIE_MONOMIAL_IS_GENERATOR(lmonpos) ?
            GeneratorN /* Don't differentiate leading generator */ : 0;
    RELATION_LIE_SUM(l) = a;
    RELATION_TO_BE_SUBSTITUTED(l) = (byte)(l < RelationN);
    ++RelationN;
#if defined(SPACE_STATISTICS)
    if(RelationN > MaxNRelation)
      MaxNRelation = RelationN;
#endif
    ReduceRelations(l);
  }
}
/*=GetScalarSum==================================================
*/
uint GetScalarSum(char **pstr)
{
  uint a, term;
  int is_par, is_negative;
  if(**pstr == '(')
  {
    is_par = YES;
    ++*pstr;
    SkipSpaces(pstr);
  }
  else
    is_par = NO;
  if(**pstr == '-')
  {
    is_negative = YES;
    ++*pstr;
    SkipSpaces(pstr);
  }
  else
    is_negative = NO;
  a = GetScalarTerm(pstr);
  if(is_negative)
    ScalarSumMinus(a);
  while(**pstr == '+' || **pstr == '-')
  {
    is_negative = (**pstr == '+') ? NO : YES;
    ++*pstr;
    SkipSpaces(pstr);
    term = GetScalarTerm(pstr);
    if(is_negative)
      ScalarSumMinus(term);
    a = ScalarSumAddition(a, term);
  }
  if(is_par)
  {
    if(**pstr != ')')
      Error(E_NO_R_PARENTHESIS);
    ++*pstr;
    SkipSpaces(pstr);
  }
  return a;
}
/*=GetScalarTerm=========================================================
 Read unsigned scalar term in Parametric regime
*/
uint GetScalarTerm(char **pstr)
{
  BIGINT nums, numh;
  int i, change_sign;
  uint m, f, a = NodeSTNew();

  /* Read numerical coefficient */

  INTEGER_STACK_NEW(nums, InputIntegerSize);
  GetInteger(nums, pstr);
  if(nums[0] == 0)
    nums[0] = nums[1] = 1;
  INTEGER_HEAP_COPY(numh, nums, i);
  SCALAR_TERM_NUMERATOR(a) = numh;

  /* Read scalar monomial */

  SkipSpaces(pstr);
  m = NIL;
  while(isalpha(**pstr) &&
        (i=FindNameInTable(*pstr,ParameterName,ParameterN)) < ParameterN)
  {
    f = NodeSFNew();
    SCALAR_FACTOR_PARAMETER(f) = i;
    SkipName(pstr);
    SkipSpaces(pstr);
    if(**pstr == '^')
    {
      ++*pstr;
      SkipSpaces(pstr);
      if(isdigit(**pstr))
      {
        i = (byte)GetUInteger(pstr);                   /* Read degree */
        SkipSpaces(pstr);
      }
      else
        Error(E_NO_GENERAL_POWER);
    }
    else
      i = 1;
    SCALAR_FACTOR_DEGREE(f) = i;                /* Degree of parameter */
    m = ScalarMonomialMultiplication(&change_sign, m, f);
    if(change_sign)
      SCALAR_TERM_MINUS(a);
  }
  SCALAR_TERM_MONOMIAL(a) = m;
  return a;
}
/*=GetWeight=================================================
 Read single generator weight from description string
*/
void GetWeight(char *str)
{
  if(GeneratorN == GeneratorMaxN)
    Error(E_TOO_MUCH_INPUT_WEIGHTS);
  if(!isdigit(*str))
    Error(E_NON_NUM_INPUT_WEIGHT);
  LIE_MONOMIAL_WEIGHT(GeneratorN++) = (byte)GetUInteger(&str);
}
/*=KeyBoardBytesToString=====================
  Returns 0 for the last input (ended by '.')
  and last position otherwise
*/
int KeyBoardBytesToString(char *str)
{
  char c;
  int inspace = YES, i = -1;
  while((c = fgetc(stdin)) != '\n')
  {
    if(c == ' ')
    {
      if(inspace)
        continue;     /* Skip extra blanks */
      inspace = YES;
    }
    else
      inspace = NO;
    PutCharacter(c);
    if(c == '\b')
    {
#if !defined(GAP)
       fseek(SessionFile, -2, SEEK_CUR);
#endif
       --i;
    }
    else
    {
      if(++i == InputStringSize)
        Error(E_INPUT_STRING_SIZE);
      str[i] = c;
    }
  }
  putchar('\n');
  if(i < 0)        /* Empty input */
    return 0;
  if(str[i] != ';' && str[i] != '.')
  {          /* Add semicolon if necessary */
    if(++i == InputStringSize)
      Error(E_INPUT_STRING_SIZE);
    str[i] = ';';
  }
  if(++i == InputStringSize)
    Error(E_INPUT_STRING_SIZE);
  str[i] = '\n';         /* Go to new line */
  return i;
}
/*=KeyBoardStringToFile==========================================
 Add string from keyboard to file between prefix and semicolon
*/
int KeyBoardStringToFile(int i_m, char * prefix, char * str, FILE *file)
{
  short itop;
  PutMessage(i_m);
  itop = KeyBoardBytesToString(str);
  if(itop)
  {
    short i = 0;
    while(*prefix)
      fputc(*(prefix++), file); /* Copy prefix to file */
    while(i <= itop)
      fputc(str[i++], file); /* Copy entered string to file */
  }
  return itop;
}
/*=ReadAndProcessStringsFromFile==================================
 Read array of strings separated by `sep' and ended with `end'
 Remove comments and unnecessary spaces
 Add ending '\0'   Process strings by (*proc_func)
*/
void ReadAndProcessStringsFromFile(void (*proc_func)(char *str), FILE *inf,
                                                     char sep, char end)
{
  char *str, *wstr;
  char line_break;
  short c;
  int i;
  str = (char *)alloca(InputStringSize);
  if(str == NULL)
    Error(E_A_STACK_INPUT_STRING);
  line_break = (sep == ' ') ? '\n' : '\0';
  wstr = str;
  i = 0;       /* Count number of characters */
  while(YES)
  {
    if((c = fgetc(inf)) == sep || c == line_break || c == end)
    {
      if(wstr != str && wstr[-1] == ' ')
        wstr[-1] = '\0';         /* Kill ending blank */
      else
        *wstr = '\0';
      (*proc_func)(str);         /* Process string */
      if(c == end)               /* Last */
        break;
      wstr = str;                /* Intermediate */
      i = 0;                     /* Count number of characters */
      while((c = SkipSpacesInFile(inf)) == LEFT_COMMENT)
        SkipCommentInFile(inf);
      if(c == EOF)
        break;
      continue;
    }
    if(isspace(c))
    {
      if(++i == InputStringSize)
        Error(E_INPUT_STRING_SIZE);
      *(wstr++) = ' ';
      SkipSpacesInFile(inf);
      continue;
    }
    if(c == LEFT_COMMENT)
    {
      if(*wstr != ' ')
      {
        if(++i == InputStringSize)
          Error(E_INPUT_STRING_SIZE);
        *(wstr++) = ' ';
      }
      do
        SkipCommentInFile(inf);
      while((c = SkipSpacesInFile(inf)) == LEFT_COMMENT);
      continue;
    }
    if(++i == InputStringSize)
      Error(E_INPUT_STRING_SIZE);
    *(wstr++) = (char)c;
  }
}
/*=ReadBooleanFromFile========================
  Read boolean constant from file
*/
int ReadBooleanFromFile(FILE *file)
{
  short c;
  int boolvar;
  c = fgetc(file);
  switch(c)
  {
    case 'Y': case 'y':
    boolvar = YES;
    break;
    case 'N': case 'n':
    boolvar = NO;
    break;
  }
  while(!isspace(c = fgetc(file)) && c != EOF)
    ;
  ungetc(c, file);
  return boolvar;
}
/*=ReadCaseFromFile=====================================
*/
int ReadCaseFromFile(FILE * file, char * case_str[], int n_cases)
{
  char file_str[CASE_STRING_SIZE];
  char * w_file_str, *w_case_str;
  int i_case;
  short c;
  while(SkipSpacesInFile(file) == LEFT_COMMENT)
    SkipCommentInFile(file);

  /* Read string ending with : from file */

  w_file_str = file_str;
  do
  {
    if((c = fgetc(file)) == EOF)
    {
      i_case = c;
      goto out;
    }
    if(c == ':')
      c = '\0';
    *(w_file_str++) = (char)c;
  }while(c);

  /* Compare strings */

  i_case = 0;
  do
  {
    w_file_str = file_str;
    w_case_str = case_str[i_case];
    do
    {
      if(*w_case_str == '\0')  /* Case is found */
      {
        while(SkipSpacesInFile(file) == LEFT_COMMENT)
          SkipCommentInFile(file);
        goto out;
      }
      if(*w_file_str == '\0')
        break;
    }while(*(w_file_str++) == *(w_case_str++));
  }while(++i_case < n_cases);
  out:
  return i_case;
}
/*=ReadDecimalFromFile===================
  Read unsigned decimal integer from file
*/
uint ReadDecimalFromFile(FILE *file)
{
  short c;
  uint i = 0;
  while(isdigit(c = fgetc(file)))
    i = i*10 + c - '0';
  ungetc(c, file);
  return i;
}
/*=ReadStringFromFile====================
*/
short ReadStringFromFile(char * str, FILE *file)
{
  short c;
  while((c = fgetc(file)) != '\n')
  {
    if(c == LEFT_COMMENT)
    {
      ungetc(c, file);
      break;
    }
    if(c == EOF)
      break;
    if(!isspace(c))
      *(str++) = (char)c;
  }
  return c;
}
/*=SkipCommentInFile=======================
*/
short SkipCommentInFile(FILE *file)
{
  short c;
  while((c = fgetc(file)) != RIGHT_COMMENT)
    if(c == EOF)
      Error(E_UNEXPECTED_EOF);
  return c;
}
/*=SkipName===============================================
*/
void SkipName(char **pstr)
{
  while(isalnum(**pstr) || **pstr == SUBSCRIPT_INPUT_SIGN)
    ++*pstr;
}
/*=SkipSpaces==================
 Skip spaces to right in string
*/
void SkipSpaces(char **pstr)
{
  while(isspace(**pstr))
    ++*pstr;
}
/*=SkipSpacesInFile==============
  Returns first non-space symbol
*/
short SkipSpacesInFile(FILE *file)
{
  short c;
  while(isspace(c = fgetc(file)))
    ;
  ungetc(c, file);
  return c;
}

/*_6_9          Output functions========================================*/

/*=AddSymbolToOutLine!===============
 Add symbol to output line OutLine
*/
void AddSymbolToOutLine(char c, int position)
{
  if(position >= OutLineSize)
    Error(E_OUT_LINE_SIZE);
  OutLine[position] = c;
}
/*=InLineLevel==============================
 Install level in OutLine
*/
void InLineLevel(int level)
{
  if(level != CurrentLevel)
  {
    AddSymbolToOutLine(LEVEL, ++PosOutLine);
    AddSymbolToOutLine((char)level, ++PosOutLine);
    CurrentLevel = level;
    if(level > MaxLevel)
      MaxLevel = level;
    else if(level < MinLevel)
      MinLevel = level;
  }
}
/*=InLineNumberInBrackets=====
 (With 2 blanks after)
*/
void InLineNumberInBrackets(uint n)
{
  InLineSymbol('(');
  InLineString(UToString(n));
  InLineString(")  ");
}
/*=InLineString!================================
 Add string to output line OutLine for 2D output
*/
void InLineString(char *str)
{
  while(*str)
    InLineSymbol(*(str++));
}
/*=InLineSubscript================
 Add symbolic subscript to OutLine
*/
void InLineSubscript(char * s)
{
  int level = CurrentLevel;
  InLineLevel(level - 1);
  InLineString(s);
  InLineLevel(level);
}
/*=InLineSymbol!================================
 Add symbol to output line OutLine for 2D output
*/
void InLineSymbol(char symbol)
{
  AddSymbolToOutLine(symbol, ++PosOutLine);
  ++LastItemEnd;
}
/*=InLineTableName===================
 Possibly subscripted
*/
void InLineTableName(char * name)
{
  PreviousEnd = LastItemEnd;
  while(*name)
  {
    if(*name == SUBSCRIPT_INPUT_SIGN)
    {
      InLineSubscript(++name);
      break;
    }
    InLineSymbol(*(name++));
  }
}
/*=UToString!====================================
 Transform unsigned long number to decimal string
*/
char * UToString(uint n)
{                            /*12345678910*/
  static char decimal_string[] = "4294967295";
  char * first = decimal_string + 10;
  do
  {
    *--first = '0' + (char)(n % 10);
    n /= 10;
  }while(n);
  return first;
}
/*=PutBasis============================================
 Set index numbers and print basis elements
*/
void PutBasis(void)
{
  int pos, ord;
  uint  i = 0;      /* Index number of basis element */
  TIME_OFF;
  if(BasisElementsPut)
    PutMessage(H_BASIS_ELEMENTS);
  for(ord = 0; ord < LieMonomialN; ord++)
  {
    pos = LIE_MONOMIAL_POSITION(ord);
    if(LIE_MONOMIAL_IS_BASIS(pos))
    {
      LIE_MONOMIAL_INDEX(pos) = ++i;
      if(BasisElementsPut)
      {
        PutStart();
        InLineNumberInBrackets(i);
        PutLieBasisElement(pos);
        InLineString(" = ");
        IN_LINE_MARGIN;
        (*PutLieMonomial)(pos);
        PutEnd();
      }
    }
  }
  if(BasisElementsPut)
    if(IncompletedBasis)
      PutDots();
  TIME_ON;
}
#if defined(GAP)
/*=PutBasisGAP=============================================

 Print basis elements into GAP file
*/
void PutBasisGAP(void)
{
  int pos, ord, tord;
  tord = LieMonomialN - 1;
  while(YES)
  {
    if(LIE_MONOMIAL_IS_BASIS(LIE_MONOMIAL_POSITION(tord)))
      break;
    --tord;
  }
  PutStringStandard(GAPBasisName);
  PutStringStandard(":=[\n");
  for(ord = 0; ord <= tord; ord++)
  {
    pos = LIE_MONOMIAL_POSITION(ord);
    if(LIE_MONOMIAL_IS_BASIS(pos))
    {
      PosOutLine = -1;
      PutLieMonomialGAP(pos);
      if(ord < tord)
        PutStringGAP(",\n");
    }
  }
  PutStringStandard("\n];\n");
}
#endif
/*=PutBlock==============================================================
 Print block of 2D output
*/
void PutBlock(void)
{
  if(!PrintEnd && (LastItemEnd > LineLength || PreviousEnd==LastItemEnd))
    PrintEnd = PreviousEnd;
  if(PrintEnd && CurrentLevel <= MAIN_LEVEL)
  {
    int xp, yp = MaxLevel, i, prlvl;
    while(yp >= MinLevel)
    {
      for(xp = 1; xp <= Margin; xp++)
        PutCharacter(' ');
      i = 0;
      while(xp <= PrintEnd)
      {
        switch(OutLine[i])
        {
          case LEVEL:
            prlvl = OutLine[++i];
            break;
          case MARGIN:
            NewMargin = xp - 1;
            break;
          default:
            PutCharacter((char)((prlvl == yp) ? OutLine[i] : ' '));
            xp++;
        }
        i++;
      }
      PutCharacter('\n');
      yp--;
    }
    PutCharacter('\n'); /* To next line */
    Margin = NewMargin;
    LastItemEnd += Margin - PrintEnd;
    PrintEnd = 0;
    xp = i;
    while(OutLine[xp] != LEVEL)
      xp--;
    *OutLine = LEVEL;
    OutLine[1] = MaxLevel = MinLevel = OutLine[xp + 1];
    xp = PosOutLine;
    PosOutLine = 1;
    while(i <= xp)
    {
      OutLine[++PosOutLine] = OutLine[i];
      if(OutLine[i++] == LEVEL)
      {
        if(OutLine[i] > MaxLevel)
          MaxLevel = OutLine[i];
        else if(OutLine[i] < MinLevel)
          MinLevel = OutLine[i];
      }
    }
  }
}
/*=PutCharacter!=============
   Echoed output of character
*/
void PutCharacter(char c)
{
#if defined(ECHO_TO_SCREEN)
  putc(c, stdout);
#endif
#if !defined(GAP)
  putc(c, SessionFile);
#endif
}
#if defined(GAP)
/*=PutCharacterGAP!======================
   Echoed output of character in GAP file
*/
void PutCharacterGAP(char c)
{
#if defined(ECHO_TO_SCREEN)
  putc(c, stdout);
#endif
/*  putc(c, SessionFile);  */
  PosOutLine++;
}
#endif
/*=PutCoefficientTable=================================
 Print list of non-zero coefficients
*/
void PutCoefficientTable(void)
{
  int i, j = 0, first = YES;
  TIME_OFF;
  if(CoeffParamTable != NULL)
  {
    for(i = FIRST_GENUINE_PARAMETER; i < ParameterN; i++)
      if(CoeffParamTable[i])
      {
        if(first)
        {
          first = NO;
          PutMessage(H_NON_ZERO_COEFFICIENTS);
        }
        PutStart();
        InLineNumberInBrackets(++j);
        InLineTableName(ParameterName + i*NameLength1);
        PutEnd();
      }
    free(CoeffParamTable);
  }
  if(CoeffSumTableN)
  {
    if(first)
    {
      first = NO;
      PutMessage(H_NON_ZERO_COEFFICIENTS);
    }
    for(i = 0; i < CoeffSumTableN; i++)
    {
      PutStart();
      InLineNumberInBrackets(++j);
      IN_LINE_MARGIN;
      PutScalarSum(CoeffSumTable[i]);
      PutEnd();
      ScalarSumKill(CoeffSumTable[i]);
    }
    free(CoeffSumTable);
  }
  TIME_ON;
}
/*=PutCommutators========================================================
 Compute and print commutators of basis elements [E(O),E(O)] -> E(O)
*/
void PutCommutators(void)
{
  int oi, oj, i, j, was_comm, set_minus;
  uint icomm, rpart;
  byte weight_j;
  TIME_OFF;
  PutMessage(H_COMMUTATORS);
  TIME_ON;
  icomm = 0;
  for(oj = 0; oj < LieMonomialN; oj++)
    if(LIE_MONOMIAL_IS_BASIS(j = LIE_MONOMIAL_POSITION(oj)))
    {
      weight_j = LIE_MONOMIAL_WEIGHT(j);
      was_comm = NO;
      set_minus = LIE_MONOMIAL_IS_EVEN(j);
      for(oi = 0; oi < oj ||
                  (oi == oj &&
                   LIE_MONOMIAL_IS_ODD(LIE_MONOMIAL_POSITION(oi))); oi++)
        if(weight_j + LIE_MONOMIAL_WEIGHT(i = LIE_MONOMIAL_POSITION(oi))
           <= LimitingWeight)
        {
          if(LIE_MONOMIAL_IS_BASIS(i))
          {
            if(LIE_MONOMIAL_IS_EVEN(i))
              set_minus = YES;
            rpart = (*PairMonomialMonomial)(j, i);
            if(rpart != NIL)
            {
              TIME_OFF;
              PutStart();
              InLineNumberInBrackets(++icomm);
              InLineSymbol('[');
              PutLieBasisElement(i);
              InLineSymbol(',');
              PutLieBasisElement(j);
              InLineString("] = ");
              IN_LINE_MARGIN;
              if(set_minus)
                (*LieSumMinus)(rpart);
              PutLieSum(PutLieBasisElement, rpart);
              LieSumKill(rpart);
              PutEnd();
              TIME_ON;
            }
          }
          was_comm = YES;
        }
        else
        {
          if(was_comm)
          {
            while(oi < oj ||
                 (oi == oj &&
                  LIE_MONOMIAL_IS_ODD(LIE_MONOMIAL_POSITION(oi))))
            {
              if(LIE_MONOMIAL_IS_BASIS(LIE_MONOMIAL_POSITION(oi)))
                ++icomm;
              ++oi;
            }
            TIME_OFF;
            PutDots();
            PutCharacter('\n');
            TIME_ON;
          }
          break;
        }
    }
  TIME_OFF;
  if(icomm == 0)
    PutMessage(H_NO_PUT_COMMUTATORS);
  TIME_ON;
}
#if defined(GAP)
/*=PutCommutatorsGAP==================================
 Transform commutators of ordinary finite-dimensional
 parameter-free Lie algebra to GAP input format
*/
void PutCommutatorsGAP(void)
{
  int oi, oj, otop, i, j;
  BIGINT num;
  uint rpart, a, b;
  otop = LieMonomialN - 1;
  while(YES)
  {
    if(LIE_MONOMIAL_IS_BASIS(LIE_MONOMIAL_POSITION(otop)))
      break;
    --otop;
  }
  PutStringStandard(GAPAlgebraName);
  PutStringStandard(":=[\n");   /* Begin main list */
  for(oj = 0; oj <= otop; oj++)
    if(LIE_MONOMIAL_IS_BASIS(j = LIE_MONOMIAL_POSITION(oj)))
    {
      PosOutLine = -1;
      PutCharacterGAP('[');
      for(oi = 0; oi <= otop; oi++)
        if(LIE_MONOMIAL_IS_BASIS(i = LIE_MONOMIAL_POSITION(oi)))
        {
          if(oi == oj)
            rpart = NIL;
          else if(oj < oi)
          {
            rpart = (*PairMonomialMonomial)(i, j);
            (*LieSumMinus)(rpart);
          }
          else
            rpart = (*PairMonomialMonomial)(j, i);
          if(rpart != NIL)
          {
            a = LIE_TERM_R(rpart);     /* Invert list */
            LIE_TERM_R(rpart) = NIL;
            while(a != NIL)
            {
              b = a;
              a = LIE_TERM_R(a);
              LIE_TERM_R(b) = rpart;
              rpart = b;
            }
            PutStringGAP("[[");
            a = rpart;      /* Put indices of basis elements */
            while(YES)
            {
              PutStringGAP(UToString(LIE_MONOMIAL_INDEX(
                                       LIE_TERM_MONOMIAL(a))));
              a = LIE_TERM_R(a);
              if(a == NIL)
                break;
              PutStringGAP(",");
            }
            PutStringGAP("],[");
            a = rpart;                  /* Put coefficients */
            while(YES)
            {
              num = LIE_TERM_NUMERATOR_INTEGER(a);
              if(INTEGER_IS_NEGATIVE(num))
                PutStringGAP("-");
              PutIntegerUnsignedGAP(num);
              if((num = LIE_TERM_DENOMINATOR_INTEGER(a)) != NULL)
              {
                PutStringGAP("/");
                PutIntegerUnsignedGAP(num);
              }
              a = LIE_TERM_R(a);
              if(a == NIL)
                break;
              PutStringGAP(",");
            }
            PutStringGAP("]]");
            LieSumKill(rpart);
          }
          else
            PutStringGAP("[[],[]]");
          if(oi != otop)
            PutStringGAP(",");
        }
      PutStringGAP("],");
      PutCharacter('\n');
    }
  PutStringStandard("-1,0];\n");  /* Close main list */
}
#endif
/*=PutDegree==========================
 Print (positive) integer degree
*/
void PutDegree(uint deg)
{
  if(deg != 1)
  {
    int level = CurrentLevel;
    InLineLevel(level + 1);
    InLineString(UToString((uint)deg));
    InLineLevel(level);
    PutBlock();
  }
}
/*=PutDimensions================================================
  Print dimensions of homogeneous components of algebra
*/
void PutDimensions(void)
{
  int next, ord, pos;
  byte   curwt;
  uint   dim;
  TIME_OFF;
  PutMessage(H_HILBERT_SERIES);
  PutStart();
  InLineString("H(t) = ");
  IN_LINE_MARGIN;
  next = NO;
  curwt = LIE_MONOMIAL_WEIGHT(0);     /* Start initial weight */
  dim = 0;
  for(ord = 0; ord < LieMonomialN; ord++)
  {
    pos = LIE_MONOMIAL_POSITION(ord);
    if(LIE_MONOMIAL_IS_BASIS(pos))
    {
      if(LIE_MONOMIAL_WEIGHT(pos) != curwt)
      {
        if(dim != 0)
        {
          if(next)
            PutString(" + ");
          else
            next = YES;
          if(dim != 1)
          {
            PutString(UToString(dim));
            PutSymbol(' ');
          }
          PutSymbol('t');
          PutDegree(curwt);
        }
        curwt = LIE_MONOMIAL_WEIGHT(pos); /* Start new weight */
        dim = 1;
      }
      else
        dim++;
    }
  }
  if(dim != 0)   /* Print last element */
  {
    if(next)
      PutString(" + ");
    if(dim != 1)
    {
      PutString(UToString(dim));
      PutSymbol(' ');
    }
    PutSymbol('t');
    PutDegree(curwt);
  }
  if(IncompletedBasis)
    PutString(" + ...");
  PutEnd();
  TIME_ON;
}
/*=PutDots=============================
 Print vertical dots
*/
void PutDots(void)
{
#if defined(ECHO_TO_SCREEN)
  printf(" .\n .\n .\n");
#endif
#if !defined(GAP)
  fprintf(SessionFile, " .\n .\n .\n");
#endif
}
/*=PutEnd=======================
   Print last block of 2D output
*/
void PutEnd(void)
{
  PreviousEnd = LastItemEnd;
  PutBlock();
}
/*=PutFormattedU=================
*/
void PutFormattedU(char * format, uint i)
{
#if !defined(GAP)
  fprintf(SessionFile, format, i);
#endif
#if defined(ECHO_TO_SCREEN)
  printf(format, i);
#endif
}
/*=PutLieBareTerm=======================================================
 put_lie_mon == PutLieMonomial     -> in terms of generators,
 put_lie_mon == PutLieBasisElement -> in terms of basis elements.
*/
void PutLieBareTerm(void (*put_lie_mon)(int a), uint a)
{
  int put_mult_sign = NO;
  if(IsParametric)
  {
    uint num = LIE_TERM_NUMERATOR_SCALAR_SUM(a),
      den = LIE_TERM_DENOMINATOR_SCALAR_SUM(a);

    /* Put numerator */

    if(SCALAR_TERM_R(num) != NIL)
    {
      PutSymbol('(');              /* Sum */
      PutScalarSum(num);
      PutSymbol(')');
      put_mult_sign = YES;
    }
    else if(SCALAR_TERM_MONOMIAL(num) != NIL ||
            INTEGER_IS_NOT_UNIT_ABS(SCALAR_TERM_NUMERATOR(num)) ||
            den != NIL)
    {
      PutScalarBareTerm(num);     /* Single term */
      put_mult_sign = YES;
    }

    /* Put denominator */

    if(den != NIL)
    {
      PutSymbol('/');
      if(SCALAR_TERM_R(den) == NIL &&  /* Single factor denominator */
         (SCALAR_TERM_MONOMIAL(den) == NIL ||
          INTEGER_IS_UNIT(SCALAR_TERM_NUMERATOR(den))))
        PutScalarBareTerm(den);
      else
      {
        PutSymbol('(');            /* Sum or multifactor denominator */
        PutScalarSum(den);
        PutSymbol(')');
      }
    }
  }
  else
  {
    BIGINT num = LIE_TERM_NUMERATOR_INTEGER(a),
        den = LIE_TERM_DENOMINATOR_INTEGER(a);
    if(den != NULL)
    {
      PutIntegerUnsigned(num);
      PutSymbol('/');
      PutIntegerUnsigned(den);
      put_mult_sign = YES;
    }
    else if(INTEGER_IS_NOT_UNIT_ABS(num))
    {
      PutIntegerUnsigned(num);
      put_mult_sign = YES;
    }
  }
  if(put_mult_sign)
    PutSymbol(' ');
  (*put_lie_mon)(LIE_TERM_MONOMIAL(a));
}
/*=PutLieBasisElement======================================
 Print Lie monomial in form E  or O
                             i     i
*/
void PutLieBasisElement(int pos)
{
  InLineSymbol((char)(LIE_MONOMIAL_PARITY(pos) ? BasisSymbolOdd :
                                          BasisSymbolEven));
  InLineSubscript(UToString(LIE_MONOMIAL_INDEX(pos)));
  PutBlock();
}
/*=PutLieMonomialLeftNormed=====================================
 Put Lie monomial in terms of generators in left normed notation
*/
void PutLieMonomialLeftNormed(int pos)
{
  if(LIE_MONOMIAL_IS_GENERATOR(pos))
  {
    InLineTableName(GeneratorName + pos*NameLength1);
    PutBlock();
  }
  else
  {
    uint deg = 1;
    int posi = LIE_MONOMIAL_LEFT(pos),
      posj = LIE_MONOMIAL_RIGHT(pos), posw;
    while(LIE_MONOMIAL_IS_COMMUTATOR(posi))
    {
      pos = posi;
      posi = LIE_MONOMIAL_LEFT(pos);
      posw = LIE_MONOMIAL_RIGHT(pos);
      if(posj == posw)
        ++deg;
      else
      {
        posi = pos;
        break;
      }
    }
    if(posi != posj)
        PutLieMonomialLeftNormed(posi);
    else
      ++deg;
    if(LIE_MONOMIAL_IS_COMMUTATOR(posj))
    {
      PutSymbol('(');
      PutLieMonomialLeftNormed(posj);
      PutSymbol(')');
    }
    else
      PutLieMonomialLeftNormed(posj);
    PutDegree(deg);
  }
}
/*=PutLieMonomialStandard============================================
 Put Lie monomial in terms of generators in standard bracket notation
*/
void PutLieMonomialStandard(int pos)
{
  if(LIE_MONOMIAL_IS_GENERATOR(pos))
  {
    InLineTableName(GeneratorName + pos*NameLength1);
    PutBlock();
  }
  else
  {
    PutSymbol('[');
    PutLieMonomialStandard(LIE_MONOMIAL_LEFT(pos));
    PutSymbol(',');
    PutLieMonomialStandard(LIE_MONOMIAL_RIGHT(pos));
    PutSymbol(']');
  }
}
#if defined(GAP)
/*=PutLieMonomialGAP============================================
 Put Lie monomial in terms of generators in GAP bracket notation
*/
void PutLieMonomialGAP(int pos)
{
  if(LIE_MONOMIAL_IS_GENERATOR(pos))
    PutStringGAP(UToString(pos+1));
  else
  {
    PutStringGAP("[");
    PutLieMonomialGAP(LIE_MONOMIAL_LEFT(pos));
    PutStringGAP(",");
    PutLieMonomialGAP(LIE_MONOMIAL_RIGHT(pos));
    PutStringGAP("]");
  }
}
#endif
/*=PutLieSum=====================================================
 put_lie_mon == PutLieMonomial     -> in terms of generators,
 put_lie_mon == PutLieBasisElement -> in terms of basis elements
*/
void PutLieSum(void (*put_lie_mon)(int a), uint a)
{
  if(a == NIL)
    PutSymbol('0');
  else
  {
    uint na;
    int is_negative = NO;
    if(IsParametric)
    {
      na = LIE_TERM_NUMERATOR_SCALAR_SUM(a);
      if(SCALAR_TERM_R(na) == NIL &&
         INTEGER_IS_NEGATIVE(SCALAR_TERM_NUMERATOR(na)))
        is_negative = YES;
    }
    else if(INTEGER_IS_NEGATIVE(LIE_TERM_NUMERATOR_INTEGER(a)))
      is_negative = YES;
    if(is_negative)
      PutString("- ");
    PutLieBareTerm(put_lie_mon, a);
    while((a = LIE_TERM_R(a)) != NIL)
    {
      is_negative = NO;
      if(IsParametric)
      {
        na = LIE_TERM_NUMERATOR_SCALAR_SUM(a);
        if(SCALAR_TERM_R(na) == NIL &&
           INTEGER_IS_NEGATIVE(SCALAR_TERM_NUMERATOR(na)))
          is_negative = YES;
      }
      else if(INTEGER_IS_NEGATIVE(LIE_TERM_NUMERATOR_INTEGER(a)))
        is_negative = YES;
      PutString(is_negative ? " - " : " + ");
      PutLieBareTerm(put_lie_mon, a);
    }
  }
}
/*=PutMessage!=========================
 Put message from MessageFile
*/
void PutMessage(int i_message)
{
#if !defined(GAP)
  static int current_message;
  short c;
  if(i_message < current_message)
  {
    rewind(MessageFile);
    current_message = 0;
  }
  while(i_message > current_message)
    if((c = fgetc(MessageFile)) == '$')
      ++current_message;
    else if(c == EOF)
      Error(E_MESSAGE);
  while(YES)
  {
    c = fgetc(MessageFile);
    if(c == '$')
    {
      ++current_message;
      return;
    }
    PutCharacter((char)c);
  }
#endif
}
/*=PutRelations====================================
 Print list of relations
*/
void PutRelations(int msg)
{
  int i;
  TIME_OFF;
  PutMessage(msg);
  for(i = 0; i < RelationN; i++)
  {
    PutStart();
    InLineNumberInBrackets(i+1);
    IN_LINE_MARGIN;
    PutLieSum(PutLieMonomial, RELATION_LIE_SUM(i));
    InLineString(" = 0");
    PutEnd();
  }
  if(IncompletedRelations)
    PutDots();
  TIME_ON;
}
#if defined(GAP)
/*=PutRelationsGAP=================================
*/
void PutRelationsGAP(void)
{
  int i;
  BIGINT num;
  uint a;
  PutStringStandard(GAPRelationsName);
  PutStringStandard(":=[\n");   /* Begin main list */
  for(i = 0; i < RelationN; i++)
  {
    a = RELATION_LIE_SUM(i);
    PosOutLine = -1;
    PutCharacterGAP('[');
    while(YES)
    {
      PutLieMonomialGAP(LIE_TERM_MONOMIAL(a));
      PutStringGAP(",");
      num = LIE_TERM_NUMERATOR_INTEGER(a);
      if(INTEGER_IS_NEGATIVE(num))
        PutStringGAP("-");
      PutIntegerUnsignedGAP(num);
      if((a = LIE_TERM_R(a)) == NIL)
        break;
      PutStringGAP(",");
    }
    PutStringGAP("]");
    if(i < RelationN - 1)
      PutStringGAP(",\n");
  }
  PutStringStandard("\n];\n");
}
#endif
/*=PutScalarBareTerm===============================
 Intermediate print of unsigned scalar term
*/
void PutScalarBareTerm(uint a)
{
  int put_1 = YES;

  /* Put integer coefficient */

  if(INTEGER_IS_NOT_UNIT_ABS(SCALAR_TERM_NUMERATOR(a)))
  {
    PutIntegerUnsigned(SCALAR_TERM_NUMERATOR(a));
    put_1 = NO;
  }

  /* Put scalar monomial */

  a = SCALAR_TERM_MONOMIAL(a); /* uint - uint mixing */
  if(a != NIL)
  {
    if(put_1 == NO)
        PutSymbol(' ');
    PutScalarFactor(a);
    while((a = SCALAR_FACTOR_R(a)) != NIL)
    {
      PutSymbol(' ');
      PutScalarFactor(a);
    }
    put_1 = NO;
  }
  if(put_1)         /* Nothing's been put before */
    PutSymbol('1');
}
/*=PutIntegerUnsigned=============================================
 Print big integer
*/
void PutIntegerUnsigned(BIGINT bn)
{
  BIGINT bnw;
  char * decstr;
  LIMB res;
  uint lw;
  int i,
    n = INTEGER_N_LIMBS(bn),
    nw = n;
  bnw = (BIGINT)alloca(sizeof(LIMB)*n);
  if(bnw == NULL)
    Error(E_A_STACK_INTEGER); /* LIMB contains 5 decimal digits */
  decstr = (char *)alloca(5*n + 1);
  if(decstr == NULL)
    Error(E_A_STACK_INTEGER_DECIMAL_STRING);
  decstr += 5*n;            /* Go to last byte of string */
  for(i = 0; i < n; i++)
    bnw[i] = *(++bn);       /* Copy body of big number */
  /* Transform big number array to decimal string */
  decstr[0] = '\0';
  do
  { /* Divide bnw number in array form by 10 on spot */
    res = 0;
    if(nw)       /* Otherwise 0/n -> 0, 0 mod n -> 0 */
    {
      i = nw;
      do
      {
        lw = (uint)res * BASE_LIMB + (uint)bnw[--i];
        res = (LIMB)(lw % 10);
        bnw[i] = (LIMB)(lw / 10);
      }while(i);
      if(bnw[nw-1] == 0)
        --nw;
    }
    *--decstr = '0' + res;
  }while(nw);
  if(n < 3)
    PutString(decstr);      /* Don't cut short number */
  else
    do
      PutSymbol(*decstr);
    while(*(++decstr));
}
#if defined(GAP)
/*=PutIntegerUnsignedGAP==========================================
 Print big integer into GAP file
*/
void PutIntegerUnsignedGAP(BIGINT bn)
{
  BIGINT bnw;
  char * decstr;
  LIMB res;
  uint lw;
  int i,
    n = INTEGER_N_LIMBS(bn),
    nw = n;
  bnw = (BIGINT)alloca(sizeof(LIMB)*n);
  if(bnw == NULL)
    Error(E_A_STACK_INTEGER); /* LIMB contains 5 decimal digits */
  decstr = (char *)alloca(5*n + 1);
  if(decstr == NULL)
    Error(E_A_STACK_INTEGER_DECIMAL_STRING);
  decstr += 5*n;            /* Go to last byte of string */
  for(i = 0; i < n; i++)
    bnw[i] = *(++bn);       /* Copy body of big number */
  /* Transform big number array to decimal string */
  decstr[0] = '\0';
  do
  { /* Divide bnw number in array form by 10 on spot */
    res = 0;
    if(nw)       /* Otherwise 0/n -> 0, 0 mod n -> 0 */
    {
      i = nw;
      do
      {
        lw = (uint)res * BASE_LIMB + (uint)bnw[--i];
        res = (LIMB)(lw % 10);
        bnw[i] = (LIMB)(lw / 10);
      }while(i);
      if(bnw[nw-1] == 0)
        --nw;
    }
    *--decstr = '0' + res;
  }while(nw);
  PutStringGAP(decstr);
}
#endif
/*=PutScalarFactor========================================================
*/
void PutScalarFactor(uint a)
{
  InLineTableName(ParameterName + SCALAR_FACTOR_PARAMETER(a)*NameLength1);
  PutDegree(SCALAR_FACTOR_DEGREE(a));
  PutBlock();
}
/*=PutScalarSum================================================
 Intermediate print of scalar sum
*/
void PutScalarSum(uint a)
{
  if(a == NIL)
    PutSymbol('0');
  else
  {
    if(INTEGER_IS_NEGATIVE(SCALAR_TERM_NUMERATOR(a)))
      PutString("- ");
    PutScalarBareTerm(a);
    while((a = SCALAR_TERM_R(a)) != NIL)
    {
      PutString(INTEGER_IS_NEGATIVE(SCALAR_TERM_NUMERATOR(a)) ?
                " - " : " + ");
      PutScalarBareTerm(a);
    }
  }
}
/*=PutStart=======================================
 Start of 2dimensional output
*/
void PutStart(void)
{
  LastItemEnd = PrintEnd = Margin = NewMargin = 0;
  AddSymbolToOutLine(LEVEL, 0);
  AddSymbolToOutLine(MAIN_LEVEL, 1);
  PosOutLine = 1;
  CurrentLevel = MaxLevel = MinLevel = MAIN_LEVEL;
}
/*=PutStatistics=========================================================
 Put time and space statistics
*/
void PutStatistics(void)
{
  uint sec, min, sec_100,
     time = CrudeTime ? TimeC : (uint)(((double)TimeC/CLOCKS_PER_SEC)*100);
  PutStringStandard("Time: ");
  if(!CrudeTime)
  {
    sec_100 = (uint)(time%100);
    time /= 100;                /* In seconds */
  }
  sec = (uint)(time%60);
  time /= 60;                   /* In minutes */
  min = (uint)(time%60);
  time /= 60;                   /* In hours */
  if(time)
    PutFormattedU("%lu h ", time);
  if(min || time)
    PutFormattedU("%lu min ", min);
  if(CrudeTime)
    PutFormattedU("%lu sec\n", sec);
  else
  {
    PutFormattedU("%u.", sec);
    if(sec_100 < 10 && sec_100 > 0)
      PutCharacter('0');
    PutFormattedU("%u sec\n", sec_100);
  }
#if defined(SPACE_STATISTICS)
  PutFormattedU("Number of relations: %14u\n", MaxNRelation);
  PutFormattedU("Number of Lie monomials: %10u\n", LieMonomialMaxN);
  PutFormattedU("Number of Lie terms: %14u\n", NodeLTTopMax - 1);
  if(IsParametric)
  {
    PutFormattedU("Number of scalar terms: %11u\n", NodeSTTopMax - 1);
    min = NodeSFTopMax - 1;
    PutFormattedU("Number of scalar factors:%10u\n", NodeSFTopMax - 1);
  }
#endif
#if defined(INTEGER_MAX_SIZE)
  PutFormattedU("Maximum size of integer in limbs: %5u\n", IntegerMaxSize);
#endif
#if defined(DEBUG)
  PutDebugU("Current Debug", Debug);
#endif
}
/*=PutString================
 2D output of string
*/
void PutString(char *str)
{
  PreviousEnd = LastItemEnd;
  InLineString(str);
  PutBlock();
}
#if defined(GAP)
/*=PutStringGAP=======================================================
 Put string in GAP file
*/
void PutStringGAP(char *str)
{
  char c = 0;
  while(*str)
  {
     if(PosOutLine == GAP_WIDTH - 2)
       if(isdigit(*str))         /* Going to write in last position */
       {
         if(isdigit(c))
           PutCharacter('\\');
         PutCharacter('\n');
         PosOutLine = -1;
         PutCharacterGAP(*str); /* Continue number in the next line */
       }
       else
       {
         PutCharacter(*str);
         PutCharacter('\n');
         PosOutLine = -1;                     /* Ready to next line */
       }
     else
       PutCharacterGAP(*str);
     c = *str;                          /* Remember previous symbol */
     ++str;
  }
}
#endif
/*=PutStringStandard==============
*/
void PutStringStandard(char *str)
{
#if defined(ECHO_TO_SCREEN)
  printf("%s", str);
#endif
#if !defined(GAP)
  fprintf(SessionFile, "%s", str);
#endif
}
/*=PutSymbol!===============
 2D output of symbol
*/
void PutSymbol(char c)
{
  PreviousEnd = LastItemEnd;
  InLineSymbol(c);
  PutBlock();
}

/*_6_10   Debugging functions===========================================*/

#if defined(DEBUG)
/*=PutDebugHeader====================================================
  Put header of tracing output.
*/
void PutDebugHeader(uint debug, char * f_name, char * in_out)
{
#if !defined(GAP)
  fprintf(SessionFile,"\nDebug==%lu %s %s:\n", debug, f_name, in_out);
#endif
  printf("\nDebug==%lu %s %s\n", debug, f_name, in_out);
}
/*=PutDebugInteger==============
 Put name and signed big number
*/
void PutDebugInteger(char * name, BIGINT u)
{
  PutStart();
  InLineString(name);
  InLineString(": ");
  if(INTEGER_IS_NEGATIVE(u))
    InLineSymbol('-');
  PutIntegerUnsigned(u);
  PutEnd();
}
/*=PutDebugLieMonomial==============
  Put name and Lie sum
*/
void PutDebugLieMonomial(char * name, int a)
{
  PutStart();
  InLineString(name);
  InLineString(": ");
  (*PutLieMonomial)(a);
  PutEnd();
}
/*=PutDebugLieMonomialTable=================================================
*/
void PutDebugLieMonomialTable(int newmon)
{
  int i, j, count;
  PutStringStandard("LieMonomial table:\n"
               "ORDER POSITION LEFT RIGHT PARITY INDEX WEIGHT MONOMIAL\n");
  i = count = 0;
  while(count < LieMonomialN)
  {
    if(LIE_MONOMIAL_IS_OCCUPIED(i))
    {
      count++;
      PutFormattedU("%5d", LIE_MONOMIAL_ORDER(i));
      PutFormattedU("%9d", i);
      PutFormattedU("%5d", LIE_MONOMIAL_LEFT(i));
      PutFormattedU("%6d", LIE_MONOMIAL_RIGHT(i));
      PutFormattedU("%4d", LIE_MONOMIAL_PARITY(i));
      if((j = LIE_MONOMIAL_INDEX(i)) < 0 )
        PutFormattedU(" Rel %-4d", ~j);
      else
        PutFormattedU("%9d", j);
      PutFormattedU("%7d ", LIE_MONOMIAL_WEIGHT(i));
      if(newmon == i)
        PutStringStandard("New! ");
      PutStart();
      (*PutLieMonomial)(i);
      PutEnd();
    }
    else
      PutFormattedU("Position %d is free!\n", i);
    i++;
  }
}
/*=PutDebugLieSum==============
  Put name and Lie sum
*/
void PutDebugLieSum(char * name, uint a)
{
  PutStart();
  InLineString(name);
  InLineString(": ");
  PutLieSum(PutLieMonomial, a);
  PutEnd();
}
/*=PutDebugLieTerm==============================================
*/
void PutDebugLieTerm(char * name, uint a)
{
  PutStart();
  InLineString(name);
  InLineString(": ");
  if(a == NIL)
    PutSymbol('0');
  else
  {
    uint na;
    int is_negative = NO;
    if(IsParametric)
    {
      na = LIE_TERM_NUMERATOR_SCALAR_SUM(a);
      if(SCALAR_TERM_R(na) == NIL &&
         INTEGER_IS_NEGATIVE(SCALAR_TERM_NUMERATOR(na)))
        is_negative = YES;
    }
    else if(INTEGER_IS_NEGATIVE(LIE_TERM_NUMERATOR_INTEGER(a)))
      is_negative = YES;
    if(is_negative)
      PutString("- ");
    PutLieBareTerm(PutLieMonomial, a);
  }
  PutEnd();
}
/*=PutDebugU==================================
  Put name and long unsigned integer
*/
void PutDebugU(char * name, uint i)
{
  printf("\n%s==%lu\n", name, i);
#if !defined(GAP)
  fprintf(SessionFile, "\n%s==%lu\n", name, i);
#endif
}
#if defined(D_PUT_RELATIONS)
/*=PutDebugRelations======================================
 Put Relation table with structure fields
*/
void PutDebugRelations(void)
{
  int i, mg;
  char * hformat = "Relations:\n  N SUB MIN_GEN  RELATION\n";
  printf(hformat);
#if !defined(GAP)
  fprintf(SessionFile, hformat);
#endif
  for(i = 0; i < RelationN; i++)
  {
    PutStart();
    if(i < 100)
      InLineSymbol(' ');
    if(i < 10)
      InLineSymbol(' ');
    InLineString(UToString(i));
    InLineString("  ");
    InLineSymbol(RELATION_TO_BE_SUBSTITUTED(i) ? 'Y':'N');
    InLineString(" ");
    mg = RELATION_MIN_GENERATOR(i);
    if(mg < 10)
      InLineSymbol(' ');
    InLineString(UToString(mg));
    InLineString(" ");
    if(mg < GeneratorN)
      PutLieMonomial(mg);
    else
      InLineString("done");
    InLineString("    ");
    IN_LINE_MARGIN;
    PutLieSum(PutLieMonomial, RELATION_LIE_SUM(i));
    PutEnd();
  }
}
#endif
/*=PutDebugScalarSum==============
  Put name and scalar sum
*/
void PutDebugScalarSum(char * name, uint a)
{
  PutStart();
  InLineString(name);
  InLineString(": ");
  PutScalarSum(a);
  PutEnd();
}
/*=PutDebugString================================
  Put name of string and string
*/
void PutDebugString(char * strname, char * str)
{
  printf("%s: %s\n", strname, str);
#if !defined(GAP)
  fprintf(SessionFile, "%s: %s\n", strname, str);
#endif
}
#endif
#if defined(MEMORY)
/*=AddLieSumNs=============================================

*/
void AddLieSumNs(uint a, int minus_or_plus,
              int *pn_lt, int *pn_int, int *pn_st, int *pn_sf)
{
  int dn_lt, dn_int;
  dn_lt = dn_int = 0;
  while(a != NIL)
  {
    ++dn_lt;
    if(IsParametric)
    {
      AddScalarSumNs(LIE_TERM_NUMERATOR_SCALAR_SUM(a),
                     minus_or_plus, pn_int, pn_st, pn_sf);
      AddScalarSumNs(LIE_TERM_DENOMINATOR_SCALAR_SUM(a),
                     minus_or_plus, pn_int, pn_st, pn_sf);
    }
    else
    {
      ++dn_int;             /* Numerator is obligatory */
      if(LIE_TERM_DENOMINATOR_INTEGER(a) != NULL)
        ++dn_int;
    }
    a = LIE_TERM_R(a);
  }
  if(minus_or_plus == PLUS)
  {
    *pn_lt += dn_lt;
    if(!IsParametric)
      *pn_int += dn_int;
  }
  else
  {
    *pn_lt -= dn_lt;
    if(!IsParametric)
      *pn_int -= dn_int;
  }
}
/*=AddScalarSumNs========================================================
*/
void AddScalarSumNs(uint a, int minus_or_plus, int *pn_int, int *pn_st, int *pn_sf)
{
  int dn_int, dn_st, dn_sf;
  uint b;
  dn_int = dn_st = dn_sf = 0;
  while(a != NIL)
  {
    ++dn_st;
    ++dn_int;    /* Numerator is obligatory */
    b = SCALAR_TERM_MONOMIAL(a);
    while(b != NIL)
    {
      ++dn_sf;
      b = SCALAR_FACTOR_R(b);
    }
    a = SCALAR_TERM_R(a);
  }
  if(minus_or_plus == PLUS)
  {
    *pn_int += dn_int;
    *pn_st += dn_st;
    *pn_sf += dn_sf;
  }
  else
  {
    *pn_int -= dn_int;
    *pn_st -= dn_st;
    *pn_sf -= dn_sf;
  }
}
/*=PutIntegerBalance===================================================
*/
void PutIntegerBalance(char * fname, int dn)
{
  PutStringStandard("\nHeap integer balance violation in function:\n");
  PutStringStandard(fname);
  if(dn > 0)
  {
#if !defined(GAP)
    fprintf(SessionFile, "\n*** %ld INTs gone to garbage\n", dn);
#endif
    printf("\n*** %ld INTs gone to garbage\n", dn);
  }
  else
  {
    dn = -dn;
#if !defined(GAP)
    fprintf(SessionFile, "\n*** %ld INTs appeared from nothing\n", dn);
#endif
    printf("\n*** %ld INTs appeared from nothing\n", dn);
  }
}
/*=PutNodeBalance======================================================
*/
void PutNodeBalance(char * type, char * fname, int dn)
{
  PutStringStandard(type);
  PutStringStandard(" balance violation in function:\n");
  PutStringStandard(fname);
  if(dn > 0)
  {
#if !defined(GAP)
    fprintf(SessionFile, "\n*** %ld nodes gone to garbage\n", dn);
#endif
    printf("\n*** %ld nodes gone to garbage", dn);
  }
  else
  {
    dn = -dn;
#if !defined(GAP)
    fprintf(SessionFile, "\n***%ld nodes appeared from nothing\n", dn);
#endif
    printf("\n***%ld nodes appeared from nothing\n", dn);
  }
}
#endif