Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.
Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.
| Download
GAP 4.8.9 installation with standard packages -- copy to your CoCalc project to get it
Project: cocalc-sagemath-dev-slelievre
Views: 418346/************************************************************************/1/* File fplsa4.c Last touch January 5, 1999 */2/* Calculation of finitely presented Lie superalgebras */3/* Version 4 of April 15, 1997. */4/* e-mail: [email protected], [email protected] */5/* Contents */6/*_0 Choice of compilation */7/*_1 System constants */8/*_2 Type definitions */9/*_3 Constants and enumerations */10/*_4 Macrodefinitions */11/*_5 Global variables and arrays */12/*_6 Function descriptions */13/*_6_0 Main and top level functions */14/*_6_1 Pairing functions */15/*_6_2 Substitution (replacing) functions */16/*_6_3 Lie and scalar algebra functions */17/*_6_4 Scalar polynomial algebraic functions */18/*_6_5 Big number functions */19/*_6_6 Copy and delete functions */20/*_6_7 Technical functions */21/*_6_8 Input functions */22/*_6_9 Output functions */23/*_6_10 Debugging functions */24/************************************************************************/252627/*_0 Choice of compilation============================================*/2829#define RATIONAL_FIELD /* Working over the field R ??30otherwise over the ring Z */31#define ECHO_TO_SCREEN /* Echo session to screen ?? */32#define RELATION_N_TO_SCREEN /* Watch increasing RelationN ?? */33//#define SPP_2000 /* Big ending computer ?? */34#define SPACE_STATISTICS /* Space statistics ?? */35//#define INTEGER_MAX_SIZE /* Multiprecision number maximum size ?? */36//#define INTEGER_ALLOCATION_CHECK /* Control of integer allocations ?? */37//#define POLY_ARRAY_ALLOCATION_CHECK /* Control of allocations of ?? plynomial arrays in stack */3839/* GAP output ?? */40/* Avoid message file, session file, */41/* compulsory suffix '.in' for input file, */42/* and printing comments to the screen ? */43#define GAP44/**/4546//#define DEBUG /* Debugging ?? */47//#define MEMORY /* Check memory balance ?? */4849/* Include files */5051#include <stdio.h>52#include <time.h>53#include <stdlib.h>54#include <ctype.h>55#include <string.h>56#if defined(SPP_2000)57#include <alloca.h> /* This file the genuine SPP compiler requires */58#endif596061#if defined(__GNUC__) || defined(__clang__)62#define ATTRIBUTE_NORETURN __attribute__ ((noreturn))63#else64#define ATTRIBUTE_NORETURN65#endif666768#if defined(DEBUG) /* Debug definitions ===============================*/6970/* Set condition for debug output ??71*/72#define D_CONDITION if(Debug>=5873)73/* Examples of D_CONDITION74`empty' if(Debug>=11951211) current!!75if(Debug>=25283LU)*e7* *e70*76if(Debug%10 == 0) if(Debug>=17566) */77/* Set condition to stop debugging ??78*/79#define D_EXIT if(Debug > 5882) EXIT;80/* Examples of D_EXIT81`empty' if(Debug > 21420LU) EXIT;82if(Debug > 26969LU) EXIT;*e7* if(Debug > 21420LU) EXIT;*e70*83if(Debug > 30881LU) EXIT; if(Debug >= 18399) EXIT;84*/8586/* Switches for debugging particular functions ?? */8788#define D_CHECK_EXACTNESS_OF_DIVISION89//#define D_ADD_PAIR_TO_LIE_MONOMIAL90#define D_GENERATE_RELATIONS91//#if defined(D_GENERATE_RELATIONS) /* Put tables ?? */92//#define D_PUT_RELATIONS93//#define D_PUT_LIE_MONOMIAL94//#define D_GET_LIE_MONOMIAL95//#define D_GET_LIE_SUM96//#define D_GET_LIE_TERM97#define D_INTEGER_CANCELLATION98#define D_INTEGER_GCD99#define D_INTEGER_PRODUCT100#define D_INTEGER_QUOTIENT101#define D_INTEGER_SUM102#define D_LIE_SUM_ADDITION103#define D_LIE_SUM_DIV_INTEGER104//#define D_LIE_SUM_DIV_SCALAR_SUM105#define D_LIE_SUM_MULT_INTEGER106//#define D_LIE_SUM_MULT_SCALAR_SUM107#define D_MAKE_RELATION_RHS108//#define D_NEW_JACOBI_RELATION109#define D_NORMALIZE_RELATION110//#define D_PAIR_MONOMIAL_MONOMIAL111//#define D_PAIR_MONOMIAL_SUM112//#define D_PAIR_SUM_MONOMIAL113//#define D_PAIR_SUM_SUM114//#define D_POLY_CONTENT115//#define D_POLY_GCD116//#define D_POLY_QUOTIENT117//#define D_POLY_PSEUDO_REMAINDER118//#define D_SCALAR_SUM_CANCELLATION119#define D_SUBSTITUTE_RELATION_IN_RELATION120#define D_IN_SET uint debug=Debug;\121D_CONDITION\122PutDebugHeader(debug, f_name, "in"),123#define D_IN_CLOSE Debug++; D_EXIT124125#define D_OUT_OPEN D_CONDITION PutDebugHeader(debug, f_name, "out"),126127#else /* Empty definitions for DEBUG off */128129#define D_IN_SET130#define D_IN_CLOSE131#define D_OUT_OPEN132133#endif134135/* Check memory balance definitions */136137#if defined(MEMORY) /* `n_...' are unique */138#define IN_SET_N_LT int n_lt = -CurrentNLT;139#define IN_SET_N_INT int n_int = -CurrentNINT;140#define IN_SET_N_ST int n_st = -CurrentNST;141#define IN_SET_N_SF int n_sf = -CurrentNSF;142143#define OUT_SET_N_LT n_lt += CurrentNLT;144#define OUT_SET_N_INT n_int += CurrentNINT;145#define OUT_SET_N_ST n_st += CurrentNST;146#define OUT_SET_N_SF n_sf += CurrentNSF;147148#define ADD_LIE_SUM_NS(a) \149AddLieSumNs(a, PLUS, &n_lt, &n_int, &n_st, &n_sf);150151#define SUBTRACT_LIE_SUM_NS(a) \152AddLieSumNs(a, MINUS, &n_lt, &n_int, &n_st, &n_sf);153154#define ADD_SCALAR_SUM_NS(a) \155AddScalarSumNs(a, PLUS, &n_int, &n_st, &n_sf);156157#define SUBTRACT_SCALAR_SUM_NS(a) \158AddScalarSumNs(a, MINUS, &n_int, &n_st, &n_sf);159160#define CHECK_NS \161if(n_lt != 0) {PutNodeBalance("\nNodeLT (Lie terms)",\162f_name, n_lt); EXIT;}\163if(n_int != 0) {PutIntegerBalance(f_name, n_int); EXIT;}\164if(n_st != 0) {PutNodeBalance("\nNodeST (scalar terms)",\165f_name, n_st); EXIT;}\166if(n_sf != 0) {PutNodeBalance("\nNodeSF (scalar factors)",\167f_name, n_sf); EXIT;}168169/* Particular functions */170/*----------------------------------------------------------*/171#define M_OUT_GET_LIE_MONOMIAL \172OUT_SET_NS\173SUBTRACT_LIE_SUM_NS(a)\174CHECK_NS175/*----------------------------------------------------------*/176#define M_OUT_GET_LIE_SUM \177OUT_SET_NS\178SUBTRACT_LIE_SUM_NS(lsum)\179CHECK_NS180/*----------------------------------------------------------*/181#define M_OUT_GET_LIE_TERM \182OUT_SET_NS\183SUBTRACT_LIE_SUM_NS(lterm)\184CHECK_NS185/*----------------------------------------------------------*/186#define M_IN_LIE_SUM_ADDITION \187ADD_LIE_SUM_NS(a)\188ADD_LIE_SUM_NS(b)189#define M_OUT_LIE_SUM_ADDITION \190OUT_SET_NS\191SUBTRACT_LIE_SUM_NS(sum)\192CHECK_NS193/*----------------------------------------------------------*/194/* No sense to print debug information, memory check only */195#define IN_LIE_SUM_COPY \196IN_SET_FUNCTION_NAME("LieSumCopy...")\197IN_SET_NS198#define OUT_LIE_SUM_COPY \199OUT_SET_NS\200SUBTRACT_LIE_SUM_NS(ca)\201CHECK_NS202/*----------------------------------------------------------*/203/* No sense to print debug information, memory check only */204#define IN_LIE_SUM_KILL \205IN_SET_FUNCTION_NAME("LieSumKill...")\206IN_SET_NS\207ADD_LIE_SUM_NS(a)208#define OUT_LIE_SUM_KILL \209OUT_SET_NS\210CHECK_NS211/*----------------------------------------------------------*/212#define M_IN_LIE_SUM_DIV_INTEGER \213ADD_LIE_SUM_NS(b)214#define M_OUT_LIE_SUM_DIV_INTEGER \215OUT_SET_NS\216SUBTRACT_LIE_SUM_NS(b)\217CHECK_NS218/*----------------------------------------------------------*/219#define M_IN_LIE_SUM_DIV_SCALAR_SUM \220ADD_LIE_SUM_NS(b)\221ADD_SCALAR_SUM_NS(den)222#define M_OUT_LIE_SUM_DIV_SCALAR_SUM \223OUT_SET_NS\224SUBTRACT_LIE_SUM_NS(b)\225CHECK_NS226/*----------------------------------------------------------*/227#define M_IN_LIE_SUM_MULT_SCALAR_SUM \228ADD_LIE_SUM_NS(b)\229ADD_SCALAR_SUM_NS(num)230#define M_OUT_LIE_SUM_MULT_SCALAR_SUM \231OUT_SET_NS\232SUBTRACT_LIE_SUM_NS(b)\233CHECK_NS234/*----------------------------------------------------------*/235#define M_IN_LIE_SUM_MULT_INTEGER \236ADD_LIE_SUM_NS(b)237#define M_OUT_LIE_SUM_MULT_INTEGER M_OUT_LIE_SUM_DIV_INTEGER238/*----------------------------------------------------------*/239#define M_OUT_MAKE_RELATION_RHS M_OUT_GET_LIE_MONOMIAL240/*----------------------------------------------------------*/241#define M_OUT_NEW_JACOBI_RELATION M_OUT_GET_LIE_MONOMIAL242/*----------------------------------------------------------*/243#define M_IN_NORMALIZE_RELATION \244ADD_LIE_SUM_NS(a)245#define M_OUT_NORMALIZE_RELATION M_OUT_GET_LIE_MONOMIAL246/*----------------------------------------------------------*/247#define M_OUT_PAIR_MONOMIAL_MONOMIAL M_OUT_GET_LIE_MONOMIAL248/*----------------------------------------------------------*/249#define M_IN_PAIR_MONOMIAL_SUM \250ADD_LIE_SUM_NS(a)251#define M_OUT_PAIR_MONOMIAL_SUM \252OUT_SET_NS\253SUBTRACT_LIE_SUM_NS(s)\254CHECK_NS255/*----------------------------------------------------------*/256#define M_IN_PAIR_SUM_MONOMIAL \257ADD_LIE_SUM_NS(a)258#define M_OUT_PAIR_SUM_MONOMIAL \259OUT_SET_NS\260SUBTRACT_LIE_SUM_NS(s)\261CHECK_NS262/*----------------------------------------------------------*/263#define M_IN_PAIR_SUM_SUM \264ADD_LIE_SUM_NS(a)\265ADD_LIE_SUM_NS(b)266#define M_OUT_PAIR_SUM_SUM \267OUT_SET_NS\268SUBTRACT_LIE_SUM_NS(a)\269CHECK_NS270/*----------------------------------------------------------*/271#define M_IN_POLY_CONTENT272#define M_OUT_POLY_CONTENT \273OUT_SET_NS\274SUBTRACT_SCALAR_SUM_NS(b)\275CHECK_NS276/*----------------------------------------------------------*/277#define M_IN_POLY_GCD278#define M_OUT_POLY_GCD \279OUT_SET_NS\280SUBTRACT_SCALAR_SUM_NS(b)\281CHECK_NS282/*----------------------------------------------------------*/283#define M_IN_POLY_QUOTIENT \284ADD_SCALAR_SUM_NS(a)\285ADD_SCALAR_SUM_NS(b)286#define M_OUT_POLY_QUOTIENT \287OUT_SET_NS\288SUBTRACT_SCALAR_SUM_NS(c)\289CHECK_NS290/*----------------------------------------------------------*/291#define M_IN_POLY_PSEUDO_REMAINDER \292ADD_SCALAR_SUM_NS(a)\293ADD_SCALAR_SUM_NS(b)294#define M_OUT_POLY_PSEUDO_REMAINDER \295OUT_SET_NS\296SUBTRACT_SCALAR_SUM_NS(a)\297CHECK_NS298/*----------------------------------------------------------*/299/* No sense to print debug information, memory check only */300#define IN_REDUCE_RELATIONS \301IN_SET_FUNCTION_NAME("ReduceRelations")\302IN_SET_NS\303{\304int o;\305for(o = 0; o < RelationN; o++)\306ADD_LIE_SUM_NS(RELATION_LIE_SUM(o))\307}308#define OUT_REDUCE_RELATIONS \309OUT_SET_NS\310{\311int o;\312for(o = 0; o < RelationN; o++)\313SUBTRACT_LIE_SUM_NS(RELATION_LIE_SUM(o))\314}\315CHECK_NS316/*----------------------------------------------------------*/317#define M_IN_SCALAR_SUM_CANCELLATION \318ADD_SCALAR_SUM_NS(*pnum)\319ADD_SCALAR_SUM_NS(*pden)320#define M_OUT_SCALAR_SUM_CANCELLATION \321OUT_SET_NS\322SUBTRACT_SCALAR_SUM_NS(*pnum)\323SUBTRACT_SCALAR_SUM_NS(*pden)\324CHECK_NS325/*----------------------------------------------------------*/326#define M_IN_SUBSTITUTE_RELATION_IN_RELATION \327ADD_LIE_SUM_NS(a)328#define M_OUT_SUBSTITUTE_RELATION_IN_RELATION \329OUT_SET_NS\330SUBTRACT_LIE_SUM_NS(a)\331CHECK_NS332333#define MM_CURRENT_N_LT ,--CurrentNLT334#define PP_CURRENT_N_LT ++CurrentNLT;335336#define MM_CURRENT_N_SF ,--CurrentNSF337#define PP_CURRENT_N_SF ++CurrentNSF;338339#define MM_CURRENT_N_ST ,--CurrentNST340#define PP_CURRENT_N_ST ++CurrentNST;341342#define MM_CURRENT_N_INT ,--CurrentNINT343#define PP_CURRENT_N_INT ;++CurrentNINT344345#else /* MEMORY off */346347#define MM_CURRENT_N_LT348#define PP_CURRENT_N_LT349#define MM_CURRENT_N_SF350#define PP_CURRENT_N_SF351#define MM_CURRENT_N_ST352#define PP_CURRENT_N_ST353#define MM_CURRENT_N_INT354#define PP_CURRENT_N_INT355356#define IN_SET_N_LT357#define IN_SET_N_INT358#define IN_SET_N_ST359#define IN_SET_N_SF360361#define OUT_SET_N_LT362#define OUT_SET_N_INT363#define OUT_SET_N_ST364#define OUT_SET_N_SF365366#define M_OUT_GET_LIE_MONOMIAL367#define M_OUT_GET_LIE_SUM368#define M_OUT_GET_LIE_TERM369#define M_IN_LIE_SUM_ADDITION370#define M_OUT_LIE_SUM_ADDITION371#define IN_LIE_SUM_COPY372#define OUT_LIE_SUM_COPY373#define IN_LIE_SUM_KILL374#define OUT_LIE_SUM_KILL375#define M_IN_LIE_SUM_DIV_INTEGER376#define M_OUT_LIE_SUM_DIV_INTEGER377#define M_IN_LIE_SUM_DIV_SCALAR_SUM378#define M_OUT_LIE_SUM_DIV_SCALAR_SUM379#define M_IN_LIE_SUM_MULT_SCALAR_SUM380#define M_OUT_LIE_SUM_MULT_SCALAR_SUM381#define M_IN_LIE_SUM_MULT_INTEGER382#define M_OUT_LIE_SUM_MULT_INTEGER383#define M_OUT_MAKE_RELATION_RHS384#define M_OUT_NEW_JACOBI_RELATION385#define M_IN_NORMALIZE_RELATION386#define M_OUT_NORMALIZE_RELATION387#define M_OUT_PAIR_MONOMIAL_MONOMIAL388#define M_IN_PAIR_MONOMIAL_SUM389#define M_OUT_PAIR_MONOMIAL_SUM390#define M_IN_PAIR_SUM_MONOMIAL391#define M_OUT_PAIR_SUM_MONOMIAL392#define M_IN_PAIR_SUM_SUM393#define M_OUT_PAIR_SUM_SUM394#define M_IN_POLY_CONTENT395#define M_OUT_POLY_CONTENT396#define M_IN_POLY_GCD397#define M_OUT_POLY_GCD398#define M_IN_POLY_QUOTIENT399#define M_OUT_POLY_QUOTIENT400#define M_IN_POLY_PSEUDO_REMAINDER401#define M_OUT_POLY_PSEUDO_REMAINDER402#define IN_REDUCE_RELATIONS403#define OUT_REDUCE_RELATIONS404#define M_IN_SCALAR_SUM_CANCELLATION405#define M_OUT_SCALAR_SUM_CANCELLATION406#define M_IN_SUBSTITUTE_RELATION_IN_RELATION407#define M_OUT_SUBSTITUTE_RELATION_IN_RELATION408#endif409410#define IN_SET_NS IN_SET_N_LT IN_SET_N_INT IN_SET_N_ST IN_SET_N_SF411#define OUT_SET_NS OUT_SET_N_LT OUT_SET_N_INT OUT_SET_N_ST OUT_SET_N_SF412413414#if defined(DEBUG) || defined(MEMORY) /* `f_name' is unique */415#define IN_SET_FUNCTION_NAME(fname) char * f_name = fname;416#else417#define IN_SET_FUNCTION_NAME(fname)418#endif419420#if defined(D_ADD_PAIR_TO_LIE_MONOMIAL) /*------------------------------*/421/* Only debugging makes sense */422#define IN_ADD_PAIR_TO_LIE_MONOMIAL \423IN_SET_FUNCTION_NAME("AddPairToLieMonomial")\424D_IN_SET\425PutDebugLieMonomial("i", i),\426PutDebugLieMonomial("j", j),\427PutDebugLieMonomialTable(-1);\428D_IN_CLOSE429430#define OUT_ADD_PAIR_TO_LIE_MONOMIAL_OLD \431D_OUT_OPEN\432PutDebugLieMonomial("Old monomial", ijp);433#define OUT_ADD_PAIR_TO_LIE_MONOMIAL_NEW \434D_OUT_OPEN\435PutDebugLieMonomialTable(ijp);436#else437#define IN_ADD_PAIR_TO_LIE_MONOMIAL D_IN_CLOSE438#define OUT_ADD_PAIR_TO_LIE_MONOMIAL_OLD439#define OUT_ADD_PAIR_TO_LIE_MONOMIAL_NEW440#endif441#if defined(D_GENERATE_RELATIONS) /*-----------------------------------*/442/* Only debugging makes sense */443#define IN_GENERATE_RELATIONS \444IN_SET_FUNCTION_NAME("GenerateRelations")\445D_IN_SET\446IN_PUT_RELATIONS\447IN_PUT_LIE_MONOMIAL\448PutDebugLieSum("rel", a),\449PutDebugLieMonomial("gen", gen);\450D_IN_CLOSE451452#define OUT_GENERATE_RELATIONS \453D_OUT_OPEN\454PutDebugLieSum("d_rel", a);455#else456#define IN_GENERATE_RELATIONS D_IN_CLOSE457#define OUT_GENERATE_RELATIONS458#endif459#if defined(D_INTEGER_CANCELLATION) /*----- (k*n)/(k*d) -> n/d --------*/460/* Only debugging makes sense */461#define IN_INTEGER_CANCELLATION \462IN_SET_FUNCTION_NAME("IntegerCancellation")\463D_IN_SET\464PutDebugInteger("num", num),\465PutDebugInteger("den", den);\466D_IN_CLOSE467#define OUT_INTEGER_CANCELLATION \468D_OUT_OPEN\469PutDebugInteger("num", num),\470PutDebugInteger("den", den);471#else472#define IN_INTEGER_CANCELLATION D_IN_CLOSE473#define OUT_INTEGER_CANCELLATION474#endif475#if defined(D_INTEGER_GCD) /*------------------------------------------*/476/* Only debugging makes sense */477#define IN_INTEGER_GCD \478IN_SET_FUNCTION_NAME("IntegerGCD")\479D_IN_SET\480PutDebugInteger("u", u),\481PutDebugInteger("v", v);\482D_IN_CLOSE483#define OUT_INTEGER_GCD \484D_OUT_OPEN\485PutDebugInteger("GCD", u);486#else487#define IN_INTEGER_GCD D_IN_CLOSE488#define OUT_INTEGER_GCD489#endif490#if defined(D_INTEGER_PRODUCT) /*------------- n*m -------------------*/491/* Only debugging makes sense */492#define IN_INTEGER_PRODUCT \493IN_SET_FUNCTION_NAME("IntegerProduct")\494D_IN_SET\495PutDebugInteger("u", u),\496PutDebugInteger("v", v);\497D_IN_CLOSE498#define OUT_INTEGER_PRODUCT \499D_OUT_OPEN\500PutDebugInteger("u*v", w0);501#else502#define IN_INTEGER_PRODUCT D_IN_CLOSE503#define OUT_INTEGER_PRODUCT504#endif505#if defined(D_INTEGER_QUOTIENT) /*------------------------------------*/506/* Only debugging makes sense */507#define IN_INTEGER_QUOTIENT \508IN_SET_FUNCTION_NAME("IntegerQuotient")\509D_IN_SET\510PutDebugInteger("a", a),\511PutDebugInteger("b", b);\512D_IN_CLOSE513#define OUT_INTEGER_QUOTIENT \514D_OUT_OPEN\515PutDebugInteger("a/b", pm);516#else517#define IN_INTEGER_QUOTIENT D_IN_CLOSE518#define OUT_INTEGER_QUOTIENT519#endif520#if defined(D_INTEGER_SUM) /*------------------------------------*/521/* Only debugging makes sense */522#define IN_INTEGER_SUM \523IN_SET_FUNCTION_NAME("IntegerSum")\524D_IN_SET\525PutDebugInteger("a", a),\526PutDebugInteger("b", b);\527D_IN_CLOSE528#define OUT_INTEGER_SUM \529D_OUT_OPEN\530PutDebugInteger("c", c);531#else532#define IN_INTEGER_SUM D_IN_CLOSE533#define OUT_INTEGER_SUM534#endif535536#if defined(D_GET_LIE_MONOMIAL) /*-------------------------------------*/537#define D_IN_GET_LIE_MONOMIAL D_IN_SET\538PutDebugString("*pstr", *pstr);\539D_IN_CLOSE540#define D_OUT_GET_LIE_MONOMIAL D_OUT_OPEN\541PutDebugLieSum("a", a);542#else543#define D_IN_GET_LIE_MONOMIAL D_IN_CLOSE544#define D_OUT_GET_LIE_MONOMIAL545#endif546#if defined(D_GET_LIE_SUM) /*------------------------------------------*/547#define D_IN_GET_LIE_SUM D_IN_SET\548PutDebugString("*pstr", *pstr);\549D_IN_CLOSE550#define D_OUT_GET_LIE_SUM D_OUT_OPEN\551PutDebugLieSum("lsum", lsum);552#else553#define D_IN_GET_LIE_SUM D_IN_CLOSE554#define D_OUT_GET_LIE_SUM555#endif556#if defined(D_GET_LIE_TERM) /*-----------------------------------------*/557#define D_IN_GET_LIE_TERM D_IN_SET\558PutDebugString("*pstr", *pstr);\559D_IN_CLOSE560#define D_OUT_GET_LIE_TERM D_OUT_OPEN\561PutDebugLieSum("lterm", lterm);562#else563#define D_IN_GET_LIE_TERM D_IN_CLOSE564#define D_OUT_GET_LIE_TERM565#endif566#if defined(D_LIE_SUM_ADDITION) /*---------- a + b (Lie) --------------*/567#define D_IN_LIE_SUM_ADDITION D_IN_SET\568PutDebugLieSum("a", a),\569PutDebugLieSum("b", b);\570D_IN_CLOSE571#define D_OUT_LIE_SUM_ADDITION D_OUT_OPEN\572PutDebugLieSum("sum", sum);573#else574#define D_IN_LIE_SUM_ADDITION D_IN_CLOSE575#define D_OUT_LIE_SUM_ADDITION576#endif577#if defined(D_LIE_SUM_DIV_INTEGER) /*--------- uint/BIGINT ----------------*/578#define D_IN_LIE_SUM_DIV_INTEGER D_IN_SET\579PutDebugLieSum("lsum", b),\580PutDebugInteger("den", den);\581D_IN_CLOSE582#define D_OUT_LIE_SUM_DIV_INTEGER D_OUT_OPEN\583PutDebugLieSum("lsum", b);584#else585#define D_IN_LIE_SUM_DIV_INTEGER D_IN_CLOSE586#define D_OUT_LIE_SUM_DIV_INTEGER587#endif588#if defined(D_LIE_SUM_DIV_SCALAR_SUM) /*--------- uint/uint --------------*/589#define D_IN_LIE_SUM_DIV_SCALAR_SUM D_IN_SET\590PutDebugLieSum("lsum", b),\591PutDebugScalarSum("den", den);\592D_IN_CLOSE593#define D_OUT_LIE_SUM_DIV_SCALAR_SUM D_OUT_OPEN\594PutDebugLieSum("lsum", b);595#else596#define D_IN_LIE_SUM_DIV_SCALAR_SUM D_IN_CLOSE597#define D_OUT_LIE_SUM_DIV_SCALAR_SUM598#endif599#if defined(D_LIE_SUM_MULT_SCALAR_SUM) /*--------- uint*uint --------------*/600#define D_IN_LIE_SUM_MULT_SCALAR_SUM D_IN_SET\601PutDebugLieSum("lsum", b),\602PutDebugScalarSum("num", num);\603D_IN_CLOSE604#define D_OUT_LIE_SUM_MULT_SCALAR_SUM D_OUT_OPEN\605PutDebugLieSum("lsum", b);606#else607#define D_IN_LIE_SUM_MULT_SCALAR_SUM D_IN_CLOSE608#define D_OUT_LIE_SUM_MULT_SCALAR_SUM609#endif610#if defined(D_LIE_SUM_MULT_INTEGER) /*------- uint*BIGINT -----------------*/611#define D_IN_LIE_SUM_MULT_INTEGER D_IN_SET\612PutDebugLieSum("lsum", b),\613PutDebugInteger("num", num);\614D_IN_CLOSE615#define D_OUT_LIE_SUM_MULT_INTEGER D_OUT_OPEN\616PutDebugLieSum("lsum", b);617#else618#define D_IN_LIE_SUM_MULT_INTEGER D_IN_CLOSE619#define D_OUT_LIE_SUM_MULT_INTEGER620#endif621#if defined(D_MAKE_RELATION_RHS) /*----------------------------------*/622#define D_IN_MAKE_RELATION_RHS D_IN_SET\623PutDebugLieSum("rel",\624RELATION_LIE_SUM(i));\625D_IN_CLOSE626#define D_OUT_MAKE_RELATION_RHS D_OUT_OPEN\627PutDebugLieSum("r.h.s", a);628#else629#define D_IN_MAKE_RELATION_RHS D_IN_CLOSE630#define D_OUT_MAKE_RELATION_RHS631#endif632#if defined(D_NEW_JACOBI_RELATION) /*----------------------------------*/633#define D_IN_NEW_JACOBI_RELATION D_IN_SET\634PutDebugLieMonomial("l", l);\635D_IN_CLOSE636#define D_OUT_NEW_JACOBI_RELATION D_OUT_OPEN\637PutDebugLieSum("a", a);638#else639#define D_IN_NEW_JACOBI_RELATION D_IN_CLOSE640#define D_OUT_NEW_JACOBI_RELATION641#endif642#if defined(D_NORMALIZE_RELATION) /*-----------------------------------*/643#define D_IN_NORMALIZE_RELATION D_IN_SET\644PutDebugLieSum("a", a);\645D_IN_CLOSE646#define D_OUT_NORMALIZE_RELATION D_OUT_OPEN\647PutDebugLieSum("a", a);648#else649#define D_IN_NORMALIZE_RELATION D_IN_CLOSE650#define D_OUT_NORMALIZE_RELATION651#endif652#if defined(D_PAIR_MONOMIAL_MONOMIAL) /*---------[mona, monb]----------*/653#define D_IN_PAIR_MONOMIAL_MONOMIAL D_IN_SET\654PutDebugLieMonomial("i", i),\655PutDebugLieMonomial("j", j);\656D_IN_CLOSE657#define D_OUT_PAIR_MONOMIAL_MONOMIAL D_OUT_OPEN\658PutDebugLieSum("a", a);659#else660#define D_IN_PAIR_MONOMIAL_MONOMIAL D_IN_CLOSE661#define D_OUT_PAIR_MONOMIAL_MONOMIAL662#endif663#if defined(D_PAIR_MONOMIAL_SUM) /*-----------[mon, a]-----------------*/664#define D_IN_PAIR_MONOMIAL_SUM D_IN_SET\665PutDebugLieMonomial("mon", mon),\666PutDebugLieSum("a", a);\667D_IN_CLOSE668#define D_OUT_PAIR_MONOMIAL_SUM D_OUT_OPEN\669PutDebugLieSum("s", s);670#else671#define D_IN_PAIR_MONOMIAL_SUM D_IN_CLOSE672#define D_OUT_PAIR_MONOMIAL_SUM673#endif674#if defined(D_PAIR_SUM_MONOMIAL) /*-----------[a, mon]-----------------*/675#define D_IN_PAIR_SUM_MONOMIAL D_IN_SET\676PutDebugLieSum("a", a),\677PutDebugLieMonomial("mon", mon);\678D_IN_CLOSE679#define D_OUT_PAIR_SUM_MONOMIAL D_OUT_OPEN\680PutDebugLieSum("s", s);681#else682#define D_IN_PAIR_SUM_MONOMIAL D_IN_CLOSE683#define D_OUT_PAIR_SUM_MONOMIAL684#endif685#if defined(D_PAIR_SUM_SUM) /*----------------[a, b]---------------*/686#define D_IN_PAIR_SUM_SUM D_IN_SET\687PutDebugLieSum("a", a),\688PutDebugLieSum("b", b);\689D_IN_CLOSE690#define D_OUT_PAIR_SUM_SUM D_OUT_OPEN\691PutDebugLieSum("a", a);692#else693#define D_IN_PAIR_SUM_SUM D_IN_CLOSE694#define D_OUT_PAIR_SUM_SUM695#endif696#if defined(D_POLY_CONTENT) /*--------------------------------------------*/697#define D_IN_POLY_CONTENT D_IN_SET\698PutDebugScalarSum("a", a),\699PutDebugString("mp", \700ParameterName + mp*NameLength1);\701D_IN_CLOSE702#define D_OUT_POLY_CONTENT D_OUT_OPEN\703PutDebugScalarSum("b", b);704#else705#define D_IN_POLY_CONTENT D_IN_CLOSE706#define D_OUT_POLY_CONTENT707#endif708#if defined(D_POLY_GCD) /*--------------------------------------------*/709#define D_IN_POLY_GCD D_IN_SET\710PutDebugScalarSum("a", a),\711PutDebugScalarSum("b", b);\712D_IN_CLOSE713#define D_OUT_POLY_GCD D_OUT_OPEN\714PutDebugScalarSum("b", b);715#else716#define D_IN_POLY_GCD D_IN_CLOSE717#define D_OUT_POLY_GCD718#endif719#if defined(D_POLY_QUOTIENT) /*--------------------------------------------*/720#define D_IN_POLY_QUOTIENT D_IN_SET\721PutDebugScalarSum("a", a),\722PutDebugScalarSum("b", b);\723D_IN_CLOSE724#define D_OUT_POLY_QUOTIENT D_OUT_OPEN\725PutDebugScalarSum("c", c);726#else727#define D_IN_POLY_QUOTIENT D_IN_CLOSE728#define D_OUT_POLY_QUOTIENT729#endif730#if defined(D_POLY_PSEUDO_REMAINDER) /*------------------------------------*/731#define D_IN_POLY_PSEUDO_REMAINDER D_IN_SET\732PutDebugScalarSum("a", a),\733PutDebugScalarSum("b", b);\734D_IN_CLOSE735#define D_OUT_POLY_PSEUDO_REMAINDER D_OUT_OPEN\736PutDebugScalarSum("a", a);737#else738#define D_IN_POLY_PSEUDO_REMAINDER D_IN_CLOSE739#define D_OUT_POLY_PSEUDO_REMAINDER740#endif741#if defined(D_SCALAR_SUM_CANCELLATION) /*----- (k*n)/(k*d) -> n/d ----*/742#define D_IN_SCALAR_SUM_CANCELLATION D_IN_SET\743PutDebugScalarSum("*pnum", *pnum),\744PutDebugScalarSum("*pden", *pden);\745D_IN_CLOSE746#define D_OUT_SCALAR_SUM_CANCELLATION D_OUT_OPEN\747PutDebugScalarSum("*pnum", *pnum),\748PutDebugScalarSum("*pden", *pden);749#else750#define D_IN_SCALAR_SUM_CANCELLATION D_IN_CLOSE751#define D_OUT_SCALAR_SUM_CANCELLATION752#endif753#if defined(D_SUBSTITUTE_RELATION_IN_RELATION) /*--------------------*/754#define D_IN_SUBSTITUTE_RELATION_IN_RELATION D_IN_SET\755PutDebugLieSum("r", r),\756PutDebugLieSum("a", a);\757D_IN_CLOSE758#define D_OUT_SUBSTITUTE_RELATION_IN_RELATION D_OUT_OPEN\759PutDebugLieSum("a", a);760#else761#define D_IN_SUBSTITUTE_RELATION_IN_RELATION D_IN_CLOSE762#define D_OUT_SUBSTITUTE_RELATION_IN_RELATION763#endif764765/* Set INs and OUTs for particular functions */766767/*--------------------------------------------------------------------*/768#define IN_GET_LIE_MONOMIAL \769IN_SET_FUNCTION_NAME("GetLieMonomial")\770IN_SET_NS\771D_IN_GET_LIE_MONOMIAL772#define OUT_GET_LIE_MONOMIAL \773D_OUT_GET_LIE_MONOMIAL\774M_OUT_GET_LIE_MONOMIAL775/*--------------------------------------------------------------------*/776#define IN_GET_LIE_SUM \777IN_SET_FUNCTION_NAME("GetLieSum")\778IN_SET_NS\779D_IN_GET_LIE_SUM780#define OUT_GET_LIE_SUM \781D_OUT_GET_LIE_SUM\782M_OUT_GET_LIE_SUM783/*--------------------------------------------------------------------*/784#define IN_GET_LIE_TERM \785IN_SET_FUNCTION_NAME("GetLieTerm")\786IN_SET_NS\787D_IN_GET_LIE_TERM788#define OUT_GET_LIE_TERM \789D_OUT_GET_LIE_TERM\790M_OUT_GET_LIE_TERM791/*--------------------------------------------------------------------*/792#define IN_LIE_SUM_ADDITION \793IN_SET_FUNCTION_NAME("LieSumAddition")\794IN_SET_NS\795D_IN_LIE_SUM_ADDITION\796M_IN_LIE_SUM_ADDITION797#define OUT_LIE_SUM_ADDITION \798D_OUT_LIE_SUM_ADDITION\799M_OUT_LIE_SUM_ADDITION800/*--------------------------------------------------------------------*/801#define IN_LIE_SUM_DIV_INTEGER \802LS_B_LSUM_DIV_INT\803IN_SET_FUNCTION_NAME("LieSumDivInteger")\804IN_SET_NS\805D_IN_LIE_SUM_DIV_INTEGER\806M_IN_LIE_SUM_DIV_INTEGER807#define OUT_LIE_SUM_DIV_INTEGER \808D_OUT_LIE_SUM_DIV_INTEGER\809M_OUT_LIE_SUM_DIV_INTEGER810/*--------------------------------------------------------------------*/811#define IN_LIE_SUM_DIV_SCALAR_SUM \812LS_B_LSUM_DIV_SS\813IN_SET_FUNCTION_NAME("LieSumDivScalarSum")\814IN_SET_NS\815D_IN_LIE_SUM_DIV_SCALAR_SUM\816M_IN_LIE_SUM_DIV_SCALAR_SUM817#define OUT_LIE_SUM_DIV_SCALAR_SUM \818D_OUT_LIE_SUM_DIV_SCALAR_SUM\819M_OUT_LIE_SUM_DIV_SCALAR_SUM820/*--------------------------------------------------------------------*/821#define IN_LIE_SUM_MULT_SCALAR_SUM \822LS_B_LSUM_MULT_SS\823IN_SET_FUNCTION_NAME("LieSumMultScalarSum")\824IN_SET_NS\825D_IN_LIE_SUM_MULT_SCALAR_SUM\826M_IN_LIE_SUM_MULT_SCALAR_SUM827#define OUT_LIE_SUM_MULT_SCALAR_SUM \828D_OUT_LIE_SUM_MULT_SCALAR_SUM\829M_OUT_LIE_SUM_MULT_SCALAR_SUM830/*--------------------------------------------------------------------*/831#define IN_LIE_SUM_MULT_INTEGER \832LS_B_LSUM_MULT_INT\833IN_SET_FUNCTION_NAME("LieSumMultInteger")\834IN_SET_NS\835D_IN_LIE_SUM_MULT_INTEGER\836M_IN_LIE_SUM_MULT_INTEGER837#define OUT_LIE_SUM_MULT_INTEGER \838D_OUT_LIE_SUM_MULT_INTEGER\839M_OUT_LIE_SUM_MULT_INTEGER840/*--------------------------------------------------------------------*/841#define IN_MAKE_RELATION_RHS \842IN_SET_FUNCTION_NAME("MakeRelationRHS...")\843IN_SET_NS\844D_IN_MAKE_RELATION_RHS845#define OUT_MAKE_RELATION_RHS \846D_OUT_MAKE_RELATION_RHS\847M_OUT_MAKE_RELATION_RHS848/*--------------------------------------------------------------------*/849#define IN_NEW_JACOBI_RELATION \850IN_SET_FUNCTION_NAME("NewJacobiRelation")\851IN_SET_NS\852D_IN_NEW_JACOBI_RELATION853#define OUT_NEW_JACOBI_RELATION \854D_OUT_NEW_JACOBI_RELATION\855M_OUT_NEW_JACOBI_RELATION856/*--------------------------------------------------------------------*/857#define IN_NORMALIZE_RELATION \858IN_SET_FUNCTION_NAME("NormalizeRelation...")\859IN_SET_NS\860D_IN_NORMALIZE_RELATION\861M_IN_NORMALIZE_RELATION862#define OUT_NORMALIZE_RELATION \863D_OUT_NORMALIZE_RELATION\864M_OUT_NORMALIZE_RELATION865/*--------------------------------------------------------------------*/866#define IN_PAIR_MONOMIAL_MONOMIAL \867IN_SET_FUNCTION_NAME("PairMonomialMonomial...")\868IN_SET_NS\869D_IN_PAIR_MONOMIAL_MONOMIAL870#define OUT_PAIR_MONOMIAL_MONOMIAL \871D_OUT_PAIR_MONOMIAL_MONOMIAL\872M_OUT_PAIR_MONOMIAL_MONOMIAL873/*--------------------------------------------------------------------*/874#define IN_PAIR_MONOMIAL_SUM \875IN_SET_FUNCTION_NAME("PairMonomialSum...")\876IN_SET_NS\877D_IN_PAIR_MONOMIAL_SUM\878M_IN_PAIR_MONOMIAL_SUM879#define OUT_PAIR_MONOMIAL_SUM \880D_OUT_PAIR_MONOMIAL_SUM\881M_OUT_PAIR_MONOMIAL_SUM882/*--------------------------------------------------------------------*/883#define IN_PAIR_SUM_MONOMIAL \884IN_SET_FUNCTION_NAME("PairSumMonomial...")\885IN_SET_NS\886D_IN_PAIR_SUM_MONOMIAL\887M_IN_PAIR_SUM_MONOMIAL888#define OUT_PAIR_SUM_MONOMIAL \889D_OUT_PAIR_SUM_MONOMIAL\890M_OUT_PAIR_SUM_MONOMIAL891/*--------------------------------------------------------------------*/892#define IN_PAIR_SUM_SUM \893IN_SET_FUNCTION_NAME("PairSumSum...")\894IN_SET_NS\895D_IN_PAIR_SUM_SUM\896M_IN_PAIR_SUM_SUM897#define OUT_PAIR_SUM_SUM \898D_OUT_PAIR_SUM_SUM\899M_OUT_PAIR_SUM_SUM900/*--------------------------------------------------------------------*/901#define IN_POLY_CONTENT \902IN_SET_FUNCTION_NAME("PolyContent")\903IN_SET_NS\904D_IN_POLY_CONTENT\905M_IN_POLY_CONTENT906#define OUT_POLY_CONTENT \907D_OUT_POLY_CONTENT\908M_OUT_POLY_CONTENT909/*--------------------------------------------------------------------*/910#define IN_POLY_GCD \911IN_SET_FUNCTION_NAME("PolyGCD")\912IN_SET_NS\913D_IN_POLY_GCD\914M_IN_POLY_GCD915#define OUT_POLY_GCD \916D_OUT_POLY_GCD\917M_OUT_POLY_GCD918/*--------------------------------------------------------------------*/919#define IN_POLY_QUOTIENT \920IN_SET_FUNCTION_NAME("PolyQuotient")\921IN_SET_NS\922D_IN_POLY_QUOTIENT\923M_IN_POLY_QUOTIENT924#define OUT_POLY_QUOTIENT \925D_OUT_POLY_QUOTIENT\926M_OUT_POLY_QUOTIENT927/*--------------------------------------------------------------------*/928#define IN_POLY_PSEUDO_REMAINDER \929IN_SET_FUNCTION_NAME("PolyPseudoRemainder")\930IN_SET_NS\931D_IN_POLY_PSEUDO_REMAINDER\932M_IN_POLY_PSEUDO_REMAINDER933#define OUT_POLY_PSEUDO_REMAINDER \934D_OUT_POLY_PSEUDO_REMAINDER\935M_OUT_POLY_PSEUDO_REMAINDER936/*--------------------------------------------------------------------*/937#define IN_SCALAR_SUM_CANCELLATION \938IN_SET_FUNCTION_NAME("ScalarSumCancellation")\939IN_SET_NS\940D_IN_SCALAR_SUM_CANCELLATION\941M_IN_SCALAR_SUM_CANCELLATION942#define OUT_SCALAR_SUM_CANCELLATION \943D_OUT_SCALAR_SUM_CANCELLATION\944M_OUT_SCALAR_SUM_CANCELLATION945/*--------------------------------------------------------------------*/946#define IN_SUBSTITUTE_RELATION_IN_RELATION \947IN_SET_FUNCTION_NAME("SubstituteRelationInRelation...")\948IN_SET_NS\949D_IN_SUBSTITUTE_RELATION_IN_RELATION\950M_IN_SUBSTITUTE_RELATION_IN_RELATION951#define OUT_SUBSTITUTE_RELATION_IN_RELATION \952D_OUT_SUBSTITUTE_RELATION_IN_RELATION\953M_OUT_SUBSTITUTE_RELATION_IN_RELATION954/* Conditional print of Relation and Monomial tables */955956#if defined(D_PUT_RELATIONS) /*-------------------------------------*/957#define IN_PUT_RELATIONS PutDebugRelations(),958#define OUT_PUT_RELATIONS D_CONDITION PutDebugRelations();959#else960#define IN_PUT_RELATIONS961#define OUT_PUT_RELATIONS962#endif963964#if defined(D_PUT_LIE_MONOMIAL) /*-------------------------------------*/965#define IN_PUT_LIE_MONOMIAL PutDebugLieMonomialTable(-1),966#define OUT_PUT_LIE_MONOMIAL D_CONDITION PutDebugLieMonomialTable(-1);967#else968#define IN_PUT_LIE_MONOMIAL969#define OUT_PUT_LIE_MONOMIAL970#endif971972#if defined(D_LIE_SUM_DIV_INTEGER) || defined(MEMORY)973#define LS_B_LSUM_DIV_INT uint b = lsum; /* For "LieSumDivInteger" */974#else975#define LS_B_LSUM_DIV_INT976#endif977#if defined(D_LIE_SUM_MULT_INTEGER) || defined(MEMORY)978#define LS_B_LSUM_MULT_INT uint b = lsum; /* For "LieSumMultInteger" */979#else980#define LS_B_LSUM_MULT_INT981#endif982#if defined(D_LIE_SUM_DIV_SCALAR_SUM) || defined(MEMORY)983#define LS_B_LSUM_DIV_SS uint b = lsum; /* For "LieSumDivScalarSum" */984#else985#define LS_B_LSUM_DIV_SS986#endif987#if defined(D_LIE_SUM_MULT_SCALAR_SUM) || defined(MEMORY)988#define LS_B_LSUM_MULT_SS uint b = lsum; /* For "LieSumMultScalarSum" */989#else990#define LS_B_LSUM_MULT_SS991#endif992993/* Test of char functions ?? */994//#define TEST_FUNCTION995996#if defined(INTEGER_ALLOCATION_CHECK)997#define INTEGER_IN_STACK(n) ;if((n)==NULL) Error(E_A_STACK_INTEGER)998#define INTEGER_IN_HEAP(n) ;if((n)==NULL) Error(E_A_HEAP_INTEGER) \999PP_CURRENT_N_INT1000#else1001#define INTEGER_IN_STACK(n)1002#define INTEGER_IN_HEAP(n) PP_CURRENT_N_INT1003#endif10041005#if defined(POLY_ARRAY_ALLOCATION_CHECK)1006#define POLY_ARRAY_IN_STACK(a) ;if((a)==NULL) Error(E_A_STACK_POLY_ARRAY)1007#else1008#define POLY_ARRAY_IN_STACK(a)1009#endif10101011/*_1 System constants=================================================*/10121013#define NIL (0u)1014#define NOTHING (~NIL)10151016#define INTEGER_SIGN_MASK ((LIMB)0x8000)1017/* 1 << 15 or 0000 0000 0000 0001 */1018#define INTEGER_N_LIMBS_MASK ((LIMB)(~INTEGER_SIGN_MASK))1019/* 1111 1111 1111 1110 */1020#define BITS_PER_LIMB (16)1021#define BASE_LIMB (0x10000lu)1022#define MAX_LIMB (0xFFFFu)1023#define FIRST_GENUINE_PARAMETER 1 /* Skip i number */10241025/*_2 Type definitions================================================*/10261027typedef unsigned char byte;1028typedef unsigned short LIMB; /* Limb of big integer */1029typedef LIMB * BIGINT; /* (Pointer to) big integer array */1030typedef unsigned int uint;10311032/* Element of LieMonomial table */10331034typedef struct1035{1036int order; /* Order of this position element */1037int position; /* Position of this order */1038int left; /* (Position of) left submonomial of Lie bracket */1039/* or ~(index of generator) if not commutator */1040int right; /* (Position of) right submonomial of Lie bracket */1041/* or 1 for 1st generator and 0 for nexts */1042struct1043{1044uint parity : 1;1045int index : 23; /* Index if basis element or */1046/* ~(index of relation with leading ordinal) */1047/* Number of relations < 4194303 */1048uint weight : 8; /* Weight of Lie monomial < 256 */1049} info;1050} LIE_MON;10511052/* Element of NodeLT pool for Lie terms */10531054typedef struct1055{1056int monomial; /* (Position of) Lie monomial in LieMonomial table */1057union1058{1059uint scalar_sum;1060BIGINT integer;1061} numerator;1062union1063{1064uint scalar_sum;1065BIGINT integer;1066} denominator; /* Pointer to next Lie term */1067uint rptr;1068} NODE_LT;10691070/* Element of NodeSF pool for scalar factors */10711072typedef struct1073{1074#if defined(SPP_2000)1075byte parameter; /* Parameter ordinal */1076byte degree; /* Degree of parameter */1077#else1078byte degree; /* Degree of parameter */1079byte parameter; /* Parameter ordinal */1080#endif1081uint rptr; /* Pointer to next parametric factor */1082} NODE_SF;10831084/* Element of NodeST pool for scalar terms */10851086typedef struct1087{1088uint monomial; /* Scalar monomial */1089BIGINT numerator; /* Integer coefficient */1090uint rptr; /* Pointer to next parametric term */1091} NODE_ST;10921093/* Element of Relation table */10941095typedef struct1096{1097uint lie_sum; /* Expression of relation (Lie sum) */1098byte min_generator; /* Minimal generator for differentiation */1099byte to_be_substituted; /* YES if relation must be substituted */1100} REL; /* into higher relations */11011102/*_3 Constants and enumerations======================================*/11031104/* Enumeration constants */11051106enum boolean {NO = 0, YES = 1};1107enum signs {PLUS = 0, MINUS = 1};1108enum parity {EVEN = 0, ODD = 1};1109enum orders {ORDER12 = -1, ORDER11 = 0, ORDER21 = 1};1110enum scalar_types /* Types of scalar factors */1111{1112I_NUMBER = 01113};1114enum init_file_cases /* Initiating file fplsa4.ini (fplsa416.ini) */1115{1116COEFFICIENT_SUM_TABLE_SIZE = 0,1117CRUDE_TIME = 1,1118ECHO_INPUT_FILE = 2,1119EVEN_BASIS_SYMBOL = 3,1120GAP_ALGEBRA_NAME = 4,1121GAP_BASIS_NAME = 5,1122GAP_RELATIONS_NAME = 6,1123GAP_OUTPUT_BASIS = 7,1124GAP_OUTPUT_COMMUTATORS = 8,1125GAP_OUTPUT_RELATIONS = 9,1126GENERATOR_MAX_N = 10,1127INPUT_DIRECTORY = 11,1128INPUT_INTEGER_SIZE = 12,1129INPUT_STRING_SIZE = 13,1130LEFT_NORMED_OUTPUT = 14,1131LIE_MONOMIAL_SIZE = 15,1132LINE_LENGTH = 16,1133NAME_LENGTH = 17,1134NODE_LT_SIZE = 18,1135NODE_SF_SIZE = 19,1136NODE_ST_SIZE = 20,1137ODD_BASIS_SYMBOL = 21,1138OUT_LINE_SIZE = 22,1139PARAMETER_MAX_N = 23,1140PUT_BASIS_ELEMENTS = 24,1141PUT_COMMUTATORS = 25,1142PUT_HILBERT_SERIES = 26,1143PUT_INITIAL_RELATIONS = 27,1144PUT_NON_ZERO_COEFFICIENTS = 28,1145PUT_PROGRAM_HEADING = 29,1146PUT_REDUCED_RELATIONS = 30,1147PUT_STATISTICS = 31,1148RELATION_SIZE = 32,1149N_INIT_CASES = 331150};1151enum input_file_cases /* Items of input files */1152{1153GENERATORS = 0,1154LIMITING_WEIGHT = 1,1155PARAMETERS = 2,1156RELATIONS = 3,1157WEIGHTS = 4,1158N_INPUT_CASES = 51159};1160enum messages1161{1162/* Head messages */11631164H_PROGRAM = 0,1165H_ENTER_FILE = 1,1166H_INPUT_FILE = 2,1167H_CREATE_NEW_FILE = 3,1168H_ENTER_GENERATORS = 4,1169H_ENTER_WEIGHTS_IN_FILE = 5,1170H_ENTER_LIMITING_WEIGHT = 6,1171H_ENTER_PARAMETERS = 7,1172H_ENTER_RELATIONS = 8,1173H_SHOW_INPUT = 9,1174H_IN_RELATIONS = 10,1175H_NON_ZERO_COEFFICIENTS = 11,1176H_REDUCED_RELATIONS = 12,1177H_BASIS_ELEMENTS = 13,1178H_HILBERT_SERIES = 14,1179H_COMMUTATORS = 15,1180H_NO_PUT_COMMUTATORS = 16,11811182/* Error messages */11831184ERROR = 17,1185E_WRONG_INI_CASE = 18,1186E_WRONG_INPUT_CASE = 19,1187E_CANCEL_PROGRAM = 20,1188E_UNEXPECTED_EOF = 21,1189E_A_GENERATOR_NAME = 22,1190E_A_PARAMETER_NAME = 23,1191E_A_OUT_LINE = 24,1192E_A_RELATION = 25,1193E_A_LIE_MONOMIAL = 26,1194E_A_NODE_LT = 27,1195E_A_NODE_ST = 28,1196E_A_NODE_SF = 29,1197E_A_HEAP_INTEGER = 30,1198E_A_COEFF_PARA_TABLE = 31,1199E_A_COEFF_SUM_TABLE = 32,1200E_ALLOC = 33,1201E_A_STACK_INPUT_STRING = 34,1202E_A_STACK_INTEGER = 35,1203E_A_STACK_INTEGER_DECIMAL_STRING = 36,1204E_A_STACK_POLY_ARRAY = 37,1205E_INPUT_STRING_SIZE = 38,1206E_OUT_LINE_SIZE = 39,1207E_LIE_MONOMIAL_SIZE = 40,1208E_RELATION_SIZE = 41,1209E_NODE_LT_SIZE = 42,1210E_NODE_SF_SIZE = 43,1211E_NODE_ST_SIZE = 44,1212E_COEFF_SUM_TABLE_SIZE = 45,1213E_GENERATOR_MAX_N = 46,1214E_PARAMETER_MAX_N = 47,1215E_TOO_MUCH_INPUT_WEIGHTS = 48,1216E_NON_NUM_INPUT_WEIGHT = 49,1217E_NO_R_PARENTHESIS = 50,1218E_NO_GENERAL_POWER = 51,1219E_UNDECLARED_GENERATOR = 52,1220E_NO_COMMUTATOR_COMMA = 53,1221E_NO_COMMUTATOR_BRACKET = 54,1222E_INVALID_CHARACTER = 55,1223E_MESSAGE = 561224};12251226/* Constants for input and output */12271228#define LEFT_COMMENT '<' /* In *.ini and input files */1229#define RIGHT_COMMENT '>'1230#define SUBSCRIPT_INPUT_SIGN '_'1231#define ODD_GENERATOR_INPUT_SIGN '-' /* At input */12321233#define LEVEL '\xFF'1234#define MAIN_LEVEL '\x2'1235#define MARGIN (LEVEL-1)12361237#define CASE_STRING_SIZE 256 /* Size of string to match case */1238#define GAP_NAME_SIZE 64 /* Including ending '\0' */1239#define GAP_WIDTH 79 /* Width of GAP page */12401241/*_4 Macrodefinitions==================================================*/1242#define COUNT_LEADING_ZERO_BITS_IN_LIMB(n,w) n=CountLeadingZeroBitsInLimb((LIMB)(w))1243#if 0 /* ?? Exclude the macro: the above function is slightly faster */1244#define COUNT_LEADING_ZERO_BITS_IN_LIMB(n, w) (n) = (w);\1245if((n) >= 0x100) \1246if((n) >= 0x1000) \1247if((n) >= 0x4000) \1248if((n) >= 0x8000) \1249(n) = 0; \1250else \1251(n) = 1; \1252else \1253if((n) >= 0x2000) \1254(n) = 2; \1255else \1256(n) = 3; \1257else \1258if((n) >= 0x400) \1259if((n) >= 0x800) \1260(n) = 4; \1261else \1262(n) = 5; \1263else \1264if((n) >= 0x200) \1265(n) = 6; \1266else \1267(n) = 7; \1268else \1269if((n) >= 0x10) \1270if((n) >= 0x40) \1271if((n) >= 0x80) \1272(n) = 8; \1273else \1274(n) = 9; \1275else \1276if((n) >= 0x20) \1277(n) = 10; \1278else \1279(n) = 11; \1280else \1281if((n) >= 0x4) \1282if((n) >= 0x8) \1283(n) = 12; \1284else \1285(n) = 13; \1286else \1287if((n) >= 0x2) \1288(n) = 14; \1289else \1290if((n)) \1291(n) = 15; \1292else \1293(n) = 161294#endif12951296#define CUT_ARRAY(arr, type, n) (arr)=(type *)realloc(arr,sizeof(type)*(n))12971298#define EXIT do { TIME_OFF; PutStatistics(); exit(1); } while(0)12991300#define IN_LINE_MARGIN OutLine[++PosOutLine]=MARGIN13011302#define INTEGER_MINUS(ia) do { if(INTEGER_IS_NEGATIVE(ia))\1303(ia)[0] &= INTEGER_N_LIMBS_MASK;\1304else\1305(ia)[0] |= INTEGER_SIGN_MASK; } while(0)1306#define INTEGER_IS_NEGATIVE(ia) (((ia)[0]&INTEGER_SIGN_MASK)!=0)1307#define INTEGER_IS_POSITIVE(ia) (((ia)[0]&INTEGER_SIGN_MASK)==0)13081309#define INTEGER_IS_UNIT(ia) (((ia)[0]==1) && ((ia)[1]==1))1310#define INTEGER_IS_UNIT_ABS(ia) ((((ia)[0]&INTEGER_N_LIMBS_MASK)==1) &&\1311((ia)[1]==1))1312#define INTEGER_IS_NOT_UNIT(ia) (((ia)[0]!=1) || ((ia)[1]!=1))1313#define INTEGER_IS_NOT_UNIT_ABS(ia) ((((ia)[0]&INTEGER_N_LIMBS_MASK)!=1) ||\1314((ia)[1]!=1))1315#define INTEGER_N_LIMBS(ia) ((ia)[0]&INTEGER_N_LIMBS_MASK)13161317#define INTEGER_SET_MINUS(ia) (ia)[0] |= INTEGER_SIGN_MASK1318#define INTEGER_SET_PLUS(ia) (ia)[0] &= INTEGER_N_LIMBS_MASK1319#define INTEGER_SIGN(ia) ((ia)[0]&INTEGER_SIGN_MASK)13201321#define INTEGER_HEAP_NEW(n,i) (n)=(BIGINT)malloc(sizeof(LIMB)*(i))\1322INTEGER_IN_HEAP(n)1323#define INTEGER_HEAP_COPY(n,o,i) (i)=INTEGER_N_LIMBS(o);\1324INTEGER_HEAP_NEW(n,++(i));\1325do{(i)--; (n)[i] = (o)[i];}while(i)1326#define INTEGER_HEAP_COPY_DOUBLE_1(n1,n2,o,i) \1327(i)=INTEGER_N_LIMBS(o);\1328INTEGER_HEAP_NEW(n1,++(i)+1);\1329INTEGER_HEAP_NEW(n2,i);\1330do{(i)--; (n1)[i] = (n2)[i] = (o)[i];}\1331while(i)1332#define INTEGER_KILL(bn) free(bn) MM_CURRENT_N_INT133313341335#define INTEGER_STACK_NEW(n,i) (n)=(BIGINT)alloca(sizeof(LIMB)*(i))\1336INTEGER_IN_STACK(n)13371338#define INTEGER_STACK_COPY(n,o,i) (i)=INTEGER_N_LIMBS(o);\1339INTEGER_STACK_NEW(n,++(i));\1340do{(i)--; (n)[i] = (o)[i];}while(i)13411342#define INTEGER_STACK_COPY_1(n,o,i) (i)=INTEGER_N_LIMBS(o);\1343INTEGER_STACK_NEW(n,++(i)+1);\1344do{(i)--; (n)[i] = (o)[i];}while(i)13451346#define LIE_MONOMIAL_I_RELATION(pos) (~LIE_MONOMIAL_INDEX(pos))1347#define LIE_MONOMIAL_INDEX_BY_ORDER(ord) \1348LIE_MONOMIAL_INDEX(LIE_MONOMIAL_POSITION(ord))1349#define LIE_MONOMIAL_IS_EVEN(pos) (LIE_MONOMIAL_PARITY(pos)==EVEN)1350#define LIE_MONOMIAL_IS_LEADING_BY_ORDER(ord) \1351LIE_MONOMIAL_IS_LEADING(LIE_MONOMIAL_POSITION(ord))1352#define LIE_MONOMIAL_IS_OCCUPIED(pos) LIE_MONOMIAL_WEIGHT(pos)1353#define LIE_MONOMIAL_IS_GENERATOR(pos) ((pos) < GeneratorN)1354#define LIE_MONOMIAL_IS_GENERATOR_BY_ORDER(ord) \1355(LIE_MONOMIAL_POSITION(ord) < GeneratorN)1356#define LIE_MONOMIAL_IS_COMMUTATOR(pos) ((pos) >= GeneratorN)1357#define LIE_MONOMIAL_IS_ODD(pos) LIE_MONOMIAL_PARITY(pos)1358#define LIE_MONOMIAL_IS_SQUARE(pos) (LIE_MONOMIAL_LEFT(pos)==\1359LIE_MONOMIAL_RIGHT(pos))1360#define LIE_MONOMIAL_IS_NOT_SQUARE(pos) (LIE_MONOMIAL_LEFT(pos)!=\1361LIE_MONOMIAL_RIGHT(pos))1362#define LIE_MONOMIAL_LEFT(pos) (LieMonomial[pos].left)1363#define LIE_MONOMIAL_RIGHT(pos) (LieMonomial[pos].right)1364#define LIE_MONOMIAL_ORDER(pos) (LieMonomial[pos].order)1365#define LIE_MONOMIAL_LEFT_ORDER(pos) \1366LIE_MONOMIAL_ORDER(LIE_MONOMIAL_LEFT(pos))1367#define LIE_MONOMIAL_RIGHT_ORDER(pos) \1368LIE_MONOMIAL_ORDER(LIE_MONOMIAL_RIGHT(pos))1369#define LIE_MONOMIAL_POSITION(ord) (LieMonomial[ord].position)1370#define LIE_MONOMIAL_WEIGHT_BY_ORDER(ord) \1371LIE_MONOMIAL_WEIGHT(LIE_MONOMIAL_POSITION(ord))1372#define LIE_MONOMIAL_INDEX(pos) (LieMonomial[pos].info.index)1373#define LIE_MONOMIAL_IS_BASIS(pos) (LieMonomial[pos].info.index >= 0)1374#define LIE_MONOMIAL_IS_LEADING(pos) (LieMonomial[pos].info.index < 0)1375#define LIE_MONOMIAL_PARITY(pos) (LieMonomial[pos].info.parity)1376#define LIE_MONOMIAL_WEIGHT(pos) (LieMonomial[pos].info.weight)13771378#define LIE_TERM_MONOMIAL(a) (NodeLT[a].monomial)1379#define LIE_TERM_MONOMIAL_ORDER(a) LIE_MONOMIAL_ORDER(\1380(NodeLT[a].monomial))1381#define LIE_TERM_R(a) (NodeLT[a].rptr)1382#define LIE_TERM_NUMERATOR_INTEGER(a) (NodeLT[a].numerator.integer)1383#define LIE_TERM_MINUS_INTEGER(a) \1384INTEGER_MINUS(NodeLT[a].numerator.integer)1385#define LIE_TERM_NUMERATOR_SCALAR_SUM(a) (NodeLT[a].numerator.scalar_sum)1386#define LIE_TERM_DENOMINATOR_INTEGER(a) (NodeLT[a].denominator.integer)1387#define LIE_TERM_DENOMINATOR_SCALAR_SUM(a) (NodeLT[a].denominator.scalar_sum)13881389#define MAX(a,b) (((a) > (b)) ? (a) : (b))13901391#define NODE_LT_KILL(a) LIE_TERM_R(a)=NodeLTTop,NodeLTTop=(a)\1392MM_CURRENT_N_LT1393#define NODE_SF_KILL(a) SCALAR_FACTOR_R(a)=NodeSFTop,NodeSFTop=(a)\1394MM_CURRENT_N_SF1395#define NODE_ST_KILL(a) SCALAR_TERM_R(a)=NodeSTTop,NodeSTTop=(a)\1396MM_CURRENT_N_ST13971398#define RELATION_LIE_SUM(i) (Relation[i].lie_sum)1399#define RELATION_MIN_GENERATOR(i) (Relation[i].min_generator)1400#define RELATION_TO_BE_SUBSTITUTED(i) (Relation[i].to_be_substituted)14011402#define SCALAR_SUM_IS_UNIT(a) (SCALAR_TERM_MONOMIAL(a)==NIL&&\1403INTEGER_IS_UNIT(SCALAR_TERM_NUMERATOR(a)))14041405#define SCALAR_SUM_IS_UNIT_ABS(a) (SCALAR_TERM_MONOMIAL(a)==NIL&&\1406INTEGER_IS_UNIT_ABS(SCALAR_TERM_NUMERATOR(a)))14071408#define SCALAR_SUM_IS_NOT_UNIT(a) (SCALAR_TERM_MONOMIAL(a)!=NIL||\1409INTEGER_IS_NOT_UNIT(\1410SCALAR_TERM_NUMERATOR(a)))14111412#define SCALAR_TERM_MONOMIAL(a) (NodeST[a].monomial)1413#define SCALAR_TERM_NUMERATOR(a) (NodeST[a].numerator)1414#define SCALAR_TERM_R(a) (NodeST[a].rptr)14151416#define SCALAR_FACTOR_PARAMETER(a) (NodeSF[a].parameter)1417#define SCALAR_FACTOR_IS_I_NUMBER(a) (NodeSF[a].parameter==I_NUMBER)1418#define SCALAR_FACTOR_DEGREE(a) (NodeSF[a].degree)1419#define SCALAR_FACTOR_WORD(a) (*(unsigned short *)(NodeSF+(a)))1420#define SCALAR_FACTOR_R(a) (NodeSF[a].rptr)14211422#define SCALAR_TERM_MINUS(a) INTEGER_MINUS(NodeST[a].numerator)1423#define SCALAR_TERM_MAIN_PARAMETER(a) \1424SCALAR_FACTOR_PARAMETER(SCALAR_TERM_MONOMIAL(a))1425#define SCALAR_TERM_MAIN_PARAMETER_WORD(a) \1426SCALAR_FACTOR_WORD(SCALAR_TERM_MONOMIAL(a))1427#define SCALAR_TERM_MAIN_DEGREE(a) \1428SCALAR_FACTOR_DEGREE(SCALAR_TERM_MONOMIAL(a))1429#define POLY_MAIN_PARAMETER(a) ((SCALAR_TERM_MONOMIAL(a)==NIL) ? -1 :\1430SCALAR_FACTOR_PARAMETER(SCALAR_TERM_MONOMIAL(a)))1431#define POLY_ARRAY_STACK_NEW(a,n) (a)=(uint*)alloca(sizeof(uint)*(n))\1432POLY_ARRAY_IN_STACK(a)14331434#define TIME_OFF TimeC += (CrudeTime ? time(NULL) : clock()) - TimeA1435#define TIME_ON TimeA = (CrudeTime ? time(NULL) : clock())14361437/*_5 Global variables and arrays=====================================*/14381439/* Files */14401441#if !defined(GAP)1442FILE *MessageFile;1443FILE *SessionFile;1444#endif14451446/* Single variables */14471448int IncompletedBasis;1449int IncompletedRelations;1450int IsParametric;1451int LieMonomialIsNew;1452int SubstitutionIsDone;1453uint LimitingWeight;1454int GeneratorN;1455int ParameterN;14561457/* Arrays */14581459LIE_MON *LieMonomial; /* Set of Lie monomials */1460int LieMonomialSize;1461int LieMonomialN;1462int LieMonomialFreePosition; /* Start search of free position */1463#if defined(SPACE_STATISTICS)1464int LieMonomialMaxN;1465#endif146614671468NODE_LT *NodeLT; /* Pool of nodes for Lie terms */1469uint NodeLTSize;1470uint NodeLTTop = 1;1471#if defined(SPACE_STATISTICS)1472uint NodeLTTopMax;1473#endif14741475NODE_SF *NodeSF; /* Pool of nodes for scalar factors */1476uint NodeSFSize;1477uint NodeSFTop = 1;1478#if defined(SPACE_STATISTICS)1479uint NodeSFTopMax;1480#endif14811482NODE_ST *NodeST; /* Pool of nodes for scalar terms */1483uint NodeSTSize;1484uint NodeSTTop = 1;1485#if defined(SPACE_STATISTICS)1486uint NodeSTTopMax;1487#endif14881489REL *Relation; /* Ordered set of relations */1490int RelationSize;1491int RelationN;1492#if defined(SPACE_STATISTICS)1493int MaxNRelation;1494#endif1495#if defined(INTEGER_MAX_SIZE)1496LIMB IntegerMaxSize;1497#endif14981499int CoeffSumTableSize; /* Non-zero coefficient tables variables */1500int CoeffSumTableN;1501uint *CoeffSumTable; /* Non-zero parametric sums */1502int *CoeffParamTable; /* Table for memorizing single */1503/* non-zero parameters */15041505/* Input and output variables */15061507char BasisSymbolEven;1508char BasisSymbolOdd;1509int CurrentLevel;1510char * GeneratorName; /* Input names */1511char * ParameterName;1512int GeneratorMaxN; /* Maximum number of input generators */1513uint InputIntegerSize; /* Maximum size of input integer in LIMBs */1514int InputStringSize; /* Size of string for reading input */1515int LastItemEnd;1516int LineLength;1517int Margin;1518int MaxLevel;1519int MinLevel;1520int NameLength1; /* Maximum length of input name (with ending '\0') */1521int NewMargin;1522int ParameterMaxN = 1; /* Maximum number of input parameters (i at least) */1523int BasisElementsPut;1524int CommutatorsPut;1525int CrudeTime; /* Prevent time variable wrapping for large tasks */1526int EchoInput;1527int GAPOutputCommutators;1528int GAPOutputBasis;1529int GAPOutputRelations;1530char GAPAlgebraName[GAP_NAME_SIZE];1531char GAPBasisName[GAP_NAME_SIZE];1532char GAPRelationsName[GAP_NAME_SIZE];1533int HeadingPut;1534int HilbertSeriesPut;1535int InitialRelationsPut;1536int NonZeroCoefficientsPut;1537int ReducedRelationsPut;1538int StatisticsPut;1539char * OutLine; /* String for preparation of output block */1540int OutLineSize;1541int PosOutLine;1542int PreviousEnd;1543int PrintEnd;15441545uint TimeA, TimeC;15461547#if defined(DEBUG)1548uint Debug;1549#endif1550#if defined(MEMORY)1551int CurrentNLT;1552int CurrentNSF;1553int CurrentNST;1554int CurrentNINT;15551556#endif15571558/*_6 Function descriptions===========================================*/15591560/*_6_0 Main and top level functions============================*/15611562void ConstructFreeAlgebraBasis(void); /* 1 call! */1563int FindNewPositionInRelation(int lmo);1564void GenerateRelations(void); /* 1 call! */1565uint NewJacobiRelation(int l); /* 1 call! */1566int ReduceRelations(int i);15671568/*_6_1 Pairing functions=======================================*/15691570int AddPairToLieMonomial(int i, int j);1571uint MakeRelationRHSInteger(int i); /* 1 call!! */1572uint MakeRelationRHSParametric(int i); /* 1 call!! */1573uint PairMonomialMonomialInteger(int i, int j);1574uint PairMonomialMonomialParametric(int i, int j);1575uint PairMonomialSumInteger(int mon, uint a);1576uint PairMonomialSumParametric(int mon, uint a);1577uint PairSumMonomialInteger(uint a, int mon);1578uint PairSumMonomialParametric(uint a, int mon);1579uint PairSumSumInteger(uint a, uint b);1580uint PairSumSumParametric(uint a, uint b);15811582/*_6_2 Substitution (replacing) functions======================*/15831584int IsMonomialInMonomial(int submon, int mon);1585uint SubstituteRelationInRelationInteger(uint r, uint a);1586uint SubstituteRelationInRelationParametric(uint r, uint a);1587uint SubstituteRHSInMonomialInteger(int mon, int lmonr, uint r);1588uint SubstituteRHSInMonomialParametric(int mon, int lmonr, uint r);15891590/*_6_3 Lie and scalar algebra functions========================*/15911592int LieLikeTermsCollectionInteger(uint a, uint b);1593int LieLikeTermsCollectionParametric(uint a, uint b);1594uint LieSumAddition(uint a, uint b);1595void LieSumDivInteger(uint lsum, BIGINT den);1596void LieSumDivScalarSum(uint lsum, uint den);1597void LieSumMinusInteger(uint a);1598void LieSumMinusParametric(uint a);1599void LieSumMultInteger(uint lsum, BIGINT num);1600#if defined(RATIONAL_FIELD)1601void LieSumMultRationalInteger(int a, BIGINT num, BIGINT den);1602#endif1603void LieSumMultScalarSum(uint lsum, uint num);1604void NormalizeRelationInteger(uint a);1605void NormalizeRelationParametric(uint a);1606uint ScalarMonomialMultiplication(int *pchange_sign, uint ma, uint mb);1607uint ScalarSumAddition(uint a, uint b);1608void ScalarSumCancellation(uint *pnum, uint *pden);1609void ScalarSumMinus(uint a);1610uint ScalarSumMultiplication(uint a, uint b);1611void ScalarTermMultiplication(uint a, uint b); /* 1 call! */16121613/*_6_4 Scalar polynomial algebraic functions===================*/16141615uint ContentOfScalarSum(uint cont, uint a);1616void InCoeffParamTable(uint cont);1617void InCoeffSumTable(uint sum);1618void InCoeffTable(uint coe);1619uint PolyCoeffAtMainParameter(uint *pa, int mp);1620uint PolyContent(uint a, int mp);1621uint PolyGCD(uint a, uint b);1622uint PolyMainParameterTerm(uint *pa, int mp, int mpdeg);1623int PolynomialsAreEqual(uint a, uint b);1624uint PolyPseudoRemainder(uint a, uint b, int mp); /* 1 call! */1625uint PolyTermGCD(uint a, uint b);1626void PolyTermQuotient(uint a, uint b);1627uint PolyQuotient(uint a, uint b);16281629/*_6_5 Big number functions====================================*/16301631int BigNMinusBigN(BIGINT a, int na, BIGINT b, int nb);1632LIMB BigNShiftLeft(BIGINT bign, int n, int cnt);1633int BigNShiftRight(BIGINT bign, int n, int cnt);1634int CountLeadingZeroBitsInLimb(LIMB w);1635void IntegerCancellation(BIGINT num, BIGINT den);1636BIGINT IntegerGCD(BIGINT u, BIGINT v);1637void IntegerProduct(BIGINT w, BIGINT u, BIGINT v);1638void IntegerQuotient(BIGINT c, BIGINT a, BIGINT b);1639void IntegerSum(BIGINT c, BIGINT a, BIGINT b);16401641/*_6_6 Copy and delete functions===============================*/16421643uint LieSumCopyInteger(uint a);1644uint LieSumCopyIntegerNegative(uint a);1645uint LieSumCopyParametric(uint a);1646void LieSumKillInteger(uint a);1647void LieSumKillParametric(uint a);1648uint LieTermFromMonomialInteger(int mon);1649uint LieTermFromMonomialParametric(int mon);1650uint ScalarSumCopy(uint a);1651void ScalarSumKill(uint a);1652uint ScalarTermCopy(uint a);16531654/*_6_7 Technical functions=====================================*/16551656void Error(int i_message) ATTRIBUTE_NORETURN;1657void Initialization(void);1658void *NewArray(uint n, uint size, int i_message);1659uint NodeLTNew(void);1660uint NodeSFNew(void);1661uint NodeSTNew(void);1662FILE *OpenFile(char * file_name, char * file_type);16631664/*_6_8 Input functions=========================================*/16651666int BinaryQuestion(int i_message);1667int FindNameInTable(char * name, char * nametab, int n_nametab);1668void GetGenerator(char * str);1669void GetInput(int n, char * fin);1670void GetInteger(BIGINT a, char **pstr);1671uint GetLieMonomial(char **pstr);1672uint GetLieSum(char **pstr);1673uint GetLieTerm(char **pstr);1674uint GetUInteger(char **pstr);1675void GetParameter(char * str);1676void GetRelation(char * str);1677uint GetScalarSum(char **pstr);1678uint GetScalarTerm(char **pstr);1679void GetWeight(char * str);1680int KeyBoardBytesToString(char * str);1681int KeyBoardStringToFile(int i_m, char * prefix, char * str, FILE *file);1682void ReadAndProcessStringsFromFile(void (*proc_func)(char * str), FILE *inf,1683char sep, char end);1684int ReadBooleanFromFile(FILE *file);1685int ReadCaseFromFile(FILE * file, char * case_str[], int n_cases);1686uint ReadDecimalFromFile(FILE *file);1687short ReadStringFromFile(char * str, FILE *file);1688short SkipCommentInFile(FILE *file);1689void SkipName(char **pstr);1690void SkipSpaces(char **pstr);1691short SkipSpacesInFile(FILE *file);16921693/*_6_9 Output functions========================================*/16941695void AddSymbolToOutLine(char c, int position);1696void InLineLevel(int level);1697void InLineNumberInBrackets(uint n);1698void InLineString(char * str);1699void InLineSubscript(char * s);1700void InLineSymbol(char symbol);1701void InLineTableName(char * name);1702char * UToString(uint n);1703void PutBasis(void);1704#if defined(GAP)1705void PutBasisGAP(void);1706#endif1707void PutBlock(void);1708void PutCharacter(char c);1709#if defined(GAP)1710void PutCharacterGAP(char c);1711#endif1712void PutCoefficientTable(void);1713void PutCommutators(void);1714#if defined(GAP)1715void PutCommutatorsGAP(void);1716#endif1717void PutDegree(uint deg);1718void PutDimensions(void);1719void PutDots(void);1720void PutEnd(void);1721void PutFormattedU(char * format, uint i);1722void PutIntegerUnsigned(BIGINT bn);1723#if defined(GAP)1724void PutIntegerUnsignedGAP(BIGINT bn);1725#endif1726void PutLieBareTerm(void (*put_lie_mon)(int a), uint a);1727void PutLieBasisElement(int pos);1728void PutLieMonomialLeftNormed(int pos);1729void PutLieMonomialStandard(int pos);1730#if defined(GAP)1731void PutLieMonomialGAP(int pos);1732#endif1733void PutLieSum(void (*put_lie_mon)(int a), uint a);1734void PutMessage(int i_message);1735void PutRelations(int i);1736#if defined(GAP)1737void PutRelationsGAP(void);1738#endif1739void PutScalarBareTerm(uint a);1740void PutScalarFactor(uint a);1741void PutScalarSum(uint a);1742void PutStart(void);1743void PutStatistics(void);1744void PutString(char * str);1745#if defined(GAP)1746void PutStringGAP(char * str);1747#endif1748void PutStringStandard(char * str);1749void PutSymbol(char c);17501751/* Global function variables */17521753int (*LieLikeTermsCollection)(uint a, uint b) = LieLikeTermsCollectionInteger;1754uint (*LieSumCopy)(uint a) = LieSumCopyInteger;1755void (*LieSumKill)(uint a) = LieSumKillInteger;1756void (*LieSumMinus)(uint a) = LieSumMinusInteger;1757uint (*LieTermFromMonomial)(int mon) = LieTermFromMonomialInteger;1758void (*NormalizeRelation)(uint a) = NormalizeRelationInteger;1759uint (*PairMonomialMonomial)(int i, int j) = PairMonomialMonomialInteger;1760uint (*PairMonomialSum)(int mon, uint a) = PairMonomialSumInteger;1761uint (*PairSumMonomial)(uint a, int mon) = PairSumMonomialInteger;1762uint (*PairSumSum)(uint a, uint b) = PairSumSumInteger;1763void (*PutLieMonomial)(int pos) = PutLieMonomialStandard;1764uint (*SubstituteRelationInRelation)(uint r, uint a) =1765SubstituteRelationInRelationInteger;17661767/*_6_10 Debugging functions===========================================*/17681769#if defined(DEBUG)1770void PutDebugHeader(uint debug, char * f_name, char * in_out);1771void PutDebugInteger(char * name, BIGINT u);1772void PutDebugLieMonomial(char * name, int a);1773void PutDebugLieMonomialTable(int newmon);1774void PutDebugLieSum(char * name, uint a);1775void PutDebugLieTerm(char * name, uint a);1776void PutDebugU(char * name, uint i);1777#if defined(D_PUT_RELATIONS)1778void PutDebugRelations(void);1779#endif1780void PutDebugScalarSum(char * name, uint a);1781void PutDebugString(char * strname, char * str);1782#endif1783#if defined(MEMORY)1784void AddLieSumNs(uint a, int minus_or_plus,1785int *pn_lt, int *pn_int, int *pn_st, int *pn_sf);1786void AddScalarSumNs(uint a, int minus_or_plus, int *pn_int, int *pn_st, int *pn_sf);1787void PutIntegerBalance(char * fname, int dn);1788void PutNodeBalance(char * type, char * fname, int dn);1789#endif17901791/*_6_0 Main and top level functions============================*/17921793#if !defined(TEST_FUNCTION)1794/*=main======================================1795*/1796int main(int narg, char ** fin)1797{1798Initialization();1799GetInput(narg, fin[1]);1800if(RelationN)1801{1802if(InitialRelationsPut)1803PutRelations(H_IN_RELATIONS);1804GenerateRelations();1805if(NonZeroCoefficientsPut)1806PutCoefficientTable();1807}1808if(RelationN)1809{1810if(ReducedRelationsPut)1811PutRelations(H_REDUCED_RELATIONS);1812}1813else /* Free algebra */1814ConstructFreeAlgebraBasis();1815PutBasis();1816if(HilbertSeriesPut)1817PutDimensions();1818if(CommutatorsPut)1819PutCommutators();1820#if defined(DEBUG)1821PutDebugU("Top Debug", Debug);1822#endif1823if(StatisticsPut)1824{1825TIME_OFF;1826PutStatistics();1827}1828#if defined(GAP)1829if(!IsParametric && !IncompletedRelations && !IncompletedBasis)1830{1831/* if(GAPOutputCommutators || GAPOutputBasis || GAPOutputRelations) */1832/* { */1833/* fclose(SessionFile); */1834/*#if defined(SPP_2000) */1835/* SessionFile = OpenFile("fplsa4.gap", "w"); */1836/*#else */1837/* SessionFile = OpenFile("fplsa4.gap", "wt"); */1838/*#endif */1839/* } */1840if(GAPOutputCommutators)1841PutCommutatorsGAP();1842if(GAPOutputBasis)1843PutBasisGAP();1844if(GAPOutputRelations)1845PutRelationsGAP();1846}1847#endif1848exit(0);1849}1850#endif1851/*=ConstructFreeAlgebraBasis==========================================1852Make all regular monomials up to LimitingWeight1853*/1854void ConstructFreeAlgebraBasis(void)1855{1856int i = 0, j, moni, monj;1857while(YES)1858{1859moni = LIE_MONOMIAL_POSITION(i);1860if(LIE_MONOMIAL_IS_NOT_SQUARE(moni)) /* ?? Not so for non-standard */1861{1862j = 0;1863while(j < i)1864{1865monj = LIE_MONOMIAL_POSITION(j);1866if(LIE_MONOMIAL_IS_NOT_SQUARE(monj))1867if(LIE_MONOMIAL_IS_GENERATOR(moni) ||1868j >= LIE_MONOMIAL_RIGHT_ORDER(moni))1869{1870if(LIE_MONOMIAL_WEIGHT(moni) + LIE_MONOMIAL_WEIGHT(0) >1871LimitingWeight) /* Out of weight for all consequent i & j */1872goto incompleted;1873if(LIE_MONOMIAL_WEIGHT(moni) + LIE_MONOMIAL_WEIGHT(monj) >1874LimitingWeight) /* Out of weight for all consequent j */1875goto next_i;1876AddPairToLieMonomial(moni, monj);1877}1878j++;1879}1880if(LIE_MONOMIAL_IS_ODD(moni))1881{ /* Add square */1882if(LIE_MONOMIAL_WEIGHT(moni) + LIE_MONOMIAL_WEIGHT(0) >1883LimitingWeight) /* Out of weight for all consequent i & j */1884goto incompleted;1885if(2*LIE_MONOMIAL_WEIGHT(moni) > LimitingWeight)1886goto next_i; /* Out of weight for all consequent j */1887AddPairToLieMonomial(moni, moni);1888}1889}1890next_i:1891if(++i >= LieMonomialN)1892return; /* One generator case */1893}1894incompleted:1895IncompletedBasis = YES;1896}1897/*=FindNewPositionInRelation======================================1898Find position of first relation with leading monomial order > lmo1899among 0, 1,..., RelationN - 11900*/1901int FindNewPositionInRelation(int lmo)1902{1903int left = 0;1904if(RelationN) /* Binary search */1905{ /* `right' must be of signed type */1906int m, right = RelationN-1;1907do1908{1909m = (left + right)/2;1910if(lmo < LIE_TERM_MONOMIAL_ORDER(RELATION_LIE_SUM(m)))1911right = --m;1912else1913left = ++m;1914}while(left <= right);1915}1916return left;1917}1918/*=GenerateRelations=======================================================1919Generate and process new relations1920*/1921void GenerateRelations(void)1922{1923int i, k, l, mon, mona,1924gen; /* LieMonomial table position of differentiating generator */1925uint a, lim_weight_i;1926i = 0;1927while(i < RelationN) /* Differentiation loop */1928if(RELATION_MIN_GENERATOR(i) < GeneratorN)1929if(LIE_MONOMIAL_IS_BASIS(gen = RELATION_MIN_GENERATOR(i)))1930{1931/* Program assures the ban of differentiation OF leading generators1932in the process of new relation adding by setting negation of the1933first IF predicate */1934a = RELATION_LIE_SUM(i);1935mona = LIE_TERM_MONOMIAL(a); /* Irregular triple criterion */1936if(LIE_MONOMIAL_IS_SQUARE(mona) ||1937LIE_MONOMIAL_ORDER(gen) < LIE_MONOMIAL_RIGHT_ORDER(mona))1938{1939IN_GENERATE_RELATIONS /*------------------------------------------------*/1940if(LIE_MONOMIAL_WEIGHT(gen) + /* Out of weight */1941LIE_MONOMIAL_WEIGHT(LIE_TERM_MONOMIAL(a)) > LimitingWeight)1942{1943RELATION_MIN_GENERATOR(i++) = GeneratorN;1944continue; /* There might be lower weight next generators */1945}1946RELATION_MIN_GENERATOR(i)++; /* For next differentiation */1947if((a = (*PairSumMonomial)((*LieSumCopy)(a), gen)) != NIL)1948{1949add_new_relation:1950if(RelationN == RelationSize)1951Error(E_RELATION_SIZE);1952gen = LIE_TERM_MONOMIAL(a);1953l = FindNewPositionInRelation(k = LIE_MONOMIAL_ORDER(gen));1954#if defined(SPACE_STATISTICS)1955if(RelationN >= MaxNRelation)1956MaxNRelation = RelationN + 1;1957#endif1958(*NormalizeRelation)(a);1959OUT_GENERATE_RELATIONS /*------------------------------------------------*/1960LIE_MONOMIAL_INDEX(gen) = ~l; /* Set position of relation */19611962/* Shift positions of higher relations in LieMonomial */19631964while(++k < LieMonomialN)1965if(LIE_MONOMIAL_IS_LEADING_BY_ORDER(k))1966--LIE_MONOMIAL_INDEX_BY_ORDER(k);19671968/* Make room for new relation */19691970for(k = RelationN; k > l; k--)1971Relation[k] = Relation[k-1];19721973if(LIE_MONOMIAL_IS_GENERATOR(gen)) /* Ban differentiating */1974RELATION_MIN_GENERATOR(l) = GeneratorN; /* lead. generat. */1975else1976{1977RELATION_MIN_GENERATOR(l) = 0;1978if(l <= i) /* Shift min. diff. index */1979i = l;1980}1981RELATION_LIE_SUM(l) = a; /* Set new relation */1982if(LIE_MONOMIAL_WEIGHT(gen) + 1 > LimitingWeight)1983RELATION_MIN_GENERATOR(l) = GeneratorN;1984if(l == RelationN++)1985RELATION_TO_BE_SUBSTITUTED(l) = NO;1986else1987{ /* New relation inside the table */1988RELATION_TO_BE_SUBSTITUTED(l) = YES;1989l = ReduceRelations(l);1990if(l <= i)1991i = l;1992}1993/* #if defined(RELATION_N_TO_SCREEN)1994TIME_OFF;1995printf("\n%10d", RelationN);1996TIME_ON;1997#endif*/1998OUT_PUT_LIE_MONOMIAL /*--------------------------------------------------*/1999OUT_PUT_RELATIONS /*-----------------------------------------------------*/2000}2001}2002else /* Any next generator >= right of leading */2003RELATION_MIN_GENERATOR(i++) = GeneratorN;2004}2005else /* Skip differentiation BY leading generator */2006RELATION_MIN_GENERATOR(i)++;2007else2008i++; /* Skip completely differentiated relation */2009IncompletedBasis = /* Limiting weight is reached */2010IncompletedRelations =2011(LIE_MONOMIAL_WEIGHT(LIE_TERM_MONOMIAL(RELATION_LIE_SUM(RelationN-1)))2012>= LimitingWeight);20132014/*??#if 0 vvvvvvvvvvvvvvvvv Off checking vvvvvvvvvvvvvvvvvvvvvvvvvvvvv*/2015/* Check regular pairs */20162017k = 0;2018while(k < LieMonomialN)2019{2020gen = LIE_MONOMIAL_POSITION(k);2021if(LIE_MONOMIAL_IS_BASIS(gen))2022{2023#if 0 /*?? Experiment with individual basis elements vvvvvvvvvvvvvvvv */2024if(LIE_MONOMIAL_IS_COMMUTATOR(gen))2025{2026/* Old pairs */2027a = NewJacobiRelation(gen);2028if(a != NIL)2029/*??*/{2030/*??*/PutDebugLieSum("\n***New Jacobi Relation from Old Basis", a);2031goto add_new_relation;2032/*??*/}2033}2034#endif /*?? Experiment ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */2035if(LIE_MONOMIAL_WEIGHT(gen) + LIE_MONOMIAL_WEIGHT(0)2036> LimitingWeight)2037{2038IncompletedBasis = YES;2039goto loops_out;2040}2041if(LIE_MONOMIAL_IS_NOT_SQUARE(gen)) /* ?? non-standard square */2042{2043lim_weight_i = LimitingWeight - LIE_MONOMIAL_WEIGHT(gen);2044i = 0;2045do2046{2047mon = LIE_MONOMIAL_POSITION(i);2048if(LIE_MONOMIAL_IS_BASIS(mon) &&2049LIE_MONOMIAL_IS_NOT_SQUARE(mon) &&2050(i != k || LIE_MONOMIAL_IS_ODD(mon)) &&2051(LIE_MONOMIAL_IS_GENERATOR(gen) ||2052i >= LIE_MONOMIAL_RIGHT_ORDER(gen)))2053{2054if(LIE_MONOMIAL_WEIGHT(mon) > lim_weight_i)2055break; /* Stop considering next right monomials */2056mon = AddPairToLieMonomial(gen, mon);2057if(LieMonomialIsNew)2058{ /* New pairs */2059#if 02060/*??*/PutDebugLieMonomial("\n***New Lie Monomial", mon);2061#endif2062a = NewJacobiRelation(mon);2063if(a != NIL)2064#if 02065/*??*/{2066/*??*/PutDebugLieSum("\n***New Jacobi Relation", a);2067#endif2068goto add_new_relation;2069#if 02070/*??*/}2071#endif2072}2073}2074}while(++i <= k);2075}2076}2077++k;2078}2079loops_out: ;2080/*??#endif ^^^^^^^^^^^^^^^^^^^^ Off checking ^^^^^^^^^^^^^^^^^^^^*/2081}2082/*=NewJacobiRelation===================================================2083Construct independent Jacobi relation containing leading commutator L2084*/2085uint NewJacobiRelation(int l)2086{2087uint a, b;2088int x, y, z;2089IN_NEW_JACOBI_RELATION /*--------------------------------------------*/2090/* Try left pair in l = [[x,y],z] to make relation: */2091/* p(x)p(y) */2092/* [[x,y],z] - [x,[y,z]] + (-1) [y,[x,z]] = 0 */20932094x = LIE_MONOMIAL_LEFT(l);2095if(LIE_MONOMIAL_IS_COMMUTATOR(x))2096{2097z = LIE_MONOMIAL_RIGHT(l);2098y = LIE_MONOMIAL_RIGHT(x);2099x = LIE_MONOMIAL_LEFT(x);21002101/* 1st of triple */21022103a = (*LieTermFromMonomial)(l);21042105/* 2nd of triple */21062107if((b = (*PairMonomialMonomial)(z, y)) != NIL)2108{2109if(LIE_MONOMIAL_IS_ODD(y) && LIE_MONOMIAL_IS_ODD(z))2110b = (*PairSumMonomial)(b, x);2111else2112b = (*PairMonomialSum)(x, b);2113a = LieSumAddition(a, b);2114}21152116/* 3d of triple */21172118if(LIE_MONOMIAL_ORDER(x) > LIE_MONOMIAL_ORDER(z))2119{2120if((b = (*PairMonomialMonomial)(x, z)) != NIL)2121{2122if(LIE_MONOMIAL_IS_ODD(x) && LIE_MONOMIAL_IS_ODD(y))2123{2124b = (*PairSumMonomial)(b, y);2125if(LIE_MONOMIAL_IS_EVEN(z))2126(*LieSumMinus)(b);2127}2128else2129b = (*PairMonomialSum)(y, b);2130a = LieSumAddition(a, b);2131}2132}2133else2134{2135if((b = (*PairMonomialMonomial)(z, x)) != NIL)2136{2137if(LIE_MONOMIAL_IS_ODD(z) &&2138LIE_MONOMIAL_PARITY(x) != LIE_MONOMIAL_PARITY(y))2139{2140b = (*PairMonomialSum)(y, b);2141if(LIE_MONOMIAL_IS_EVEN(x))2142(*LieSumMinus)(b);2143}2144else2145b = (*PairSumMonomial)(b, y);2146a = LieSumAddition(a, b);2147}2148}21492150if(a != NIL)2151goto out;2152}21532154/* Try right pair in l = [x,[y,z]] to make relation: */2155/* p(x)p(y) */2156/* [x,[y,z]] - [[x,y],z] - (-1) [y,[x,z]] = 0 */21572158y = LIE_MONOMIAL_RIGHT(l);2159if(LIE_MONOMIAL_IS_COMMUTATOR(y))2160{2161x = LIE_MONOMIAL_LEFT(l);2162z = LIE_MONOMIAL_RIGHT(y);2163y = LIE_MONOMIAL_LEFT(y);21642165/* 1st of triple */21662167a = (*LieTermFromMonomial)(l);21682169/* 2nd of triple */21702171if((b = (*PairMonomialMonomial)(x, y)) != NIL)2172{2173b = (*PairSumMonomial)(b, z);2174(*LieSumMinus)(b);2175a = LieSumAddition(a, b);2176}21772178/* 3d of triple */21792180if((b = (*PairMonomialMonomial)(x, z)) != NIL)2181{2182b = (*PairMonomialSum)(y, b);2183if(LIE_MONOMIAL_IS_EVEN(x) || LIE_MONOMIAL_IS_EVEN(y))2184(*LieSumMinus)(b);2185a = LieSumAddition(a, b);2186}21872188goto out;2189}2190a = NIL;2191out:2192OUT_NEW_JACOBI_RELATION /*-------------------------------------------*/2193return a;2194}2195/*=ReduceRelations=====================================================2196Reduce the system of relations starting from Ith one. For further2197differentiations returns lowest new positon (or starting one int).2198*/2199int ReduceRelations(int i)2200{2201int i_min, j, lordj, lordl, min_gen, lmoni, k, l, m;2202uint ai, aj;2203i_min = i;2204do /* While relations with new leading monomialsls arise */2205if(RELATION_TO_BE_SUBSTITUTED(i))2206{2207IN_REDUCE_RELATIONS /*----------------------------------------------*/2208new_i:2209ai = RELATION_LIE_SUM(i);2210lmoni = LIE_TERM_MONOMIAL(ai);2211j = i + 1;2212while(j < RelationN)2213{2214aj = RELATION_LIE_SUM(j);2215lordj = LIE_TERM_MONOMIAL_ORDER(aj);2216aj = (*SubstituteRelationInRelation)(ai, aj);2217test_substitution_result:2218if(SubstitutionIsDone)2219{2220if(aj == NIL) /* Killed relation */2221{2222--RelationN;2223k = lordj;2224l = LIE_MONOMIAL_POSITION(k);2225LIE_MONOMIAL_INDEX(l) = 0; /* To ban usage */2226while(++k < LieMonomialN) /* Shift int's of relations */2227{2228l = LIE_MONOMIAL_POSITION(k);2229if(LIE_MONOMIAL_IS_LEADING(l)) /* LieMonomial */2230++LIE_MONOMIAL_INDEX(l);2231}2232for(k = j; k < RelationN; k++) /* Remove gap shifting */2233Relation[k] = Relation[k+1]; /* down top relations */2234continue;2235}2236else /* Non-killing substitution has been done */2237{2238(*NormalizeRelation)(aj);2239lordl = LIE_TERM_MONOMIAL_ORDER(aj);2240if(lordl < lordj) /* Change of leading ordinal */2241{2242l = 0;2243k = j - 1; /* Binary search of new position */2244while(l <= k)2245{2246m = (l + k)/2;2247if(lordl < LIE_TERM_MONOMIAL_ORDER(RELATION_LIE_SUM(m)))2248k = --m;2249else2250l = ++m;2251}2252if(l &&2253LIE_TERM_MONOMIAL_ORDER(RELATION_LIE_SUM(l-1)) == lordl)2254{ /* Substitute once more to avoid collision */2255aj = (*SubstituteRelationInRelation)2256(RELATION_LIE_SUM(l-1), aj);2257goto test_substitution_result;2258}2259/* Set index of dropped relation in LieMonomial */2260LIE_MONOMIAL_INDEX_BY_ORDER(lordl) = ~l;2261/* Set zero to ban using old leading monomial2262of dropped relation */2263LIE_MONOMIAL_INDEX_BY_ORDER(lordj) = 0;2264min_gen = RELATION_MIN_GENERATOR(j);2265/* Shift upper relations and their indices */2266k = (j+1 == RelationN) ? LieMonomialN :2267LIE_TERM_MONOMIAL_ORDER(RELATION_LIE_SUM(j+1));2268while(--k > lordl)2269if(LIE_MONOMIAL_IS_LEADING_BY_ORDER(k))2270--LIE_MONOMIAL_INDEX_BY_ORDER(k); /* Means ++ */2271for(k = j; k > l; k--)2272Relation[k] = Relation[k-1];2273RELATION_MIN_GENERATOR(l) =2274LIE_MONOMIAL_IS_GENERATOR_BY_ORDER(lordl) ?2275GeneratorN : /* Don't differentiate leading generator */2276min_gen; /* Avoiding redifferentiation (or set 0 ??) */2277RELATION_TO_BE_SUBSTITUTED(l) = (byte)(l < RelationN-1);2278RELATION_LIE_SUM(l) = aj;2279if(l <= i)2280{ /* New substituted relation with lesser ordinal */2281i = l;2282if(i < i_min)2283i_min = i;2284goto new_i;2285}2286}2287else /* No change of leading ordinal */2288RELATION_LIE_SUM(j) = aj;2289}2290}2291j++;2292}2293/* Remove consequences of leading monomial of substituted relation2294from LieMonomial table */2295lordj = LieMonomialN-1;2296m = LIE_MONOMIAL_WEIGHT(lmoni);2297while(LIE_MONOMIAL_WEIGHT_BY_ORDER(lordj) > (uint)m)2298{2299l = LIE_MONOMIAL_POSITION(lordj);2300if(IsMonomialInMonomial(lmoni, l))2301{2302if(l < LieMonomialFreePosition) /* Set possibly lowest */2303LieMonomialFreePosition = l; /* free position */2304LIE_MONOMIAL_IS_OCCUPIED(l) = NO; /* Mark free position */2305LieMonomialN--;2306lordl = k = lordj;2307while(k < LieMonomialN)2308{2309l = LIE_MONOMIAL_POSITION(++lordl); /* Shift positions */2310LIE_MONOMIAL_POSITION(k++) = l; /* of upper orders */2311--LIE_MONOMIAL_ORDER(l); /* Decrease orders of upper monomials */2312}2313}2314lordj--;2315}2316RELATION_TO_BE_SUBSTITUTED(i) = NO;2317OUT_REDUCE_RELATIONS /*----------------------------------------------*/2318}2319while(++i < RelationN);2320return i_min;2321}23222323/*_6_1 Pairing functions=======================================*/23242325/*=AddPairToLieMonomial=================================================2326Find position in LieMonomial table for regular pair [i,j].2327LieMonomialIsNew is global boolean signal:2328LieMonomialIsNew == YES if pair has to be added in table,2329if so, add the pair to the table,2330LieMonomialIsNew == NO if pair exists already in table.2331Return position for [i,j]2332*/2333int AddPairToLieMonomial(int i, int j)2334{2335uint wt = LIE_MONOMIAL_WEIGHT(i) + LIE_MONOMIAL_WEIGHT(j);2336int ijo /* left */, r /* right */, m /* middle */, ijp;2337IN_ADD_PAIR_TO_LIE_MONOMIAL /*------------------------------------------*/2338ijo = 0;2339r = LieMonomialN - 1;2340do2341{2342m = (ijo + r)/2;2343ijp = LIE_MONOMIAL_POSITION(m);23442345/* Compare wrt weights */23462347if(wt > LIE_MONOMIAL_WEIGHT(ijp))2348goto shift_left;2349if(wt < LIE_MONOMIAL_WEIGHT(ijp))2350goto shift_right;23512352/* Equal weights: compare lexicographically */23532354if(LIE_MONOMIAL_ORDER(i) > LIE_MONOMIAL_LEFT_ORDER(ijp))2355goto shift_left;2356if(LIE_MONOMIAL_ORDER(i) < LIE_MONOMIAL_LEFT_ORDER(ijp))2357goto shift_right;2358if(LIE_MONOMIAL_ORDER(j) > LIE_MONOMIAL_RIGHT_ORDER(ijp))2359goto shift_left;2360if(LIE_MONOMIAL_ORDER(j) < LIE_MONOMIAL_RIGHT_ORDER(ijp))2361goto shift_right;2362LieMonomialIsNew = NO;2363OUT_ADD_PAIR_TO_LIE_MONOMIAL_OLD /*-------------------------------------*/2364return ijp;2365shift_left:2366ijo = ++m;2367continue;2368shift_right:2369r = --m;2370}while(ijo <= r);23712372/* Add new monomial to table */23732374LieMonomialIsNew = YES;2375if(LieMonomialN >= LieMonomialSize)2376{2377TIME_OFF;2378PutStatistics(); /* No room for new element */2379Error(E_LIE_MONOMIAL_SIZE);2380}2381m = r = LieMonomialN++;2382#if defined(SPACE_STATISTICS)2383if(LieMonomialN > LieMonomialMaxN)2384LieMonomialMaxN = LieMonomialN;2385#endif2386while(m > ijo)2387{2388ijp = LIE_MONOMIAL_POSITION(--r); /* Shift positions of upper orders */2389LIE_MONOMIAL_POSITION(m--) = ijp;2390++LIE_MONOMIAL_ORDER(ijp); /* Increase orders of upper elements */2391}2392while(LIE_MONOMIAL_IS_OCCUPIED(LieMonomialFreePosition))2393LieMonomialFreePosition++; /* Search free position */2394ijp = LieMonomialFreePosition++;2395LIE_MONOMIAL_ORDER(ijp) = ijo;2396LIE_MONOMIAL_POSITION(ijo) = ijp; /* Parity is 1bit field: + is mod 2 */2397LIE_MONOMIAL_PARITY(ijp) = LIE_MONOMIAL_PARITY(i) + LIE_MONOMIAL_PARITY(j);2398LIE_MONOMIAL_LEFT(ijp) = i;2399LIE_MONOMIAL_RIGHT(ijp) = j;2400LIE_MONOMIAL_WEIGHT(ijp) = wt;2401LIE_MONOMIAL_INDEX(ijp) = 0; /* Set type of basis element */2402OUT_ADD_PAIR_TO_LIE_MONOMIAL_NEW /*-------------------------------------*/2403return ijp;2404}2405/*=MakeRelationRHSInteger============================================2406Make r.h.s. of i-th relation (Integer regime)2407Relation is in normalized form (no denominators, no content)2408*/2409uint MakeRelationRHSInteger(int i)2410{2411uint c;2412if((c = LIE_TERM_R(RELATION_LIE_SUM(i))) != NIL)2413{2414BIGINT n, cn;2415#if !defined(RATIONAL_FIELD)2416BIGINT ln = LIE_TERM_NUMERATOR_INTEGER(RELATION_LIE_SUM(i));2417#endif2418uint a, ea;2419IN_MAKE_RELATION_RHS /*-------------------------------------------*/2420/* Copy r.h.s setting negations for numerators */2421a = ea = NodeLTNew();2422LIE_TERM_MONOMIAL(ea) = LIE_TERM_MONOMIAL(c);2423n = LIE_TERM_NUMERATOR_INTEGER(c);2424INTEGER_HEAP_COPY(cn, n, i);2425INTEGER_MINUS(cn);2426LIE_TERM_NUMERATOR_INTEGER(ea) = cn;2427if((n = LIE_TERM_DENOMINATOR_INTEGER(c)) != NULL)2428{2429INTEGER_HEAP_COPY(cn, n, i);2430LIE_TERM_DENOMINATOR_INTEGER(ea) = cn;2431}2432else2433LIE_TERM_DENOMINATOR_INTEGER(ea) = NULL;2434while((c = LIE_TERM_R(c)) != NIL)2435{2436LIE_TERM_R(ea) = NodeLTNew();2437ea = LIE_TERM_R(ea);2438LIE_TERM_MONOMIAL(ea) = LIE_TERM_MONOMIAL(c);2439n = LIE_TERM_NUMERATOR_INTEGER(c);2440INTEGER_HEAP_COPY(cn, n, i);2441INTEGER_MINUS(cn);2442LIE_TERM_NUMERATOR_INTEGER(ea)= cn;2443if((n = LIE_TERM_DENOMINATOR_INTEGER(c)) != NULL)2444{2445INTEGER_HEAP_COPY(cn, n, i);2446LIE_TERM_DENOMINATOR_INTEGER(ea) = cn;2447}2448else2449LIE_TERM_DENOMINATOR_INTEGER(ea) = NULL;2450}2451#if !defined(RATIONAL_FIELD)2452/* Divide by leading coefficient */2453if(INTEGER_IS_NOT_UNIT(ln))2454{2455INTEGER_STACK_COPY(cn, ln, i);2456LieSumDivInteger(a, cn);2457}2458#endif2459OUT_MAKE_RELATION_RHS /*-------------------------------------------*/2460return a;2461}2462return NIL;2463}2464/*=MakeRelationRHSParametric============================================2465Make r.h.s. of i-th relation (Parametric regime)2466Relation is in normalized form (no denominators, no content)2467*/2468uint MakeRelationRHSParametric(int i)2469{2470uint a;2471IN_MAKE_RELATION_RHS /*----------------------------------------------*/2472if((a = LieSumCopyParametric(LIE_TERM_R(RELATION_LIE_SUM(i)))) != NIL)2473{2474uint lc;24752476/* Negation */24772478LieSumMinusParametric(a);24792480/* Divide by leading coefficient */24812482lc = LIE_TERM_NUMERATOR_SCALAR_SUM(RELATION_LIE_SUM(i));2483if(SCALAR_SUM_IS_NOT_UNIT(lc))2484{2485lc = ScalarSumCopy(lc);2486LieSumDivScalarSum(a, lc);2487}2488}2489OUT_MAKE_RELATION_RHS /*----------------------------------------------*/2490return a;2491}2492/*=PairMonomialMonomialInteger==============================2493Make regular expression from two monomials (Integer regime)2494Caller ensures ORDER(i) >= ORDER(j)2495*/2496uint PairMonomialMonomialInteger(int i, int j)2497{2498uint a;2499int k;2500IN_PAIR_MONOMIAL_MONOMIAL /*-----------------------------*/2501if(LIE_MONOMIAL_IS_SQUARE(j))2502{ /* [i,[j,j]] = 2 [[i,j],j] */2503j = LIE_MONOMIAL_LEFT(j);2504a = PairMonomialMonomialInteger(i, j);2505a = PairSumMonomialInteger(a, j);2506if(a != NIL)2507{2508LIMB two[2] = {1, 2};2509LieSumMultInteger(a, two);2510}2511}2512else if(LIE_MONOMIAL_IS_SQUARE(i))2513{2514LIMB two[2] = {1, 2};2515i = LIE_MONOMIAL_LEFT(i);2516if(i == j)2517a = NIL; /* [[i,i],i] = 0 */2518else2519{2520if(LIE_MONOMIAL_ORDER(i) < LIE_MONOMIAL_ORDER(j))2521a = PairMonomialMonomialInteger(j, i);2522/* [[i,i],j] = - 2 [[j,i],i] */2523else2524{2525a = PairMonomialMonomialInteger(i, j);2526if(LIE_MONOMIAL_IS_EVEN(j))2527goto last_pairing; /* [[i,i],j] = 2 [[i,j],i] */2528} /* IS_ODD(j) -> [[i,i],j] = - 2 [[i,j],i] */2529INTEGER_SET_MINUS(two); /* 2 -> -2 */2530last_pairing:2531a = PairSumMonomialInteger(a, i);2532if(a != NIL)2533LieSumMultInteger(a, two);2534}2535}2536else if(LIE_MONOMIAL_IS_COMMUTATOR(i) && /* Irregular triple */2537LIE_MONOMIAL_ORDER(j) < LIE_MONOMIAL_RIGHT_ORDER(i))2538{2539uint b; /* i > k > j => [[i,k],j]] = */2540k = LIE_MONOMIAL_RIGHT(i);2541i = LIE_MONOMIAL_LEFT(i);2542a = PairMonomialMonomialInteger(i, j);2543a = PairSumMonomialInteger(a, k); /* + [[i,j],k] */2544if(LIE_MONOMIAL_IS_ODD(j) && LIE_MONOMIAL_IS_ODD(k))2545{2546uint c = a; /* - [[i,j],k] */2547while(c != NIL)2548{2549LIE_TERM_MINUS_INTEGER(c);2550c = LIE_TERM_R(c);2551}2552}2553b = PairMonomialMonomialInteger(k, j);2554b = PairSumMonomialInteger(b, i); /* + [[k,j],i] */2555if(LIE_MONOMIAL_IS_EVEN(i) ||2556LIE_MONOMIAL_PARITY(j) == LIE_MONOMIAL_PARITY(k))2557{2558uint c = b; /* - [[k,j],i] */2559while(c != NIL)2560{2561LIE_TERM_MINUS_INTEGER(c);2562c = LIE_TERM_R(c);2563}2564}2565a = LieSumAddition(a, b);2566}2567else if(i == j && LIE_MONOMIAL_IS_EVEN(i))2568a = NIL; /* [i,i] = 0 for even i */2569else2570{ /* Regular pair */2571k = AddPairToLieMonomial(i, j);2572if(LIE_MONOMIAL_IS_LEADING(k))2573a = MakeRelationRHSInteger(LIE_MONOMIAL_I_RELATION(k));2574else2575{2576BIGINT n;2577a = NodeLTNew();2578LIE_TERM_MONOMIAL(a) = k;2579INTEGER_HEAP_NEW(n, 2);2580n[0] = n[1] = 1;2581LIE_TERM_NUMERATOR_INTEGER(a) = n;2582LIE_TERM_DENOMINATOR_INTEGER(a) = NULL;2583}2584}2585OUT_PAIR_MONOMIAL_MONOMIAL /*----------------------------*/2586return a;2587}2588/*=PairMonomialMonomialParametric==============================2589Make regular expression from two monomials (Parametric regime)2590Caller ensures i >= j by Shirshov2591*/2592uint PairMonomialMonomialParametric(int i, int j)2593{2594uint a;2595int k;2596IN_PAIR_MONOMIAL_MONOMIAL /*--------------------------------*/2597if(LIE_MONOMIAL_IS_SQUARE(j))2598{ /* [i,[j,j]] = 2 [[i,j],j] */2599j = LIE_MONOMIAL_LEFT(j);2600a = PairMonomialMonomialParametric(i, j);2601a = PairSumMonomialParametric(a, j);2602if(a != NIL)2603{2604uint two;2605BIGINT n2;2606INTEGER_HEAP_NEW(n2, 2);2607n2[0] = 1;2608n2[1] = 2;2609two = NodeSTNew();2610SCALAR_TERM_MONOMIAL(two) = NIL;2611SCALAR_TERM_NUMERATOR(two) = n2;2612LieSumMultScalarSum(a, two);2613}2614}2615else if(LIE_MONOMIAL_IS_SQUARE(i))2616{2617BIGINT n2;2618i = LIE_MONOMIAL_LEFT(i);2619if(i == j)2620a = NIL; /* [[i,i],i] = 0 */2621else2622{2623INTEGER_HEAP_NEW(n2, 2);2624n2[0] = 1;2625n2[1] = 2;2626if(LIE_MONOMIAL_ORDER(i) < LIE_MONOMIAL_ORDER(j))2627a = PairMonomialMonomialParametric(j, i);2628/* [[i,i],j] = - 2 [[j,i],i] */2629else2630{2631a = PairMonomialMonomialParametric(i, j);2632if(LIE_MONOMIAL_IS_EVEN(j))2633goto last_pairing; /* [[i,i],j] = 2 [[i,j],i] */2634} /* IS_ODD(j) -> [[i,i],j] = - 2 [[i,j],i] */2635INTEGER_SET_MINUS(n2); /* 2 -> -2 */2636last_pairing:2637a = PairSumMonomialParametric(a, i);2638if(a != NIL)2639{2640uint two = NodeSTNew();2641SCALAR_TERM_MONOMIAL(two) = NIL;2642SCALAR_TERM_NUMERATOR(two) = n2;2643LieSumMultScalarSum(a, two);2644}2645else2646INTEGER_KILL(n2);2647}2648}2649else if(LIE_MONOMIAL_IS_COMMUTATOR(i) && /* Irregular triple */2650LIE_MONOMIAL_ORDER(j) < LIE_MONOMIAL_RIGHT_ORDER(i))2651{2652uint b; /* i > k > j => [[i,k],j]] = */2653k = LIE_MONOMIAL_RIGHT(i);2654i = LIE_MONOMIAL_LEFT(i);2655a = PairMonomialMonomialParametric(i, j);2656a = PairSumMonomialParametric(a, k); /* + [[i,j],k] */2657if(LIE_MONOMIAL_IS_ODD(j) && LIE_MONOMIAL_IS_ODD(k))2658LieSumMinusParametric(a); /* - [[i,j],k] */2659b = PairMonomialMonomialParametric(k, j);2660b = PairSumMonomialParametric(b, i); /* + [[k,j],i] */2661if(LIE_MONOMIAL_IS_EVEN(i) ||2662LIE_MONOMIAL_PARITY(j) == LIE_MONOMIAL_PARITY(k))2663LieSumMinusParametric(b); /* - [[k,j],i] */2664a = LieSumAddition(a, b);2665}2666else if(i == j && LIE_MONOMIAL_IS_EVEN(i))2667a = NIL; /* [i,i] = 0 for even i */2668else2669{ /* Regular pair */2670k = AddPairToLieMonomial(i, j);2671if(LIE_MONOMIAL_IS_LEADING(k))2672a = MakeRelationRHSParametric(LIE_MONOMIAL_I_RELATION(k));2673else2674{2675BIGINT n;2676uint c;2677a = NodeLTNew();2678LIE_TERM_MONOMIAL(a) = k;2679INTEGER_HEAP_NEW(n, 2);2680n[0] = n[1] = 1;2681c = NodeSTNew();2682SCALAR_TERM_MONOMIAL(c) = NIL;2683SCALAR_TERM_NUMERATOR(c) = n;2684LIE_TERM_NUMERATOR_SCALAR_SUM(a) = c;2685LIE_TERM_DENOMINATOR_SCALAR_SUM(a) = NIL;2686}2687}2688OUT_PAIR_MONOMIAL_MONOMIAL /*--------------------------------*/2689return a;2690}2691/*=PairMonomialSumInteger==========================================26922693Make commutator of the form [mon, Lie_sum] (Integer regime)2694*/2695uint PairMonomialSumInteger(int mon, uint a)2696{2697uint b, s;2698BIGINT nb, db; /* Sum has definite parity */2699int monb, change_sign = (LIE_MONOMIAL_IS_EVEN(mon) ||2700LIE_MONOMIAL_IS_EVEN(LIE_TERM_MONOMIAL(a)));2701IN_PAIR_MONOMIAL_SUM /*----------------------------------------*/2702s = NIL;2703while(a != NIL)2704{2705b = a;2706a = LIE_TERM_R(a);2707monb = LIE_TERM_MONOMIAL(b); /* Take full info from `b' */2708nb = LIE_TERM_NUMERATOR_INTEGER(b);2709db = LIE_TERM_DENOMINATOR_INTEGER(b);2710NODE_LT_KILL(b);2711if(mon != monb || LIE_MONOMIAL_IS_ODD(mon))2712{ /* [mon, monb] != 0 */2713if(LIE_MONOMIAL_ORDER(mon) < LIE_MONOMIAL_ORDER(monb))2714{ /* Swap monomials */2715if(change_sign)2716INTEGER_MINUS(nb);2717b = PairMonomialMonomialInteger(monb, mon);2718}2719else2720b = PairMonomialMonomialInteger(mon, monb);2721if(INTEGER_IS_NOT_UNIT(nb))2722LieSumMultInteger(b, nb);2723if(db != NULL)2724LieSumDivInteger(b, db);2725s = LieSumAddition(s, b);2726}2727INTEGER_KILL(nb);2728if(db != NULL)2729INTEGER_KILL(db);2730}2731OUT_PAIR_MONOMIAL_SUM /*---------------------------------------*/2732return s;2733}2734/*=PairMonomialSumParametric======================================2735Make commutator of the form [mon, Lie_sum] (Parametric regime)2736*/2737uint PairMonomialSumParametric(int mon, uint a)2738{2739uint b, s, nb, db; /* Sum has definite parity */2740int monb, change_sign = (LIE_MONOMIAL_IS_EVEN(mon) ||2741LIE_MONOMIAL_IS_EVEN(LIE_TERM_MONOMIAL(a)));2742IN_PAIR_MONOMIAL_SUM /*---------------------------------------*/2743s = NIL;2744while(a != NIL)2745{2746b = a;2747a = LIE_TERM_R(a);2748monb = LIE_TERM_MONOMIAL(b); /* Take full info from `b' */2749nb = LIE_TERM_NUMERATOR_SCALAR_SUM(b);2750db = LIE_TERM_DENOMINATOR_SCALAR_SUM(b);2751NODE_LT_KILL(b);2752if(mon != monb || LIE_MONOMIAL_IS_ODD(mon))2753{ /* [mon, monb] != 0 */2754if(LIE_MONOMIAL_ORDER(mon) < LIE_MONOMIAL_ORDER(monb))2755{ /* Swap monomials */2756if(change_sign)2757ScalarSumMinus(nb);2758b = PairMonomialMonomialParametric(monb, mon);2759}2760else2761b = PairMonomialMonomialParametric(mon, monb);2762if(SCALAR_SUM_IS_NOT_UNIT(nb))2763LieSumMultScalarSum(b, nb);2764else2765ScalarSumKill(nb);2766if(db != NIL)2767LieSumDivScalarSum(b, db);2768s = LieSumAddition(s, b);2769}2770else2771{2772ScalarSumKill(nb);2773ScalarSumKill(db);2774}2775}2776OUT_PAIR_MONOMIAL_SUM /*---------------------------------------*/2777return s;2778}2779/*=PairSumMonomialInteger==========================================27802781Make commutator of the form [Lie_sum, mon] (Integer regime)2782*/2783uint PairSumMonomialInteger(uint a, int mon)2784{2785uint b, s;2786BIGINT nb, db; /* Sum has definite parity */2787int monb, change_sign = (LIE_MONOMIAL_IS_EVEN(mon) ||2788LIE_MONOMIAL_IS_EVEN(LIE_TERM_MONOMIAL(a)));2789IN_PAIR_SUM_MONOMIAL /*---------------------------------------*/2790s = NIL;2791while(a != NIL)2792{2793b = a;2794a = LIE_TERM_R(a);2795monb = LIE_TERM_MONOMIAL(b); /* Take full info from `b' */2796nb = LIE_TERM_NUMERATOR_INTEGER(b);2797db = LIE_TERM_DENOMINATOR_INTEGER(b);2798NODE_LT_KILL(b);2799if(mon != monb || LIE_MONOMIAL_IS_ODD(mon))2800{ /* [monb, mon] != 0 */2801if(LIE_MONOMIAL_ORDER(monb) < LIE_MONOMIAL_ORDER(mon))2802{ /* Swap monomials */2803if(change_sign)2804INTEGER_MINUS(nb);2805b = PairMonomialMonomialInteger(mon, monb);2806}2807else2808b = PairMonomialMonomialInteger(monb, mon);2809if(INTEGER_IS_NOT_UNIT(nb))2810LieSumMultInteger(b, nb);2811if(db != NULL)2812LieSumDivInteger(b, db);2813s = LieSumAddition(s, b);2814}2815INTEGER_KILL(nb);2816if(db != NULL)2817INTEGER_KILL(db);2818}2819OUT_PAIR_SUM_MONOMIAL /*---------------------------------------*/2820return s;2821}2822/*=PairSumMonomialParametric======================================2823Make commutator of the form [Lie_sum, mon] (Parametric regime)2824*/2825uint PairSumMonomialParametric(uint a, int mon)2826{2827uint b, s, nb, db; /* Sum has definite parity */2828int monb, change_sign = (LIE_MONOMIAL_IS_EVEN(mon) ||2829LIE_MONOMIAL_IS_EVEN(LIE_TERM_MONOMIAL(a)));2830IN_PAIR_SUM_MONOMIAL /*---------------------------------------*/2831s = NIL;2832while(a != NIL)2833{2834b = a;2835a = LIE_TERM_R(a);2836monb = LIE_TERM_MONOMIAL(b); /* Take full info from `b' */2837nb = LIE_TERM_NUMERATOR_SCALAR_SUM(b);2838db = LIE_TERM_DENOMINATOR_SCALAR_SUM(b);2839NODE_LT_KILL(b);2840if(mon != monb || LIE_MONOMIAL_IS_ODD(mon))2841{ /* [monb, mon] != 0 */2842if(LIE_MONOMIAL_ORDER(monb) < LIE_MONOMIAL_ORDER(mon))2843{ /* Swap monomials */2844if(change_sign)2845ScalarSumMinus(nb);2846b = PairMonomialMonomialParametric(mon, monb);2847}2848else2849b = PairMonomialMonomialParametric(monb, mon);2850if(SCALAR_SUM_IS_NOT_UNIT(nb))2851LieSumMultScalarSum(b, nb);2852else2853ScalarSumKill(nb);2854if(db != NIL)2855LieSumDivScalarSum(b, db);2856s = LieSumAddition(s, b);2857}2858}2859OUT_PAIR_SUM_MONOMIAL /*---------------------------------------*/2860return s;2861}2862/*=PairSumSumInteger======================================2863Commutator of the form [Lie_sum,Lie_sum] (Integer regime)2864*/2865uint PairSumSumInteger(uint a, uint b)2866{2867IN_PAIR_SUM_SUM /*-------------------------------------*/2868if(a == NIL)2869LieSumKillInteger(b);2870else if(b == NIL)2871{2872LieSumKillInteger(a);2873a = b;2874}2875else2876{2877BIGINT num, den;2878uint c, d, s;2879s = NIL;2880do2881{2882c = a;2883a = LIE_TERM_R(c);2884d = (a != NIL) ? LieSumCopyInteger(b) : b;2885d = PairMonomialSumInteger(LIE_TERM_MONOMIAL(c), d);2886num = LIE_TERM_NUMERATOR_INTEGER(c);2887den = LIE_TERM_DENOMINATOR_INTEGER(c);2888NODE_LT_KILL(c);2889if(INTEGER_IS_NOT_UNIT(num))2890LieSumMultInteger(d, num);2891INTEGER_KILL(num);2892if(den != NULL)2893{2894LieSumDivInteger(d, den);2895INTEGER_KILL(den);2896}2897s = LieSumAddition(s, d);2898}while(a != NIL);2899a = s;2900}2901OUT_PAIR_SUM_SUM /*-------------------------------------*/2902return a;2903}2904/*=PairSumSumParametric======================================2905Commutator of the form [Lie_sum,Lie_sum] (Parametric regime)2906*/2907uint PairSumSumParametric(uint a, uint b)2908{2909IN_PAIR_SUM_SUM /*----------------------------------------*/2910if(a == NIL)2911LieSumKillParametric(b);2912else if(b == NIL)2913{2914LieSumKillParametric(a);2915a = b;2916}2917else2918{2919uint num, den, c, d, s;2920s = NIL;2921do2922{2923c = a;2924a = LIE_TERM_R(c);2925d = (a != NIL) ? LieSumCopyParametric(b) : b;2926d = PairMonomialSumParametric(LIE_TERM_MONOMIAL(c), d);2927num = LIE_TERM_NUMERATOR_SCALAR_SUM(c);2928den = LIE_TERM_DENOMINATOR_SCALAR_SUM(c);2929NODE_LT_KILL(c);2930if(SCALAR_SUM_IS_NOT_UNIT(num))2931LieSumMultScalarSum(d, num);2932else2933ScalarSumKill(num);2934if(den != NIL)2935LieSumDivScalarSum(d, den);2936s = LieSumAddition(s, d);2937}while(a != NIL);2938a = s;2939}2940OUT_PAIR_SUM_SUM /*----------------------------------------*/2941return a;2942}29432944/*_6_2 Substitution (replacing) functions======================*/29452946/*=IsMonomialInMonomial==========================================29472948Check whether `mon' contains `submon'2949*/2950int IsMonomialInMonomial(int submon, int mon)2951{2952if(LIE_MONOMIAL_ORDER(submon) > LIE_MONOMIAL_ORDER(mon))2953return NO;2954if(submon == mon)2955return YES;2956if(LIE_MONOMIAL_IS_GENERATOR(mon))2957return NO;2958return IsMonomialInMonomial(submon, LIE_MONOMIAL_LEFT(mon)) ||2959IsMonomialInMonomial(submon, LIE_MONOMIAL_RIGHT(mon));2960}2961/*=SubstituteRelationInRelationInteger================================2962R is donor and unchanged, A is acceptor and changed. Integer regime2963*/2964uint SubstituteRelationInRelationInteger(uint r, uint a)2965{2966uint b, bl, bf, rhs;2967BIGINT nb;2968#if defined(RATIONAL_FIELD)2969BIGINT db;2970#endif2971int lmon, monb, lord;2972IN_SUBSTITUTE_RELATION_IN_RELATION /*-------------------------------*/2973lmon = LIE_TERM_MONOMIAL(r);2974lord = LIE_MONOMIAL_ORDER(lmon);2975SubstitutionIsDone = NO;2976bl = NIL;2977bf = b = a;2978do2979{2980monb = LIE_TERM_MONOMIAL(b);2981if(LIE_MONOMIAL_ORDER(monb) < lord)2982{2983if(SubstitutionIsDone)2984{2985a = LieSumAddition(a, bf);2986LieSumKillInteger(rhs);2987}2988goto out;2989}2990if(IsMonomialInMonomial(lmon, monb))2991{2992if(SubstitutionIsDone == NO) /* First substituent term */2993{ /* Make right hand side of R */2994rhs = LieSumCopyIntegerNegative(LIE_TERM_R(r));2995#if !defined(RATIONAL_FIELD)2996if(rhs != NIL)2997{2998nb = LIE_TERM_NUMERATOR_INTEGER(r);2999if(INTEGER_IS_NOT_UNIT(nb))3000{3001BIGINT den;3002int i; /* Divide by leading coefficient */3003INTEGER_STACK_COPY(den, nb, i);3004LieSumDivInteger(rhs, den);3005}3006}3007#endif3008SubstitutionIsDone = YES;3009a = NIL;3010}3011nb = LIE_TERM_NUMERATOR_INTEGER(b);3012if((r = SubstituteRHSInMonomialInteger(monb, lmon, rhs)) != NIL)3013{3014#if defined(RATIONAL_FIELD) /* Field R case compiling */3015db = LIE_TERM_DENOMINATOR_INTEGER(b);3016if(INTEGER_IS_UNIT(nb))3017{3018if(db != NULL)3019{3020LieSumDivInteger(r, db);3021INTEGER_KILL(db);3022}3023}3024else3025if(db != NULL)3026{3027LieSumMultRationalInteger(r, nb, db);3028INTEGER_KILL(db);3029}3030else3031LieSumMultInteger(r, nb);3032#else /* Ring Z case compiling */3033if(INTEGER_IS_NOT_UNIT(nb))3034LieSumMultInteger(r, nb);3035#endif3036a = LieSumAddition(a, r);3037}3038INTEGER_KILL(nb);3039if(bl != NIL)3040{ /* There is substitution-free sublist before */3041LIE_TERM_R(bl) = NIL;3042a = LieSumAddition(a, bf);3043}3044bl = b;3045bf = b = LIE_TERM_R(b);3046NODE_LT_KILL(bl);3047bl = NIL;3048}3049else3050{ /* Skip substitution-free term */3051bl = b;3052b = LIE_TERM_R(b);3053}3054}while(b != NIL);3055if(SubstitutionIsDone)3056{ /* Append substitution-free tail */3057if(bl != NIL)3058a = LieSumAddition(a, bf);3059LieSumKillInteger(rhs);3060}3061out:3062OUT_SUBSTITUTE_RELATION_IN_RELATION /*------------------------------*/3063return a;3064}3065/*=SubstituteRelationInRelationParametric================================3066R is donor and unchanged, A is acceptor and changed. Parametric regime3067*/3068uint SubstituteRelationInRelationParametric(uint r, uint a)3069{3070uint b, bl, bf, rhs, pb;3071int lmon, monb, lord;3072IN_SUBSTITUTE_RELATION_IN_RELATION /*----------------------------------*/3073lmon = LIE_TERM_MONOMIAL(r);3074lord = LIE_MONOMIAL_ORDER(lmon);3075SubstitutionIsDone = NO;3076bl = NIL;3077bf = b = a;3078do3079{3080monb = LIE_TERM_MONOMIAL(b);3081if(LIE_MONOMIAL_ORDER(monb) < lord)3082{3083if(SubstitutionIsDone)3084{3085a = LieSumAddition(a, bf);3086LieSumKillParametric(rhs);3087}3088goto out;3089}3090if(IsMonomialInMonomial(lmon, monb))3091{3092if(SubstitutionIsDone == NO) /* First substituent term */3093{ /* Make right hand side of R */3094if((rhs = LieSumCopyParametric(LIE_TERM_R(r))) != NIL)3095{3096pb = LIE_TERM_NUMERATOR_SCALAR_SUM(r);3097if(SCALAR_SUM_IS_NOT_UNIT(pb)) /* Divide by leading coeff. */3098LieSumDivScalarSum(rhs, ScalarSumCopy(pb));3099LieSumMinusParametric(rhs); /* Negation */3100}3101SubstitutionIsDone = YES;3102a = NIL;3103}3104pb = LIE_TERM_NUMERATOR_SCALAR_SUM(b);3105if((r = SubstituteRHSInMonomialParametric(monb, lmon, rhs)) != NIL)3106{3107if(SCALAR_SUM_IS_NOT_UNIT(pb))3108LieSumMultScalarSum(r, pb);3109else3110ScalarSumKill(pb);3111a = LieSumAddition(a, r);3112}3113if(bl != NIL)3114{ /* There is substitution-free sublist before */3115LIE_TERM_R(bl) = NIL;3116a = LieSumAddition(a, bf);3117}3118bl = b;3119bf = b = LIE_TERM_R(b);3120NODE_LT_KILL(bl);3121bl = NIL;3122}3123else3124{ /* Skip substitution-free term */3125bl = b;3126b = LIE_TERM_R(b);3127}3128}while(b != NIL);3129if(SubstitutionIsDone)3130{ /* Append substitution-free tail */3131if(bl != NIL)3132a = LieSumAddition(a, bf);3133LieSumKillParametric(rhs);3134}3135out:3136OUT_SUBSTITUTE_RELATION_IN_RELATION /*---------------------------------*/3137return a;3138}3139/*=SubstituteRHSInMonomialInteger=======================================3140Insert right hand side `r' of relation with leading monomial `lmonr'3141in monomial `mon'.3142Returned NOTHING means "no substitution", NIL means 0.3143Function saves input `r'. Integer regime.3144*/3145uint SubstituteRHSInMonomialInteger(int mon, int lmonr, uint r)3146{3147/* Single monomial matching case */31483149if(mon == lmonr)3150return LieSumCopyInteger(r);31513152/* Possible substitution(s) in submonomial(s) */31533154if(LIE_MONOMIAL_ORDER(mon) > LIE_MONOMIAL_ORDER(lmonr)3155&& LIE_MONOMIAL_IS_COMMUTATOR(mon))3156{3157uint res;3158int monl = LIE_MONOMIAL_LEFT(mon);3159mon = LIE_MONOMIAL_RIGHT(mon);3160res = SubstituteRHSInMonomialInteger(monl, lmonr, r);3161if(res == NIL)3162return NIL; /* [0, x] -> 0 */3163{3164uint resr = SubstituteRHSInMonomialInteger(mon, lmonr, r);3165if(res == NOTHING)3166{3167if(resr == NOTHING)3168return NOTHING; /* No substitutions in both submonomials */3169if(resr == NIL)3170return NIL; /* [x, 0] -> 0 */3171return PairMonomialSumInteger(monl, resr); /* -> [monl, resr] */3172}3173if(resr == NOTHING)3174return PairSumMonomialInteger(res, mon); /* -> [res, mon] */3175if(resr == NIL)3176{3177LieSumKillInteger(res);3178return NIL;3179}3180return PairSumSumInteger(res, resr); /* -> [res, resr] */3181}3182}3183return NOTHING;3184}3185/*=SubstituteRHSInMonomialParametric======================================3186Insert right hand side `r' of relation with leading monomial `lmonr'3187in monomial `mon'.3188Returned NOTHING means "no substitution", NIL means 0.3189Function saves input `r'. Parametric regime.3190*/3191uint SubstituteRHSInMonomialParametric(int mon, int lmonr, uint r)3192{3193/* Single monomial matching case */31943195if(mon == lmonr)3196return LieSumCopyParametric(r);31973198/* Possible substitution(s) in submonomial(s) */31993200if(LIE_MONOMIAL_ORDER(mon) > LIE_MONOMIAL_ORDER(lmonr)3201&& LIE_MONOMIAL_IS_COMMUTATOR(mon))3202{3203uint res;3204int monl = LIE_MONOMIAL_LEFT(mon);3205mon = LIE_MONOMIAL_RIGHT(mon);3206res = SubstituteRHSInMonomialParametric(monl, lmonr, r);3207if(res == NIL)3208return NIL; /* [0, x] -> 0 */3209{3210uint resr = SubstituteRHSInMonomialParametric(mon, lmonr, r);3211if(res == NOTHING)3212{3213if(resr == NOTHING)3214return NOTHING; /* No substitutions in both submonomials */3215if(resr == NIL)3216return NIL; /* [x, 0] -> 0 */3217return PairMonomialSumParametric(monl, resr); /* -> [monl, resr] */3218}3219if(resr == NOTHING)3220return PairSumMonomialParametric(res, mon); /* -> [res, mon] */3221if(resr == NIL)3222{3223LieSumKillParametric(res);3224return NIL;3225}3226return PairSumSumParametric(res, resr); /* -> [res, resr] */3227}3228}3229return NOTHING;3230}32313232/*_6_3 Lie and scalar algebra functions========================*/32333234/*=LieLikeTermsCollectionInteger==========================================3235For Lie terms `a' and `b' sum rational integers `na/da' and `nb/db'3236destructing input terms, set nonzero result in `a', return YES for3237nonzero result, otherwise kill `a' and return NO;3238`da', `db'= NIL means 1.3239*/3240int LieLikeTermsCollectionInteger(uint a, uint b)3241{3242BIGINT na, da, nb, db;3243nb = LIE_TERM_NUMERATOR_INTEGER(b);3244db = LIE_TERM_DENOMINATOR_INTEGER(b);3245NODE_LT_KILL(b);3246na = LIE_TERM_NUMERATOR_INTEGER(a);3247da = LIE_TERM_DENOMINATOR_INTEGER(a);3248if(da != NULL) if(db != NULL) /* `da' != 1, `db' != 1 */3249{3250BIGINT g, h;3251int i;3252INTEGER_STACK_COPY(g, da, i);3253INTEGER_STACK_COPY(h, db, i);3254if((g = IntegerGCD(g, h)) != NULL) /* g = GCD(da, db) > 1 */3255{3256BIGINT k, m, daa;3257INTEGER_STACK_COPY(k, g, i); /* k = GCD(da, db)' */3258INTEGER_STACK_NEW(m, 2+INTEGER_N_LIMBS(da)-INTEGER_N_LIMBS(k));3259INTEGER_STACK_COPY_1(daa, da, i);3260INTEGER_KILL(da);3261IntegerQuotient(m, daa, k); /* m = da/GCD(da, db)' */3262INTEGER_STACK_NEW(h, 1+INTEGER_N_LIMBS(nb)+INTEGER_N_LIMBS(m));3263IntegerProduct(h, nb, m); /* h = nb*da'/GCD(da, db)' */3264INTEGER_KILL(nb);3265INTEGER_STACK_COPY_1(k, db, i); /* k = db' */3266INTEGER_STACK_COPY(da, g, i); /* da = GCD(da, db)' */3267INTEGER_STACK_NEW(nb, 2+INTEGER_N_LIMBS(k)-INTEGER_N_LIMBS(da));3268IntegerQuotient(nb, k, da); /* nb = db'/GCD(da, db)' */3269INTEGER_STACK_NEW(k, 1+INTEGER_N_LIMBS(na)+INTEGER_N_LIMBS(nb));3270IntegerProduct(k, na, nb); /* k = na*db'/GCD(da, db)' */3271INTEGER_KILL(na);3272INTEGER_STACK_NEW(na, 3+MAX(INTEGER_N_LIMBS(h),INTEGER_N_LIMBS(k)));3273IntegerSum(na, h, k); /* na = h + k */3274if(INTEGER_N_LIMBS(na) != 0)3275{3276INTEGER_STACK_COPY(h, na, i); /* Numerator */3277if((h = IntegerGCD(h, g)) != NULL)3278{3279INTEGER_STACK_COPY(g, h, i);3280INTEGER_HEAP_NEW(k, 2+INTEGER_N_LIMBS(na)-INTEGER_N_LIMBS(g));3281IntegerQuotient(k, na, g);3282LIE_TERM_NUMERATOR_INTEGER(a) = k;3283INTEGER_STACK_NEW(k, 2+INTEGER_N_LIMBS(db)-INTEGER_N_LIMBS(h));3284INTEGER_STACK_COPY_1(daa, db, i);3285INTEGER_KILL(db);3286IntegerQuotient(k, daa, h);3287INTEGER_HEAP_NEW(db, 1+INTEGER_N_LIMBS(m)+INTEGER_N_LIMBS(k));3288IntegerProduct(db, m, k);3289if(INTEGER_IS_UNIT(db))3290{3291INTEGER_KILL(db);3292db = NULL; /* Standard convention for unit denominators */3293}3294LIE_TERM_DENOMINATOR_INTEGER(a) = db;3295}3296else3297{3298INTEGER_HEAP_COPY(h, na, i);3299LIE_TERM_NUMERATOR_INTEGER(a) = h;3300INTEGER_HEAP_NEW(k, 1+INTEGER_N_LIMBS(m)+INTEGER_N_LIMBS(db));3301IntegerProduct(k, m, db);3302INTEGER_KILL(db);3303if(INTEGER_IS_UNIT(k))3304{3305INTEGER_KILL(k);3306k = NULL; /* Standard convention for unit denominators */3307}3308LIE_TERM_DENOMINATOR_INTEGER(a) = k;3309}3310goto non_zero;3311}3312else3313{3314INTEGER_KILL(db);3315goto zero;3316}3317}3318else /* Mutually prime da, db */3319{3320INTEGER_STACK_NEW(g, 1+INTEGER_N_LIMBS(nb)+INTEGER_N_LIMBS(da));3321IntegerProduct(g, nb, da); /* g = nb*da */3322INTEGER_KILL(nb);3323INTEGER_STACK_NEW(nb, 1+INTEGER_N_LIMBS(na)+INTEGER_N_LIMBS(db));3324IntegerProduct(nb, na, db); /* nb = na*db */3325INTEGER_KILL(na);3326INTEGER_HEAP_NEW(na, 2+MAX(INTEGER_N_LIMBS(nb),INTEGER_N_LIMBS(g)));3327IntegerSum(na, nb, g); /* na = nb*da + na*db */3328if(INTEGER_N_LIMBS(na) != 0)3329{3330LIE_TERM_NUMERATOR_INTEGER(a) = na;3331INTEGER_HEAP_NEW(nb, 1+INTEGER_N_LIMBS(da)+INTEGER_N_LIMBS(db));3332IntegerProduct(nb, da, db); /* nb = da*db */3333LIE_TERM_DENOMINATOR_INTEGER(a) = nb;3334INTEGER_KILL(da);3335INTEGER_KILL(db);3336goto non_zero;3337}3338else3339{3340INTEGER_KILL(na);3341INTEGER_KILL(da);3342INTEGER_KILL(db);3343goto zero;3344}3345}3346}3347else /* `da' != 1, `db' = 1 */3348{3349INTEGER_STACK_NEW(db, 1+INTEGER_N_LIMBS(nb)+INTEGER_N_LIMBS(da));3350IntegerProduct(db, nb, da);3351INTEGER_KILL(nb);3352INTEGER_HEAP_NEW(nb, 2+MAX(INTEGER_N_LIMBS(db),INTEGER_N_LIMBS(na)));3353IntegerSum(nb, db, na);3354INTEGER_KILL(na);3355LIE_TERM_NUMERATOR_INTEGER(a) = nb;3356goto non_zero;3357}3358else if(db != NULL) /* `da' = 1, `db' != 1 */3359{3360INTEGER_STACK_NEW(da, 1+INTEGER_N_LIMBS(na)+INTEGER_N_LIMBS(db));3361IntegerProduct(da, na, db);3362INTEGER_KILL(na);3363INTEGER_HEAP_NEW(na, 2+MAX(INTEGER_N_LIMBS(da),INTEGER_N_LIMBS(nb)));3364IntegerSum(na, da, nb);3365INTEGER_KILL(nb);3366LIE_TERM_NUMERATOR_INTEGER(a) = na;3367LIE_TERM_DENOMINATOR_INTEGER(a) = db;3368goto non_zero;3369}3370else /* `da' = `db' = 1 */3371{3372INTEGER_HEAP_NEW(da, 2+MAX(INTEGER_N_LIMBS(na),INTEGER_N_LIMBS(nb)));3373IntegerSum(da, na, nb);3374INTEGER_KILL(na);3375INTEGER_KILL(nb);3376if(INTEGER_N_LIMBS(da) != 0)3377{3378LIE_TERM_NUMERATOR_INTEGER(a) = da;3379goto non_zero;3380}3381else3382{3383INTEGER_KILL(da);3384goto zero;3385}3386}3387non_zero:3388return YES;3389zero:3390NODE_LT_KILL(a); /* `na' + `nb' = 0 */3391return NO;3392}3393/*=LieLikeTermsCollectionParametric=======================================3394For Lie terms `a' and `b' sum rational functions `na/da' and `nb/db'3395destructing input terms, set nonzero result in `a', return YES for3396nonzero result, otherwise kill `a' and return NO;3397`da', `db'= NIL means 1.3398*/3399int LieLikeTermsCollectionParametric(uint a, uint b)3400{3401uint na, da, nb, db;3402nb = LIE_TERM_NUMERATOR_SCALAR_SUM(b);3403db = LIE_TERM_DENOMINATOR_SCALAR_SUM(b);3404NODE_LT_KILL(b);3405na = LIE_TERM_NUMERATOR_SCALAR_SUM(a);3406da = LIE_TERM_DENOMINATOR_SCALAR_SUM(a);3407if(da != NIL) if(db != NIL) /* `da' != 1 and `db' != 1 */3408{3409uint g;3410if((g = PolyGCD(da, db)) != NIL) /* g = GCD(da, db) != 1 */3411{3412da = PolyQuotient(da, g); /* da' = da/GCD(da, db) */3413/* nb' = nb*da' */3414nb = ScalarSumMultiplication(nb, ScalarSumCopy(da));3415/* na' = na*db/GCD(da, db) */3416na = ScalarSumMultiplication(na, PolyQuotient(ScalarSumCopy(db), g));3417na = ScalarSumAddition(na, nb); /* na'' = na' + nb' */3418if(na != NIL)3419{3420if((nb = PolyGCD(na, g)) != NIL)3421{ /* Set na''/GCD(na'', g) */3422LIE_TERM_NUMERATOR_SCALAR_SUM(a) = PolyQuotient(na, nb);3423db = PolyQuotient(db, nb); /* db' = db/GCD(na'', g) */3424ScalarSumKill(nb); /* Kill GCD(na'', g) */3425}3426else3427LIE_TERM_NUMERATOR_SCALAR_SUM(a) = na;3428ScalarSumKill(g); /* Kill GCD(da, db) */3429da = ScalarSumMultiplication(da, db);3430if(SCALAR_SUM_IS_UNIT(da))3431{3432ScalarSumKill(da);3433da = NIL;3434}3435LIE_TERM_DENOMINATOR_SCALAR_SUM(a) = da;3436goto non_zero;3437}3438ScalarSumKill(g); /* Nontrivial g at na == NIL branch */3439}3440else /* Mutually prime da, db */3441{ /* na' = na*db' + nb*da' */3442na = ScalarSumAddition(ScalarSumMultiplication(na, ScalarSumCopy(db)),3443ScalarSumMultiplication(nb, ScalarSumCopy(da)));3444if(na != NIL)3445{3446LIE_TERM_NUMERATOR_SCALAR_SUM(a) = na;3447LIE_TERM_DENOMINATOR_SCALAR_SUM(a) = ScalarSumMultiplication(da, db);3448goto non_zero;3449}3450}3451ScalarSumKill(da); /* na == NIL branch */3452ScalarSumKill(db);3453goto zero;3454}3455else /* `da' != 1 and `db' = 1 --> (na + nb*da')/da */3456{3457LIE_TERM_NUMERATOR_SCALAR_SUM(a) =3458ScalarSumAddition(na, ScalarSumMultiplication(nb, ScalarSumCopy(da)));3459goto non_zero;3460}3461else if(db != NIL) /* `da' = 1 and `db' != 1 --> (nb + na*db')/db */3462{3463LIE_TERM_NUMERATOR_SCALAR_SUM(a) =3464ScalarSumAddition(nb, ScalarSumMultiplication(na, ScalarSumCopy(db)));3465LIE_TERM_DENOMINATOR_SCALAR_SUM(a) = db;3466goto non_zero;3467}3468else if((na = ScalarSumAddition(na, nb)) != NIL)3469{ /* `da' = `db' = 1 --> (na + nb)/1 (!= 0) */3470LIE_TERM_NUMERATOR_SCALAR_SUM(a) = na;3471goto non_zero;3472}3473zero:3474NODE_LT_KILL(a); /* `na' + `nb' = 0 */3475return NO;3476non_zero:3477return YES;3478}3479/*=LieSumAddition=============================================3480Sum of two Lie expressions3481*/3482uint LieSumAddition(uint a, uint b)3483{3484uint sum = NIL, last, wa, wb;3485IN_LIE_SUM_ADDITION /*---------------------------------------*/3486while(YES)3487{3488next_pair:3489if(b == NIL)3490{ /* List b is ended, append rest of a */3491if(sum == NIL)3492sum = a;3493else3494LIE_TERM_R(last) = a;3495break;3496}3497if(a == NIL)3498{ /* List a is ended, append rest of b */3499if(sum == NIL)3500sum = b;3501else3502LIE_TERM_R(last) = b;3503break;3504}35053506/* Compare algebra terms */35073508if(LIE_TERM_MONOMIAL_ORDER(a) > LIE_TERM_MONOMIAL_ORDER(b))3509goto order_12;3510if(LIE_TERM_MONOMIAL_ORDER(a) < LIE_TERM_MONOMIAL_ORDER(b))3511goto order_21;35123513/* Reduce like algebra terms */35143515wa = a;3516wb = b;3517a = LIE_TERM_R(a);3518b = LIE_TERM_R(b);35193520/* Sum rational coefficients */35213522if((*LieLikeTermsCollection)(wa, wb))3523goto append_term;3524else3525goto next_pair;35263527order_12:3528wa = a;3529a = LIE_TERM_R(a);3530goto append_term;35313532order_21:3533wa = b;3534b = LIE_TERM_R(b);35353536append_term:3537if(sum == NIL)3538sum = wa;3539else3540LIE_TERM_R(last) = wa;3541last = wa;3542}3543OUT_LIE_SUM_ADDITION /*--------------------------------------*/3544return sum;3545}3546/*=LieSumDivInteger=======================================================3547Divide Lie sum by integer (of unknown nature) on spot in Integer regime3548Integer `den' is spoiled3549*/3550void LieSumDivInteger(uint lsum, BIGINT den)3551{3552if(lsum != NIL)3553{3554BIGINT d, da, dao;3555int i, n;3556uint a;3557IN_LIE_SUM_DIV_INTEGER /*----------------------------------------------*/3558n = INTEGER_N_LIMBS(den);3559INTEGER_STACK_NEW(d, 1+n); /* Space for copies input `den' */3560do3561{3562a = lsum;3563lsum = LIE_TERM_R(lsum);3564if(lsum != NIL)3565{3566i = n;3567do3568d[i] = den[i];3569while(i--);3570}3571else3572d = den;3573IntegerCancellation(LIE_TERM_NUMERATOR_INTEGER(a), d);3574if(INTEGER_IS_NOT_UNIT(d))3575{3576if((dao = LIE_TERM_DENOMINATOR_INTEGER(a)) != NULL)3577{ /* Nontrivial old denominator `dao' */3578INTEGER_HEAP_NEW(da, 1+INTEGER_N_LIMBS(d)+INTEGER_N_LIMBS(dao));3579IntegerProduct(da, d, dao);3580INTEGER_KILL(dao);3581}3582else3583{3584INTEGER_HEAP_COPY(da, d, i); /* Composite statement */3585}3586LIE_TERM_DENOMINATOR_INTEGER(a) = da;3587}3588}while(lsum != NIL);3589OUT_LIE_SUM_DIV_INTEGER /*---------------------------------------------*/3590}3591}3592/*=LieSumDivScalarSum======================================3593Divide Lie sum by scalar sum on spot in Parametric regime3594`den' is killed3595*/3596void LieSumDivScalarSum(uint lsum, uint den)3597{3598if(lsum == NIL)3599ScalarSumKill(den);3600else3601{3602uint n, d, a;3603IN_LIE_SUM_DIV_SCALAR_SUM /*-----------------------------*/3604do3605{3606a = lsum;3607lsum = LIE_TERM_R(lsum);3608d = (lsum != NIL) ? ScalarSumCopy(den) : den;3609n = LIE_TERM_NUMERATOR_SCALAR_SUM(a);3610ScalarSumCancellation(&n, &d);3611LIE_TERM_NUMERATOR_SCALAR_SUM(a) = n;3612if(SCALAR_SUM_IS_NOT_UNIT(d))3613{3614if((n = LIE_TERM_DENOMINATOR_SCALAR_SUM(a)) != NIL)3615d = ScalarSumMultiplication(d, n); /* Absorb */3616LIE_TERM_DENOMINATOR_SCALAR_SUM(a) = d;3617}3618else3619ScalarSumKill(d); /* Kill unit */3620}while(lsum != NIL);3621OUT_LIE_SUM_DIV_SCALAR_SUM /*----------------------------*/3622}3623}3624/*=LieSumMinusInteger====================3625Change signs in Lie sum (Integer regime)3626*/3627void LieSumMinusInteger(uint a)3628{3629while(a != NIL)3630{3631LIE_TERM_MINUS_INTEGER(a);3632a = LIE_TERM_R(a);3633}3634}36353636/*=LieSumMinusParametric====================3637Change signs in Lie sum (Parametric regime)3638*/3639void LieSumMinusParametric(uint a)3640{3641uint b;3642while(a != NIL)3643{3644b = LIE_TERM_NUMERATOR_SCALAR_SUM(a);3645do3646SCALAR_TERM_MINUS(b);3647while((b = SCALAR_TERM_R(b)) != NIL);3648a = LIE_TERM_R(a);3649}3650}3651/*=LieSumMultInteger=======================================================3652Multiply Lie sum by integer (of unknown nature) on spot in Integer regime3653Integer `num' is spoiled3654*/3655void LieSumMultInteger(uint lsum, BIGINT num)3656{3657if(lsum != NIL)3658{3659BIGINT nw, nao, da;3660int i, n;3661uint a;3662IN_LIE_SUM_MULT_INTEGER /*----------------------------------------------*/3663n = INTEGER_N_LIMBS(num);3664INTEGER_STACK_NEW(nw, 1+n); /* Space for copies of input `num' */3665do3666{3667a = lsum;3668lsum = LIE_TERM_R(lsum);3669nao = LIE_TERM_NUMERATOR_INTEGER(a); /* Old numerator `nao' */3670if((da = LIE_TERM_DENOMINATOR_INTEGER(a)) != NULL)3671{ /* Nontrivial denominator `da' */3672if(lsum != NIL)3673{ /* Copy if not last */3674i = n;3675do3676nw[i] = num[i];3677while(i--);3678}3679else3680nw = num;3681IntegerCancellation(nw, da);3682if(INTEGER_IS_UNIT(da))3683{ /* Trivialize unit denominator */3684INTEGER_KILL(da);3685LIE_TERM_DENOMINATOR_INTEGER(a) = NULL;3686}3687if(INTEGER_IS_NOT_UNIT(nw))3688{3689INTEGER_HEAP_NEW(da, 1+INTEGER_N_LIMBS(nw)+INTEGER_N_LIMBS(nao));3690IntegerProduct(da, nw, nao);3691goto stick_new;3692}3693}3694else3695{3696INTEGER_HEAP_NEW(da, 1+INTEGER_N_LIMBS(num)+INTEGER_N_LIMBS(nao));3697IntegerProduct(da, num, nao);3698stick_new:3699INTEGER_KILL(nao);3700LIE_TERM_NUMERATOR_INTEGER(a) = da;3701}3702}while(lsum != NIL);3703OUT_LIE_SUM_MULT_INTEGER /*----------------------------------------------*/3704}3705}3706#if defined(RATIONAL_FIELD)3707/*=LieSumMultRationalInteger=============================================3708num and den are non-NULL integers of unknown nature3709*/3710void LieSumMultRationalInteger(int a, BIGINT num, BIGINT den)3711{3712BIGINT numc, denc, numa, dena, w;3713int i, nn = INTEGER_N_LIMBS(num), nd = INTEGER_N_LIMBS(den);3714INTEGER_STACK_NEW(numc, 1+nn);3715INTEGER_STACK_NEW(denc, 1+nd);3716while(a != NIL)3717{3718for(i = 0; i <= nn; i++) /* Copy input numerator */3719numc[i] = num[i];3720if((dena = LIE_TERM_DENOMINATOR_INTEGER(a)) != NULL)3721{ /* n1/d2 */3722IntegerCancellation(numc, dena);3723if(INTEGER_IS_UNIT(dena))3724{3725INTEGER_KILL(dena);3726dena = NULL;3727}3728}3729numa = LIE_TERM_NUMERATOR_INTEGER(a);3730for(i = 0; i <= nd; i++) /* Copy input denominator */3731denc[i] = den[i];3732IntegerCancellation(numa, denc); /* n2/d1 */3733INTEGER_HEAP_NEW(w, 1+INTEGER_N_LIMBS(numc)+INTEGER_N_LIMBS(numa));3734IntegerProduct(w, numc, numa); /* n1*n2 */3735INTEGER_KILL(numa);3736LIE_TERM_NUMERATOR_INTEGER(a) = w;3737if(INTEGER_IS_UNIT(denc))3738if(dena == NULL)3739w = NULL; /* 1*1 */3740else3741{ /* Copy to avoid garbage tail in w */3742INTEGER_HEAP_COPY(w, dena, i); /* 1*d2 */3743INTEGER_KILL(dena);3744}3745else if(dena == NULL)3746{3747INTEGER_HEAP_COPY(w, denc, i); /* d1*1 */3748}3749else3750{3751INTEGER_HEAP_NEW(w, 1+INTEGER_N_LIMBS(denc)+INTEGER_N_LIMBS(dena));3752IntegerProduct(w, denc, dena); /* d1*d2 */3753INTEGER_KILL(dena);3754}3755LIE_TERM_DENOMINATOR_INTEGER(a) = w;3756a = LIE_TERM_R(a);3757}3758}3759#endif3760/*=LieSumMultScalarSum===============================================3761Multiply Lie sum by scalar sum on spot in Parametric regime3762`num' is killed3763*/3764void LieSumMultScalarSum(uint lsum, uint num)3765{3766if(lsum == NIL)3767ScalarSumKill(num);3768else3769{3770uint n, d, a;3771IN_LIE_SUM_MULT_SCALAR_SUM /*--------------------------------------*/3772do3773{3774a = lsum;3775lsum = LIE_TERM_R(lsum);3776n = (lsum != NIL) ? ScalarSumCopy(num) : num;3777if((d = LIE_TERM_DENOMINATOR_SCALAR_SUM(a)) != NIL)3778{3779ScalarSumCancellation(&n, &d);3780if(SCALAR_SUM_IS_UNIT(d))3781{3782ScalarSumKill(d); /* Kill unit */3783d = NIL;3784}3785LIE_TERM_DENOMINATOR_SCALAR_SUM(a) = d;3786}3787LIE_TERM_NUMERATOR_SCALAR_SUM(a) =3788ScalarSumMultiplication(LIE_TERM_NUMERATOR_SCALAR_SUM(a), n);3789}while(lsum != NIL);3790OUT_LIE_SUM_MULT_SCALAR_SUM /*-------------------------------------*/3791}3792}3793/*=NormalizeRelationInteger================================================3794Normalize sign, remove GCD of integer numerators,3795remove denominators for non-NIL relation3796*/3797void NormalizeRelationInteger(uint a)3798{3799uint b;3800BIGINT n2, n1 = LIE_TERM_NUMERATOR_INTEGER(a);3801IN_NORMALIZE_RELATION /*-----------------------------------------------*/3802/* Normalize sign */38033804if(INTEGER_IS_NEGATIVE(n1))3805{3806INTEGER_SET_PLUS(n1);3807b = LIE_TERM_R(a);3808while(b != NIL)3809{3810INTEGER_MINUS(LIE_TERM_NUMERATOR_INTEGER(b));3811b = LIE_TERM_R(b);3812}3813}3814#if defined(RATIONAL_FIELD) /* Field R case compiling */3815n2 = LIE_TERM_DENOMINATOR_INTEGER(a);3816if(INTEGER_IS_UNIT(n1)) /* Either 1 (nothing to do or 1/n2 */3817{3818if(n2 != NULL) /* Leading coefficient in form 1/n2 */3819{3820LIE_TERM_DENOMINATOR_INTEGER(a) = NULL;3821LieSumMultInteger(LIE_TERM_R(a), n2);3822INTEGER_KILL(n2); /* Free spoiled array n2 */3823}3824}3825else /* Leading coefficient is either n1 or n1/n2 */3826{3827if(n2 != NULL)3828{3829LIE_TERM_DENOMINATOR_INTEGER(a) = NULL;3830LieSumMultRationalInteger(LIE_TERM_R(a), n2, n1);3831INTEGER_KILL(n2); /* Free spoiled array n2 */3832}3833else /* Leading coefficient in form n1 */3834LieSumDivInteger(LIE_TERM_R(a), n1);3835n1[0] = n1[1] = 1; /* Set big number unit in old array */3836}3837#else /* Ring Z case compiling */3838{3839BIGINT gcd;3840int i;38413842/* Remove GCD of numerators */38433844if(INTEGER_IS_UNIT(n1))3845goto kill_denominators;3846INTEGER_STACK_COPY(gcd, n1, i);3847b = LIE_TERM_R(a);3848while(b != NIL)3849{3850n1 = LIE_TERM_NUMERATOR_INTEGER(b);3851if(INTEGER_IS_UNIT_ABS(n1))3852goto kill_denominators;3853INTEGER_HEAP_COPY(n2, n1, i); /* Working heap integer */3854gcd = IntegerGCD(gcd, n2);3855INTEGER_KILL(n2); /* Free heap integer */3856if(gcd == NULL)3857goto kill_denominators;3858b = LIE_TERM_R(b);3859}3860LieSumDivInteger(a, gcd);38613862/* Remove denominators */38633864kill_denominators:3865b = a;3866do3867if(LIE_TERM_DENOMINATOR_INTEGER(b) != NULL)3868{3869BIGINT n3, n4, lcm;38703871/* Make first LCM */38723873n1 = LIE_TERM_DENOMINATOR_INTEGER(b);3874INTEGER_HEAP_COPY(lcm, n1, i);3875while((b = LIE_TERM_R(b)) != NIL)3876if((n1 = LIE_TERM_DENOMINATOR_INTEGER(b)) != NULL)3877{3878/* Make copy of previous LCM */38793880INTEGER_HEAP_COPY(n3, lcm, i);38813882/* Make 2 copies of current denominator */38833884INTEGER_HEAP_COPY_DOUBLE_1(n2, n4, n1, i);38853886/* GCD of LCM and current denominator (stored in `n3') */38873888gcd = IntegerGCD(n3, n4);3889INTEGER_KILL(n4);38903891/* Divide current denominator by GCD */38923893if(gcd != NULL)3894{3895INTEGER_HEAP_NEW(n1, 2+INTEGER_N_LIMBS(n2)-INTEGER_N_LIMBS(gcd));3896IntegerQuotient(n1, n2, gcd);3897INTEGER_KILL(n2);3898}3899else3900n1 = n2;3901INTEGER_KILL(n3);39023903/* New LCM */39043905n2 = lcm;3906INTEGER_HEAP_NEW(lcm, 1+INTEGER_N_LIMBS(n1)+INTEGER_N_LIMBS(n2));3907IntegerProduct(lcm, n1, n2);3908INTEGER_KILL(n2);3909INTEGER_KILL(n1);3910}39113912/* Kill denominators */39133914LieSumMultInteger(a, lcm);3915INTEGER_KILL(lcm);3916break;3917}3918while((b = LIE_TERM_R(b)) != NIL);3919}3920#endif3921OUT_NORMALIZE_RELATION /*----------------------------------------------*/3922}3923/*=NormalizeRelationParametric===========================================3924Normalize sign, remove GCD of polynomial numerators, remove denominators3925for non-NIL relation, set in table common factor and leading coefficient3926*/3927void NormalizeRelationParametric(uint a)3928{3929uint b, c = LIE_TERM_NUMERATOR_SCALAR_SUM(a), d, e;3930IN_NORMALIZE_RELATION /*---------------------------------------------*/39313932/* Normalize sign */39333934if(INTEGER_IS_NEGATIVE(SCALAR_TERM_NUMERATOR(c)))3935LieSumMinusParametric(a);39363937/* Remove GCD of numerators */39383939if(SCALAR_SUM_IS_UNIT(c))3940goto kill_denominators;3941b = LIE_TERM_R(a);3942c = ScalarSumCopy(c);3943while(b != NIL)3944{3945d = LIE_TERM_NUMERATOR_SCALAR_SUM(b);3946if(SCALAR_SUM_IS_UNIT_ABS(d))3947{3948ScalarSumKill(c);3949goto kill_denominators;3950}3951e = c;3952c = PolyGCD(e, d);3953ScalarSumKill(e); /* Kill old GCD */3954if(c == NIL)3955goto kill_denominators;3956b = LIE_TERM_R(b);3957}3958if(NonZeroCoefficientsPut)3959InCoeffTable(ScalarSumCopy(c));3960LieSumDivScalarSum(a, c);39613962/* Remove denominators */39633964kill_denominators:3965b = a;3966do3967if(LIE_TERM_DENOMINATOR_SCALAR_SUM(b) != NIL)3968{39693970/* Make first LCM */39713972c = ScalarSumCopy(LIE_TERM_DENOMINATOR_SCALAR_SUM(b));3973while((b = LIE_TERM_R(b)) != NIL)3974if((d = LIE_TERM_DENOMINATOR_SCALAR_SUM(b)) != NIL)3975{3976/* GCD of LCM and current denominator */39773978e = PolyGCD(c, d);39793980/* Divide current denominator by GCD */39813982d = ScalarSumCopy(d);3983if(e != NIL)3984{3985d = PolyQuotient(d, e);3986ScalarSumKill(e);3987}39883989/* New LCM */39903991c = ScalarSumMultiplication(c, d);3992}39933994/* Kill denominators */39953996LieSumMultScalarSum(a, c);3997if(NonZeroCoefficientsPut)3998InCoeffTable(ScalarSumCopy(LIE_TERM_NUMERATOR_SCALAR_SUM(a)));3999break;4000}4001while((b = LIE_TERM_R(b)) != NIL);4002OUT_NORMALIZE_RELATION /*--------------------------------------------*/4003}4004/*=ScalarMonomialMultiplication============================================4005*/4006uint ScalarMonomialMultiplication(int *pchange_sign, uint ma, uint mb)4007{4008uint mc, wa, wb, last;4009*pchange_sign = NO;4010mc = NIL;4011while(YES)4012{4013next_pair:4014if(mb == NIL)4015{ /* List mb is ended, append rest of ma */4016if(mc == NIL)4017mc = ma;4018else4019SCALAR_FACTOR_R(last) = ma;4020break;4021}4022if(ma == NIL)4023{ /* List ma is ended, append rest of mb */4024if(mc == NIL)4025mc = mb;4026else4027SCALAR_FACTOR_R(last) = mb;4028break;4029}40304031/* Compare scalar factors */40324033if(SCALAR_FACTOR_PARAMETER(ma) > SCALAR_FACTOR_PARAMETER(mb))4034goto order_12;4035if(SCALAR_FACTOR_PARAMETER(ma) < SCALAR_FACTOR_PARAMETER(mb))4036goto order_21;40374038/* Reduce like factors */40394040wa = ma;4041wb = mb;4042ma = SCALAR_FACTOR_R(ma);4043mb = SCALAR_FACTOR_R(mb);4044if(SCALAR_FACTOR_IS_I_NUMBER(wa)) /* Imaginary unit i*i --> -1 */4045{4046NODE_SF_KILL(wa);4047NODE_SF_KILL(wb);4048if(*pchange_sign)4049*pchange_sign = NO; /* Convey change of sign to up */4050else4051*pchange_sign = YES;4052goto next_pair;4053}4054SCALAR_FACTOR_DEGREE(wa) += SCALAR_FACTOR_DEGREE(wb); /* Sum degrees */4055NODE_SF_KILL(wb);4056goto append_term;40574058order_12:4059wa = ma;4060ma = SCALAR_FACTOR_R(ma);4061goto append_term;40624063order_21:4064wa = mb;4065mb = SCALAR_FACTOR_R(mb);40664067append_term:4068if(mc == NIL)4069mc = wa;4070else4071SCALAR_FACTOR_R(last) = wa;4072last = wa;4073}4074return mc;4075}4076/*=ScalarSumAddition=====================================================4077Sum of two scalar (polynomial) expressions4078*/4079uint ScalarSumAddition(uint a, uint b)4080{4081uint sum = NIL, last, wa, wb, ma, mb;4082BIGINT na, nb, nc;4083while(YES)4084{4085next_pair:4086if(b == NIL)4087{ /* List b is ended, append rest of a */4088if(sum == NIL)4089sum = a;4090else4091SCALAR_TERM_R(last) = a;4092break;4093}4094if(a == NIL)4095{ /* List a is ended, append rest of b */4096if(sum == NIL)4097sum = b;4098else4099SCALAR_TERM_R(last) = b;4100break;4101}41024103/* Compare scalar terms */41044105ma = SCALAR_TERM_MONOMIAL(a);4106mb = SCALAR_TERM_MONOMIAL(b);4107while(YES)4108{4109if(ma == NIL)4110{4111if(mb != NIL)4112goto order_21; /* a-monomial is contained in b-monomial */4113break; /* a-monomial == b-monomial */4114} /* (including both are NILs) */4115if(mb == NIL)4116goto order_12; /* a-monomial contains b-monomial */4117if(SCALAR_FACTOR_WORD(ma) > SCALAR_FACTOR_WORD(mb))4118goto order_12;4119if(SCALAR_FACTOR_WORD(ma) < SCALAR_FACTOR_WORD(mb))4120goto order_21;4121ma = SCALAR_FACTOR_R(ma); /* Skip equal factors in monomials */4122mb = SCALAR_FACTOR_R(mb);4123}41244125/* Reduce like scalar terms */41264127wa = a;4128wb = b;4129a = SCALAR_TERM_R(a);4130b = SCALAR_TERM_R(b);41314132/* Sum integer coefficients */41334134na = SCALAR_TERM_NUMERATOR(wa);4135nb = SCALAR_TERM_NUMERATOR(wb);4136INTEGER_HEAP_NEW(nc, 2+MAX(INTEGER_N_LIMBS(na),INTEGER_N_LIMBS(nb)));4137IntegerSum(nc, na, nb);4138INTEGER_KILL(na);4139INTEGER_KILL(nb);41404141/* Kill head of wb and its monomial */41424143mb = SCALAR_TERM_MONOMIAL(wb);4144NODE_ST_KILL(wb);4145while(mb != NIL)4146{4147ma = mb;4148mb = SCALAR_FACTOR_R(mb);4149NODE_SF_KILL(ma);4150}41514152/* Nonzero collection */41534154if(INTEGER_N_LIMBS(nc) != 0)4155{4156SCALAR_TERM_NUMERATOR(wa) = nc;4157goto append_term;4158}41594160/* Zero collection: kill nc and wa */41614162INTEGER_KILL(nc);4163ma = SCALAR_TERM_MONOMIAL(wa);4164NODE_ST_KILL(wa);4165while(ma != NIL)4166{4167mb = ma;4168ma = SCALAR_FACTOR_R(ma);4169NODE_SF_KILL(mb);4170}4171goto next_pair;41724173order_12:4174wa = a;4175a = SCALAR_TERM_R(a);4176goto append_term;41774178order_21:4179wa = b;4180b = SCALAR_TERM_R(b);41814182append_term:4183if(sum == NIL)4184sum = wa;4185else4186SCALAR_TERM_R(last) = wa;4187last = wa;4188}4189return sum;4190}4191/*=ScalarSumCancellation===================4192*/4193void ScalarSumCancellation(uint *pnum, uint *pden)4194{4195uint g;4196IN_SCALAR_SUM_CANCELLATION /*------------*/4197if((g = PolyGCD(*pnum, *pden)) != NIL)4198{4199*pnum = PolyQuotient(*pnum, g);4200*pden = PolyQuotient(*pden, g);4201ScalarSumKill(g);4202}4203OUT_SCALAR_SUM_CANCELLATION /*-----------*/4204}4205/*=ScalarSumMinus===========================4206Change sign of scalar sum in Parametric regime4207*/4208void ScalarSumMinus(uint a)4209{4210while(a != NIL)4211{4212SCALAR_TERM_MINUS(a);4213a = SCALAR_TERM_R(a);4214}4215}4216/*=ScalarSumMultiplication===========================================4217Expanded product of two positive nonzero general scalar expressions,4218caller ensures A != NIL and B != NIL.4219*/4220uint ScalarSumMultiplication(uint a, uint b)4221{4222uint s, aw, bcur, bw, ac, bc;4223s = NIL;4224while(a != NIL)4225{4226aw = a;4227a = SCALAR_TERM_R(a);4228bcur = b;4229while(bcur != NIL)4230{4231bw = bcur;4232bcur = SCALAR_TERM_R(bcur);4233if(a == NIL)4234{4235bc = bw;4236SCALAR_TERM_R(bc) = NIL;4237}4238else4239bc = ScalarTermCopy(bw);4240if(bcur == NIL)4241{4242ac = aw;4243SCALAR_TERM_R(ac) = NIL;4244}4245else4246ac = ScalarTermCopy(aw);4247ScalarTermMultiplication(ac, bc);4248s = ScalarSumAddition(s, ac);4249}4250}4251return s;4252}4253/*=ScalarTermMultiplication===============================================4254Product of two scalar terms on place of A, B deleted4255*/4256void ScalarTermMultiplication(uint a, uint b)4257{4258BIGINT na, nb, nc;4259uint ma, mb, mc, last, aa;42604261/* Multiply integer coefficients */42624263na = SCALAR_TERM_NUMERATOR(a);4264nb = SCALAR_TERM_NUMERATOR(b);4265INTEGER_HEAP_NEW(nc, 1+INTEGER_N_LIMBS(na)+INTEGER_N_LIMBS(nb));4266IntegerProduct(nc, na, nb);4267INTEGER_KILL(na);4268INTEGER_KILL(nb);4269SCALAR_TERM_NUMERATOR(a) = nc;42704271/* Multiply monomials */42724273ma = SCALAR_TERM_MONOMIAL(a);4274mb = SCALAR_TERM_MONOMIAL(b);4275NODE_ST_KILL(b);4276mc = NIL;4277while(YES)4278{4279next_pair:4280if(mb == NIL)4281{ /* List mb is ended, append rest of ma */4282if(mc == NIL)4283mc = ma;4284else4285SCALAR_FACTOR_R(last) = ma;4286break;4287}4288if(ma == NIL)4289{ /* List ma is ended, append rest of mb */4290if(mc == NIL)4291mc = mb;4292else4293SCALAR_FACTOR_R(last) = mb;4294break;4295}42964297/* Compare scalar factors */42984299if(SCALAR_FACTOR_PARAMETER(ma) > SCALAR_FACTOR_PARAMETER(mb))4300goto order_12;4301if(SCALAR_FACTOR_PARAMETER(ma) < SCALAR_FACTOR_PARAMETER(mb))4302goto order_21;43034304/* Reduce like factors */43054306aa = ma;4307b = mb; /* uint-uint mixion */4308ma = SCALAR_FACTOR_R(ma);4309mb = SCALAR_FACTOR_R(mb);4310if(SCALAR_FACTOR_IS_I_NUMBER(aa)) /* Imaginary unit i*i --> -1 */4311{4312NODE_SF_KILL(aa);4313NODE_SF_KILL(b);4314INTEGER_MINUS(nc);4315goto next_pair;4316}4317SCALAR_FACTOR_DEGREE(aa) += SCALAR_FACTOR_DEGREE(b); /* Sum degrees */4318NODE_SF_KILL(b);4319goto append_term;43204321order_12:4322aa = ma;4323ma = SCALAR_FACTOR_R(ma);4324goto append_term;43254326order_21:4327aa = mb;4328mb = SCALAR_FACTOR_R(mb);43294330append_term:4331if(mc == NIL)4332mc = aa;4333else4334SCALAR_FACTOR_R(last) = aa;4335last = aa;4336}4337SCALAR_TERM_MONOMIAL(a) = mc;4338}43394340/*_6_4 Scalar polynomial algebraic functions===================*/43414342/*=ContentOfScalarSum=========================================4343Returns relative (or full if initial CONT == NIL) single term4344content of scalar sum. NIL corresponds to 1.4345CONT is destroyed, A remains.4346*/4347uint ContentOfScalarSum(uint cont, uint a)4348{4349if(cont == NIL)4350{4351cont = ScalarTermCopy(a);4352INTEGER_SET_PLUS(SCALAR_TERM_NUMERATOR(cont));4353if((a = SCALAR_TERM_R(a)) == NIL)4354goto out_cont;4355}4356while((cont = PolyTermGCD(cont, a)) != NIL4357&& (a = SCALAR_TERM_R(a)) != NIL)4358;4359out_cont:4360return cont;4361}4362/*=InCoeffParamTable===============================================4363Set in CoeffParamTable parameters of scalar term `cont' killing it4364*/4365void InCoeffParamTable(uint cont)4366{4367uint a = SCALAR_TERM_MONOMIAL(cont);4368INTEGER_KILL(SCALAR_TERM_NUMERATOR(cont));4369NODE_ST_KILL(cont);43704371if(a != NIL)4372{4373if(SCALAR_FACTOR_IS_I_NUMBER(a))4374{4375NODE_SF_KILL(a);4376return;4377}4378if(CoeffParamTable == NULL)4379CoeffParamTable = (int*)NewArray(ParameterN, sizeof(int),4380E_A_COEFF_PARA_TABLE);4381do4382{4383cont = a;4384a = SCALAR_FACTOR_R(a);4385CoeffParamTable[SCALAR_FACTOR_PARAMETER(cont)] = YES;4386NODE_SF_KILL(cont);4387}while(a != NIL);4388}4389}4390/*=InCoeffSumTable========================================================4391Insert parametric content-free SUM in table or delete if already exists4392*/4393void InCoeffSumTable(uint sum)4394{4395if(SCALAR_FACTOR_IS_I_NUMBER(SCALAR_TERM_MONOMIAL(sum)))4396ScalarSumKill(sum); /* Kill complex number a*i + b */4397else4398{4399int i;4400uint gcd, quocoe, quosum;4401for(i = 0; i < CoeffSumTableN; i++)4402if(PolynomialsAreEqual(sum, CoeffSumTable[i]))4403{4404ScalarSumKill(sum);4405return;4406}4407else4408{4409gcd = PolyGCD(sum, CoeffSumTable[i]);4410if(gcd != NIL)4411{4412quocoe = PolyQuotient(CoeffSumTable[i], gcd);4413quosum = PolyQuotient(sum, gcd);4414--CoeffSumTableN; /* Remove gap in ith position */4415while(i < CoeffSumTableN)4416{4417CoeffSumTable[i] = CoeffSumTable[i+1];4418++i;4419}4420InCoeffSumTable(gcd);4421if(SCALAR_TERM_R(quocoe) == NIL) /* Either sum or 1 certainly */4422{4423INTEGER_KILL(SCALAR_TERM_NUMERATOR(quocoe)); /* Kill 1 */4424NODE_ST_KILL(quocoe);4425}4426else4427InCoeffSumTable(quocoe);4428if(SCALAR_TERM_R(quosum) == NIL) /* Either sum or 1 certainly */4429{4430INTEGER_KILL(SCALAR_TERM_NUMERATOR(quosum)); /* Kill 1 */4431NODE_ST_KILL(quosum);4432}4433else4434InCoeffSumTable(quosum);4435return;4436}4437}4438if(CoeffSumTableN >= CoeffSumTableSize)4439Error(E_COEFF_SUM_TABLE_SIZE);4440if(CoeffSumTable == NULL)4441CoeffSumTable = (uint *)NewArray(CoeffSumTableSize, sizeof(uint),4442E_A_COEFF_SUM_TABLE);4443CoeffSumTable[CoeffSumTableN++] = sum;4444}4445}4446/*=InCoeffTable===========================================4447Set in tables components of non-NIL parametric polynomial4448*/4449void InCoeffTable(uint coe)4450{4451if(SCALAR_TERM_R(coe) != NIL)4452{4453uint cont;4454if((cont = ContentOfScalarSum(NIL, coe)) != NIL)4455{4456coe = PolyQuotient(coe, cont);4457InCoeffParamTable(cont);4458}4459InCoeffSumTable(coe);4460}4461else4462InCoeffParamTable(coe);4463}4464/*=PolyCoeffAtMainParameter===========================================4465Polynomial coefficient with normalized sign at the current degree of4466the main parameter. *PA initially points to the start of list, at the4467end of work it points to tail of list. NIL corresponds to 1.4468Initial polynomial remains.4469*/4470uint PolyCoeffAtMainParameter(uint *pa, int mp)4471{4472uint a;4473int isnegative = INTEGER_IS_NEGATIVE(SCALAR_TERM_NUMERATOR(*pa));4474if(POLY_MAIN_PARAMETER(*pa) < mp)4475{ /* Free term */4476a = ScalarSumCopy(*pa);4477if(isnegative) /* Normalize sign */4478{4479*pa = a;4480do4481SCALAR_TERM_MINUS(*pa);4482while((*pa = SCALAR_TERM_R(*pa)) != NIL);4483}4484else4485*pa = NIL;4486}4487else4488{4489uint b;4490uint mf;4491int mppow = SCALAR_TERM_MAIN_DEGREE(*pa);4492a = b = ScalarTermCopy(*pa);4493if(isnegative)4494SCALAR_TERM_MINUS(a);4495mf = SCALAR_TERM_MONOMIAL(a); /* Strike out main parameter */4496SCALAR_TERM_MONOMIAL(a) = SCALAR_FACTOR_R(mf);4497NODE_SF_KILL(mf);4498while((*pa = SCALAR_TERM_R(*pa)) != NIL4499&& POLY_MAIN_PARAMETER(*pa) == mp4500&& SCALAR_TERM_MAIN_DEGREE(*pa) == mppow)4501{4502SCALAR_TERM_R(b) = ScalarTermCopy(*pa);4503b = SCALAR_TERM_R(b);4504if(isnegative)4505{4506SCALAR_TERM_MINUS(b);4507}4508mf = SCALAR_TERM_MONOMIAL(b); /* Strike out main parameter */4509SCALAR_TERM_MONOMIAL(b) = SCALAR_FACTOR_R(mf);4510NODE_SF_KILL(mf);4511}4512}4513if(SCALAR_SUM_IS_UNIT(a))4514{4515INTEGER_KILL(SCALAR_TERM_NUMERATOR(a));4516NODE_ST_KILL(a);4517a = NIL;4518}4519return a;4520}4521/*=PolyContent=============================================4522Polynomial content of polynomial w.r.t. main parameter MP.4523A remains unchanged.4524*/4525uint PolyContent(uint a, int mp)4526{4527uint b;4528IN_POLY_CONTENT /*---------------------------------------*/4529if((b = PolyCoeffAtMainParameter(&a, mp)) != NIL)4530{4531uint c, d;4532while(a != NIL)4533{4534if((c = PolyCoeffAtMainParameter(&a, mp)) == NIL)4535{4536ScalarSumKill(b);4537b = NIL;4538break;4539}4540d = b;4541if((b = PolyGCD(b, c)) == NIL)4542{4543ScalarSumKill(d);4544ScalarSumKill(c);4545break;4546}4547}4548}4549OUT_POLY_CONTENT /*--------------------------------------*/4550return b;4551}4552/*=PolyGCD==============================================================4553Returns Greatest Common Divisor of two multivariate polynomials in4554the form GCD(PP(A), PP(B)) * GCD(CONT(A), CONT(B)).4555A, B unchanged.4556Returned NIL means trivial GCD = 14557*/4558uint PolyGCD(uint a, uint b)4559{4560uint c;4561IN_POLY_GCD /*--------------------------------------------------------*/4562if(SCALAR_TERM_R(a) == NIL || SCALAR_TERM_R(b) == NIL)4563{ /* At least one of the polynomials is not a sum */4564if(SCALAR_TERM_R(a) != NIL)4565{ /* Set A to be a single term */4566c = a;4567a = b;4568b = c;4569}4570c = ScalarSumCopy(a);4571INTEGER_SET_PLUS(SCALAR_TERM_NUMERATOR(c));4572b = ContentOfScalarSum(c, b);4573}4574else /* Both are polynomials really */4575{4576uint conta, contb;4577int mp, mpb; /* Main parameters */4578mp = SCALAR_FACTOR_PARAMETER(SCALAR_TERM_MONOMIAL(a));4579mpb = SCALAR_FACTOR_PARAMETER(SCALAR_TERM_MONOMIAL(b));4580if(mpb > mp ||4581(mpb == mp &&4582SCALAR_TERM_MAIN_DEGREE(b) > SCALAR_TERM_MAIN_DEGREE(a)))4583{ /* Parameters go in DECREASING order! */4584c = a; /* Swap polynomials */4585a = b;4586b = c;4587mp = mpb;4588}4589a = ScalarSumCopy(a);4590b = ScalarSumCopy(b);4591contb = PolyContent(b, mp);4592if((conta = PolyContent(a, mp)) != NIL)4593{ /* Make primitive parts and GCD of contents */4594if(contb != NIL)4595{4596c = PolyGCD(conta, contb);4597b = PolyQuotient(b, contb); /* Primitive part */4598ScalarSumKill(contb);4599}4600else4601c = NIL;4602a = PolyQuotient(a, conta); /* Primitive part */4603ScalarSumKill(conta);4604}4605else4606{4607if(contb != NIL)4608{4609b = PolyQuotient(b, contb); /* Primitive part */4610ScalarSumKill(contb);4611}4612c = NIL;4613}4614while((conta = PolyPseudoRemainder(a, ScalarSumCopy(b), mp)) != NIL)4615{4616if(SCALAR_TERM_MONOMIAL(conta) == NIL || /* Pure number */4617SCALAR_FACTOR_PARAMETER(SCALAR_TERM_MONOMIAL(conta)) != mp)4618{ /* Zero degree with respect to main parameter */4619ScalarSumKill(b);4620ScalarSumKill(conta);4621b = c; /* char is content ?? */4622goto out;4623}4624a = b;4625if((contb = ContentOfScalarSum(NIL, conta)) != NIL)4626{ /* Term content */4627conta = PolyQuotient(conta, contb);4628ScalarSumKill(contb);4629}4630if((contb = PolyContent(conta, mp)) != NIL)4631{ /* Polynomial content */4632b = PolyQuotient(conta, contb); /* Primitive part */4633ScalarSumKill(contb);4634}4635else4636b = conta;4637}4638if(c != NIL)4639b = ScalarSumMultiplication(b, c);4640if(INTEGER_IS_NEGATIVE(SCALAR_TERM_NUMERATOR(b)))4641{ /* Standardize sign of GCD */4642c = b;4643do4644SCALAR_TERM_MINUS(c);4645while((c = SCALAR_TERM_R(c)) != NIL);4646}4647if(SCALAR_SUM_IS_UNIT(b))4648{ /* Standardize trivial GCD */4649INTEGER_KILL(SCALAR_TERM_NUMERATOR(b));4650NODE_ST_KILL(b);4651b = NIL;4652}4653}4654out:4655OUT_POLY_GCD /*-------------------------------------------------------*/4656return b;4657}4658/*=PolyMainParameterTerm================================================4659Take sublist of terms, containing main parameter MP of given degree4660MPDEG, *PA points to the next term after end of A.4661This function is applied in succession starting from top degree.4662Initial expression *PA is destructed.4663*/4664uint PolyMainParameterTerm(uint *pa, int mp, int mpdeg)4665{4666uint a;4667if(mpdeg)4668{4669unsigned short w; /* Word combining degree and index of parameter (MPDEG,MP) */4670w = mpdeg;4671#if defined(SPP_2000)4672*((byte*)&w) = mp;4673#else4674*((byte*)&w+1) = mp;4675#endif4676if(SCALAR_TERM_MONOMIAL(*pa) != NIL &&4677SCALAR_TERM_MAIN_PARAMETER_WORD(*pa) == w)4678{ /* There is such degree */4679uint b;4680a = *pa;4681while(YES)4682{4683b = *pa; /* Remember previous term */4684*pa = SCALAR_TERM_R(*pa);4685if(*pa == NIL)4686break; /* Whole expression is homogeneous */4687if(SCALAR_TERM_MONOMIAL(*pa) == NIL ||4688SCALAR_TERM_MAIN_PARAMETER_WORD(*pa) != w)4689{ /* End of homogeneous part is found */4690SCALAR_TERM_R(b) = NIL; /* Set end of sublist */4691break;4692}4693}4694}4695else /* No such degree */4696a = NIL;4697}4698else4699{ /* Free term */4700a = *pa;4701*pa = NIL;4702}4703return a;4704}4705/*=PolynomialsAreEqual====================================4706*/4707int PolynomialsAreEqual(uint a, uint b)4708{4709uint ma, mb;4710BIGINT na, nb;4711while(YES)4712{4713/* Compare monomials */47144715ma = SCALAR_TERM_MONOMIAL(a);4716mb = SCALAR_TERM_MONOMIAL(b);4717while(YES)4718{4719if(ma == NIL)4720if(mb == NIL)4721break;4722else4723goto no;4724else if(mb == NIL)4725goto no;4726if(SCALAR_FACTOR_WORD(ma) != SCALAR_FACTOR_WORD(mb))4727goto no;4728ma = SCALAR_FACTOR_R(ma);4729mb = SCALAR_FACTOR_R(mb);4730}4731/* Compare numerators */47324733na = SCALAR_TERM_NUMERATOR(a);4734nb = SCALAR_TERM_NUMERATOR(b);4735if(na[0] != nb[0])4736goto no;4737ma = INTEGER_N_LIMBS(na);4738do4739if(na[ma] != nb[ma])4740goto no;4741while(ma-- > 0);47424743a = SCALAR_TERM_R(a);4744b = SCALAR_TERM_R(b);4745if(a == NIL)4746if(b == NIL)4747return YES;4748else4749goto no;4750else if(b == NIL)4751goto no;4752}4753no:4754return NO;4755}4756/*=PolyPseudoRemainder==================================================4757Returns pseudo-remainder of two polynomials. MP is main parameter.4758main_degree(A) >= main_degree(B). A, B destructed.4759*/4760uint PolyPseudoRemainder(uint a, uint b, int mp)4761{4762IN_POLY_PSEUDO_REMAINDER /*-------------------------------------------*/4763if(SCALAR_TERM_MAIN_PARAMETER(b) != mp)4764{ /* B doesn't contain MP => return 0 (NIL) */4765ScalarSumKill(a);4766ScalarSumKill(b);4767a = NIL;4768}4769else4770{4771uint *u, *v, vn, c, w;4772int m, n, j, k;4773m = SCALAR_TERM_MAIN_DEGREE(a);4774n = SCALAR_TERM_MAIN_DEGREE(b);4775POLY_ARRAY_STACK_NEW(u, m+1);4776POLY_ARRAY_STACK_NEW(v, n);4777for(j = m; j >= 0; j--) /* j */4778u[j] = PolyMainParameterTerm(&a, mp, j); /* u[j] = mp u etc. */4779vn = PolyMainParameterTerm(&b, mp, n); /* j */4780for(j = n - 1; j >= 0; j--)4781v[j] = PolyMainParameterTerm(&b, mp, j);4782for(k = m - n; k >= 0; k--)4783{4784j = n + k - 1;4785while(j >= k) /* n+j */4786{ /* mp (u = v u ) */4787if(u[j] != NIL) /* j n j */4788u[j] = ScalarSumMultiplication(ScalarSumCopy(vn), u[j]);4789if(u[n+k] != NIL)4790if(v[j-k] != NIL)4791{4792a = ScalarSumMultiplication(ScalarSumCopy(v[j-k]),4793ScalarSumCopy(u[n+k]));4794b = a;4795do4796SCALAR_TERM_MINUS(b);4797while((b = SCALAR_TERM_R(b)) != NIL); /* n+j */4798u[j] = ScalarSumAddition(u[j], a); /* - mp v u */4799} /* j-k n+k */4800if(u[j] != NIL)4801{4802c = u[j]; /* Drop degree of main parameter */4803do4804if((SCALAR_TERM_MAIN_DEGREE(c) -= n) == 0)4805{ /* Strike out main parameter */4806w = SCALAR_TERM_MONOMIAL(c);4807SCALAR_TERM_MONOMIAL(c) = SCALAR_FACTOR_R(w);4808NODE_SF_KILL(w);4809}4810while((c = SCALAR_TERM_R(c)) != NIL);4811}4812j--;4813}4814if(u[n+k] != NIL)4815ScalarSumKill(u[n+k]);4816while(j >= 0)4817{4818if(u[j] != NIL)4819{4820u[j] = ScalarSumMultiplication(ScalarSumCopy(vn), u[j]);4821c = u[j]; /* Drop degree of main parameter */4822do4823if((SCALAR_TERM_MAIN_DEGREE(c) -= n) == 0)4824{ /* Strike out main parameter */4825w = SCALAR_TERM_MONOMIAL(c);4826SCALAR_TERM_MONOMIAL(c) = SCALAR_FACTOR_R(w);4827NODE_SF_KILL(w);4828}4829while((c = SCALAR_TERM_R(c)) != NIL);4830}4831j--;4832}4833}4834ScalarSumKill(vn);4835for(j = n - 1; j >= 0; j--)4836if(v[j] != NIL)4837ScalarSumKill(v[j]);4838j = n - 1;4839while(j >= 0 && u[j] == NIL)4840j--; /* Search first nonzero term u[j] */4841if(j >= 0)4842{ /* Concatenate pseudoremainder from array u[] */4843a = u[j--];4844b = a;4845while(SCALAR_TERM_R(b) != NIL)4846b = SCALAR_TERM_R(b);4847while(j >= 0)4848{4849if(u[j] != NIL)4850{4851SCALAR_TERM_R(b) = u[j];4852while(SCALAR_TERM_R(b) != NIL)4853b = SCALAR_TERM_R(b);4854}4855j--;4856}4857}4858else /* All u[j] are zeros */4859a = NIL;4860}4861OUT_POLY_PSEUDO_REMAINDER /*------------------------------------------*/4862return a;4863}4864/*=PolyTermGCD==========================================================4865GCD of two single (non-NIL) terms, A is destroyed, B remains,4866caller ensures A is positive, returned NIL corresponds to 14867*/4868uint PolyTermGCD(uint a, uint b)4869{4870BIGINT na, naa, nb, nbb;4871int i;4872uint ma, maa;48734874/* Do integer coefficients */48754876na = SCALAR_TERM_NUMERATOR(a);4877nb = SCALAR_TERM_NUMERATOR(b);4878naa = na; /* Anyway it will be killed */4879INTEGER_STACK_COPY(nbb, nb, i);4880if((naa = IntegerGCD(naa, nbb)) != NULL) /* naa = GCD(na, nb) > 1 */4881{4882INTEGER_HEAP_COPY(nb, naa, i);4883SCALAR_TERM_NUMERATOR(a) = nb;4884}4885INTEGER_KILL(na);48864887/* Do parametric monomials (parameters go in DECREASING order) */48884889maa = NIL;4890if((ma = SCALAR_TERM_MONOMIAL(a)) != NIL)4891{4892uint maw, mb, mal;4893mb = SCALAR_TERM_MONOMIAL(b);4894while(YES)4895{4896if(mb == NIL)4897{ /* MB is ended - kill tail of MA and break loop */4898while(ma != NIL)4899{4900maw = ma;4901ma = SCALAR_FACTOR_R(ma);4902NODE_SF_KILL(maw);4903}4904if(maa != NIL)4905SCALAR_FACTOR_R(mal) = NIL;4906break;4907}4908if(SCALAR_FACTOR_PARAMETER(ma) > SCALAR_FACTOR_PARAMETER(mb))4909{4910maw = ma; /* No match for MA - kill it */4911ma = SCALAR_FACTOR_R(ma);4912NODE_SF_KILL(maw);4913if(ma == NIL)4914{4915if(maa != NIL)4916SCALAR_FACTOR_R(mal) = NIL;4917break;4918}4919}4920else if(SCALAR_FACTOR_PARAMETER(ma) < SCALAR_FACTOR_PARAMETER(mb))4921mb = SCALAR_FACTOR_R(mb); /* No match for MB - shift it */4922else4923{ /* Match - set minimum degree */4924if(SCALAR_FACTOR_DEGREE(mb) < SCALAR_FACTOR_DEGREE(ma))4925SCALAR_FACTOR_DEGREE(ma) = SCALAR_FACTOR_DEGREE(mb);4926if(maa == NIL)4927maa = ma;4928else /* Append to last MAA */4929SCALAR_FACTOR_R(mal) = ma;4930mal = ma;4931ma = SCALAR_FACTOR_R(ma);4932if(ma == NIL)4933break;4934}4935}4936}4937SCALAR_TERM_MONOMIAL(a) = maa; /* Set constructed monomial */4938if(naa == NULL)4939{4940if(maa == NIL)4941{ /* Trivial GCD */4942NODE_ST_KILL(a);4943a = NIL;4944}4945else /* Make standard integer coefficient */4946{4947INTEGER_HEAP_NEW(na, 2);4948na[0] = na[1] = 1;4949SCALAR_TERM_NUMERATOR(a) = na;4950}4951}4952return a;4953}4954/*=PolyTermQuotient===================================================4955Exact division of term A by term B: A = char*B on place of A, B remains.4956Parameters go in decreasing order.4957*/4958void PolyTermQuotient(uint a, uint b)4959{4960BIGINT na, nb, naa, nbb, nc;4961int i;4962uint mb;49634964/* Divide integer numerator */49654966na = SCALAR_TERM_NUMERATOR(a);4967nb = SCALAR_TERM_NUMERATOR(b);4968INTEGER_STACK_COPY(nbb, nb, i);4969INTEGER_HEAP_NEW(nc, 2+INTEGER_N_LIMBS(na)-INTEGER_N_LIMBS(nbb));4970INTEGER_STACK_COPY_1(naa, na, i);4971INTEGER_KILL(na);4972IntegerQuotient(nc, naa, nbb);4973SCALAR_TERM_NUMERATOR(a) = nc;49744975/* Divide parametric monomial */49764977if((mb = SCALAR_TERM_MONOMIAL(b)) != NIL)4978{4979uint ma, maa, mae, maw;4980ma = SCALAR_TERM_MONOMIAL(a);4981maa = NIL;4982do4983{4984while(SCALAR_FACTOR_PARAMETER(ma) > SCALAR_FACTOR_PARAMETER(mb))4985{4986if(maa == NIL)4987maa = ma;4988else4989SCALAR_FACTOR_R(mae) = ma;4990mae = ma;4991ma = SCALAR_FACTOR_R(ma);4992}4993if((SCALAR_FACTOR_DEGREE(ma) -= SCALAR_FACTOR_DEGREE(mb)) != 0)4994{4995if(maa == NIL)4996maa = ma;4997else4998SCALAR_FACTOR_R(mae) = ma;4999mae = ma;5000ma = SCALAR_FACTOR_R(ma);5001}5002else5003{5004maw = ma;5005ma = SCALAR_FACTOR_R(ma);5006NODE_SF_KILL(maw);5007}5008}while((mb = SCALAR_FACTOR_R(mb)) != NIL);5009if(maa == NIL) /* Append tail of numerator */5010maa = ma;5011else5012SCALAR_FACTOR_R(mae) = ma;5013SCALAR_TERM_MONOMIAL(a) = maa;5014}5015}5016/*=PolyQuotient====================================================5017Exact division of polynomial A by polynomial B: A = char*B, return char.5018Caller ensures A, B != NIL, B is positive.5019A is destructed, B remains unchanged.5020*/5021uint PolyQuotient(uint a, uint b)5022{5023uint c;5024IN_POLY_QUOTIENT /*----------------------------------------------*/5025if(SCALAR_TERM_R(b) == NIL) /* Division by single term B */5026{5027c = a;5028if(SCALAR_SUM_IS_NOT_UNIT(b)) /* Nontrivial term */5029do5030PolyTermQuotient(a, b);5031while((a = SCALAR_TERM_R(a)) != NIL);5032}5033else /* Division by polynomial B */5034{5035uint aw, bw, cw;5036BIGINT n;5037bw = SCALAR_TERM_R(b);5038c = NIL;5039do5040{5041aw = a;5042a = SCALAR_TERM_R(a);5043SCALAR_TERM_R(aw) = NIL;5044PolyTermQuotient(aw, b);5045cw = ScalarTermCopy(aw);5046n = SCALAR_TERM_NUMERATOR(aw);5047INTEGER_MINUS(n);5048aw = ScalarSumMultiplication(aw, ScalarSumCopy(bw));5049a = ScalarSumAddition(a, aw); /* Remainder of A */5050c = ScalarSumAddition(c, cw); /* Quotient char */5051}while(a != NIL);5052}5053OUT_POLY_QUOTIENT /*---------------------------------------------*/5054return c;5055}50565057/*_6_5 Big number functions====================================*/50585059/*=BigNMinusBigN================================================5060Subtract two big numbers on the place of first one: `a' -= `b',5061caller provides `a' > `b', returns new size of `a'5062*/5063int BigNMinusBigN(BIGINT a, int na, BIGINT b, int nb)5064{5065uint lw;5066LIMB k = 1;5067int i = 0;5068while(i < nb) /* Common part */5069{5070lw = MAX_LIMB + (uint)k + (uint)a[i] - (uint)b[i];5071k = (lw > MAX_LIMB);5072a[i++] = (LIMB)lw;5073}5074while(i < na)5075{5076lw = MAX_LIMB + (uint)k + (uint)a[i];5077k = (lw > MAX_LIMB);5078a[i++] = (LIMB)lw;5079}5080while(--i >= 0)5081if(a[i] != 0)5082break;5083return (++i);5084}5085/*=BigNShiftLeft==============================================5086Add on spot 0 <= `cnt' < BITS_PER_LIMB lowest zero bits to5087`bign' of size `n', return the bits shifted out from the most5088significant LIMB digit5089*/5090LIMB BigNShiftLeft(BIGINT bign, int n, int cnt)5091{5092if(cnt)5093{5094int cocnt = BITS_PER_LIMB - cnt;5095LIMB low_limb,5096high_limb = bign[--n],5097pushed_out = high_limb >> cocnt;5098while(--n >= 0)5099{5100low_limb = bign[n];5101bign[n+1] = (high_limb << cnt) | (low_limb >> cocnt);5102high_limb = low_limb;5103}5104bign[n+1] = high_limb << cnt;5105return pushed_out;5106}5107return 0;5108}5109/*=BigNShiftRight==========================================5110Remove on spot 0 <= `cnt' < BITS_PER_LIMB lowest bits from5111`bign' of size `n', return size of the result5112*/5113int BigNShiftRight(BIGINT bign, int n, int cnt)5114{5115if(cnt)5116{5117BIGINT bigni;5118int high_limb, low_limb, i, cocnt = BITS_PER_LIMB - cnt;5119low_limb = *bign;5120bigni = bign;5121bigni++;5122for(i = n-1; i > 0; i--)5123{5124high_limb = *(bigni++);5125*(bign++) = (low_limb >> cnt) | (high_limb << cocnt);5126low_limb = high_limb;5127}5128low_limb >>= cnt;5129if(low_limb != 0)5130*bign = low_limb;5131else5132n--;5133#if 05134LIMB high_limb, low_limb;5135int i, cocnt = BITS_PER_LIMB - cnt;5136low_limb = bign[0];5137for(i = 1; i < n; i++)5138{5139high_limb = bign[i];5140bign[i-1] = (low_limb >> cnt) | (high_limb << cocnt);5141low_limb = high_limb;5142}5143low_limb >>= cnt;5144if(low_limb != 0)5145bign[i-1] = low_limb;5146else5147n--;5148#endif5149}5150return n;5151}5152/*=CountLeadingZeroBitsInLimb=========================5153Count number of leading zero bits in LIMB word5154*/5155int CountLeadingZeroBitsInLimb(LIMB w)5156{5157if(w >= 0x100u) /* [0, 7] */5158if(w >= 0x1000u) /* [0, 3] */5159if(w >= 0x4000u) /* [0, 1] */5160if(w >= 0x8000u)5161return 0;5162else5163return 1;5164else /* [2, 3] */5165if(w >= 0x2000u)5166return 2;5167else5168return 3;5169else /* [4, 7] */5170if(w >= 0x400u) /* [4, 5] */5171if(w >= 0x800u)5172return 4;5173else5174return 5;5175else /* [6, 7] */5176if(w >= 0x200u)5177return 6;5178else5179return 7;5180else /* [8, 16] */5181if(w >= 0x10u) /* [ 8, 11] */5182if(w >= 0x40u) /* [ 8, 9] */5183if(w >= 0x80u)5184return 8;5185else5186return 9;5187else /* [10, 11] */5188if(w >= 0x20u)5189return 10;5190else5191return 11;5192else /* [12, 16] */5193if(w >= 0x4u) /* [12, 13] */5194if(w >= 0x8u)5195return 12;5196else5197return 13;5198else /* [14, 16] */5199if(w >= 0x2u)5200return 14;5201else /* [15, 16] */5202if(w)5203return 15;5204else5205return 16;5206}5207/*=IntegerCancellation========================================5208Results are placed in `num' and `den' arrays5209*/5210void IntegerCancellation(BIGINT num, BIGINT den)5211{5212BIGINT n, d, g;5213int i;5214IN_INTEGER_CANCELLATION /*----------------------*/5215INTEGER_STACK_COPY_1(n, num, i);5216INTEGER_STACK_COPY_1(d, den, i);5217if((g = IntegerGCD(num, den)) != NULL)5218{5219BIGINT gg;5220INTEGER_STACK_COPY(gg, g, i); /* `g' in `num' */52215222/* Cancel `den' */52235224IntegerQuotient(den, d, g);52255226/* Cancel `num' */52275228IntegerQuotient(num, n, gg);5229}5230else /* Lozh' vzad */5231{5232i = INTEGER_N_LIMBS(n);5233do5234num[i] = n[i];5235while(i--);5236i = d[0]; /* Denominators and GCDs are positive always */5237do5238den[i] = d[i];5239while(i--);5240}5241OUT_INTEGER_CANCELLATION /*----------------------*/5242}5243/*=IntegerGCD=========================================================5244Binary algorithm is used,5245Returns pointer to array of greatest common divisor5246or NULL interpreted as 1 by caller,5247the function spoils both arrays `u' and `v',5248result is placed in `u' array5249*/5250BIGINT IntegerGCD(BIGINT u, BIGINT v)5251{5252int i, nu, nv, bcnt, w_bcnt;5253BIGINT u0;5254LIMB carry_digit;5255IN_INTEGER_GCD /*--------------------------------------------------*/5256nu = INTEGER_N_LIMBS(u);5257nv = INTEGER_N_LIMBS(v);5258u0 = ++u; /* Skip size information limbs and memorize begin */5259++v;5260i = 0; /* Shift down uint to make it odd */5261while(*u == 0)5262{ /* Skip zero limbs */5263++i;5264++u;5265}5266COUNT_LEADING_ZERO_BITS_IN_LIMB(bcnt, *u & -*u);5267bcnt = BITS_PER_LIMB - 1 - bcnt;5268nu = BigNShiftRight(u, nu - i, bcnt);5269bcnt += i * BITS_PER_LIMB;5270w_bcnt = bcnt;5271i = 0; /* Shift down void to make it odd */5272while(*v == 0)5273{ /* Skip zero limbs */5274++i;5275++v;5276}5277COUNT_LEADING_ZERO_BITS_IN_LIMB(bcnt, *v & -*v);5278bcnt = BITS_PER_LIMB - 1 - bcnt;5279nv = BigNShiftRight(v, nv - i, bcnt);5280bcnt += i * BITS_PER_LIMB;5281if(bcnt < w_bcnt)5282w_bcnt = bcnt; /* Number of common 2 factors. */5283while(YES)5284{5285if(nu > nv)5286goto u_greater_v;5287if(nu < nv)5288goto v_greater_u;5289i = nu;5290while(--i >= 0)5291{5292if(u[i] > v[i])5293goto u_greater_v;5294if(v[i] > u[i])5295goto v_greater_u;5296}5297break; /* If uint and void have become equal, we have found the GCD */5298u_greater_v: /* Replace uint by (uint - void) >> cnt making uint odd again */5299nu = BigNMinusBigN(u, nu, v, nv);5300while(*u == 0)5301{5302--nu;5303++u;5304}5305COUNT_LEADING_ZERO_BITS_IN_LIMB(bcnt, *u & -*u);5306bcnt = BITS_PER_LIMB - 1 - bcnt;5307nu = BigNShiftRight(u, nu, bcnt);5308continue;5309v_greater_u: /* Replace void by (void - uint) >> cnt making void odd again */5310nv = BigNMinusBigN(v, nv, u, nu);5311while(*v == 0)5312{5313--nv;5314++v;5315}5316COUNT_LEADING_ZERO_BITS_IN_LIMB(bcnt, *v & -*v);5317bcnt = BITS_PER_LIMB - 1 - bcnt;5318nv = BigNShiftRight(v, nv, bcnt);5319}5320/* GCD(U_IN, V_IN) now is uint * 2**W_BCNT. */5321carry_digit = BigNShiftLeft(u, nu, w_bcnt % BITS_PER_LIMB);5322i = w_bcnt / BITS_PER_LIMB;5323u -= i;5324i += nu;5325if(carry_digit != 0)5326{5327if(u > u0)5328{5329u0 = u--; /* Shift to left by 1 limb to make room for carry */5330for(nu = 0; nu < i; nu++)5331u[nu] = u0[nu];5332}5333u[i++] = carry_digit;5334}5335if(i == 1 && u[0] == 1)5336return NULL;5337--u;5338u[0] = i; /* PLUS == 0 assumed */5339OUT_INTEGER_GCD /*--------------------------------------------------*/5340return u;5341}5342/*=IntegerProduct======================================5343Traditional multiplication of two signed non-zero5344big numbers uint and void, result in W, W[] != uint[] or void[]5345*/5346void IntegerProduct(BIGINT w, BIGINT u, BIGINT v)5347{5348LIMB carry;5349uint luw;5350BIGINT w0;5351int set_minus, i, j, n, m;5352IN_INTEGER_PRODUCT /*-------------------------------*/5353n = INTEGER_N_LIMBS(u);5354m = INTEGER_N_LIMBS(v);5355set_minus = (INTEGER_SIGN(u) != INTEGER_SIGN(v));5356u++;5357v++;5358w0 = w;5359w++;5360i = j = 0;5361do5362w[i] = 0;5363while(++i < n);5364do5365{5366i = carry = 0;5367do5368{5369luw = (uint)u[i]*(uint)v[j] + (uint)w[i+j] + (uint)carry;5370w[i+j] = (LIMB)luw;5371carry = (LIMB)(luw/BASE_LIMB);5372}while(++i < n);5373w[j+n] = carry;5374}while(++j < m);5375n += m - (carry == 0);5376#if defined(INTEGER_MAX_SIZE)5377if(n > IntegerMaxSize)5378IntegerMaxSize = n;5379#endif5380w0[0] = n;5381if(set_minus)5382INTEGER_SET_MINUS(w0);5383OUT_INTEGER_PRODUCT /*-------------------------------*/5384}5385/*=IntegerQuotient=====================================================5386Exact division of big numbers:5387Quotient in char[*PM - N + 1 or *PM - N stored in *PM] = A[M] / B[N],5388char != A, char != B.5389Function is spoiling input A and B.5390Array for A should have 1 additional LIMB at the top5391for increasing A at normalizing B.5392*/5393void IntegerQuotient(BIGINT c, BIGINT a, BIGINT b)5394{5395uint lw;5396LIMB q;5397int i, n, set_minus = (INTEGER_SIGN(a) != INTEGER_SIGN(b));5398#if defined(D_CHECK_EXACTNESS_OF_DIVISION)5399int nr;5400#endif5401BIGINT pm = c;5402IN_INTEGER_QUOTIENT /*----------------------------------------------*/5403*pm = INTEGER_N_LIMBS(a);5404n = INTEGER_N_LIMBS(b);5405a++;5406b++;5407c++;5408if(n == 1)5409{ /* Division by short number */5410i = *pm - 1;5411q = a[i] % *b; /* Carried residue */5412if((c[i] = a[i] / *b) == 0)5413--*pm;5414while(i)5415{5416lw = (uint)(q)*BASE_LIMB + (uint)a[--i];5417q = (LIMB)(lw % *b);5418c[i] = (LIMB)(lw / *b);5419}5420#if defined(D_CHECK_EXACTNESS_OF_DIVISION)5421nr = (q != 0);5422#endif5423}5424else5425{ /* Division by big number: n > 1 */5426BIGINT aw, bq;5427int k, j, shift, n1 = n + 1;5428LIMB carry, aj, aj1, aj2, bn1, bn2;5429INTEGER_STACK_NEW(bq, n1); /* For B[] * Q */5430COUNT_LEADING_ZERO_BITS_IN_LIMB(shift, b[n-1]);5431if(shift) /* Normalize to make b[n-1] >= BASE_LIMB/2 */5432{5433a[*pm] = BigNShiftLeft(a, *pm, shift);5434BigNShiftLeft(b, n, shift);5435}5436else5437a[*pm] = 0;5438bn1 = b[n-1];5439bn2 = b[n-2];5440j = *pm;5441k = *pm - n; /* Top digit char[M-N] may be zero */5442*pm = k + 1; /* Number of char[K] getting iterations */5443aw = a + k; /* Start of current subarray of A of the length N+1 */5444do5445{5446aj = a[j];5447aj1 = a[j-1];5448aj2 = a[j-2];5449lw = (uint)aj * BASE_LIMB + (uint)aj1;5450q = (aj == bn1) ? MAX_LIMB : (LIMB)(lw/bn1);5451lw -= (uint)q*(uint)bn1;5452if(lw < BASE_LIMB && (uint)bn2*(uint)q > lw*BASE_LIMB + (uint)aj2)5453{ /* Knuth's criterion shows Q is too big */5454q--;5455lw += (uint)bn1;5456if(lw < BASE_LIMB && (uint)bn2*(uint)q > lw*BASE_LIMB + (uint)aj2)5457q--; /* Q was still too big */5458}5459if(q)5460{ /* Multiply and subtract */5461i = carry = 0; /* Make copy of product B by Q */5462do5463{5464lw = (uint)q * (uint)b[i] + (uint)carry;5465carry = (LIMB)(lw/BASE_LIMB);5466bq[i] = (LIMB)lw;5467}while(++i < n);5468bq[i] = carry; /* BQ[] padded by zero if needs */5469i = n;5470do5471if(aw[i] != bq[i])5472{ /* AW[] - BQ[] != 0 */5473if(aw[i] < bq[i])5474{ /* AW[] - BQ[] is negative */5475q--; /* Additional correction */5476BigNMinusBigN(bq, n1, b, n);5477}5478break;5479}5480while(--i >= 0);5481#if defined(D_CHECK_EXACTNESS_OF_DIVISION)5482nr =5483#endif5484BigNMinusBigN(aw, n1, bq, n1); /* AW[] - BQ[] on spot */5485}5486c[k--] = q; /* Set current LIMB of quotient */5487--aw; /* Shift subarray of A digits */5488}while(--j >= n);5489if(c[*pm-1] == 0)5490--*pm; /* Real quotient size is M - N, not M - N + 1 */5491#if defined(D_CHECK_EXACTNESS_OF_DIVISION)5492if(nr && shift)5493nr = BigNShiftRight(a, nr, shift); /* Normalize remainder */5494#endif5495}5496if(set_minus)5497INTEGER_SET_MINUS(pm);5498OUT_INTEGER_QUOTIENT /*----------------------------------------------*/5499#if defined(D_CHECK_EXACTNESS_OF_DIVISION)5500if(nr)5501{5502PutFormattedU("\n***Division violation at Debug==%u\n", Debug);5503EXIT;5504}5505#endif5506}5507/*=IntegerSum=========================================5508Sum of two signed big numbers A and B, result in char5509*/5510void IntegerSum(BIGINT c, BIGINT a, BIGINT b)5511{5512int set_minus, i, na, nb;5513uint lw;5514LIMB carry;5515IN_INTEGER_SUM /*-----------------------------------*/5516if(INTEGER_N_LIMBS(a) < INTEGER_N_LIMBS(b))5517{5518BIGINT w = a; /* Swap input numbers if necessary */5519a = b;5520b = w;5521}5522na = INTEGER_N_LIMBS(a);5523nb = INTEGER_N_LIMBS(b);5524if(INTEGER_SIGN(a) == INTEGER_SIGN(b))5525{ /* The same signs: addition */5526set_minus = INTEGER_IS_NEGATIVE(a);5527i = 1;5528carry = 0;5529while(i <= nb) /* Common part */5530{5531lw = (uint)carry + (uint)a[i] + (uint)b[i];5532carry = (lw > MAX_LIMB);5533c[i++] = (LIMB)lw;5534}5535while(i <= na) /* Tail */5536{5537lw = (uint)carry + (uint)a[i];5538carry = (lw > MAX_LIMB);5539c[i++] = (LIMB)lw;5540}5541if(carry)5542c[i] = 1;5543else5544i--;5545#if defined(INTEGER_MAX_SIZE)5546if(i > IntegerMaxSize)5547IntegerMaxSize = i;5548#endif5549}5550else /* Different signs: subtraction */5551{5552if(na == nb)5553{5554i = na;5555while(i > 0)5556{5557if(a[i] < b[i])5558{ /* Swap numbers */5559BIGINT w = a;5560a = b;5561b = w;5562goto subtract;5563}5564else if(a[i] != b[i])5565goto subtract;5566i--;5567}5568c[0] = 0; /* Zero result of subtraction */5569goto out;5570}5571subtract:5572set_minus = INTEGER_IS_NEGATIVE(a);5573i = carry = 1;5574while(i <= nb) /* Common part */5575{5576lw = MAX_LIMB + (uint)carry + (uint)a[i] - (uint)b[i];5577carry = (lw > MAX_LIMB);5578c[i++] = (LIMB)lw;5579}5580while(i <= na)5581{5582lw = MAX_LIMB + (uint)carry + (uint)a[i];5583carry = (lw > MAX_LIMB);5584c[i++] = (LIMB)lw;5585}5586while(--i > 0)5587if(c[i] != 0)5588break;5589}5590c[0] = i; /* Number of LIMBs in sum */5591if(set_minus)5592INTEGER_SET_MINUS(c);5593out:5594OUT_INTEGER_SUM /*----------------------------------*/5595return;5596}55975598/*_6_6 Copy and delete functions===============================*/55995600/*=LieSumCopyInteger===================================5601Integer regime5602*/5603uint LieSumCopyInteger(uint a)5604{5605if(a != NIL)5606{5607uint ca, eca;5608BIGINT n, cn;5609int i;5610IN_LIE_SUM_COPY /*----------------------------------*/5611ca = eca = NodeLTNew();5612LIE_TERM_MONOMIAL(eca) = LIE_TERM_MONOMIAL(a);5613n = LIE_TERM_NUMERATOR_INTEGER(a);5614INTEGER_HEAP_COPY(cn, n, i);5615LIE_TERM_NUMERATOR_INTEGER(eca) = cn;5616if((n = LIE_TERM_DENOMINATOR_INTEGER(a)) != NULL)5617{5618INTEGER_HEAP_COPY(cn, n, i);5619LIE_TERM_DENOMINATOR_INTEGER(eca) = cn;5620}5621else5622LIE_TERM_DENOMINATOR_INTEGER(eca) = NULL;5623while((a = LIE_TERM_R(a)) != NIL)5624{5625LIE_TERM_R(eca) = NodeLTNew();5626eca = LIE_TERM_R(eca);5627LIE_TERM_MONOMIAL(eca) = LIE_TERM_MONOMIAL(a);5628n = LIE_TERM_NUMERATOR_INTEGER(a);5629INTEGER_HEAP_COPY(cn, n, i);5630LIE_TERM_NUMERATOR_INTEGER(eca)= cn;5631if((n = LIE_TERM_DENOMINATOR_INTEGER(a)) != NULL)5632{5633INTEGER_HEAP_COPY(cn, n, i);5634LIE_TERM_DENOMINATOR_INTEGER(eca) = cn;5635}5636else5637LIE_TERM_DENOMINATOR_INTEGER(eca) = NULL;5638}5639OUT_LIE_SUM_COPY /*----------------------------------*/5640return ca;5641}5642return NIL;5643}5644/*=LieSumCopyIntegerNegative===================================5645Copy changing sign. Integer regime5646*/5647uint LieSumCopyIntegerNegative(uint a)5648{5649if(a != NIL)5650{5651uint ca, eca;5652BIGINT n, cn;5653int i;5654ca = eca = NodeLTNew();5655LIE_TERM_MONOMIAL(eca) = LIE_TERM_MONOMIAL(a);5656n = LIE_TERM_NUMERATOR_INTEGER(a);5657INTEGER_HEAP_COPY(cn, n, i);5658INTEGER_MINUS(cn);5659LIE_TERM_NUMERATOR_INTEGER(eca) = cn;5660if((n = LIE_TERM_DENOMINATOR_INTEGER(a)) != NULL)5661{5662INTEGER_HEAP_COPY(cn, n, i);5663LIE_TERM_DENOMINATOR_INTEGER(eca) = cn;5664}5665else5666LIE_TERM_DENOMINATOR_INTEGER(eca) = NULL;5667while((a = LIE_TERM_R(a)) != NIL)5668{5669LIE_TERM_R(eca) = NodeLTNew();5670eca = LIE_TERM_R(eca);5671LIE_TERM_MONOMIAL(eca) = LIE_TERM_MONOMIAL(a);5672n = LIE_TERM_NUMERATOR_INTEGER(a);5673INTEGER_HEAP_COPY(cn, n, i);5674INTEGER_MINUS(cn);5675LIE_TERM_NUMERATOR_INTEGER(eca)= cn;5676if((n = LIE_TERM_DENOMINATOR_INTEGER(a)) != NULL)5677{5678INTEGER_HEAP_COPY(cn, n, i);5679LIE_TERM_DENOMINATOR_INTEGER(eca) = cn;5680}5681else5682LIE_TERM_DENOMINATOR_INTEGER(eca) = NULL;5683}5684return ca;5685}5686return NIL;5687}5688/*=LieSumCopyParametric=====================================5689Parametric regime5690*/5691uint LieSumCopyParametric(uint a)5692{5693if(a != NIL)5694{5695uint ca, eca;5696IN_LIE_SUM_COPY /*---------------------------------------*/5697ca = eca = NodeLTNew();5698LIE_TERM_MONOMIAL(eca) = LIE_TERM_MONOMIAL(a);5699LIE_TERM_NUMERATOR_SCALAR_SUM(eca) =57005701ScalarSumCopy(LIE_TERM_NUMERATOR_SCALAR_SUM(a));5702if(LIE_TERM_DENOMINATOR_SCALAR_SUM(a) != NIL)5703LIE_TERM_DENOMINATOR_SCALAR_SUM(eca) =5704ScalarSumCopy(LIE_TERM_DENOMINATOR_SCALAR_SUM(a));5705else5706LIE_TERM_DENOMINATOR_SCALAR_SUM(eca) = NIL;5707while((a = LIE_TERM_R(a)) != NIL)5708{5709LIE_TERM_R(eca) = NodeLTNew();5710eca = LIE_TERM_R(eca);5711LIE_TERM_MONOMIAL(eca) = LIE_TERM_MONOMIAL(a);5712LIE_TERM_NUMERATOR_SCALAR_SUM(eca) =5713ScalarSumCopy(LIE_TERM_NUMERATOR_SCALAR_SUM(a));5714if(LIE_TERM_DENOMINATOR_SCALAR_SUM(a) != NIL)5715LIE_TERM_DENOMINATOR_SCALAR_SUM(eca) =5716ScalarSumCopy(LIE_TERM_DENOMINATOR_SCALAR_SUM(a));5717else5718LIE_TERM_DENOMINATOR_SCALAR_SUM(eca) = NIL;5719}5720OUT_LIE_SUM_COPY /*---------------------------------------*/5721return ca;5722}5723return NIL;5724}5725/*=LieSumKillInteger================================5726a == NIL is admitted (Integer regime)5727*/5728void LieSumKillInteger(uint a)5729{5730uint b;5731BIGINT d;5732IN_LIE_SUM_KILL /*-------------------------------*/5733while(a != NIL)5734{5735b = a;5736a = LIE_TERM_R(a);5737INTEGER_KILL(LIE_TERM_NUMERATOR_INTEGER(b));5738if((d = LIE_TERM_DENOMINATOR_INTEGER(b)) != NULL)5739INTEGER_KILL(d);5740NODE_LT_KILL(b);5741}5742OUT_LIE_SUM_KILL /*------------------------------*/5743}5744/*=LieSumKillParametric===============================5745a == NIL is admitted (Parametric regime)5746*/5747void LieSumKillParametric(uint a)5748{5749uint b;5750IN_LIE_SUM_KILL /*---------------------------------*/5751while(a != NIL)5752{5753b = a;5754a = LIE_TERM_R(a);5755ScalarSumKill(LIE_TERM_NUMERATOR_SCALAR_SUM(b));5756ScalarSumKill(LIE_TERM_DENOMINATOR_SCALAR_SUM(b));5757NODE_LT_KILL(b);5758}5759OUT_LIE_SUM_KILL /*--------------------------------*/5760}5761/*=LieTermFromMonomialInteger============5762*/5763uint LieTermFromMonomialInteger(int mon)5764{5765BIGINT num;5766uint a = NodeLTNew();57675768LIE_TERM_MONOMIAL(a) = mon;57695770INTEGER_HEAP_NEW(num, 2);5771num[0] = num[1] = 1;5772LIE_TERM_NUMERATOR_INTEGER(a) = num;57735774LIE_TERM_DENOMINATOR_INTEGER(a) = NULL;57755776return a;5777}5778/*=LieTermFromMonomialParametric===========5779*/5780uint LieTermFromMonomialParametric(int mon)5781{5782BIGINT num;5783uint c = NodeSTNew(),5784a = NodeLTNew();57855786LIE_TERM_MONOMIAL(a) = mon;57875788SCALAR_TERM_MONOMIAL(c) = NIL;5789INTEGER_HEAP_NEW(num, 2);5790num[0] = num[1] = 1;5791SCALAR_TERM_NUMERATOR(c) = num;57925793LIE_TERM_NUMERATOR_SCALAR_SUM(a) = c;57945795LIE_TERM_DENOMINATOR_SCALAR_SUM(a) = NIL;57965797return a;5798}5799/*=ScalarSumCopy=======================================5800Caller ensures a != NIL5801*/5802uint ScalarSumCopy(uint a)5803{5804int i;5805BIGINT n, o;5806uint ca, bca, b, cb;5807bca = ca = NodeSTNew();5808while(YES)5809{5810/* Copy integer coefficient */58115812o = SCALAR_TERM_NUMERATOR(a);5813INTEGER_HEAP_COPY(n, o, i);5814SCALAR_TERM_NUMERATOR(ca) = n;58155816/* Copy scalar monomial */58175818b = SCALAR_TERM_MONOMIAL(a);5819if(b != NIL)5820{5821SCALAR_TERM_MONOMIAL(ca) = cb = NodeSFNew();5822SCALAR_FACTOR_WORD(cb) = SCALAR_FACTOR_WORD(b);5823while((b = SCALAR_FACTOR_R(b)) != NIL)5824{5825SCALAR_FACTOR_R(cb) = NodeSFNew();5826cb = SCALAR_FACTOR_R(cb);5827SCALAR_FACTOR_WORD(cb) = SCALAR_FACTOR_WORD(b);5828}5829}5830else5831SCALAR_TERM_MONOMIAL(ca) = NIL;58325833if((a = SCALAR_TERM_R(a)) == NIL)5834break;5835SCALAR_TERM_R(ca) = NodeSTNew();5836ca = SCALAR_TERM_R(ca);5837}5838return bca;5839}5840/*=ScalarSumKill=================================================5841Only at IsParametric == YES, a == NIL is admitted5842*/5843void ScalarSumKill(uint a)5844{5845uint b, c;5846while(a != NIL)5847{5848b = a;5849a = SCALAR_TERM_R(a);5850INTEGER_KILL(SCALAR_TERM_NUMERATOR(b));5851c = SCALAR_TERM_MONOMIAL(b);5852NODE_ST_KILL(b);5853while(c != NIL) /* Scalar monomial may be NIL, uint-uint mix */5854{5855b = c;5856c = SCALAR_FACTOR_R(c);5857NODE_SF_KILL(b);5858}5859}5860}5861/*=ScalarTermCopy====================================5862Caller ensures a != NIL5863*/5864uint ScalarTermCopy(uint a)5865{5866int i;5867BIGINT cn, n;5868uint m, ca = NodeSTNew();58695870/* Copy integer coefficient */58715872n = SCALAR_TERM_NUMERATOR(a);5873INTEGER_HEAP_COPY(cn, n, i);5874SCALAR_TERM_NUMERATOR(ca) = cn;58755876/* Copy monomial */58775878if((m = SCALAR_TERM_MONOMIAL(a)) != NIL)5879{5880uint cm;5881SCALAR_TERM_MONOMIAL(ca) = cm = NodeSFNew();5882SCALAR_FACTOR_WORD(cm) = SCALAR_FACTOR_WORD(m);5883while((m = SCALAR_FACTOR_R(m)) != NIL)5884{5885SCALAR_FACTOR_R(cm) = NodeSFNew();5886cm = SCALAR_FACTOR_R(cm);5887SCALAR_FACTOR_WORD(cm) = SCALAR_FACTOR_WORD(m);5888}5889}5890else5891SCALAR_TERM_MONOMIAL(ca) = NIL;5892return ca;5893}58945895/*_6_7 Technical functions=====================================*/58965897/*=Error!===============5898*/5899void Error(int i_message)5900{5901PutMessage(ERROR);5902PutMessage(i_message);5903EXIT;5904}5905/*=Initialization==========================================================5906*/5907void Initialization(void)5908{5909FILE *inif;5910short c;5911uint i, j;5912char * init_case[N_INIT_CASES];59135914/* Set cases strings */59155916init_case[COEFFICIENT_SUM_TABLE_SIZE] = "Co";5917init_case[CRUDE_TIME] = "Cr";5918init_case[ECHO_INPUT_FILE] = "Ec";5919init_case[EVEN_BASIS_SYMBOL] = "Ev";5920init_case[GAP_ALGEBRA_NAME] = "GAP a";5921init_case[GAP_BASIS_NAME] = "GAP b";5922init_case[GAP_RELATIONS_NAME] = "GAP r";5923init_case[GAP_OUTPUT_BASIS] = "GAP output b";5924init_case[GAP_OUTPUT_COMMUTATORS] = "GAP output c";5925init_case[GAP_OUTPUT_RELATIONS] = "GAP output r";5926init_case[GENERATOR_MAX_N] = "Ge";5927init_case[INPUT_DIRECTORY] = "Input d";5928init_case[INPUT_INTEGER_SIZE] = "Input i";5929init_case[INPUT_STRING_SIZE] = "Input s";5930init_case[LEFT_NORMED_OUTPUT] = "Le";5931init_case[LIE_MONOMIAL_SIZE] = "Lie";5932init_case[LINE_LENGTH] = "Lin";5933init_case[NAME_LENGTH] = "Na";5934init_case[NODE_LT_SIZE] = "Node L";5935init_case[NODE_SF_SIZE] = "Node scalar f";5936init_case[NODE_ST_SIZE] = "Node scalar t";5937init_case[ODD_BASIS_SYMBOL] = "Od";5938init_case[OUT_LINE_SIZE] = "Ou";5939init_case[PARAMETER_MAX_N] = "Pa";5940init_case[PUT_BASIS_ELEMENTS] = "Put b";5941init_case[PUT_COMMUTATORS] = "Put c";5942init_case[PUT_HILBERT_SERIES] = "Put H";5943init_case[PUT_INITIAL_RELATIONS] = "Put i";5944init_case[PUT_NON_ZERO_COEFFICIENTS] = "Put n";5945init_case[PUT_PROGRAM_HEADING] = "Put p";5946init_case[PUT_REDUCED_RELATIONS] = "Put r";5947init_case[PUT_STATISTICS] = "Put s";5948init_case[RELATION_SIZE] = "Re";59495950#if !defined(GAP)5951#if defined(SPP_2000)5952MessageFile = OpenFile("fplsa4.msg", "r");5953SessionFile = OpenFile("fplsa4.ses", "w");5954inif = OpenFile("fplsa4.ini", "r");5955#else5956MessageFile = OpenFile("fplsa4.msg", "rt");5957SessionFile = OpenFile("fplsa4.ses", "wt");5958inif = OpenFile("fplsa4.ini", "rt");5959#endif5960#else5961#if defined(SPP_2000)5962inif = OpenFile("fplsa4.ini", "r");5963#else5964inif = OpenFile("fplsa4.ini", "rt");5965#endif5966#endif5967while(YES)5968switch(ReadCaseFromFile(inif, init_case, N_INIT_CASES))5969{5970case COEFFICIENT_SUM_TABLE_SIZE:5971CoeffSumTableSize = ReadDecimalFromFile(inif);5972break;5973case CRUDE_TIME:5974CrudeTime = ReadBooleanFromFile(inif);5975break;5976case ECHO_INPUT_FILE:5977EchoInput = ReadBooleanFromFile(inif);5978break;5979case EVEN_BASIS_SYMBOL:5980BasisSymbolEven = fgetc(inif);5981break;5982case GAP_ALGEBRA_NAME:5983if(ReadStringFromFile(GAPAlgebraName, inif) == EOF)5984goto out;5985break;5986case GAP_BASIS_NAME:5987if(ReadStringFromFile(GAPBasisName, inif) == EOF)5988goto out;5989break;5990case GAP_RELATIONS_NAME:5991if(ReadStringFromFile(GAPRelationsName, inif) == EOF)5992goto out;5993break;5994case GAP_OUTPUT_BASIS:5995GAPOutputBasis = ReadBooleanFromFile(inif);5996break;5997case GAP_OUTPUT_COMMUTATORS:5998GAPOutputCommutators = ReadBooleanFromFile(inif);5999break;6000case GAP_OUTPUT_RELATIONS:6001GAPOutputRelations = ReadBooleanFromFile(inif);6002break;6003case GENERATOR_MAX_N:6004GeneratorMaxN = ReadDecimalFromFile(inif);6005break;6006case INPUT_DIRECTORY:6007while((c = fgetc(inif)) != '\n')6008{6009if(c == LEFT_COMMENT)6010{6011ungetc(c, inif);6012break;6013}6014if(c == EOF)6015goto out;6016if(!isspace(c))6017OutLine[PosOutLine++] = (char)c;6018}6019break;6020case INPUT_INTEGER_SIZE:6021InputIntegerSize = (uint)ReadDecimalFromFile(inif);6022InputIntegerSize++; /* For head */6023break;6024case INPUT_STRING_SIZE:6025InputStringSize = (uint)ReadDecimalFromFile(inif);6026break;6027case LEFT_NORMED_OUTPUT:6028if(ReadBooleanFromFile(inif))6029PutLieMonomial = PutLieMonomialLeftNormed;6030break;6031case LINE_LENGTH:6032LineLength = (uint)ReadDecimalFromFile(inif);6033break;6034case LIE_MONOMIAL_SIZE:6035LieMonomialSize = (int)ReadDecimalFromFile(inif);6036LieMonomial = (LIE_MON*)NewArray(LieMonomialSize, sizeof(LIE_MON),6037E_A_LIE_MONOMIAL);6038break;6039case NAME_LENGTH:6040NameLength1 = ReadDecimalFromFile(inif);6041NameLength1++;6042break;6043case NODE_LT_SIZE:6044NodeLTSize = (uint)ReadDecimalFromFile(inif);6045NodeLT =6046(NODE_LT*)NewArray(NodeLTSize, sizeof(NODE_LT), E_A_NODE_LT);6047i = 1;6048j = 2;6049while(j < NodeLTSize) /* Install NodeLT for Lie terms */6050LIE_TERM_R(i++) = j++;6051#if defined(SPACE_STATISTICS)6052LIE_TERM_R(i) = NOTHING;6053#else6054LIE_TERM_R(i) = NIL;6055#endif6056break;6057case NODE_SF_SIZE:6058NodeSFSize = (uint)ReadDecimalFromFile(inif);6059break;6060case NODE_ST_SIZE:6061NodeSTSize = (uint)ReadDecimalFromFile(inif);6062break;6063case ODD_BASIS_SYMBOL:6064BasisSymbolOdd = fgetc(inif);6065break;6066case OUT_LINE_SIZE:6067OutLineSize = (uint)ReadDecimalFromFile(inif);6068OutLine = (char *)NewArray(OutLineSize, 1, E_A_OUT_LINE);6069break;6070case PARAMETER_MAX_N:6071ParameterMaxN = ReadDecimalFromFile(inif);6072break;6073case PUT_BASIS_ELEMENTS:6074BasisElementsPut = ReadBooleanFromFile(inif);6075break;6076case PUT_COMMUTATORS:6077CommutatorsPut = ReadBooleanFromFile(inif);6078break;6079case PUT_HILBERT_SERIES:6080HilbertSeriesPut = ReadBooleanFromFile(inif);6081break;6082case PUT_INITIAL_RELATIONS:6083InitialRelationsPut = ReadBooleanFromFile(inif);6084break;6085case PUT_NON_ZERO_COEFFICIENTS:6086NonZeroCoefficientsPut = ReadBooleanFromFile(inif);6087break;6088case PUT_PROGRAM_HEADING:6089HeadingPut = ReadBooleanFromFile(inif);6090break;6091case PUT_REDUCED_RELATIONS:6092ReducedRelationsPut = ReadBooleanFromFile(inif);6093break;6094case PUT_STATISTICS:6095StatisticsPut = ReadBooleanFromFile(inif);6096break;6097case RELATION_SIZE:6098RelationSize = (int)ReadDecimalFromFile(inif);6099Relation = (REL*)NewArray(RelationSize, sizeof(REL), E_A_RELATION);6100break;6101case EOF:6102goto out;6103default:6104Error(E_WRONG_INI_CASE);6105}6106out:6107fclose(inif);6108if(HeadingPut)6109PutMessage(H_PROGRAM);6110GeneratorName = (char *)NewArray(GeneratorMaxN*NameLength1, sizeof(char),6111E_A_GENERATOR_NAME);6112}6113/*=NewArray=========================================6114*/6115void * NewArray(uint n, uint size, int i_message)6116{6117void * new_pointer = (void *)calloc(n, size);6118if(new_pointer == NULL && n != 0)6119{6120const char * format = "%u elements of size %u\n%u bytes\n";6121PutMessage(E_ALLOC);6122PutMessage(i_message);6123#if defined(ECHO_TO_SCREEN)6124printf(format, n, size, n*size);6125#endif6126#if !defined(GAP)6127fprintf(SessionFile, format, n, size, n*size);6128#endif6129EXIT;6130}6131return new_pointer;6132}6133/*=NodeLTNew===================6134Get node from NodeLT pool.6135*/6136uint NodeLTNew(void)6137{6138uint a = NodeLTTop;6139#if !defined(SPACE_STATISTICS)6140if(a == NIL)6141{6142TIME_OFF;6143PutStatistics();6144Error(E_NODE_LT_SIZE);6145}6146#endif6147NodeLTTop = LIE_TERM_R(a);6148#if defined(SPACE_STATISTICS)6149if(NodeLTTopMax < NodeLTTop)6150{6151if(NodeLTTop > NodeLTSize)6152{6153TIME_OFF;6154PutStatistics();6155Error(E_NODE_LT_SIZE);6156}6157NodeLTTopMax = NodeLTTop;6158}6159#endif6160LIE_TERM_R(a) = NIL;6161PP_CURRENT_N_LT /* MEMORY */6162return a;6163}6164/*=NodeSFNew=====================6165Get node from NodeSF pool.6166*/6167uint NodeSFNew(void)6168{6169uint a = NodeSFTop;6170#if !defined(SPACE_STATISTICS)6171if(a == NIL)6172{6173TIME_OFF;6174PutStatistics();6175Error(E_NODE_SF_SIZE);6176}6177#endif6178NodeSFTop = SCALAR_FACTOR_R(a);6179#if defined(SPACE_STATISTICS)6180if(NodeSFTopMax < NodeSFTop)6181{6182if(NodeSFTop > NodeSFSize)6183{6184TIME_OFF;6185PutStatistics();6186Error(E_NODE_SF_SIZE);6187}6188NodeSFTopMax = NodeSFTop;6189}6190#endif6191SCALAR_FACTOR_R(a) = NIL;6192PP_CURRENT_N_SF /* MEMORY */6193return a;6194}6195/*=NodeSTNew===================6196Get node from NodeST pool.6197*/6198uint NodeSTNew(void)6199{6200uint a = NodeSTTop;6201#if !defined(SPACE_STATISTICS)6202if(a == NIL)6203{6204TIME_OFF;6205PutStatistics();6206Error(E_NODE_ST_SIZE);6207}6208#endif6209NodeSTTop = SCALAR_TERM_R(a);6210#if defined(SPACE_STATISTICS)6211if(NodeSTTopMax < NodeSTTop)6212{6213if(NodeSTTop > NodeSTSize)6214{6215TIME_OFF;6216PutStatistics();6217Error(E_NODE_ST_SIZE);6218}6219NodeSTTopMax = NodeSTTop;6220}6221#endif6222SCALAR_TERM_R(a) = NIL;6223PP_CURRENT_N_ST /* MEMORY */6224return a;6225}6226/*=OpenFile================================6227*/6228FILE *OpenFile(char * file_name, char * file_type)6229{6230FILE *file = fopen(file_name, file_type);6231if(file == NULL)6232{6233printf("\nNo file: %s", file_name);6234exit(1);6235}6236return file;6237}62386239/*_6_8 Input functions=========================================*/62406241/*=BinaryQuestion!================6242*/6243int BinaryQuestion(int i_message)6244{6245char c[2];6246get_symbol:6247PutMessage(i_message);6248scanf("%1s", c);6249#if !defined(GAP)6250fputc(c[0], SessionFile);6251#endif6252switch(c[0])6253{6254case 'y': case 'Y': case '\n':6255return YES;6256case 'n': case 'N':6257break;6258case 'c': case 'C':6259Error(E_CANCEL_PROGRAM);6260default:6261goto get_symbol;6262}6263return NO;6264}6265/*=FindNameInTable==============================================6266Find name from string in table ...NameIn6267*/6268int FindNameInTable(char * name, char * nametab, int n_nametab)6269{6270char *w_nametab, *w_name;6271int j = 0;6272while(j < n_nametab)6273{6274w_nametab = nametab;6275w_name = name;6276while(YES)6277{6278if(*w_nametab == '\0') /* Table name ended */6279{6280if(!isalnum(*w_name) && *w_name != SUBSCRIPT_INPUT_SIGN)6281goto out; /* String name ended */6282break;6283} /* Table name goes on */6284if(!isalnum(*w_name) && *w_name != SUBSCRIPT_INPUT_SIGN)6285break; /* String name ended */6286if(*w_nametab != *w_name)6287break; /* Different names */6288w_nametab++;6289w_name++;6290}6291nametab += NameLength1;6292++j;6293}6294out:6295return j;6296}6297/*=GetGenerator===========================================================6298Read single generator from description string6299*/6300void GetGenerator(char * str)6301{6302char * name = GeneratorName + GeneratorN*NameLength1;6303if(GeneratorN == GeneratorMaxN)6304Error(E_GENERATOR_MAX_N);6305do6306{6307if(*str == ODD_GENERATOR_INPUT_SIGN) /* EVEN == 0 is assumed */6308LIE_MONOMIAL_PARITY(GeneratorN) = ODD;6309else6310*(name++) = *str;6311}while(*(str++));6312LIE_MONOMIAL_ORDER(GeneratorN) = /* Standard settings */6313LIE_MONOMIAL_POSITION(GeneratorN) = GeneratorN;6314LIE_MONOMIAL_RIGHT(GeneratorN) = 1; /* To avoid SQUARE interpretation */6315LIE_MONOMIAL_WEIGHT(GeneratorN++) = 1;6316}6317/*=GetInput===============================================================6318*/6319void GetInput(int n, char * fin)6320{6321char *instr, *sfname, *in_case[N_INPUT_CASES];6322FILE *inf;6323uint i, j;63246325in_case[GENERATORS] = "G";6326in_case[LIMITING_WEIGHT] = "L";6327in_case[PARAMETERS] = "P";6328in_case[RELATIONS] = "R";6329in_case[WEIGHTS] = "W";63306331sfname = OutLine + PosOutLine;6332if(n == 1) /* No input file at call */6333{6334PutMessage(H_ENTER_FILE);6335if (fgets(sfname, OutLineSize - PosOutLine, stdin)) {6336// success6337sfname[strcspn(sfname, "\n")] = '\0';6338}6339else {6340// failure6341exit(1);6342}6343}6344else6345do6346*(sfname++) = *fin;6347while(*(fin++));6348sfname = OutLine;6349while(YES)6350{6351if(*sfname == '.')6352break;6353if(*sfname == '\0')6354{6355*sfname = '.';6356*++sfname = 'i';6357*++sfname = 'n';6358*++sfname = '\0';6359break;6360}6361++sfname;6362}6363/*6364PutMessage(H_INPUT_FILE);6365PutStringStandard(OutLine);6366*/6367#if defined(SPP_2000)6368if((inf = fopen(OutLine, "r")) == NULL)6369#else6370if((inf = fopen(OutLine, "rt")) == NULL)6371#endif6372{ /* New file */6373if(!BinaryQuestion(H_CREATE_NEW_FILE))6374exit(1);6375instr = (char *)alloca(InputStringSize);6376if(instr == NULL)6377Error(E_A_STACK_INPUT_STRING);6378fgetc(stdin);6379#if defined(SPP_2000)6380inf = OpenFile(OutLine, "w");6381#else6382inf = OpenFile(OutLine, "wt");6383#endif6384KeyBoardStringToFile(H_ENTER_GENERATORS, "Generators: ", instr, inf);6385KeyBoardStringToFile(H_ENTER_WEIGHTS_IN_FILE, "Weights: ", instr, inf);6386KeyBoardStringToFile(H_ENTER_LIMITING_WEIGHT, "Limiting weight: ",6387instr, inf);6388KeyBoardStringToFile(H_ENTER_PARAMETERS, "Parameters: ", instr, inf);6389if(KeyBoardStringToFile(H_ENTER_RELATIONS, "Relations:\n", instr, inf))6390{6391n = YES; /* Now non-last */6392while(n && (j = KeyBoardBytesToString(instr)) > 0)6393for(i = 0; i <= j; i++)6394{6395if(instr[i] == '.')6396n = NO; /* Last input */6397fputc(instr[i], inf); /* Copy entered string to file */6398}6399}6400fclose(inf);6401#if defined(SPP_2000)6402inf = OpenFile(OutLine, "r");6403#else6404inf = OpenFile(OutLine, "rt");6405#endif6406}6407if(EchoInput)6408{6409short c;6410PutMessage(H_SHOW_INPUT);6411while((c = fgetc(inf)) != EOF)6412PutCharacter((char)c);6413}6414rewind(inf);6415while(YES)6416switch(ReadCaseFromFile(inf, in_case, N_INPUT_CASES))6417{6418case GENERATORS:6419ReadAndProcessStringsFromFile(GetGenerator, inf, ' ', ';');6420CUT_ARRAY(GeneratorName, char, NameLength1*GeneratorN);6421LieMonomialFreePosition = LieMonomialN = GeneratorN;6422#if defined(SPACE_STATISTICS)6423LieMonomialMaxN = LieMonomialN;6424#endif6425break;6426case WEIGHTS: /* Should come after generators */6427GeneratorMaxN = GeneratorN;6428GeneratorN = 0;6429ReadAndProcessStringsFromFile(GetWeight, inf, ' ', ';');6430GeneratorN = GeneratorMaxN;64316432/* Reorder generators in accordance with new weights */64336434i = 1;6435while(i < (uint)GeneratorN)6436{6437for(j = GeneratorN - 1; j >= i; j--)6438if(LIE_MONOMIAL_WEIGHT(j-1) > LIE_MONOMIAL_WEIGHT(j))6439{6440byte wt; /* To save swapped walues */64416442/* Swap generator names */64436444fin = GeneratorName + j*NameLength1; /* Next name */6445instr = fin - NameLength1; /* Previous name */6446for(n = 1; n < NameLength1; n++)6447{6448wt = *instr;6449*(instr++) = *fin;6450*(fin++) = wt;6451}64526453/* Swap weights */64546455wt = LIE_MONOMIAL_WEIGHT(j-1);6456LIE_MONOMIAL_WEIGHT(j-1) = LIE_MONOMIAL_WEIGHT(j);6457LIE_MONOMIAL_WEIGHT(j) = wt;64586459/* Swap parities */64606461wt = LIE_MONOMIAL_PARITY(j-1);6462LIE_MONOMIAL_PARITY(j-1) = LIE_MONOMIAL_PARITY(j);6463LIE_MONOMIAL_PARITY(j) = wt;6464}6465i++;6466}6467break;6468case PARAMETERS:6469IsParametric = YES;6470LieLikeTermsCollection = LieLikeTermsCollectionParametric;6471LieSumCopy = LieSumCopyParametric;6472LieSumKill = LieSumKillParametric;6473LieSumMinus = LieSumMinusParametric;6474NormalizeRelation = NormalizeRelationParametric;6475LieTermFromMonomial = LieTermFromMonomialParametric;6476PairMonomialMonomial = PairMonomialMonomialParametric;6477PairMonomialSum = PairMonomialSumParametric;6478PairSumMonomial = PairSumMonomialParametric;6479PairSumSum = PairSumSumParametric;6480SubstituteRelationInRelation =6481SubstituteRelationInRelationParametric;6482ParameterName = (char *)NewArray(NameLength1*ParameterMaxN,6483sizeof(char), E_A_PARAMETER_NAME);6484ParameterName[0] = 'i'; /* Obligatory imaginary unit */6485ParameterN = 1;6486ReadAndProcessStringsFromFile(GetParameter, inf, ' ', ';');6487CUT_ARRAY(ParameterName, char, NameLength1*ParameterN);6488NodeST =6489(NODE_ST*)NewArray(NodeSTSize, sizeof(NODE_ST), E_A_NODE_ST);6490i = 1;6491j = 2;6492while(j < NodeSTSize) /* Install NodeST for scalar terms */6493SCALAR_TERM_R(i++) = j++;6494#if defined(SPACE_STATISTICS)6495SCALAR_TERM_R(i) = NOTHING;6496#else6497SCALAR_TERM_R(i) = NIL;6498#endif6499NodeSF =6500(NODE_SF*)NewArray(NodeSFSize, sizeof(NODE_SF), E_A_NODE_SF);6501i = 1;6502j = 2;6503while(j < NodeSFSize) /* Install NodeSF for scalar factors */6504SCALAR_FACTOR_R(i++) = j++;6505#if defined(SPACE_STATISTICS)6506SCALAR_FACTOR_R(i) = NOTHING;6507#else6508SCALAR_FACTOR_R(i) = NIL;6509#endif6510break;6511case RELATIONS:6512if(!IsParametric)6513NonZeroCoefficientsPut = NO;6514ReadAndProcessStringsFromFile(GetRelation, inf, ';', '.');6515break;6516case LIMITING_WEIGHT:6517LimitingWeight = ReadDecimalFromFile(inf);6518while(fgetc(inf) != '\n') /* Go to next line */6519;6520break;6521case EOF:6522goto out;6523default:6524Error(E_WRONG_INPUT_CASE);6525}6526out:6527fclose(inf);6528if(LimitingWeight == 0)6529{6530PutMessage(H_ENTER_LIMITING_WEIGHT);6531scanf("%d", &LimitingWeight);6532#if !defined(GAP)6533fprintf(SessionFile, "%d\n", LimitingWeight);6534#endif6535}6536TIME_ON; /* First start of time */6537}6538/*=GetInteger===========================================================6539Read big integer with shift in string.6540A is already allocated array of LIMBs A[].6541*/6542void GetInteger(BIGINT a, char **pstr)6543{6544BIGINT w;6545int i;6546LIMB digit[2], ten[2];6547digit[0] = ten[0] = 1;6548ten[1] = 10;6549INTEGER_STACK_NEW(w, InputIntegerSize);6550while(**pstr == '0') /* Skip leading zeros */6551++*pstr;6552if(isdigit(**pstr))6553{ /* First digit: IntegerProduct does not eat 0 */6554a[0] = 1;6555a[1] = (LIMB)(**pstr - '0');6556while(isdigit(*(++*pstr)))6557{6558IntegerProduct(w, a, ten); /* a*10 */6559i = w[0];6560while(i >= 0) {6561a[i] = w[i]; /* Copy w to a */6562i--;6563}6564if((digit[1] = (LIMB)(**pstr - '0')) != 0)6565IntegerSum(a, a, digit); /* a*10 + digit */6566}6567}6568else /* Caller interprets absence as */6569a[0] = 0; /* 1 (or 0) depending on context */6570}6571/*=GetLieMonomial=========================================================6572Read monomial from string with transformations and substitutions6573*/6574uint GetLieMonomial(char **pstr)6575{6576int mon;6577uint a;6578IN_GET_LIE_MONOMIAL /*------------------------------------------------*/6579if(isalpha(**pstr))6580if((mon=FindNameInTable(*pstr,GeneratorName,GeneratorN)) < GeneratorN)6581{6582BIGINT num;6583SkipName(pstr);6584SkipSpaces(pstr);6585a = NodeLTNew();6586LIE_TERM_MONOMIAL(a) = mon;6587INTEGER_HEAP_NEW(num, 2); /* Make integer 1 */6588num[0] = num[1] = 1;6589if(IsParametric)6590{6591uint st = NodeSTNew();6592SCALAR_TERM_MONOMIAL(st) = NIL;6593SCALAR_TERM_NUMERATOR(st) = num;6594LIE_TERM_NUMERATOR_SCALAR_SUM(a) = st;6595LIE_TERM_DENOMINATOR_SCALAR_SUM(a) = NIL;6596}6597else6598{6599LIE_TERM_NUMERATOR_INTEGER(a) = num;6600LIE_TERM_DENOMINATOR_INTEGER(a) = NULL;6601}6602}6603else6604Error(E_UNDECLARED_GENERATOR);6605else if(**pstr == '[')6606{6607uint b;6608SkipSpaces(pstr);6609++*pstr;6610a = GetLieMonomial(pstr);6611SkipSpaces(pstr);6612if(**pstr != ',')6613Error(E_NO_COMMUTATOR_COMMA);6614++*pstr;6615SkipSpaces(pstr);6616b = GetLieMonomial(pstr);6617SkipSpaces(pstr);6618if(**pstr != ']')6619Error(E_NO_COMMUTATOR_BRACKET);6620++*pstr;6621a = (*PairSumSum)(a, b);6622}6623else6624Error(E_INVALID_CHARACTER);6625OUT_GET_LIE_MONOMIAL /*------------------------------------------------*/6626return a;6627}6628/*=GetLieSum=====================================================6629Read Lie expression from string and make internal representation6630*/6631uint GetLieSum(char **pstr)6632{6633uint lsum, term;6634int sign = PLUS;6635IN_GET_LIE_SUM /*--------------------------------------------*/6636SkipSpaces(pstr);6637if(**pstr == '-')6638{6639sign = MINUS;6640++*pstr;6641SkipSpaces(pstr);6642}6643lsum = GetLieTerm(pstr);6644SkipSpaces(pstr);6645if(sign)6646LieSumMinus(lsum);6647while(**pstr == '+' || **pstr == '-')6648{6649sign = (**pstr == '+') ? PLUS : MINUS;6650++*pstr;6651SkipSpaces(pstr);6652term = GetLieTerm(pstr);6653SkipSpaces(pstr);6654if(sign)6655LieSumMinus(term);6656lsum = LieSumAddition(lsum, term);6657}6658OUT_GET_LIE_SUM /*--------------------------------------------*/6659return lsum;6660}6661/*=GetLieTerm===========================================================6662*/6663uint GetLieTerm(char **pstr)6664{6665uint lterm;6666IN_GET_LIE_TERM /*--------------------------------------------------*/6667if(IsParametric)6668{6669uint num = GetScalarSum(pstr),6670den = NIL;6671SkipSpaces(pstr);6672if(**pstr == '/')6673{6674++*pstr;6675SkipSpaces(pstr);6676den = GetScalarSum(pstr);6677ScalarSumCancellation(&num, &den);6678}6679lterm = GetLieMonomial(pstr); /* May be sum with generic coeffs */6680if(den != NIL)6681LieSumDivScalarSum(lterm, den);6682LieSumMultScalarSum(lterm, num);6683}6684else6685{6686BIGINT num, den;6687INTEGER_STACK_NEW(num, InputIntegerSize);6688INTEGER_STACK_NEW(den, InputIntegerSize);6689GetInteger(num, pstr);6690den[0] = 0;6691SkipSpaces(pstr);6692if(**pstr == '/')6693{6694++*pstr;6695SkipSpaces(pstr);6696GetInteger(den, pstr);6697SkipSpaces(pstr);6698IntegerCancellation(num, den);6699}6700lterm = GetLieMonomial(pstr); /* May be sum with generic coeffs */6701if(den[0] != 0 && (den[0] != 1 || den[1] != 1))6702LieSumDivInteger(lterm, den);6703if(num[0] != 0 &&6704(INTEGER_IS_NEGATIVE(num) || num[0] != 1 || num[1] != 1))6705LieSumMultInteger(lterm, num);6706}6707OUT_GET_LIE_TERM /*--------------------------------------------------*/6708return lterm;6709}6710/*=GetUInteger====================================6711Read with shift long unsigned integer from string6712*/6713uint GetUInteger(char **pstr)6714{6715uint i = 0;6716while(isdigit(**pstr))6717{6718i = i*10 + **pstr - '0';6719++*pstr;6720}6721return i;6722}6723/*=GetParameter=====================================================6724Read single parameter from description string6725*/6726void GetParameter(char * str)6727{6728if(str[0] != 'i' || str[1] != '\0') /* Skip already settled `i' */6729{6730char * name = ParameterName + ParameterN*NameLength1;6731if(ParameterN == ParameterMaxN)6732Error(E_PARAMETER_MAX_N);6733do6734*(name++) = *str;6735while(*(str++));6736++ParameterN;6737}6738}6739/*=GetRelation===============================================================6740Read single relation from string, reduce and set in array6741*/6742void GetRelation(char * str)6743{6744uint a;6745if(str[0] != '\0' && (a = GetLieSum(&str)) != NIL)6746{6747int lmonpos, pos, i, l;6748if(RelationN + 1 == RelationSize)6749Error(E_RELATION_SIZE);6750(*NormalizeRelation)(a);6751lmonpos = LIE_TERM_MONOMIAL(a);6752i = LIE_MONOMIAL_ORDER(lmonpos);6753l = FindNewPositionInRelation(i);6754LIE_MONOMIAL_INDEX(lmonpos) = ~l; /* Set int of relation in LieMonomial */6755while(++i < LieMonomialN) /* Shift int's of relations in LieMonomial */6756if(LIE_MONOMIAL_IS_LEADING(pos = LIE_MONOMIAL_POSITION(i)))6757--LIE_MONOMIAL_INDEX(pos);6758for(i = RelationN; i > l; i--) /* Make room moving Relation structures */6759Relation[i] = Relation[i-1];6760RELATION_MIN_GENERATOR(l) = LIE_MONOMIAL_IS_GENERATOR(lmonpos) ?6761GeneratorN /* Don't differentiate leading generator */ : 0;6762RELATION_LIE_SUM(l) = a;6763RELATION_TO_BE_SUBSTITUTED(l) = (byte)(l < RelationN);6764++RelationN;6765#if defined(SPACE_STATISTICS)6766if(RelationN > MaxNRelation)6767MaxNRelation = RelationN;6768#endif6769ReduceRelations(l);6770}6771}6772/*=GetScalarSum==================================================6773*/6774uint GetScalarSum(char **pstr)6775{6776uint a, term;6777int is_par, is_negative;6778if(**pstr == '(')6779{6780is_par = YES;6781++*pstr;6782SkipSpaces(pstr);6783}6784else6785is_par = NO;6786if(**pstr == '-')6787{6788is_negative = YES;6789++*pstr;6790SkipSpaces(pstr);6791}6792else6793is_negative = NO;6794a = GetScalarTerm(pstr);6795if(is_negative)6796ScalarSumMinus(a);6797while(**pstr == '+' || **pstr == '-')6798{6799is_negative = (**pstr == '+') ? NO : YES;6800++*pstr;6801SkipSpaces(pstr);6802term = GetScalarTerm(pstr);6803if(is_negative)6804ScalarSumMinus(term);6805a = ScalarSumAddition(a, term);6806}6807if(is_par)6808{6809if(**pstr != ')')6810Error(E_NO_R_PARENTHESIS);6811++*pstr;6812SkipSpaces(pstr);6813}6814return a;6815}6816/*=GetScalarTerm=========================================================6817Read unsigned scalar term in Parametric regime6818*/6819uint GetScalarTerm(char **pstr)6820{6821BIGINT nums, numh;6822int i, change_sign;6823uint m, f, a = NodeSTNew();68246825/* Read numerical coefficient */68266827INTEGER_STACK_NEW(nums, InputIntegerSize);6828GetInteger(nums, pstr);6829if(nums[0] == 0)6830nums[0] = nums[1] = 1;6831INTEGER_HEAP_COPY(numh, nums, i);6832SCALAR_TERM_NUMERATOR(a) = numh;68336834/* Read scalar monomial */68356836SkipSpaces(pstr);6837m = NIL;6838while(isalpha(**pstr) &&6839(i=FindNameInTable(*pstr,ParameterName,ParameterN)) < ParameterN)6840{6841f = NodeSFNew();6842SCALAR_FACTOR_PARAMETER(f) = i;6843SkipName(pstr);6844SkipSpaces(pstr);6845if(**pstr == '^')6846{6847++*pstr;6848SkipSpaces(pstr);6849if(isdigit(**pstr))6850{6851i = (byte)GetUInteger(pstr); /* Read degree */6852SkipSpaces(pstr);6853}6854else6855Error(E_NO_GENERAL_POWER);6856}6857else6858i = 1;6859SCALAR_FACTOR_DEGREE(f) = i; /* Degree of parameter */6860m = ScalarMonomialMultiplication(&change_sign, m, f);6861if(change_sign)6862SCALAR_TERM_MINUS(a);6863}6864SCALAR_TERM_MONOMIAL(a) = m;6865return a;6866}6867/*=GetWeight=================================================6868Read single generator weight from description string6869*/6870void GetWeight(char *str)6871{6872if(GeneratorN == GeneratorMaxN)6873Error(E_TOO_MUCH_INPUT_WEIGHTS);6874if(!isdigit(*str))6875Error(E_NON_NUM_INPUT_WEIGHT);6876LIE_MONOMIAL_WEIGHT(GeneratorN++) = (byte)GetUInteger(&str);6877}6878/*=KeyBoardBytesToString=====================6879Returns 0 for the last input (ended by '.')6880and last position otherwise6881*/6882int KeyBoardBytesToString(char *str)6883{6884char c;6885int inspace = YES, i = -1;6886while((c = fgetc(stdin)) != '\n')6887{6888if(c == ' ')6889{6890if(inspace)6891continue; /* Skip extra blanks */6892inspace = YES;6893}6894else6895inspace = NO;6896PutCharacter(c);6897if(c == '\b')6898{6899#if !defined(GAP)6900fseek(SessionFile, -2, SEEK_CUR);6901#endif6902--i;6903}6904else6905{6906if(++i == InputStringSize)6907Error(E_INPUT_STRING_SIZE);6908str[i] = c;6909}6910}6911putchar('\n');6912if(i < 0) /* Empty input */6913return 0;6914if(str[i] != ';' && str[i] != '.')6915{ /* Add semicolon if necessary */6916if(++i == InputStringSize)6917Error(E_INPUT_STRING_SIZE);6918str[i] = ';';6919}6920if(++i == InputStringSize)6921Error(E_INPUT_STRING_SIZE);6922str[i] = '\n'; /* Go to new line */6923return i;6924}6925/*=KeyBoardStringToFile==========================================6926Add string from keyboard to file between prefix and semicolon6927*/6928int KeyBoardStringToFile(int i_m, char * prefix, char * str, FILE *file)6929{6930short itop;6931PutMessage(i_m);6932itop = KeyBoardBytesToString(str);6933if(itop)6934{6935short i = 0;6936while(*prefix)6937fputc(*(prefix++), file); /* Copy prefix to file */6938while(i <= itop)6939fputc(str[i++], file); /* Copy entered string to file */6940}6941return itop;6942}6943/*=ReadAndProcessStringsFromFile==================================6944Read array of strings separated by `sep' and ended with `end'6945Remove comments and unnecessary spaces6946Add ending '\0' Process strings by (*proc_func)6947*/6948void ReadAndProcessStringsFromFile(void (*proc_func)(char *str), FILE *inf,6949char sep, char end)6950{6951char *str, *wstr;6952char line_break;6953short c;6954int i;6955str = (char *)alloca(InputStringSize);6956if(str == NULL)6957Error(E_A_STACK_INPUT_STRING);6958line_break = (sep == ' ') ? '\n' : '\0';6959wstr = str;6960i = 0; /* Count number of characters */6961while(YES)6962{6963if((c = fgetc(inf)) == sep || c == line_break || c == end)6964{6965if(wstr != str && wstr[-1] == ' ')6966wstr[-1] = '\0'; /* Kill ending blank */6967else6968*wstr = '\0';6969(*proc_func)(str); /* Process string */6970if(c == end) /* Last */6971break;6972wstr = str; /* Intermediate */6973i = 0; /* Count number of characters */6974while((c = SkipSpacesInFile(inf)) == LEFT_COMMENT)6975SkipCommentInFile(inf);6976if(c == EOF)6977break;6978continue;6979}6980if(isspace(c))6981{6982if(++i == InputStringSize)6983Error(E_INPUT_STRING_SIZE);6984*(wstr++) = ' ';6985SkipSpacesInFile(inf);6986continue;6987}6988if(c == LEFT_COMMENT)6989{6990if(*wstr != ' ')6991{6992if(++i == InputStringSize)6993Error(E_INPUT_STRING_SIZE);6994*(wstr++) = ' ';6995}6996do6997SkipCommentInFile(inf);6998while((c = SkipSpacesInFile(inf)) == LEFT_COMMENT);6999continue;7000}7001if(++i == InputStringSize)7002Error(E_INPUT_STRING_SIZE);7003*(wstr++) = (char)c;7004}7005}7006/*=ReadBooleanFromFile========================7007Read boolean constant from file7008*/7009int ReadBooleanFromFile(FILE *file)7010{7011short c;7012int bool;7013c = fgetc(file);7014switch(c)7015{7016case 'Y': case 'y':7017bool = YES;7018break;7019case 'N': case 'n':7020bool = NO;7021break;7022}7023while(!isspace(c = fgetc(file)) && c != EOF)7024;7025ungetc(c, file);7026return bool;7027}7028/*=ReadCaseFromFile=====================================7029*/7030int ReadCaseFromFile(FILE * file, char * case_str[], int n_cases)7031{7032char file_str[CASE_STRING_SIZE];7033char * w_file_str, *w_case_str;7034int i_case;7035short c;7036while(SkipSpacesInFile(file) == LEFT_COMMENT)7037SkipCommentInFile(file);70387039/* Read string ending with : from file */70407041w_file_str = file_str;7042do7043{7044if((c = fgetc(file)) == EOF)7045{7046i_case = c;7047goto out;7048}7049if(c == ':')7050c = '\0';7051*(w_file_str++) = (char)c;7052}while(c);70537054/* Compare strings */70557056i_case = 0;7057do7058{7059w_file_str = file_str;7060w_case_str = case_str[i_case];7061do7062{7063if(*w_case_str == '\0') /* Case is found */7064{7065while(SkipSpacesInFile(file) == LEFT_COMMENT)7066SkipCommentInFile(file);7067goto out;7068}7069if(*w_file_str == '\0')7070break;7071}while(*(w_file_str++) == *(w_case_str++));7072}while(++i_case < n_cases);7073out:7074return i_case;7075}7076/*=ReadDecimalFromFile===================7077Read unsigned decimal integer from file7078*/7079uint ReadDecimalFromFile(FILE *file)7080{7081short c;7082uint i = 0;7083while(isdigit(c = fgetc(file)))7084i = i*10 + c - '0';7085ungetc(c, file);7086return i;7087}7088/*=ReadStringFromFile====================7089*/7090short ReadStringFromFile(char * str, FILE *file)7091{7092short c;7093while((c = fgetc(file)) != '\n')7094{7095if(c == LEFT_COMMENT)7096{7097ungetc(c, file);7098break;7099}7100if(c == EOF)7101break;7102if(!isspace(c))7103*(str++) = (char)c;7104}7105return c;7106}7107/*=SkipCommentInFile=======================7108*/7109short SkipCommentInFile(FILE *file)7110{7111short c;7112while((c = fgetc(file)) != RIGHT_COMMENT)7113if(c == EOF)7114Error(E_UNEXPECTED_EOF);7115return c;7116}7117/*=SkipName===============================================7118*/7119void SkipName(char **pstr)7120{7121while(isalnum(**pstr) || **pstr == SUBSCRIPT_INPUT_SIGN)7122++*pstr;7123}7124/*=SkipSpaces==================7125Skip spaces to right in string7126*/7127void SkipSpaces(char **pstr)7128{7129while(isspace(**pstr))7130++*pstr;7131}7132/*=SkipSpacesInFile==============7133Returns first non-space symbol7134*/7135short SkipSpacesInFile(FILE *file)7136{7137short c;7138while(isspace(c = fgetc(file)))7139;7140ungetc(c, file);7141return c;7142}71437144/*_6_9 Output functions========================================*/71457146/*=AddSymbolToOutLine!===============7147Add symbol to output line OutLine7148*/7149void AddSymbolToOutLine(char c, int position)7150{7151if(position >= OutLineSize)7152Error(E_OUT_LINE_SIZE);7153OutLine[position] = c;7154}7155/*=InLineLevel==============================7156Install level in OutLine7157*/7158void InLineLevel(int level)7159{7160if(level != CurrentLevel)7161{7162AddSymbolToOutLine(LEVEL, ++PosOutLine);7163AddSymbolToOutLine((char)level, ++PosOutLine);7164CurrentLevel = level;7165if(level > MaxLevel)7166MaxLevel = level;7167else if(level < MinLevel)7168MinLevel = level;7169}7170}7171/*=InLineNumberInBrackets=====7172(With 2 blanks after)7173*/7174void InLineNumberInBrackets(uint n)7175{7176InLineSymbol('(');7177InLineString(UToString(n));7178InLineString(") ");7179}7180/*=InLineString!================================7181Add string to output line OutLine for 2D output7182*/7183void InLineString(char *str)7184{7185while(*str)7186InLineSymbol(*(str++));7187}7188/*=InLineSubscript================7189Add symbolic subscript to OutLine7190*/7191void InLineSubscript(char * s)7192{7193int level = CurrentLevel;7194InLineLevel(level - 1);7195InLineString(s);7196InLineLevel(level);7197}7198/*=InLineSymbol!================================7199Add symbol to output line OutLine for 2D output7200*/7201void InLineSymbol(char symbol)7202{7203AddSymbolToOutLine(symbol, ++PosOutLine);7204++LastItemEnd;7205}7206/*=InLineTableName===================7207Possibly subscripted7208*/7209void InLineTableName(char * name)7210{7211PreviousEnd = LastItemEnd;7212while(*name)7213{7214if(*name == SUBSCRIPT_INPUT_SIGN)7215{7216InLineSubscript(++name);7217break;7218}7219InLineSymbol(*(name++));7220}7221}7222/*=UToString!====================================7223Transform unsigned long number to decimal string7224*/7225char * UToString(uint n)7226{ /*12345678910*/7227static char decimal_string[] = "4294967295";7228char * first = decimal_string + 10;7229do7230{7231*--first = '0' + (char)(n % 10);7232n /= 10;7233}while(n);7234return first;7235}7236/*=PutBasis============================================7237Set index numbers and print basis elements7238*/7239void PutBasis(void)7240{7241int pos, ord;7242uint i = 0; /* Index number of basis element */7243TIME_OFF;7244if(BasisElementsPut)7245PutMessage(H_BASIS_ELEMENTS);7246for(ord = 0; ord < LieMonomialN; ord++)7247{7248pos = LIE_MONOMIAL_POSITION(ord);7249if(LIE_MONOMIAL_IS_BASIS(pos))7250{7251LIE_MONOMIAL_INDEX(pos) = ++i;7252if(BasisElementsPut)7253{7254PutStart();7255InLineNumberInBrackets(i);7256PutLieBasisElement(pos);7257InLineString(" = ");7258IN_LINE_MARGIN;7259(*PutLieMonomial)(pos);7260PutEnd();7261}7262}7263}7264if(BasisElementsPut)7265if(IncompletedBasis)7266PutDots();7267TIME_ON;7268}7269#if defined(GAP)7270/*=PutBasisGAP=============================================72717272Print basis elements into GAP file7273*/7274void PutBasisGAP(void)7275{7276int pos, ord, tord;7277tord = LieMonomialN - 1;7278while(YES)7279{7280if(LIE_MONOMIAL_IS_BASIS(LIE_MONOMIAL_POSITION(tord)))7281break;7282--tord;7283}7284PutStringStandard(GAPBasisName);7285PutStringStandard(":=[\n");7286for(ord = 0; ord <= tord; ord++)7287{7288pos = LIE_MONOMIAL_POSITION(ord);7289if(LIE_MONOMIAL_IS_BASIS(pos))7290{7291PosOutLine = -1;7292PutLieMonomialGAP(pos);7293if(ord < tord)7294PutStringGAP(",\n");7295}7296}7297PutStringStandard("\n];\n");7298}7299#endif7300/*=PutBlock==============================================================7301Print block of 2D output7302*/7303void PutBlock(void)7304{7305if(!PrintEnd && (LastItemEnd > LineLength || PreviousEnd==LastItemEnd))7306PrintEnd = PreviousEnd;7307if(PrintEnd && CurrentLevel <= MAIN_LEVEL)7308{7309int xp, yp = MaxLevel, i, prlvl;7310while(yp >= MinLevel)7311{7312for(xp = 1; xp <= Margin; xp++)7313PutCharacter(' ');7314i = 0;7315while(xp <= PrintEnd)7316{7317switch(OutLine[i])7318{7319case LEVEL:7320prlvl = OutLine[++i];7321break;7322case MARGIN:7323NewMargin = xp - 1;7324break;7325default:7326PutCharacter((char)((prlvl == yp) ? OutLine[i] : ' '));7327xp++;7328}7329i++;7330}7331PutCharacter('\n');7332yp--;7333}7334PutCharacter('\n'); /* To next line */7335Margin = NewMargin;7336LastItemEnd += Margin - PrintEnd;7337PrintEnd = 0;7338xp = i;7339while(OutLine[xp] != LEVEL)7340xp--;7341*OutLine = LEVEL;7342OutLine[1] = MaxLevel = MinLevel = OutLine[xp + 1];7343xp = PosOutLine;7344PosOutLine = 1;7345while(i <= xp)7346{7347OutLine[++PosOutLine] = OutLine[i];7348if(OutLine[i++] == LEVEL)7349{7350if(OutLine[i] > MaxLevel)7351MaxLevel = OutLine[i];7352else if(OutLine[i] < MinLevel)7353MinLevel = OutLine[i];7354}7355}7356}7357}7358/*=PutCharacter!=============7359Echoed output of character7360*/7361void PutCharacter(char c)7362{7363#if defined(ECHO_TO_SCREEN)7364putc(c, stdout);7365#endif7366#if !defined(GAP)7367putc(c, SessionFile);7368#endif7369}7370#if defined(GAP)7371/*=PutCharacterGAP!======================7372Echoed output of character in GAP file7373*/7374void PutCharacterGAP(char c)7375{7376#if defined(ECHO_TO_SCREEN)7377putc(c, stdout);7378#endif7379/* putc(c, SessionFile); */7380PosOutLine++;7381}7382#endif7383/*=PutCoefficientTable=================================7384Print list of non-zero coefficients7385*/7386void PutCoefficientTable(void)7387{7388int i, j = 0, first = YES;7389TIME_OFF;7390if(CoeffParamTable != NULL)7391{7392for(i = FIRST_GENUINE_PARAMETER; i < ParameterN; i++)7393if(CoeffParamTable[i])7394{7395if(first)7396{7397first = NO;7398PutMessage(H_NON_ZERO_COEFFICIENTS);7399}7400PutStart();7401InLineNumberInBrackets(++j);7402InLineTableName(ParameterName + i*NameLength1);7403PutEnd();7404}7405free(CoeffParamTable);7406}7407if(CoeffSumTableN)7408{7409if(first)7410{7411first = NO;7412PutMessage(H_NON_ZERO_COEFFICIENTS);7413}7414for(i = 0; i < CoeffSumTableN; i++)7415{7416PutStart();7417InLineNumberInBrackets(++j);7418IN_LINE_MARGIN;7419PutScalarSum(CoeffSumTable[i]);7420PutEnd();7421ScalarSumKill(CoeffSumTable[i]);7422}7423free(CoeffSumTable);7424}7425TIME_ON;7426}7427/*=PutCommutators========================================================7428Compute and print commutators of basis elements [E(O),E(O)] -> E(O)7429*/7430void PutCommutators(void)7431{7432int oi, oj, i, j, was_comm, set_minus;7433uint icomm, rpart;7434byte weight_j;7435TIME_OFF;7436PutMessage(H_COMMUTATORS);7437TIME_ON;7438icomm = 0;7439for(oj = 0; oj < LieMonomialN; oj++)7440if(LIE_MONOMIAL_IS_BASIS(j = LIE_MONOMIAL_POSITION(oj)))7441{7442weight_j = LIE_MONOMIAL_WEIGHT(j);7443was_comm = NO;7444set_minus = LIE_MONOMIAL_IS_EVEN(j);7445for(oi = 0; oi < oj ||7446(oi == oj &&7447LIE_MONOMIAL_IS_ODD(LIE_MONOMIAL_POSITION(oi))); oi++)7448if(weight_j + LIE_MONOMIAL_WEIGHT(i = LIE_MONOMIAL_POSITION(oi))7449<= LimitingWeight)7450{7451if(LIE_MONOMIAL_IS_BASIS(i))7452{7453if(LIE_MONOMIAL_IS_EVEN(i))7454set_minus = YES;7455rpart = (*PairMonomialMonomial)(j, i);7456if(rpart != NIL)7457{7458TIME_OFF;7459PutStart();7460InLineNumberInBrackets(++icomm);7461InLineSymbol('[');7462PutLieBasisElement(i);7463InLineSymbol(',');7464PutLieBasisElement(j);7465InLineString("] = ");7466IN_LINE_MARGIN;7467if(set_minus)7468(*LieSumMinus)(rpart);7469PutLieSum(PutLieBasisElement, rpart);7470LieSumKill(rpart);7471PutEnd();7472TIME_ON;7473}7474}7475was_comm = YES;7476}7477else7478{7479if(was_comm)7480{7481while(oi < oj ||7482(oi == oj &&7483LIE_MONOMIAL_IS_ODD(LIE_MONOMIAL_POSITION(oi))))7484{7485if(LIE_MONOMIAL_IS_BASIS(LIE_MONOMIAL_POSITION(oi)))7486++icomm;7487++oi;7488}7489TIME_OFF;7490PutDots();7491PutCharacter('\n');7492TIME_ON;7493}7494break;7495}7496}7497TIME_OFF;7498if(icomm == 0)7499PutMessage(H_NO_PUT_COMMUTATORS);7500TIME_ON;7501}7502#if defined(GAP)7503/*=PutCommutatorsGAP==================================7504Transform commutators of ordinary finite-dimensional7505parameter-free Lie algebra to GAP input format7506*/7507void PutCommutatorsGAP(void)7508{7509int oi, oj, otop, i, j;7510BIGINT num;7511uint rpart, a, b;7512otop = LieMonomialN - 1;7513while(YES)7514{7515if(LIE_MONOMIAL_IS_BASIS(LIE_MONOMIAL_POSITION(otop)))7516break;7517--otop;7518}7519PutStringStandard(GAPAlgebraName);7520PutStringStandard(":=[\n"); /* Begin main list */7521for(oj = 0; oj <= otop; oj++)7522if(LIE_MONOMIAL_IS_BASIS(j = LIE_MONOMIAL_POSITION(oj)))7523{7524PosOutLine = -1;7525PutCharacterGAP('[');7526for(oi = 0; oi <= otop; oi++)7527if(LIE_MONOMIAL_IS_BASIS(i = LIE_MONOMIAL_POSITION(oi)))7528{7529if(oi == oj)7530rpart = NIL;7531else if(oj < oi)7532{7533rpart = (*PairMonomialMonomial)(i, j);7534(*LieSumMinus)(rpart);7535}7536else7537rpart = (*PairMonomialMonomial)(j, i);7538if(rpart != NIL)7539{7540a = LIE_TERM_R(rpart); /* Invert list */7541LIE_TERM_R(rpart) = NIL;7542while(a != NIL)7543{7544b = a;7545a = LIE_TERM_R(a);7546LIE_TERM_R(b) = rpart;7547rpart = b;7548}7549PutStringGAP("[[");7550a = rpart; /* Put indices of basis elements */7551while(YES)7552{7553PutStringGAP(UToString(LIE_MONOMIAL_INDEX(7554LIE_TERM_MONOMIAL(a))));7555a = LIE_TERM_R(a);7556if(a == NIL)7557break;7558PutStringGAP(",");7559}7560PutStringGAP("],[");7561a = rpart; /* Put coefficients */7562while(YES)7563{7564num = LIE_TERM_NUMERATOR_INTEGER(a);7565if(INTEGER_IS_NEGATIVE(num))7566PutStringGAP("-");7567PutIntegerUnsignedGAP(num);7568if((num = LIE_TERM_DENOMINATOR_INTEGER(a)) != NULL)7569{7570PutStringGAP("/");7571PutIntegerUnsignedGAP(num);7572}7573a = LIE_TERM_R(a);7574if(a == NIL)7575break;7576PutStringGAP(",");7577}7578PutStringGAP("]]");7579LieSumKill(rpart);7580}7581else7582PutStringGAP("[[],[]]");7583if(oi != otop)7584PutStringGAP(",");7585}7586PutStringGAP("],");7587PutCharacter('\n');7588}7589PutStringStandard("-1,0];\n"); /* Close main list */7590}7591#endif7592/*=PutDegree==========================7593Print (positive) integer degree7594*/7595void PutDegree(uint deg)7596{7597if(deg != 1)7598{7599int level = CurrentLevel;7600InLineLevel(level + 1);7601InLineString(UToString((uint)deg));7602InLineLevel(level);7603PutBlock();7604}7605}7606/*=PutDimensions================================================7607Print dimensions of homogeneous components of algebra7608*/7609void PutDimensions(void)7610{7611int next, ord, pos;7612byte curwt;7613uint dim;7614TIME_OFF;7615PutMessage(H_HILBERT_SERIES);7616PutStart();7617InLineString("H(t) = ");7618IN_LINE_MARGIN;7619next = NO;7620curwt = LIE_MONOMIAL_WEIGHT(0); /* Start initial weight */7621dim = 0;7622for(ord = 0; ord < LieMonomialN; ord++)7623{7624pos = LIE_MONOMIAL_POSITION(ord);7625if(LIE_MONOMIAL_IS_BASIS(pos))7626{7627if(LIE_MONOMIAL_WEIGHT(pos) != curwt)7628{7629if(dim != 0)7630{7631if(next)7632PutString(" + ");7633else7634next = YES;7635if(dim != 1)7636{7637PutString(UToString(dim));7638PutSymbol(' ');7639}7640PutSymbol('t');7641PutDegree(curwt);7642}7643curwt = LIE_MONOMIAL_WEIGHT(pos); /* Start new weight */7644dim = 1;7645}7646else7647dim++;7648}7649}7650if(dim != 0) /* Print last element */7651{7652if(next)7653PutString(" + ");7654if(dim != 1)7655{7656PutString(UToString(dim));7657PutSymbol(' ');7658}7659PutSymbol('t');7660PutDegree(curwt);7661}7662if(IncompletedBasis)7663PutString(" + ...");7664PutEnd();7665TIME_ON;7666}7667/*=PutDots=============================7668Print vertical dots7669*/7670void PutDots(void)7671{7672#if defined(ECHO_TO_SCREEN)7673printf(" .\n .\n .\n");7674#endif7675#if !defined(GAP)7676fprintf(SessionFile, " .\n .\n .\n");7677#endif7678}7679/*=PutEnd=======================7680Print last block of 2D output7681*/7682void PutEnd(void)7683{7684PreviousEnd = LastItemEnd;7685PutBlock();7686}7687/*=PutFormattedU=================7688*/7689void PutFormattedU(char * format, uint i)7690{7691#if !defined(GAP)7692fprintf(SessionFile, format, i);7693#endif7694#if defined(ECHO_TO_SCREEN)7695printf(format, i);7696#endif7697}7698/*=PutLieBareTerm=======================================================7699put_lie_mon == PutLieMonomial -> in terms of generators,7700put_lie_mon == PutLieBasisElement -> in terms of basis elements.7701*/7702void PutLieBareTerm(void (*put_lie_mon)(int a), uint a)7703{7704int put_mult_sign = NO;7705if(IsParametric)7706{7707uint num = LIE_TERM_NUMERATOR_SCALAR_SUM(a),7708den = LIE_TERM_DENOMINATOR_SCALAR_SUM(a);77097710/* Put numerator */77117712if(SCALAR_TERM_R(num) != NIL)7713{7714PutSymbol('('); /* Sum */7715PutScalarSum(num);7716PutSymbol(')');7717put_mult_sign = YES;7718}7719else if(SCALAR_TERM_MONOMIAL(num) != NIL ||7720INTEGER_IS_NOT_UNIT_ABS(SCALAR_TERM_NUMERATOR(num)) ||7721den != NIL)7722{7723PutScalarBareTerm(num); /* Single term */7724put_mult_sign = YES;7725}77267727/* Put denominator */77287729if(den != NIL)7730{7731PutSymbol('/');7732if(SCALAR_TERM_R(den) == NIL && /* Single factor denominator */7733(SCALAR_TERM_MONOMIAL(den) == NIL ||7734INTEGER_IS_UNIT(SCALAR_TERM_NUMERATOR(den))))7735PutScalarBareTerm(den);7736else7737{7738PutSymbol('('); /* Sum or multifactor denominator */7739PutScalarSum(den);7740PutSymbol(')');7741}7742}7743}7744else7745{7746BIGINT num = LIE_TERM_NUMERATOR_INTEGER(a),7747den = LIE_TERM_DENOMINATOR_INTEGER(a);7748if(den != NULL)7749{7750PutIntegerUnsigned(num);7751PutSymbol('/');7752PutIntegerUnsigned(den);7753put_mult_sign = YES;7754}7755else if(INTEGER_IS_NOT_UNIT_ABS(num))7756{7757PutIntegerUnsigned(num);7758put_mult_sign = YES;7759}7760}7761if(put_mult_sign)7762PutSymbol(' ');7763(*put_lie_mon)(LIE_TERM_MONOMIAL(a));7764}7765/*=PutLieBasisElement======================================7766Print Lie monomial in form E or O7767i i7768*/7769void PutLieBasisElement(int pos)7770{7771InLineSymbol((char)(LIE_MONOMIAL_PARITY(pos) ? BasisSymbolOdd :7772BasisSymbolEven));7773InLineSubscript(UToString(LIE_MONOMIAL_INDEX(pos)));7774PutBlock();7775}7776/*=PutLieMonomialLeftNormed=====================================7777Put Lie monomial in terms of generators in left normed notation7778*/7779void PutLieMonomialLeftNormed(int pos)7780{7781if(LIE_MONOMIAL_IS_GENERATOR(pos))7782{7783InLineTableName(GeneratorName + pos*NameLength1);7784PutBlock();7785}7786else7787{7788uint deg = 1;7789int posi = LIE_MONOMIAL_LEFT(pos),7790posj = LIE_MONOMIAL_RIGHT(pos), posw;7791while(LIE_MONOMIAL_IS_COMMUTATOR(posi))7792{7793pos = posi;7794posi = LIE_MONOMIAL_LEFT(pos);7795posw = LIE_MONOMIAL_RIGHT(pos);7796if(posj == posw)7797++deg;7798else7799{7800posi = pos;7801break;7802}7803}7804if(posi != posj)7805PutLieMonomialLeftNormed(posi);7806else7807++deg;7808if(LIE_MONOMIAL_IS_COMMUTATOR(posj))7809{7810PutSymbol('(');7811PutLieMonomialLeftNormed(posj);7812PutSymbol(')');7813}7814else7815PutLieMonomialLeftNormed(posj);7816PutDegree(deg);7817}7818}7819/*=PutLieMonomialStandard============================================7820Put Lie monomial in terms of generators in standard bracket notation7821*/7822void PutLieMonomialStandard(int pos)7823{7824if(LIE_MONOMIAL_IS_GENERATOR(pos))7825{7826InLineTableName(GeneratorName + pos*NameLength1);7827PutBlock();7828}7829else7830{7831PutSymbol('[');7832PutLieMonomialStandard(LIE_MONOMIAL_LEFT(pos));7833PutSymbol(',');7834PutLieMonomialStandard(LIE_MONOMIAL_RIGHT(pos));7835PutSymbol(']');7836}7837}7838#if defined(GAP)7839/*=PutLieMonomialGAP============================================7840Put Lie monomial in terms of generators in GAP bracket notation7841*/7842void PutLieMonomialGAP(int pos)7843{7844if(LIE_MONOMIAL_IS_GENERATOR(pos))7845PutStringGAP(UToString(pos+1));7846else7847{7848PutStringGAP("[");7849PutLieMonomialGAP(LIE_MONOMIAL_LEFT(pos));7850PutStringGAP(",");7851PutLieMonomialGAP(LIE_MONOMIAL_RIGHT(pos));7852PutStringGAP("]");7853}7854}7855#endif7856/*=PutLieSum=====================================================7857put_lie_mon == PutLieMonomial -> in terms of generators,7858put_lie_mon == PutLieBasisElement -> in terms of basis elements7859*/7860void PutLieSum(void (*put_lie_mon)(int a), uint a)7861{7862if(a == NIL)7863PutSymbol('0');7864else7865{7866uint na;7867int is_negative = NO;7868if(IsParametric)7869{7870na = LIE_TERM_NUMERATOR_SCALAR_SUM(a);7871if(SCALAR_TERM_R(na) == NIL &&7872INTEGER_IS_NEGATIVE(SCALAR_TERM_NUMERATOR(na)))7873is_negative = YES;7874}7875else if(INTEGER_IS_NEGATIVE(LIE_TERM_NUMERATOR_INTEGER(a)))7876is_negative = YES;7877if(is_negative)7878PutString("- ");7879PutLieBareTerm(put_lie_mon, a);7880while((a = LIE_TERM_R(a)) != NIL)7881{7882is_negative = NO;7883if(IsParametric)7884{7885na = LIE_TERM_NUMERATOR_SCALAR_SUM(a);7886if(SCALAR_TERM_R(na) == NIL &&7887INTEGER_IS_NEGATIVE(SCALAR_TERM_NUMERATOR(na)))7888is_negative = YES;7889}7890else if(INTEGER_IS_NEGATIVE(LIE_TERM_NUMERATOR_INTEGER(a)))7891is_negative = YES;7892PutString(is_negative ? " - " : " + ");7893PutLieBareTerm(put_lie_mon, a);7894}7895}7896}7897/*=PutMessage!=========================7898Put message from MessageFile7899*/7900void PutMessage(int i_message)7901{7902#if !defined(GAP)7903static int current_message;7904short c;7905if(i_message < current_message)7906{7907rewind(MessageFile);7908current_message = 0;7909}7910while(i_message > current_message)7911if((c = fgetc(MessageFile)) == '$')7912++current_message;7913else if(c == EOF)7914Error(E_MESSAGE);7915while(YES)7916{7917c = fgetc(MessageFile);7918if(c == '$')7919{7920++current_message;7921return;7922}7923PutCharacter((char)c);7924}7925#endif7926}7927/*=PutRelations====================================7928Print list of relations7929*/7930void PutRelations(int msg)7931{7932int i;7933TIME_OFF;7934PutMessage(msg);7935for(i = 0; i < RelationN; i++)7936{7937PutStart();7938InLineNumberInBrackets(i+1);7939IN_LINE_MARGIN;7940PutLieSum(PutLieMonomial, RELATION_LIE_SUM(i));7941InLineString(" = 0");7942PutEnd();7943}7944if(IncompletedRelations)7945PutDots();7946TIME_ON;7947}7948#if defined(GAP)7949/*=PutRelationsGAP=================================7950*/7951void PutRelationsGAP(void)7952{7953int i;7954BIGINT num;7955uint a;7956PutStringStandard(GAPRelationsName);7957PutStringStandard(":=[\n"); /* Begin main list */7958for(i = 0; i < RelationN; i++)7959{7960a = RELATION_LIE_SUM(i);7961PosOutLine = -1;7962PutCharacterGAP('[');7963while(YES)7964{7965PutLieMonomialGAP(LIE_TERM_MONOMIAL(a));7966PutStringGAP(",");7967num = LIE_TERM_NUMERATOR_INTEGER(a);7968if(INTEGER_IS_NEGATIVE(num))7969PutStringGAP("-");7970PutIntegerUnsignedGAP(num);7971if((a = LIE_TERM_R(a)) == NIL)7972break;7973PutStringGAP(",");7974}7975PutStringGAP("]");7976if(i < RelationN - 1)7977PutStringGAP(",\n");7978}7979PutStringStandard("\n];\n");7980}7981#endif7982/*=PutScalarBareTerm===============================7983Intermediate print of unsigned scalar term7984*/7985void PutScalarBareTerm(uint a)7986{7987int put_1 = YES;79887989/* Put integer coefficient */79907991if(INTEGER_IS_NOT_UNIT_ABS(SCALAR_TERM_NUMERATOR(a)))7992{7993PutIntegerUnsigned(SCALAR_TERM_NUMERATOR(a));7994put_1 = NO;7995}79967997/* Put scalar monomial */79987999a = SCALAR_TERM_MONOMIAL(a); /* uint - uint mixing */8000if(a != NIL)8001{8002if(put_1 == NO)8003PutSymbol(' ');8004PutScalarFactor(a);8005while((a = SCALAR_FACTOR_R(a)) != NIL)8006{8007PutSymbol(' ');8008PutScalarFactor(a);8009}8010put_1 = NO;8011}8012if(put_1) /* Nothing's been put before */8013PutSymbol('1');8014}8015/*=PutIntegerUnsigned=============================================8016Print big integer8017*/8018void PutIntegerUnsigned(BIGINT bn)8019{8020BIGINT bnw;8021char * decstr;8022LIMB res;8023uint lw;8024int i,8025n = INTEGER_N_LIMBS(bn),8026nw = n;8027bnw = (BIGINT)alloca(sizeof(LIMB)*n);8028if(bnw == NULL)8029Error(E_A_STACK_INTEGER); /* LIMB contains 5 decimal digits */8030decstr = (char *)alloca(5*n + 1);8031if(decstr == NULL)8032Error(E_A_STACK_INTEGER_DECIMAL_STRING);8033decstr += 5*n; /* Go to last byte of string */8034for(i = 0; i < n; i++)8035bnw[i] = *(++bn); /* Copy body of big number */8036/* Transform big number array to decimal string */8037decstr[0] = '\0';8038do8039{ /* Divide bnw number in array form by 10 on spot */8040res = 0;8041if(nw) /* Otherwise 0/n -> 0, 0 mod n -> 0 */8042{8043i = nw;8044do8045{8046lw = (uint)res * BASE_LIMB + (uint)bnw[--i];8047res = (LIMB)(lw % 10);8048bnw[i] = (LIMB)(lw / 10);8049}while(i);8050if(bnw[nw-1] == 0)8051--nw;8052}8053*--decstr = '0' + res;8054}while(nw);8055if(n < 3)8056PutString(decstr); /* Don't cut short number */8057else8058do8059PutSymbol(*decstr);8060while(*(++decstr));8061}8062#if defined(GAP)8063/*=PutIntegerUnsignedGAP==========================================8064Print big integer into GAP file8065*/8066void PutIntegerUnsignedGAP(BIGINT bn)8067{8068BIGINT bnw;8069char * decstr;8070LIMB res;8071uint lw;8072int i,8073n = INTEGER_N_LIMBS(bn),8074nw = n;8075bnw = (BIGINT)alloca(sizeof(LIMB)*n);8076if(bnw == NULL)8077Error(E_A_STACK_INTEGER); /* LIMB contains 5 decimal digits */8078decstr = (char *)alloca(5*n + 1);8079if(decstr == NULL)8080Error(E_A_STACK_INTEGER_DECIMAL_STRING);8081decstr += 5*n; /* Go to last byte of string */8082for(i = 0; i < n; i++)8083bnw[i] = *(++bn); /* Copy body of big number */8084/* Transform big number array to decimal string */8085decstr[0] = '\0';8086do8087{ /* Divide bnw number in array form by 10 on spot */8088res = 0;8089if(nw) /* Otherwise 0/n -> 0, 0 mod n -> 0 */8090{8091i = nw;8092do8093{8094lw = (uint)res * BASE_LIMB + (uint)bnw[--i];8095res = (LIMB)(lw % 10);8096bnw[i] = (LIMB)(lw / 10);8097}while(i);8098if(bnw[nw-1] == 0)8099--nw;8100}8101*--decstr = '0' + res;8102}while(nw);8103PutStringGAP(decstr);8104}8105#endif8106/*=PutScalarFactor========================================================8107*/8108void PutScalarFactor(uint a)8109{8110InLineTableName(ParameterName + SCALAR_FACTOR_PARAMETER(a)*NameLength1);8111PutDegree(SCALAR_FACTOR_DEGREE(a));8112PutBlock();8113}8114/*=PutScalarSum================================================8115Intermediate print of scalar sum8116*/8117void PutScalarSum(uint a)8118{8119if(a == NIL)8120PutSymbol('0');8121else8122{8123if(INTEGER_IS_NEGATIVE(SCALAR_TERM_NUMERATOR(a)))8124PutString("- ");8125PutScalarBareTerm(a);8126while((a = SCALAR_TERM_R(a)) != NIL)8127{8128PutString(INTEGER_IS_NEGATIVE(SCALAR_TERM_NUMERATOR(a)) ?8129" - " : " + ");8130PutScalarBareTerm(a);8131}8132}8133}8134/*=PutStart=======================================8135Start of 2dimensional output8136*/8137void PutStart(void)8138{8139LastItemEnd = PrintEnd = Margin = NewMargin = 0;8140AddSymbolToOutLine(LEVEL, 0);8141AddSymbolToOutLine(MAIN_LEVEL, 1);8142PosOutLine = 1;8143CurrentLevel = MaxLevel = MinLevel = MAIN_LEVEL;8144}8145/*=PutStatistics=========================================================8146Put time and space statistics8147*/8148void PutStatistics(void)8149{8150uint sec, min, sec_100,8151time = CrudeTime ? TimeC : (uint)(((double)TimeC/CLOCKS_PER_SEC)*100);8152PutStringStandard("Time: ");8153if(!CrudeTime)8154{8155sec_100 = (uint)(time%100);8156time /= 100; /* In seconds */8157}8158sec = (uint)(time%60);8159time /= 60; /* In minutes */8160min = (uint)(time%60);8161time /= 60; /* In hours */8162if(time)8163PutFormattedU("%lu h ", time);8164if(min || time)8165PutFormattedU("%lu min ", min);8166if(CrudeTime)8167PutFormattedU("%lu sec\n", sec);8168else8169{8170PutFormattedU("%u.", sec);8171if(sec_100 < 10 && sec_100 > 0)8172PutCharacter('0');8173PutFormattedU("%u sec\n", sec_100);8174}8175#if defined(SPACE_STATISTICS)8176PutFormattedU("Number of relations: %14u\n", MaxNRelation);8177PutFormattedU("Number of Lie monomials: %10u\n", LieMonomialMaxN);8178PutFormattedU("Number of Lie terms: %14u\n", NodeLTTopMax - 1);8179if(IsParametric)8180{8181PutFormattedU("Number of scalar terms: %11u\n", NodeSTTopMax - 1);8182min = NodeSFTopMax - 1;8183PutFormattedU("Number of scalar factors:%10u\n", NodeSFTopMax - 1);8184}8185#endif8186#if defined(INTEGER_MAX_SIZE)8187PutFormattedU("Maximum size of integer in limbs: %5u\n", IntegerMaxSize);8188#endif8189#if defined(DEBUG)8190PutDebugU("Current Debug", Debug);8191#endif8192}8193/*=PutString================81942D output of string8195*/8196void PutString(char *str)8197{8198PreviousEnd = LastItemEnd;8199InLineString(str);8200PutBlock();8201}8202#if defined(GAP)8203/*=PutStringGAP=======================================================8204Put string in GAP file8205*/8206void PutStringGAP(char *str)8207{8208char c = 0;8209while(*str)8210{8211if(PosOutLine == GAP_WIDTH - 2)8212if(isdigit(*str)) /* Going to write in last position */8213{8214if(isdigit(c))8215PutCharacter('\\');8216PutCharacter('\n');8217PosOutLine = -1;8218PutCharacterGAP(*str); /* Continue number in the next line */8219}8220else8221{8222PutCharacter(*str);8223PutCharacter('\n');8224PosOutLine = -1; /* Ready to next line */8225}8226else8227PutCharacterGAP(*str);8228c = *str; /* Remember previous symbol */8229++str;8230}8231}8232#endif8233/*=PutStringStandard==============8234*/8235void PutStringStandard(char *str)8236{8237#if defined(ECHO_TO_SCREEN)8238printf("%s", str);8239#endif8240#if !defined(GAP)8241fprintf(SessionFile, "%s", str);8242#endif8243}8244/*=PutSymbol!===============82452D output of symbol8246*/8247void PutSymbol(char c)8248{8249PreviousEnd = LastItemEnd;8250InLineSymbol(c);8251PutBlock();8252}82538254/*_6_10 Debugging functions===========================================*/82558256#if defined(DEBUG)8257/*=PutDebugHeader====================================================8258Put header of tracing output.8259*/8260void PutDebugHeader(uint debug, char * f_name, char * in_out)8261{8262#if !defined(GAP)8263fprintf(SessionFile,"\nDebug==%lu %s %s:\n", debug, f_name, in_out);8264#endif8265printf("\nDebug==%lu %s %s\n", debug, f_name, in_out);8266}8267/*=PutDebugInteger==============8268Put name and signed big number8269*/8270void PutDebugInteger(char * name, BIGINT u)8271{8272PutStart();8273InLineString(name);8274InLineString(": ");8275if(INTEGER_IS_NEGATIVE(u))8276InLineSymbol('-');8277PutIntegerUnsigned(u);8278PutEnd();8279}8280/*=PutDebugLieMonomial==============8281Put name and Lie sum8282*/8283void PutDebugLieMonomial(char * name, int a)8284{8285PutStart();8286InLineString(name);8287InLineString(": ");8288(*PutLieMonomial)(a);8289PutEnd();8290}8291/*=PutDebugLieMonomialTable=================================================8292*/8293void PutDebugLieMonomialTable(int newmon)8294{8295int i, j, count;8296PutStringStandard("LieMonomial table:\n"8297"ORDER POSITION LEFT RIGHT PARITY INDEX WEIGHT MONOMIAL\n");8298i = count = 0;8299while(count < LieMonomialN)8300{8301if(LIE_MONOMIAL_IS_OCCUPIED(i))8302{8303count++;8304PutFormattedU("%5d", LIE_MONOMIAL_ORDER(i));8305PutFormattedU("%9d", i);8306PutFormattedU("%5d", LIE_MONOMIAL_LEFT(i));8307PutFormattedU("%6d", LIE_MONOMIAL_RIGHT(i));8308PutFormattedU("%4d", LIE_MONOMIAL_PARITY(i));8309if((j = LIE_MONOMIAL_INDEX(i)) < 0 )8310PutFormattedU(" Rel %-4d", ~j);8311else8312PutFormattedU("%9d", j);8313PutFormattedU("%7d ", LIE_MONOMIAL_WEIGHT(i));8314if(newmon == i)8315PutStringStandard("New! ");8316PutStart();8317(*PutLieMonomial)(i);8318PutEnd();8319}8320else8321PutFormattedU("Position %d is free!\n", i);8322i++;8323}8324}8325/*=PutDebugLieSum==============8326Put name and Lie sum8327*/8328void PutDebugLieSum(char * name, uint a)8329{8330PutStart();8331InLineString(name);8332InLineString(": ");8333PutLieSum(PutLieMonomial, a);8334PutEnd();8335}8336/*=PutDebugLieTerm==============================================8337*/8338void PutDebugLieTerm(char * name, uint a)8339{8340PutStart();8341InLineString(name);8342InLineString(": ");8343if(a == NIL)8344PutSymbol('0');8345else8346{8347uint na;8348int is_negative = NO;8349if(IsParametric)8350{8351na = LIE_TERM_NUMERATOR_SCALAR_SUM(a);8352if(SCALAR_TERM_R(na) == NIL &&8353INTEGER_IS_NEGATIVE(SCALAR_TERM_NUMERATOR(na)))8354is_negative = YES;8355}8356else if(INTEGER_IS_NEGATIVE(LIE_TERM_NUMERATOR_INTEGER(a)))8357is_negative = YES;8358if(is_negative)8359PutString("- ");8360PutLieBareTerm(PutLieMonomial, a);8361}8362PutEnd();8363}8364/*=PutDebugU==================================8365Put name and long unsigned integer8366*/8367void PutDebugU(char * name, uint i)8368{8369printf("\n%s==%lu\n", name, i);8370#if !defined(GAP)8371fprintf(SessionFile, "\n%s==%lu\n", name, i);8372#endif8373}8374#if defined(D_PUT_RELATIONS)8375/*=PutDebugRelations======================================8376Put Relation table with structure fields8377*/8378void PutDebugRelations(void)8379{8380int i, mg;8381char * hformat = "Relations:\n N SUB MIN_GEN RELATION\n";8382printf(hformat);8383#if !defined(GAP)8384fprintf(SessionFile, hformat);8385#endif8386for(i = 0; i < RelationN; i++)8387{8388PutStart();8389if(i < 100)8390InLineSymbol(' ');8391if(i < 10)8392InLineSymbol(' ');8393InLineString(UToString(i));8394InLineString(" ");8395InLineSymbol(RELATION_TO_BE_SUBSTITUTED(i) ? 'Y':'N');8396InLineString(" ");8397mg = RELATION_MIN_GENERATOR(i);8398if(mg < 10)8399InLineSymbol(' ');8400InLineString(UToString(mg));8401InLineString(" ");8402if(mg < GeneratorN)8403PutLieMonomial(mg);8404else8405InLineString("done");8406InLineString(" ");8407IN_LINE_MARGIN;8408PutLieSum(PutLieMonomial, RELATION_LIE_SUM(i));8409PutEnd();8410}8411}8412#endif8413/*=PutDebugScalarSum==============8414Put name and scalar sum8415*/8416void PutDebugScalarSum(char * name, uint a)8417{8418PutStart();8419InLineString(name);8420InLineString(": ");8421PutScalarSum(a);8422PutEnd();8423}8424/*=PutDebugString================================8425Put name of string and string8426*/8427void PutDebugString(char * strname, char * str)8428{8429printf("%s: %s\n", strname, str);8430#if !defined(GAP)8431fprintf(SessionFile, "%s: %s\n", strname, str);8432#endif8433}8434#endif8435#if defined(MEMORY)8436/*=AddLieSumNs=============================================84378438*/8439void AddLieSumNs(uint a, int minus_or_plus,8440int *pn_lt, int *pn_int, int *pn_st, int *pn_sf)8441{8442int dn_lt, dn_int;8443dn_lt = dn_int = 0;8444while(a != NIL)8445{8446++dn_lt;8447if(IsParametric)8448{8449AddScalarSumNs(LIE_TERM_NUMERATOR_SCALAR_SUM(a),8450minus_or_plus, pn_int, pn_st, pn_sf);8451AddScalarSumNs(LIE_TERM_DENOMINATOR_SCALAR_SUM(a),8452minus_or_plus, pn_int, pn_st, pn_sf);8453}8454else8455{8456++dn_int; /* Numerator is obligatory */8457if(LIE_TERM_DENOMINATOR_INTEGER(a) != NULL)8458++dn_int;8459}8460a = LIE_TERM_R(a);8461}8462if(minus_or_plus == PLUS)8463{8464*pn_lt += dn_lt;8465if(!IsParametric)8466*pn_int += dn_int;8467}8468else8469{8470*pn_lt -= dn_lt;8471if(!IsParametric)8472*pn_int -= dn_int;8473}8474}8475/*=AddScalarSumNs========================================================8476*/8477void AddScalarSumNs(uint a, int minus_or_plus, int *pn_int, int *pn_st, int *pn_sf)8478{8479int dn_int, dn_st, dn_sf;8480uint b;8481dn_int = dn_st = dn_sf = 0;8482while(a != NIL)8483{8484++dn_st;8485++dn_int; /* Numerator is obligatory */8486b = SCALAR_TERM_MONOMIAL(a);8487while(b != NIL)8488{8489++dn_sf;8490b = SCALAR_FACTOR_R(b);8491}8492a = SCALAR_TERM_R(a);8493}8494if(minus_or_plus == PLUS)8495{8496*pn_int += dn_int;8497*pn_st += dn_st;8498*pn_sf += dn_sf;8499}8500else8501{8502*pn_int -= dn_int;8503*pn_st -= dn_st;8504*pn_sf -= dn_sf;8505}8506}8507/*=PutIntegerBalance===================================================8508*/8509void PutIntegerBalance(char * fname, int dn)8510{8511PutStringStandard("\nHeap integer balance violation in function:\n");8512PutStringStandard(fname);8513if(dn > 0)8514{8515#if !defined(GAP)8516fprintf(SessionFile, "\n*** %ld INTs gone to garbage\n", dn);8517#endif8518printf("\n*** %ld INTs gone to garbage\n", dn);8519}8520else8521{8522dn = -dn;8523#if !defined(GAP)8524fprintf(SessionFile, "\n*** %ld INTs appeared from nothing\n", dn);8525#endif8526printf("\n*** %ld INTs appeared from nothing\n", dn);8527}8528}8529/*=PutNodeBalance======================================================8530*/8531void PutNodeBalance(char * type, char * fname, int dn)8532{8533PutStringStandard(type);8534PutStringStandard(" balance violation in function:\n");8535PutStringStandard(fname);8536if(dn > 0)8537{8538#if !defined(GAP)8539fprintf(SessionFile, "\n*** %ld nodes gone to garbage\n", dn);8540#endif8541printf("\n*** %ld nodes gone to garbage", dn);8542}8543else8544{8545dn = -dn;8546#if !defined(GAP)8547fprintf(SessionFile, "\n***%ld nodes appeared from nothing\n", dn);8548#endif8549printf("\n***%ld nodes appeared from nothing\n", dn);8550}8551}8552#endif8553855485558556