Artifact 214be8515d465b76679f312506545ee6485a0ed133a14ca1e916491b8cfb1fff:
- File
r30/solve.red
— part of check-in
[eb17ceb7f6]
at
2020-04-21 19:40:01
on branch master
— Add Reduce 3.0 to the historical section of the archive, and some more
files relating to version sof PSL from the early 1980s. Thanks are due to
Paul McJones and Nelson Beebe for these, as well as to all the original
authors.git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/historical@5328 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 35048) [annotate] [blame] [check-ins using] [more...]
COMMENT SOLVE MODULE; %******************* Global Declarations ***************************; SYMBOLIC; FLAG('(!*SOLVEWRITE), 'SHARE); ARRAY !!CF(12), !!INTERVAL(10,2), !!EXACT(10); GLOBAL '(!!HIPOW !!GCD !*SOLVESINGULAR SM!* MP!* !*ALLBRANCH !*SOLVEWRITE !!ARBINT !*SOLVEINTERVAL !!INTERVALARRAY); !*SOLVESINGULAR := T; % Solves consistent, singular eqns (0=0) if T; !*ALLBRANCH := T; % Returns all branches of solutions if T; !*SOLVEWRITE := T; % Prints solutions if T; %!*SOLVEINTERVAL = NIL;% Attempts to isolate insoluble, real roots if T; !!INTERVALARRAY := '!!INTERVAL; % Value is the name of an array used to % pass args to RealRoot routines; !!ARBINT := 0; % Index for arbitrary constants; % !!HIPOW : SOLVECOEFF returns highest power of its arg in this % !!GCD : SOLVECOEFF returns GCD of powers of its arg in this % !!CF : Array of coeffs from SOLVECOEFF % % SM!* : List of solutions % MP!* : List of multiplicities; ALGEBRAIC MATRIX SOLN, MULTIPLICITY; %******************* Utility Functions *****************************; SYMBOLIC PROCEDURE RPLACX U; BEGIN SCALAR CARU; CARU := CAR U; RETURN RPLACD(RPLACA(U,CDR U),CARU) END; SYMBOLIC PROCEDURE UNIVARIATEP F; % F is a standard form. Non-nil iff F is univariate or a constant; DOMAINP F OR (DOMAINP LC F AND (DOMAINP RED F OR ((MVAR F = MVAR RED F) AND UNIVARIATEP RED F) )); SYMBOLIC SMACRO PROCEDURE SUBTRSQ(U,V); ADDSQ(U, NEGSQ V); SYMBOLIC SMACRO PROCEDURE VARLIS U; %U is an r-polynomial. %value is an ordered list of variables in U; VARLIS1(U,NIL); SYMBOLIC SMACRO PROCEDURE LFCTR U; COMMENT RETURNS LEFTFACTOR OF A PAIR. USED BY SUMFACTORS IN IEQN.RED; CAAR U; SYMBOLIC OPERATOR LCMD; SYMBOLIC PROCEDURE LCMD(C,D); COMMENT C and D are prefix rational numbers. Returns integer least-common-multiple of their denominators; LCM(DENR SIMP!* C, DENR SIMP!* D); SYMBOLIC PROCEDURE VARLIS1(U,V); IF DOMAINP U THEN V ELSE VARLIS1(CDR U,VARLIS1(CDAR U,ORDAS(CAAAR U,V))); SYMBOLIC PROCEDURE ORDAS(A,L); IF NULL L THEN LIST A ELSE IF A=CAR L THEN L ELSE IF ORDP(A,CAR L) THEN A . L ELSE CAR L . ORDAS(A,CDR L); SYMBOLIC PROCEDURE RATNUMP X; COMMENT Returns T iff any prefix expression X is a rational number; ATOM NUMR(X:=SIMP!* X) AND ATOM DENR X; FLAG ('(RATNUMP), 'DIRECT); SYMBOLIC PROCEDURE KARGLIS(KNAME, KLIS); COMMENT KNAME evaluates to an atom and KLIS to a list of kernels. Returns the list of kernels named KNAME in KLIS; IF NULL KLIS THEN NIL ELSE UNION(KARG1(KNAME, CAR KLIS), KARGLIS(KNAME,CDR KLIS)); SYMBOLIC PROCEDURE KARG1(KNAME, KRN); COMMENT KNAME evaluates to an atom and KRN to a kernel. Returns a list of kernels named KNAME in KRN; IF ATOM KRN THEN NIL ELSE IF CAR KRN=KNAME THEN UNION(KARGLIS(KNAME,CDR KRN), LIST(KRN)) ELSE KARGLIS(KNAME, CDR KRN); SYMBOLIC PROCEDURE ALLKERN ELST; COMMENT Returns list of all top-level kernels in the list of standard forms ELST; IF NULL ELST THEN NIL ELSE UNION(VARLIS CAR ELST, ALLKERN CDR ELST); SYMBOLIC OPERATOR FREEOFKERN; SYMBOLIC PROCEDURE FREEOFKERN(X,U); COMMENT Returns T iff any expression U is free of kernel X; IF ATOM X THEN FREEOF(U,X) ELSE FREEOF(SUBST('!!DUM,X,U),'!!DUM); FLAG('(FREEOFKERN),'DIRECT); SYMBOLIC PROCEDURE TOPKERN(EX, X); BEGIN COMMENT Returns list of toplevel kernels in the standard form EX that contain the kernel X; SCALAR ALLK, WITHX; ALLK := VARLIS EX; WHILE ALLK DO<< IF NOT FREEOFKERN(X,CAR ALLK) THEN WITHX:=CAR ALLK.WITHX; ALLK:=CDR ALLK>>; RETURN WITHX END; SYMBOLIC PROCEDURE COEFLIS(EX); % EX is a standard form. % Returns a list of the coefficients of the main variable % in ex in the form ((expon.coeff) (expon.coeff) ... ), % where the expon's occur in increasing order, and entries % do not occur of zero coefficients; BEGIN SCALAR X, ANS, OLDKORD, VAR; X := EX; IF DOMAINP(X) THEN RETURN (0 . X); VAR := MVAR(EX); WHILE (NOT DOMAINP(X)) AND MVAR(X)=VAR DO << ANS := (LDEG(X) . LC(X)) . ANS; X := RED(X) >>; IF X THEN ANS := (0 . X) . ANS; RETURN ANS END; %******************* Temporary Factoring Routine *******************; % The following square free factoring routine, based on the Reduce % function SQFRF, will eventually be replaced by the Norman-Moore % complete factorization technique.; FLUID '(!*GCD); SYMBOLIC PROCEDURE FACTLIS(EX, KLIST); % EX is a standard form. % KLIST is a list of kernels. % Returns a list of square free factors containing the elements of % KLIST in the form ((integer exponent . standard form factor) ...).; % Factors constant with respect to KLIST are discarded; BEGIN SCALAR FIRST, ANS, OLDGCD, OLDKORD; INTEGER EXPON; OLDGCD := !*GCD; !*GCD := T; % Must be on for SQFRF; OLDKORD := SETKORDER(KLIST); EX := REORDER(EX); WHILE (NOT DOMAINP(EX)) AND (MVAR(EX) MEMBER KLIST) DO << FIRST := PP(EX); IF NOT DOMAINP(FIRST) THEN << % Non-zero roots; EX := QUOTF(EX, FIRST); FIRST := SQFRF(FIRST); FOR EACH X IN FIRST DO IF NOT DOMAINP X THEN ANS := RPLACX X . ANS >> ELSE << % Zero root (possibly multiple); ANS := (LDEG(EX) . !*K2F(MVAR(EX))) . ANS; EX := QUOTF(EX, !*P2F(LPOW(EX))) >> >>; % Restore the state of the world; SETKORDER(OLDKORD); !*GCD := OLDGCD; RETURN ANS END; %******************* SOLVE Statement ******************************; SYMBOLIC PROCEDURE SIMPSOLVE ARGLIST; BEGIN INTEGER NARGS; NARGS := LENGTH(ARGLIST); RETURN !*F2Q IF NARGS=1 THEN SOLVE0(CAR ARGLIST,NIL) ELSE IF NARGS=2 THEN SOLVE0(CAR ARGLIST, CADR ARGLIST) ELSE SOLVE0(CAR ARGLIST,'LST . CDR ARGLIST) END; PUT ('SOLVE,'SIMPFN,'SIMPSOLVE); %******************* Fundamental SOLVE Procedures ******************; SYMBOLIC PROCEDURE SOLVE0(ELST, XLST); BEGIN COMMENT ELST is any prefix expression, including the kernel named LST with any number of arguments. XLST is a kernel, perhaps named LST with any number of arguments. Solves eqns in ELST for vars in XLST, putting solutions and multiplicities in SOLN and MULTIPLICITIES. Prints SOLN if !*SOLVEWRITE is non-nil. Returns number of rows in global matrix SOLN; SCALAR FLST, VARS, NONLIN; INTEGER NEQN, I; %/ MAYBELOADMATR(); ALGEBRAIC CLEAR SOLN, MULTIPLICITY; SM!* := MP!* := NIL; IF NOT ATOM ELST AND CAR ELST='LST THEN ELST:=CDR ELST ELSE ELST:=LIST ELST; NEQN:=0; WHILE ELST DO <<FLST:= NUMR SIMP!* CAR ELST . FLST; NEQN:=NEQN+1; ELST:= CDR ELST >>; % Note that ELST and XLST are reversed from the order entered; IF NULL XLST THEN <<VARS := ALLKERN FLST; WRITE "UNKNOWNS:"; MAPCAR(REVERSE VARS, FUNCTION MATHPRINT); TERPRI()>> ELSE<<IF ATOM XLST OR NOT(CAR XLST='LST)THEN XLST:=LIST(XLST) ELSE XLST:=CDR XLST; WHILE XLST DO<< VARS:=MVAR !*A2F CAR XLST.VARS; XLST:= CDR XLST>>>>; IF NOT(NEQN=LENGTH VARS) THEN REDERR "SOLVE CALLED WITH UNEQUAL NUMBER OF EXPRESSIONS AND UNKNOWNS"; IF NEQN=1 THEN IF NULL (FLST:=CAR FLST) THEN IF !*SOLVESINGULAR THEN <<!!ARBINT:=!!ARBINT+1; CONSSMMP(SIMP!* LIST('ARBCOMPLEX,!!ARBINT), 1) >> ELSE RETURN 0 ELSE <<VARS:=CAR VARS; SOLVE1(FLST./1, VARS, 1) >> COMMENT More than one equation; ELSE << SM!* := TP1(SOLVESYS(FLST, VARS)); MP!* := LIST(LIST(MK!*SQ(!*F2Q(1)))) >>; SM!* := MAPC2(SM!*, FUNCTION MK!*SQ); PUT('MULTIPLICITY, 'MATRIX, 'MAT . MP!*); PUT('SOLN, 'MATRIX, 'MAT . SM!*); IF !*SOLVEWRITE THEN MATPRI(SM!*, 'SOLN); RETURN LENGTH SM!* END; SYMBOLIC PROCEDURE CONSSMMP(S, M); BEGIN COMMENT S is a standard quotient and M is an integer. Conses (S) to global variable SM!* and (M) to global variable MP!*; SM!* := LIST(S) . SM!*; MP!* := LIST(MK!*SQ(M./1)) . MP!* END; SYMBOLIC PROCEDURE SOLVEF(F, V); % F is a standard form, V is a kernel. Returns a list of % pairs, each of which car is a standard quotient and cdr an % integer. If the integer is positive, the SQ is a zero of % F with multiplicity equal to the integer. Otherwise it is % an insoluble factor, with multiplicity the absolute value of % the integer; BEGIN SCALAR OLDSOLVEWRITE, ANS; OLDSOLVEWRITE := !*SOLVEWRITE; !*SOLVEWRITE := NIL; SOLVE0(MK!*SQ(!*F2Q(F)), V); ANS := PAIR(MAPCAR(SM!*, FUNCTION LAMBDA(X); SIMP!*(CAR(X))), MAPCAR(MP!*, FUNCTION CAR) ); !*SOLVEWRITE := OLDSOLVEWRITE; RETURN ANS END; %******************* Procedures for solving a single eqn ***********; SYMBOLIC PROCEDURE SOLVE1 (EX, X, MUL); BEGIN COMMENT Factors standard quotient EX with respect to toplevel occurrences of X and kernels containing variable X. Factors containing more than one such kernel are appended to SM!*, with a negative multiplicity indicating unsolvability, and SOLVE2 is applied to the other factors. Integer MUL is the multiplicity passed from any previous factorizations. Returns NIL; SCALAR E1, X1, TKLIST; INTEGER MU; EX := NUMR EX; TKLIST := TOPKERN(EX,X); IF NULL TKLIST THEN RETURN NIL; EX := FACTLIS(EX, TKLIST); WHILE EX DO << E1 := CDAR(EX); X1 := TOPKERN(E1, X); MU := MUL*CAAR EX; IF X1 THEN IF NULL CDR X1 THEN SOLVE2(E1,CAR X1,X,MU) ELSE IF SMEMQ('SOL, (X1:=SIMP!* LIST('SOL,MK!*SQ(E1./1), X))) THEN CONSSMMP(E1./1, -MU) ELSE SOLVE1(X1,X,MU); EX := CDR(EX) >> END; SYMBOLIC PROCEDURE SOLVE2(E1, X1, X, MU); BEGIN COMMENT E1 is a standard form, MU is an integer, X1 and X are kernels. Uses roots of unity, known inverses, together with quadratic, cubic and quartic formulas, treating other cases as unsolvable. Returns NIL; SCALAR B, C, D, F; INTEGER N; F:= ERRORSET(SOLVECOEFF(E1, X1),NIL,NIL); N:= !!GCD; COMMENT Test for single power of X1; IF ATOM(F) THEN CONSSMMP(E1./1, -MU) ELSE IF (F:=CAR F)=-1 THEN << B:= LIST('EXPT, MK!*SQ QUOTSQ(NEGSQ SIMP!* GETELV(LIST('!!CF,0)), SIMP!* GETELV(LIST('!!CF,1))), MK!*SQ(1 ./!!GCD)); FOR K := 0:N-1 DO << SETELV(LIST('!!CF,1), SIMP!* LIST('TIMES,B, MKEXP LIST('QUOTIENT,LIST('TIMES,K,2,'PI),N))); COMMENT x = b; IF X1=X THEN CONSSMMP(GETELV(LIST('!!CF, 1)), MU) COMMENT LOG(x) = b; ELSE IF CAR X1 = 'LOG THEN SOLVE1 (SUBTRSQ(SIMP!* CADR X1,SIMP!* LIST('EXPT,'E,MK!*SQ GETELV(LIST('!!CF, 1)))),X,MU) ELSE IF CAR X1 = 'EXPT THEN COMMENT c**(...) = b; IF FREEOF(CADR X1,X) THEN << IF !*ALLBRANCH THEN <<!!ARBINT:=!!ARBINT+1; C:=LIST('TIMES,2,'I,'PI,LIST('ARBINT,!!ARBINT)) >> ELSE C:=0; SOLVE1(SUBTRSQ(SIMP!* CADDR X1,QUOTSQ(ADDSQ( SIMP!* LIST('LOG,MK!*SQ GETELV(LIST('!!CF, 1))),SIMP!* C), SIMP!* LIST('LOG,CADR X1))),X,MU) >> ELSE IF FREEOF(CADDR X1,X) THEN COMMENT (...)**(m/n) = b; IF RATNUMP CADDR X1 THEN SOLVE1(SUBTRSQ( EXPTSQ(SIMP!* CADR X1,NUMR SIMP!* CADDR X1), SIMP!* LIST('EXPT,MK!*SQ GETELV(LIST('!!CF, 1)),MK!*SQ(DENR SIMP!* CADDR X1./1))),X,MU) COMMENT (...)**c = b; ELSE << IF !*ALLBRANCH THEN <<!!ARBINT:=!!ARBINT+1; C:=MKEXP LIST('TIMES,LIST ('ARBREAL,!!ARBINT)) >> ELSE C:=1; SOLVE1(SUBTRSQ(SIMP!* CADR X1,MULTSQ(SIMP!* LIST('EXPT,MK!*SQ GETELV(LIST('!!CF, 1)), MK!*SQ INVSQ SIMP!* CADDR X1),SIMP!* C)), X, MU) >> COMMENT (...)**(...) = b : transcendental; ELSE CONSSMMP(SUBTRSQ(SIMP!* X1,GETELV(LIST('!!CF, 1))), -MU) COMMENT SIN(...) = b; ELSE IF CAR X1='SIN THEN<< IF !*ALLBRANCH THEN << !!ARBINT:=!!ARBINT+1; F:=LIST('TIMES,2,'PI,LIST('ARBINT,!!ARBINT)) >> ELSE F:=0; C:=SIMP!* CADR X1; D:=LIST('ASIN,MK!*SQ GETELV(LIST('!!CF, 1))); SOLVE1(SUBTRSQ(C,SIMP!* LIST('PLUS,D,F)),X,MU); IF !*ALLBRANCH THEN SOLVE1(SUBTRSQ(C,SIMP!* LIST ('PLUS,'PI,MK!*SQ SUBTRSQ(SIMP!* F,SIMP!* D))), X, MU) >> COMMENT COS(...) = b; ELSE IF CAR X1='COS THEN<< IF !*ALLBRANCH THEN<<!!ARBINT:=!!ARBINT+1; C:=LIST('TIMES,2,'PI,LIST('ARBINT,!!ARBINT))>> ELSE C:=0; C:=SUBTRSQ(SIMP!* CADR X1,SIMP!* C); D:=SIMP!* LIST('ACOS,MK!*SQ GETELV(LIST('!!CF,1))); SOLVE1(SUBTRSQ(C,D), X, MU); IF !*ALLBRANCH THEN SOLVE1(ADDSQ(C,D), X, MU) >> COMMENT Unknown inverse; ELSE IF NULL(F:=GET(CAR X1,'INVERSE))THEN CONSSMMP(SUBTRSQ(SIMP!* X1,GETELV(LIST('!!CF,1))), -MU) COMMENT Other, known inverse; ELSE SOLVE1(SUBTRSQ(SIMP!* CADR X1,SIMP!* LIST(F,MK!*SQ GETELV(LIST('!!CF,1)))), X, MU)>> >> COMMENT Test for 2 powers of X1; ELSE IF F>=0 THEN << D:= SIMP!* GETELV(LIST('!!CF,2)); C := QUOTSQ(SIMP!* GETELV(LIST('!!CF,0)),D); D := QUOTSQ(SIMP!* GETELV(LIST('!!CF,1)),MULTSQ((2 ./1),D)); C:=SIMP!* LIST('EXPT, MK!*SQ SUBTRSQ(EXPTSQ(D,2),C), MK!*SQ(1 ./2)); D := ADDSQ(D, EXPTSQ(SIMP!* X1, N)); SOLVE1(SUBTRSQ(D,C), X, MU); SOLVE1(ADDSQ(D,C), X, MU) >> ELSE SOLVE22(E1,X1,X,MU) END; SYMBOLIC PROCEDURE SOLVE22(E1,X1,X,MU); BEGIN SCALAR B,C,D,F; INTEGER N; COMMENT Test for reciprocal equation, cubic, or quartic; F:=(!!HIPOW+1)/2; D:=EXPTSQ(SIMP!* X1,N); IF (FOR J:=0:F DO IF NOT(GETELV(LIST('!!CF,J)) =GETELV(LIST('!!CF,!!HIPOW-J)) ) THEN RETURN T) THEN IF (FOR J:=0:F DO IF NUMR ADDSQ(SIMP!* GETELV(LIST('!!CF,J)), SIMP!* GETELV(LIST('!!CF,!!HIPOW-J))) THEN RETURN T) THEN IF !!HIPOW=3 THEN SOLVECUBIC(D,X,MU,T) ELSE IF !!HIPOW=4 THEN SOLVEQUARTIC(D,X,MU) ELSE IF !*SOLVEINTERVAL AND UNIVARIATEP E1 THEN SOLVEINTERVAL(E1,MU) ELSE CONSSMMP(E1./1, -MU) COMMENT Antisymmetric reciprocal equation; ELSE << C:=ADDSQ(D,(-1 ./1)); SOLVE1(C, X, MU); E1:= QUOTSQ(E1./1, C); IF F+F = !!HIPOW THEN <<C:=ADDSQ(D,(1 ./1)); SOLVE1(C, X, MU); E1:= QUOTSQ(E1, C) >>; SOLVE1(E1, X, MU) >> COMMENT Symmetric reciprocal equation; ELSE IF F+F=!!HIPOW+1 THEN << C:=ADDSQ(D, 1 ./1); SOLVE1(C,X,MU); SOLVE1(QUOTSQ(E1./1, C), X, MU) >> ELSE << B:=SM!*; SETELV(LIST('!!CF, 0), SIMP!* 2); SETELV(LIST('!!CF, 1), SIMP!* '!!X); C:=ADDSQ(MULTSQ(SIMP!* GETELV(LIST('!!CF,F+1)), GETELV(LIST('!!CF,1))), SIMP!* GETELV(LIST('!!CF,F))); FOR J:=2:F DO << SETELV(LIST('!!CF, J), SUBTRSQ(MULTSQ(GETELV(LIST('!!CF,1)), GETELV(LIST('!!CF,J-1))), GETELV(LIST('!!CF,J-2)))); C:=ADDSQ(C,MULTSQ(GETELV(LIST('!!CF,J)), SIMP!* GETELV(LIST('!!CF,F+J)) )) >>; SOLVE1(C,'!!X,MU); C:=F:=NIL; WHILE NOT(SM!*=B) DO << C:=CAR SM!* . C; F:=CAR MP!* . F; SM!*:=CDR SM!*; MP!*:=CDR MP!* >>; WHILE C DO << SOLVE1(ADDSQ(1 ./1,MULTSQ(D,SUBTRSQ(D,CAAR C))), X, !*A2F CAAR F*MU); C:=CDR C >> >> END; SYMBOLIC PROCEDURE MKEXP U; (LAMBDA X; LIST('PLUS,LIST('COS,X),LIST('TIMES,'I,LIST('SIN,X)))) REVAL U; SYMBOLIC PROCEDURE SOLVECOEFF(EX, VAR); % EX is a standard form. % VAR is a kernel. % Puts the coefficients (as prefix standard quotients) of % VAR in EX into the elements of !!CF, with index equal % to the exponent divided by the GCD of all the % exponents. This GCD is put into !!GCD, and the % highest power divided by the GCD is put into % !!HIPOW. % Returns the lowest power if the highest is equal to 2; % -1 if the highest power is less than 2, or -1 if % the highest power is greater than 2. % This bizarre behaviour stems from the rewriting of the % Reduce COEFF function originally used by SOLVE. % Hopefully this will be rewritten someday without % the kludginess. % Note that !!CF (an array), !!GCD, and !!HIPOW are globals.; BEGIN SCALAR CLIST, X, OLDKORD; OLDKORD := SETKORDER(LIST(VAR)); CLIST := REORDER (EX); SETKORDER(OLDKORD); !!HIPOW := LDEG(CLIST); CLIST := COEFLIS(CLIST); !!GCD := 0; X := CLIST; WHILE X DO << !!GCD := GCDN(CAAR(X), !!GCD); X := CDR(X) >>; X := CLIST; FOR I := 0:(CAR(DIMENSION('!!CF))-1) DO SETELV(LIST('!!CF, I), NIL); WHILE X DO << SETELV(LIST('!!CF, CAAR(X)/!!GCD), MK!*SQ(CDAR(X) ./ 1)); X := CDR(X) >>; !!HIPOW := !!HIPOW/!!GCD; IF !!HIPOW=2 THEN RETURN CAAR(CLIST)/!!GCD ELSE IF !!HIPOW<2 THEN RETURN -1 ELSE RETURN -2 END; SYMBOLIC PROCEDURE SOLVEINTERVAL(EX, MUL); % EX is a standard form, MUL is an integer. Isolates % insoluble, real roots of EX in rational intervals, % stuffing result in the form INTERVL(Lowlim,Highlim) % into SM!* with multiplicity MUL put into MP!*.; BEGIN INTEGER I; REALROOT(PREPF EX,PREPSQ !*K2Q MVAR EX,!!INTERVALARRAY,'!!EXACT); FOR I := 1:GETELV LIST('!!EXACT,0) DO CONSSMMP(SIMP!* GETELV LIST('!!EXACT,I), MUL); FOR I := 1:GETELV LIST(!!INTERVALARRAY,0,0) DO CONSSMMP(SIMP!* LIST('INTERVL, GETELV LIST(!!INTERVALARRAY,I,1), GETELV LIST(!!INTERVALARRAY,I,2) ), MUL) END; SYMBOLIC PROCEDURE REALROOT(U,V,W,X); REDERR("Real root finding not yet implemented"); %***************** Procedures for solving Cubic and Quartic eqns ***; SYMBOLIC PROCEDURE SOLVECUBIC(X1, X, MU, CUBE3) ; BEGIN COMMENT Solves !!CF(3)*X1**3 + !!CF(2)*X1**2 ... X1 and X are kernels, M and MU are integers, CUBE3 is T or NIL. Returns NIL; SCALAR A,B,C,D; D:=SIMP!* GETELV(LIST('!!CF,3)); C:=QUOTSQ(SIMP!* GETELV(LIST('!!CF,2)),D); B:=QUOTSQ(SIMP!* GETELV(LIST('!!CF,1)),D); A:=QUOTSQ(SIMP!* GETELV(LIST('!!CF,0)),D); A:=MULTSQ(ADDSQ(MULTSQ((9 ./1),MULTSQ(C,B)), ADDSQ(MULTSQ ((-27 ./1),A),MULTSQ((-2 ./1),EXPTSQ(C,3)))),(1 ./54)); B := MULTSQ((-1 ./9),ADDSQ(EXPTSQ(C,2),MULTSQ((-3 ./1),B))); D := SIMP!* LIST('EXPT, MK!*SQ ADDSQ(EXPTSQ(B,3), EXPTSQ(A,2)), MK!*SQ(1 ./2)); D := SIMP!* LIST('EXPT, MK!*SQ ADDSQ(A,D),MK!*SQ(1 ./3)); A := NEGSQ QUOTSQ(B,D); B := ADDSQ(D, A); C := ADDSQ(X1, MULTSQ(C,(1 ./3))); SOLVE1(SUBTRSQ(C,B), X, MU); IF CUBE3 THEN <<C := ADDSQ(MULTSQ(B,(1 ./2)), C); D := MULTSQ(SIMP!* LIST('EXPT,MK!*SQ(-3 ./4),MK!*SQ (1 ./2)), SUBTRSQ(D,A)); SOLVE1(ADDSQ(C,D), X, MU); SOLVE1(SUBTRSQ(C,D), X, MU)>> END; SYMBOLIC PROCEDURE SOLVEQUARTIC(X1,X,MU) ; BEGIN COMMENT Solves !!CF(4)*X1**4 + !!CF(3)*X1**3 + .... X1 is a standard quotient, X is a kernel, MU is an integer, CUBE3 is T or NIL. Returns NIL; SCALAR A,B,C,D,F; F:=SIMP!* GETELV(LIST('!!CF,4)); A:=QUOTSQ(SIMP!* GETELV(LIST('!!CF,0)),F); B:=QUOTSQ(SIMP!* GETELV(LIST('!!CF,1)),F); C:=QUOTSQ(SIMP!* GETELV(LIST('!!CF,2)),F); D:=QUOTSQ(SIMP!* GETELV(LIST('!!CF,3)),F); F := ADDSQ(EXPTSQ(D,2), MULTSQ((-4 ./1),C)); SETELV(LIST('!!CF, 0), MK!*SQ NEGSQ ADDSQ(EXPTSQ(B,2),MULTSQ(A,F))); SETELV(LIST('!!CF, 1), MK!*SQ ADDSQ(MULTSQ(B,D),MULTSQ((-4 ./1),A))); SETELV(LIST('!!CF, 2), MK!*SQ NEGSQ C); SETELV(LIST('!!CF, 3), 1); SOLVECUBIC(SIMP!* X, X, MU, NIL); B := CAAR SM!*; SM!* := CDR SM!*; MP!*:= CDR MP!*; A := SIMP!* LIST('EXPT, MK!*SQ ADDSQ(EXPTSQ(B,2),MULTSQ(A, (-4 ./1))), MK!*SQ(1 ./2)); F := SIMP!* LIST('EXPT, MK!*SQ ADDSQ(F,MULTSQ(B,(4 ./1))), MK!*SQ(1 ./2)); SOLVE1(ADDSQ(EXPTSQ(X1,2),MULTSQ((1 ./2),ADDSQ(MULTSQ(X1,ADDSQ (D,F)), ADDSQ(B,A)))), X, MU); SOLVE1(ADDSQ(EXPTSQ(X1,2),MULTSQ((1 ./2),ADDSQ(MULTSQ(X1, SUBTRSQ(D,F)), SUBTRSQ(B,A)))), X, MU); END; %******************* Procedures for solving a system of eqns *******; SYMBOLIC PROCEDURE SOLVESYS(EXLIST,VARLIST); % EXLIST is a list of standard forms. % VARLIST is a list of kernels. % If EXLIST and VARLIST are of the same length and the % elements of VARLIST are linear in the elements of % exlist, and further the system of linear eqns so % defined is non-singular, then SOLVESYS returns a % list of standard quotients which are solutions of % the system, ordered as in VARLIST. % Otherwise an error results.; BEGIN SCALAR MTRX, RHS; % Coeffs and right side of system; SCALAR ROW, OLDKORD; IF LENGTH(EXLIST) NEQ LENGTH(VARLIST) THEN REDERR "SOLVESYS given unequal number of eqns & unknowns"; OLDKORD := SETKORDER(VARLIST); EXLIST := MAPCAR(EXLIST, 'REORDER); FOR EACH EX IN EXLIST DO << ROW := NIL; FOR EACH VAR IN VARLIST DO<< IF NOT DOMAINP(EX) AND (MVAR(EX)=VAR AND LDEG(EX)>1 OR (NOT FREEOFKERN(VAR, LC(EX))) OR (NOT FREEOFKERN(VAR, RED(EX))) ) THEN REDERR "SOLVE given system of non linear-fractional equations"; IF NOT DOMAINP(EX) AND MVAR(EX)=VAR THEN << ROW := !*F2Q(LC(EX)) . ROW; EX := RED(EX) >> ELSE ROW := !*F2Q(NIL) . ROW >>; RHS := LIST(!*F2Q(NEGF(EX))) . RHS; MTRX := ROW . MTRX >>; SETKORDER(OLDKORD); RETURN SOLVELNRSOLVE(MTRX, RHS) END; SYMBOLIC PROCEDURE SOLVELNRSOLVE(U,V); % U is a matrix canonical form, V a compatible matrix form. % Result is the solution,y, to the matrix equation U*y=V. % If !*SOLVESINGULAR is non-nil, introduces arbitrary constants % if necessary. Returns an error if the system represented is % inconsistent or if !*SOLVESINGULAR is nil and U is singular.; BEGIN INTEGER N, K; SCALAR X,!*S!*, PERM; X := !*EXP; !*EXP := T; N := LENGTH U; PERM := INDEXLIS(1, N); U := CAR NORMMAT AUGMENT(U,V); IF NOT !*SOLVESINGULAR THEN U := BAREISS U ELSE << U := SOLVEBAREISS(U, PERM); IF U THEN U := INSERTARBCONSTS(CDR(U), CAR(U)+1, FUNCTION MAKEARBCOMPLEX) >>; !*S!* := BACKSUB(U,N); U := MAPC2(RHSIDE(CAR !*S!*,N), FUNCTION (LAMBDA J; CANCEL(J . CDR !*S!*))); !*EXP := X; RETURN PERMUTE(U, PERM); END; SYMBOLIC PROCEDURE SOLVEBAREISS(U, PERM); %The 2-step integer preserving elimination method of Bareiss %based on the implementation of Lipson; %This is based on the Bareiss function in the Reduce matrix package, %modified to reduce singular matrices. If PERM is nil, behaves %as BAREISS, except a pair is returned for non-singular U, whose %cdr is the triangularized U. The car is the rank of U, which in %this case is always LENGTH(U). %Otherwise PERM is a list of the integers 1,2...length(U). %As columns are interchanged, then so are the elements of PERM. %In this case a pair is returned whose car is the rank of U and %whose cdr is the triangularized U. Note that, just as in BAREISS, the %lower triangular portion of the returned matrix standard form is only %implicitly all nils--the requisite RPLACAs are not performed. Also, %if PERM is non-nil and the rank,r, of U is less than the order of U, %only the first r rows of the upper triangular portion are explicitly %set. The all nil rows are only implicitly all nils. %U is a list of lists of standard forms (a matrix standard form) %corresponding to the appropriate augmented matrix. %If the value of procedure is NIL then U is singular, otherwise the %value is the triangularized form of U (in the same form); BEGIN SCALAR AA,C0,CI1,CI2,IK1,IJ,KK1,KJ,K1J,K1K1, UI,U1,X,K1COL,KIJ,FLG; INTEGER K,K1,COL,MAXCOL; %U1 points to K-1th row of U %UI points to Ith row of U %IJ points to U(I,J) %K1J points to U(K-1,J) %KJ points to U(K,J) %IK1 points to U(I,K-1) %KK1 points to U(K,K-1) %K1K1 points to U(K-1,K-1) %M in comments is number of rows in U %N in comments is number of columns in U; MAXCOL := LENGTH(U); AA:= 1; K:= 2; K1:=1; U1:=U; GO TO PIVOT; AGN: U1 := CDR U1; IF NULL CDR U1 OR NULL CDDR U1 THEN IF PERM AND CDR(U1) AND NULL(CAR(IJ := PNTH(NTH(U, MAXCOL), MAXCOL))) THEN << MAPC(CDR(IJ), FUNCTION LAMBDA(X); IF X THEN RETURN NIL); RETURN (MAXCOL-1).U >> ELSE RETURN MAXCOL.U; AA:=NTH(CAR U1,K); %AA := U(K,K); K:=K+2; K1:=K-1; U1:=CDR U1; PIVOT: %pivot algorithm; COL := K1; K1J:= K1K1 := PNTH(CAR U1,K1); PIV1: K1COL := PNTH(CAR(U1), COL); IF CAR K1COL THEN GO TO L2; UI:= CDR U1; %I := K; L: IF NULL UI THEN IF PERM THEN IF COL>=MAXCOL THEN RETURN (K1-1).U ELSE << COL := COL+1; GO TO PIV1 >> ELSE RETURN NIL ELSE IF NULL CAR(IJ := PNTH(CAR UI,COL)) THEN GO TO L1; L0: IF NULL IJ THEN GO TO L2; X := CAR IJ; RPLACA(IJ,NEGF CAR K1COL); RPLACA(K1COL,X); IJ:= CDR IJ; K1COL:= CDR K1COL; GO TO L0; L1: UI:= CDR UI; GO TO L; L2: SWAPCOLUMNS(U, K1, COL, PERM); COL := K; PIV2: UI:= CDR U1; %I:= K; L21: IF NULL UI THEN IF PERM THEN IF COL>=MAXCOL THEN << FLG := T; WHILE FLG AND U1 DO << IK1 := PNTH(CAR(U1), K1); IJ := PNTH(IK1, MAXCOL-K1+2); KIJ := PNTH(K1K1, MAXCOL-K1+2); WHILE FLG AND IJ DO IF ADDF(MULTF(CAR(K1K1), CAR(IJ)), MULTF(CAR(IK1), NEGF(CAR(KIJ))) ) THEN FLG := NIL ELSE IJ := CDR(IJ); U1 := CDR(U1) >>; IF FLG THEN RETURN (K-1).U ELSE RETURN NIL >> ELSE << COL := COL+1; GO TO PIV2 >> ELSE RETURN NIL; IJ:= PNTH(CAR UI,K1); C0:= ADDF(MULTF(CAR K1K1,NTH(IJ, COL-K+2)), MULTF(NTH(K1K1, COL-K+2),NEGF CAR IJ)); IF C0 THEN GO TO L3; UI:= CDR UI; %I:= I+1; GO TO L21; L3: SWAPCOLUMNS(U, K, COL, PERM); C0:= QUOTF!*(C0,AA); KK1 := KJ := PNTH(CADR U1,K1); %KK1 := U(K,K-1); IF CDR U1 AND NULL CDDR U1 THEN GO TO EV0 ELSE IF UI EQ CDR U1 THEN GO TO COMP; L31: IF NULL IJ THEN GO TO COMP; %IF I>N THEN GO TO COMP; X:= CAR IJ; RPLACA(IJ,NEGF CAR KJ); RPLACA(KJ,X); IJ:= CDR IJ; KJ:= CDR KJ; GO TO L31; %pivoting complete; COMP: IF NULL CDR U1 THEN GO TO EV; UI:= CDDR U1; %I:= K+1; COMP1: IF NULL UI THEN GO TO EV; %IF I>M THEN GO TO EV; IK1:= PNTH(CAR UI,K1); CI1:= QUOTF!*(ADDF(MULTF(CADR K1K1,CAR IK1), MULTF(CAR K1K1,NEGF CADR IK1)), AA); CI2:= QUOTF!*(ADDF(MULTF(CAR KK1,CADR IK1), MULTF(CADR KK1,NEGF CAR IK1)), AA); IF NULL CDDR K1K1 THEN GO TO COMP3;%IF J>N THEN GO TO COMP3; IJ:= CDDR IK1; %J:= K+1; KJ:= CDDR KK1; K1J:= CDDR K1K1; COMP2: IF NULL IJ THEN GO TO COMP3; RPLACA(IJ,QUOTF!*(ADDF(MULTF(CAR IJ,C0), ADDF(MULTF(CAR KJ,CI1), MULTF(CAR K1J,CI2))), AA)); IJ:= CDR IJ; KJ:= CDR KJ; K1J:= CDR K1J; GO TO COMP2; COMP3: UI:= CDR UI; GO TO COMP1; EV0:IF NULL C0 THEN RETURN; EV: KJ := CDR KK1; X := CDDR K1K1; %X := U(K-1,K+1); RPLACA(KJ,C0); EV1:KJ:= CDR KJ; IF NULL KJ THEN GO TO AGN; RPLACA(KJ,QUOTF!*(ADDF(MULTF(CAR K1K1,CAR KJ), MULTF(CAR KK1,NEGF CAR X)), AA)); X := CDR X; GO TO EV1 END; SYMBOLIC PROCEDURE SWAPCOLUMNS(MATRX, COL1, COL2, PERM); IF COL1=COL2 THEN MATRX ELSE << SWAPELEMENTS(PERM, COL1, COL2); FOR EACH U IN MATRX DO SWAPELEMENTS(U, COL1, COL2); MATRX >>; SYMBOLIC PROCEDURE SWAPELEMENTS(LST, I, J); % Swaps the Ith and Jth elements of the list LST al la % RPLACA and returns nil.; BEGIN SCALAR TEMP; IF I>J THEN << TEMP := I; I := J; J := TEMP >>; LST := PNTH(LST, I); I := J-I+1; TEMP := NTH(LST, I); RPLACA(PNTH(LST, I), CAR(LST)); RPLACA(LST, TEMP) END; SYMBOLIC PROCEDURE INDEXLIS(M, N); % M,N are integers. Returns the list (M M+1 M+2 ... N-1 N); IF M<=N THEN M . INDEXLIS(M+1,N); SYMBOLIC PROCEDURE INSERTARBCONSTS(M, ZEROROW, ARBFN); % M is a matrix standard form, representing a % matrix which has been row reduced. All elements below % the principal diagonal are implicitly nil, as are all % elements in row ZEROROW and below. It is such a form % as is returned by SOLVEBAREISS with a non-nil second % argument. It inserts approriate arbitrary constants in % the inhomogeneous portion, and 1's on the main diagonal % except for the last row, which gets the new determinant % of the square submatrix. Calls ARBFN to generate arbitrary % constants.; BEGIN SCALAR U, V, NEWDET; INTEGER N; N := LENGTH(M); IF ZEROROW<=N THEN << NEWDET := 1; U := M; FOR I := 1:(ZEROROW-1) DO << NEWDET := MULTF(NEWDET, NTH(CAR(U), I)); U := CDR(U) >>; FOR I := ZEROROW:(N-1) DO << V := PNTH(CAR(U), I); RPLACA(V, 1); V := CDR(V); FOR J := I+1:N DO << RPLACA(V, NIL); V := CDR(V) >>; WHILE V DO << RPLACA(V, !*K2F EVAL LIST ARBFN); V := CDR(V) >>; U := CDR(U) >>; V := PNTH(CAR(U), N); RPLACA(V, NEWDET); V := CDR(V); WHILE V DO << RPLACA(V, MULTF(NEWDET, !*K2F EVAL LIST ARBFN)); V := CDR(V) >> >>; RETURN M END; SYMBOLIC PROCEDURE PERMUTE(U, V); % U is a list. V is a list of the numbers 1,2,...LENGTH(U), permuted; % Returns a constructed list of the elements of U permuted by V.; IF V THEN NCONC(LIST(NTH(U,CAR(V))), PERMUTE(U, CDR(V))); SYMBOLIC PROCEDURE MAKEARBCOMPLEX(); BEGIN SCALAR ANS; ANS := NUMR(SIMP!*(LIST('ARBCOMPLEX, !!ARBINT))); !!ARBINT := !!ARBINT+1; RETURN ANS END; %******** Algebraic Let Statements and related declarations ********; PUT('ASIN, 'INVERSE, 'SIN); PUT('ACOS, 'INVERSE, 'COS); ALGEBRAIC << OPERATOR SOL, INTERVL, ARBCOMPLEX, ARBREAL, ARBINT, LST; COMMENT Supply missing argument and simplify 1/4 roots of unity; LET E**(I*PI/2) = I, E**(I*PI) = -1, E**(3*I*PI/2)=-I; FOR ALL N SUCH THAT FIXP N LET COS((N*PI)/2)= 0; LET COS(PI/2)=0; FOR ALL N SUCH THAT FIXP N LET SIN((N*PI)/2)= IF REMAINDER(ABS N,4)<2 THEN 1 ELSE -1; LET SIN(PI/2)=1; FOR ALL N SUCH THAT FIXP N LET COS((N*PI)/3)= (IF N=4 OR REMAINDER(ABS N+2,6)>3 THEN -1 ELSE 1)/2; LET COS(PI/3)=1/2; FOR ALL N SUCH THAT FIXP N LET SIN((N*PI)/3)= (IF REMAINDER(ABS N,6)<3 THEN 1 ELSE -1)*SQRT(3)/2; LET SIN(PI/3)=SQRT(3)/2; FOR ALL N SUCH THAT FIXP N LET COS((N*PI)/4)= (IF REMAINDER(ABS N+2,8)<4 THEN 1 ELSE -1)*SQRT(2)/2; LET COS(PI/4)=SQRT 2/2; FOR ALL N SUCH THAT FIXP N LET SIN((N*PI)/4)= (IF REMAINDER(ABS N,8)<4 THEN 1 ELSE -1)*SQRT(2)/2; LET SIN(PI/4)=SQRT(2)/2; FOR ALL N SUCH THAT FIXP N LET COS((N*PI)/6)= (IF REMAINDER(ABS N+2,12)<6 THEN 1 ELSE -1)*SQRT(3)/2; LET COS(PI/6)=SQRT 3/2; FOR ALL N SUCH THAT FIXP N LET SIN((N*PI)/6)= (IF REMAINDER(ABS N,12)<6 THEN 1 ELSE -1)/2; LET SIN(PI/6)=1/2; COMMENT Rules for reducing the number of distinct kernels in an equation; FOR ALL A,B,X SUCH THAT RATNUMP C AND RATNUMP D LET SOL(A**C-B**D, X) = A**(C*LCMD(C,D)) - B**(D*LCMD(C,D)); FOR ALL A,B,C,D,X SUCH THAT FREEOFKERN(X,A) AND FREEOFKERN(X,C) LET SOL(A**B-C**D, X) = E**(B*LOG A - D*LOG C); FOR ALL A,B,C,D,X SUCH THAT FREEOFKERN(X,A) AND FREEOFKERN(X,C) LET SOL(A*LOG B + C*LOG D, X) = B**A*D**C - 1, SOL(A*LOG B - C*LOG D, X) = B**A - D**C; FOR ALL A,B,C,D,F,X SUCH THAT FREEOFKERN(X,A) AND FREEOFKERN(X,C) LET SOL(A*LOG B + C*LOG D + F, X) = SOL(LOG(B**A*D**C) + F, X), SOL(A*LOG B + C*LOG D - F, X) = SOL(LOG(B**A*D**C) - F, X), SOL(A*LOG B - C*LOG D + F, X) = SOL(LOG(B**A/D**C) + F, X), SOL(A*LOG B - C*LOG D - F, X) = SOL(LOG(B**A/D**C) - F, X); FOR ALL A,B,D,F,X SUCH THAT FREEOFKERN(X,A) LET SOL(A*LOG B + LOG D + F, X) = SOL(LOG(B**A*D) + F, X), SOL(A*LOG B + LOG D - F, X) = SOL(LOG(B**A*D) - F, X), SOL(A*LOG B - LOG D + F, X) = SOL(LOG(B**A/D) + F, X), SOL(A*LOG B - LOG D - F, X) = SOL(LOG(B**A/D) - F, X), SOL(LOG D - A*LOG B + F, X) = SOL(LOG(D/B**A) + F, X), SOL(LOG D - A*LOG B - F, X) = SOL(LOG(D/B**A) - F, X); FOR ALL A,B,D,X SUCH THAT FREEOFKERN(X,A) LET SOL(A*LOG B + LOG D, X) = B**A*D - 1, SOL(A*LOG B - LOG D, X) = B**A - D, SOL(LOG D - A*LOG B, X) = D - B**A; FOR ALL A,B,C,X LET SOL(LOG A + LOG B + C, X) = SOL(LOG(A*B) + C, X), SOL(LOG A - LOG B + C, X) = SOL(LOG(A/B) + C, X), SOL(LOG A + LOG B - C, X) = SOL(LOG(A*B) - C, X), SOL(LOG A - LOG B - C, X) = SOL(LOG(A/B) - C, X); FOR ALL A,C,X SUCH THAT FREEOFKERN(X,C) LET SOL(LOG A + C, X) = A - E**C, SOL(LOG A - C, X) = A - E**(-C); FOR ALL A,B,X LET SOL(LOG A + LOG B, X) = A*B - 1, SOL(LOG A - LOG B, X) = A - B, SOL(COS A - SIN B, X) = SOL(COS A - COS(PI/2-B), X), SOL(SIN A + COS B, X) = SOL(SIN A - SIN(B-PI/2), X), SOL(SIN A - COS B, X) = SOL(SIN A - SIN(PI/2-B), X), SOL(SIN A + SIN B, X) = SOL(SIN A - SIN(-B), X), SOL(SIN A - SIN B, X) = IF !*ALLBRANCH THEN SIN((A-B)/2)* COS((A+B)/2) ELSE A-B, SOL(COS A + COS B, X) = IF !*ALLBRANCH THEN COS((A+B)/2)* COS((A-B)/2) ELSE A+B, SOL(COS A - COS B, X) = IF !*ALLBRANCH THEN SIN((A+B)/2)* SIN((A-B)/2) ELSE A-B, SOL(ASIN A - ASIN B, X) = A-B, SOL(ASIN A + ASIN B, X) = A+B, SOL(ACOS A - ACOS B, X) = A-B, SOL(ACOS A + ACOS B, X) = A+B; LET COS(PI/2)=0>>; END;