1
0
mirror of https://github.com/PDP-10/its.git synced 2026-01-13 15:27:28 +00:00
PDP-10.its/src/rat/psolve.22
Eric Swenson 85994ed770 Added files to support building and running Macsyma.
Resolves #284.

Commented out uses of time-origin in maxtul; mcldmp (init) until we
can figure out why it gives arithmetic overflows under the emulators.

Updated the expect script statements in build_macsyma_portion to not
attempt to match expected strings, but simply sleep for some time
since in some cases the matching appears not to work.
2018-03-11 13:10:19 -07:00

278 lines
7.0 KiB
Plaintext
Executable File
Raw Blame History

This file contains invisible Unicode characters

This file contains invisible Unicode characters that are indistinguishable to humans but may be processed differently by a computer. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

;; -*- Mode: Lisp; Package: Macsyma; -*-
;; (c) Copyright 1976, 1984 Massachusetts Institute of Technology
;; All Rights Reserved.
;; Enhancements (c) Copyright 1984 Symbolics Inc.
;; All Rights Reserved.
;; The data and information in the Enhancements are proprietary to, and
;; a valuable trade secret of, SYMBOLICS, INC., a Delaware corporation.
;; They are given in confidence by SYMBOLICS, pursuant to the license
;; agreement between Symbolics and their recipient, and may not be used,
;; reproduced, or copied, or distributed to any other party, in whole or
;; in part, without the prior written consent of SYMBOLICS except as
;; permitted by the license agreement.
(macsyma-module psolve)
(DECLARE (GENPREFIX PSO)
(SPECIAL MULT *ROOTS *FAILURES $SOLVEFACTORS))
(DECLARE (SPLITFILE SCUBIC))
(DEFMVAR FLAG4 NIL)
(DEFMFUN SOLVECUBIC (X)
(PROG (S1 A0 A1 A2 DISCR LCOEF ADIV3 OMEGA^2 PDIV3 QDIV-2
OMEGA Y1 U y2)
(SETQ X (CDR X))
(SETQ LCOEF (CADR X))
(SETQ ADIV3
(LIST '(MTIMES)
'((RAT) -1. 3.)
(RDIS (SETQ A2 (RATREDUCE (PTERM X 2.)
LCOEF)))))
(SETQ A1 (RATREDUCE (PTERM X 1.) LCOEF))
(SETQ A0 (RATREDUCE (PTERM X 0.) LCOEF))
(SETQ S1 '((MTIMES)
((RAT) 1. 2.)
$%I
((MEXPT) 3. ((RAT) 1. 2.))))
(SETQ OMEGA (LIST '(MPLUS)
'((RAT) -1. 2.)
S1)
OMEGA^2 (LIST '(MPLUS)
'((RAT) -1. 2.)
(LIST '(MTIMES) -1. S1)))
(SETQ PDIV3
(RDIS (RATPLUS (RATTIMES A1 '(1. . 3.) T)
(RATTIMES (RATEXPT A2 2.)
'(-1. . 9.)
T))))
(AND (NOT (EQUAL PDIV3 0.)) (GO HARDER))
(SETQ
Y1
(SIMPTIMES
(LIST
'(MTIMES)
'((RAT) 1. 3.)
(LIST '(MPLUS)
(SIMPNRT (RDIS (setq y2 (RATPLUS (RATEXPT A2 3.)
(RATTIMES '(-27. . 1.)
A0
T))))
3)
(LIST '(MTIMES) -1. (RDIS A2))))
1.
NIL))
(AND FLAG4 (RETURN (SOLVE3 Y1 MULT)))
(setq y2 (simpnrt (rdis (rattimes y2 '(1. . 27.) t))
3))
(RETURN (MAPC #'(LAMBDA (J) (SOLVE3 J MULT))
(LIST Y1
(LIST '(MPLUS)
(LIST '(MTIMES)
OMEGA
Y2)
ADIV3)
(LIST '(MPLUS)
(LIST '(MTIMES)
OMEGA^2
Y2)
ADIV3))))
HARDER
(SETQ
QDIV-2
(RDIS (RATPLUS (RATTIMES (RATPLUS (RATTIMES A1 A2 T)
(RATTIMES '(-3. . 1.)
A0
T))
'(1. . 6.)
T)
(RATTIMES (RATEXPT A2 3.)
'(-1. . 27.)
T))))
(COND ((EQUAL QDIV-2 0.)
(SETQ U (SIMPNRT PDIV3 2))
(SETQ Y1 ADIV3))
(T (SETQ DISCR (SIMPLUS (LIST '(MPLUS)
(LIST '(MEXPT)
PDIV3
3.)
(LIST '(MEXPT)
QDIV-2
2.))
1.
NIL))
(COND ((EQUAL DISCR 0.)
(SETQ U (SIMPNRT QDIV-2 3)))
(T (SETQ DISCR (SIMPNRT DISCR 2))
(AND (COMPLICATED DISCR)
(SETQ DISCR (adispline DISCR)))
(SETQ U (SIMPEXPT (LIST '(MEXPT)
(LIST '(MPLUS)
QDIV-2
DISCR)
'((RAT) 1 3)) 1 NIL))
(AND (COMPLICATED U)
(SETQ U (adispline U)))))))
(IF (EQUAL U 0) (MERROR "Arithmetic overflow - SOLVECUBIC"))
(OR Y1
(SETQ Y1 (SIMPLUS (LIST '(MPLUS)
ADIV3
U
(LIST '(MTIMES)
-1.
PDIV3
(LIST '(MEXPT) U -1)))
1.
NIL)))
(RETURN
(COND (FLAG4 (SOLVE3 Y1 MULT))
(T (MAPC
#'(LAMBDA (J) (SOLVE3 J MULT))
(LIST Y1
(LIST '(MPLUS)
ADIV3
(LIST '(MTIMES) OMEGA U)
(LIST '(MTIMES)
-1.
PDIV3
OMEGA^2
(LIST '(MEXPT) U -1)))
(LIST '(MPLUS)
ADIV3
(LIST '(MTIMES) OMEGA^2 U)
(LIST '(MTIMES)
-1.
PDIV3
OMEGA
(LIST '(MEXPT) U -1))))))))))
(DECLARE (SPLITFILE SQUART))
(DEFMVAR FLAG4 NIL)
(DEFMFUN SOLVEQUARTIC (X)
(PROG (A0 A1 A2 B1 B2 B3 B0 LCOEF Z1 R TR1 TR2 D D1 E SQB3)
(SETQ X (CDR X) LCOEF (CADR X))
(SETQ B3 (RATREDUCE (PTERM X 3) LCOEF))
(SETQ B2 (RATREDUCE (PTERM X 2) LCOEF))
(SETQ B1 (RATREDUCE (PTERM X 1) LCOEF))
(SETQ B0 (RATREDUCE (PTERM X 0) LCOEF))
(SETQ A2 (RATMINUS B2))
(SETQ A1 (RATDIF (RATTIMES B1 B3 T)
(SETQ A0 (RATTIMES B0 '(4 . 1) T))))
(SETQ A0
(RATDIF (RATDIF (RATTIMES B2 A0 T)
(RATTIMES (SETQ SQB3
(RATEXPT B3 2))
B0
T))
(RATEXPT B1 2.)))
(SETQ
TR2
(SIMPLIFY (RDIS
(RATTIMES
'(1. . 4.)
(RATDIF (RATDIF (RATTIMES B3
(RATTIMES B2
'(4. . 1.)
T)
T)
(RATTIMES '(8. . 1.) B1 T))
(RATTIMES SQB3 B3 NIL))
T))))
(SETQ Z1 (RESOLVENT A2 A1 A0))
(SETQ R
(SIMPLUS (LIST '(MPLUS)
Z1
(RDIS (RATDIF (RATTIMES SQB3
'(1. . 4.)
T)
B2)))
1.
NIL))
(AND (EQUAL R 0.) (GO L0))
(SETQ R (SIMPNRT R 2))
(AND (COMPLICATED R) (SETQ R (adispline R)))
(AND (COMPLICATED TR2) (SETQ TR2 (adispline TR2)))
(SETQ TR1
(SIMPLUS (LIST '(MPLUS)
(RDIS (RATDIF (RATTIMES SQB3
'(1. . 2.)
T)
B2))
(LIST '(MTIMES) -1. Z1))
1.
NIL))
(AND (COMPLICATED TR1) (SETQ TR1 (adispline TR1)))
(SETQ TR2 (DIV* TR2 R))
(GO LB1)
L0 (SETQ D1
(SIMPNRT (SIMPLIFY (LIST '(MPLUS)
(LIST '(MEXPT) Z1 2.)
(LIST '(MTIMES)
-4.
(RDIS B0))))
2))
(SETQ TR2 (SIMPLIFY (LIST '(MTIMES) 2. D1)))
(AND (COMPLICATED TR2) (SETQ TR2 (adispline TR2)))
(SETQ TR1
(SIMPLIFY (RDIS (RATDIF (RATTIMES SQB3 '(3. . 4.) T)
(RATTIMES B2 '(2. . 1.) T)))))
(AND (COMPLICATED TR1) (SETQ TR1 (adispline TR1)))
LB1 (SETQ D (SIMPNRT (SIMPLIFY (LIST '(MPLUS) TR1 TR2)) 2))
(SETQ E
(SIMPNRT (SIMPLIFY (LIST '(MPLUS)
TR1
(LIST '(MTIMES) -1. TR2)))
2))
(SETQ D (DIV* D 2.))
(AND (COMPLICATED D) (SETQ D (adispline D)))
(SETQ E (DIV* E 2.))
(AND (COMPLICATED E) (SETQ E (adispline E)))
(SETQ A2 (RDIS (RATTIMES B3 '(-1. . 4.) T)))
(SETQ A1 (DIV* R 2.))
(SETQ Z1
(LIST (LIST '(MPLUS) A2 A1 D)
(LIST '(MPLUS)
A2
A1
(LIST '(MTIMES) -1. D))
(LIST '(MPLUS)
A2
(LIST '(MTIMES) -1. A1)
E)
(LIST '(MPLUS)
A2
(LIST '(MTIMES) -1. A1)
(LIST '(MTIMES) -1. E))))
(RETURN (MAPC #'(LAMBDA (J) (SOLVE3 J MULT)) Z1))))
;;; Solves resolvent cubic equation generated from quartic.
(DEFUN RESOLVENT (A2 A1 A0)
(PROG (*ROOTS FLAG4 *FAILURES $solvefactors) ;undoes binding in
(SETQ FLAG4 T $solvefactors t) ;algsys
(SOLVE (SIMPLUS (LIST '(MPLUS)
(LIST '(MEXPT)
'YY
3.)
(LIST '(MTIMES)
(RDIS A2)
(LIST '(MEXPT)
'YY
2))
(LIST '(MTIMES)
(RDIS A1)
'YY)
(RDIS A0))
1
NIL)
'YY
1)
(IF (MEMBER 0 *ROOTS) (RETURN 0))
(RETURN (CADDAR (CDR (REVERSE *ROOTS))))))