Trailing-Edge
-
PDP-10 Archives
-
BB-H506D-SM_1983
-
cobol/source/expon.mac
There are 7 other files named expon.mac in the archive. Click here to see a list.
; UPD ID= 2970 on 6/25/80 at 11:48 AM by NIXON
TITLE EXPON FOR LIBOL
SUBTTL EXPONENTIATION KERRY BENSMAN/ALB/CAM/DMN
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
; OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.
;
;COPYRIGHT (C) 1974, 1981 BY DIGITAL EQUIPMENT CORPORATION
;EDIT HISTORY
;V12*******************
;NAME DATE COMMENTS
;DMN 14-APR-79 ADD DOUBLE PRECISION FUNCTIONS
;EHM 28-APR-78 [527] FIX ON SIZE ERROR FOR EXPONENTIATION
EXTERNAL OVFLO.
SEARCH LBLPRM
HISEG
SALL
;ACCUMULATOR DEFINITIONS
A= 0
B= 1
C= 4
D= 5
E= 6
F= 7
G= 10
H= 11
K= 13
L= 14
PA= 16
P= 17
DEFINE DOUBLE (A,B)<EXP A,B>
;THESE ROUTINES CALCULATE A**B
;THE GENERAL CALLING SEQUENCE IS:
; MOVE 0,A
; MOVEI 16,B
; PUSHJ 17,<ROUTINE>
;FOR THE SINGLE PRECISION CASES, THE BASE IS IN ACCUMULATOR 0
;AND THE EXPONENT IS REFERENCED BY ACCUMULATOR 16.
;FOR DOUBLE PRECISION BASES,THE HIGH ORDER WORD IS IN
;ACCUMULATOR 0 AND THE LOW ORDER WORD IS IN ACCUMULATOR 1.
;FOR DOUBLE PRECISION EXPONENTS, THE HIGH ORDER WORD IS
;REFERENCED BY 0(16) AND THE LOW ORDER WORD IS REFERENCED BY 1(16).
;NOTE ACCS 2 AND 3 CAN NOT BE USED AS SCRATCH SINCE ACC 16 MAY POINT TO THEM.
;ALL ANSWERS ARE RETURNED IN FLOATING POINT FORMAT.
SUBTTL COMP-1 RAISED TO 1-WORD COMP POWER
;THIS ROUTINE CALCULATES A COMP-1 FLOATING POINT NUMBER
;TO A FIXED POINT (1-WORD) POWER.
;THE CALCULATION IS A**B, WHERE B IS OF THE FORM
; B=Q(0) + Q(1)*2 + Q(2)*4 + ...WHERE Q(I)=0 OR 1
;THERE ARE NO RESTRICTIONS ON THE BASE OR EXPONENT
ENTRY E.C3C1
E.C3C1: MOVE C,0(PA) ;MOVE EXP TO SAFE PLACE
EC3C1: JUMPE A,ATEST ;SPECIAL CASE IF BASE = 0.
MOVSI D,(1.0) ;GET 1.0 IN AC D.
JUMPGE C, FEXP2 ;IS EXPONENT POSITIVE?
MOVMS C ;NO, MAKE IT POSITIVE
PUSHJ P, FEXP2 ;CALL MAIN PART OF PROGRAM.
MOVSI C, (1.0) ;GET 1.0 IN "B2.
FDVM C, A ;FORM 1/(A**B) FOR NEG. EXPONENT.
POPJ P, ;RETURN.
FEXP1: FMP A, A ;FORM A**N, FLOATING POINT.
JFCL 1,AINF ;IF OVER/UNDERFLOW, GO TO AINF.
LSH C, -1 ;SHIFT EXPONENT FOR NEXT BIT.
FEXP2: JFCL 1,.+1 ;CLEAR OVERFLOW FLAG
TRZE C, 1 ;IS THE BIT ON?
FMP D, A ;YES, MULTIPLY ANSWER BY A**N.
JFCL 1,AINF ;IF OVER/UNDERFLOW, GO TO AINF.
JUMPN C, FEXP1 ;UPDATE A**N UNLESS ALL THROUGH.
MOVE A, D ;PICK UP RESULT FROM D.
POPJ P,
ATEST: JUMPGE C,AZERO ;EXIT IF BASE =0,EXP>=0.
AINF: HRLOI A,377777 ;ANS.=+INFINITY
SETOM OVFLO. ;[527] SET OVERFLOW ON
POPJ 17, ;AND EXIT.
AZERO: SETZ A, ;ANSWER IS ZERO
POPJ P,
SUBTTL COMP-1 RAISED TO COMP-1 POWER
;THIS ROUTINE CALCULATES A COMP-1 FLOATING POINT NUMBER
;RAISED TO A COMP-1 FLOATING POINT POWER.
;THE CALCULATION IS
; A**B= EXP(B*LOG(A))
;IF THE EXPONENT IS AN INTEGER < 2**35 IN MAGNITUDE, THE
;RESULT WILL BE COMPUTED USING "E.C3C1" AND THE ANSWER
;WILL HAVE THE CORRECT SIGN. (REMEMBER THAT THE "INTEGER"
;HAS ONLY 27 SIGNIFCANT BITS.)
;SINCE NEGATIVE NUMBERS RAISED TO NON-INTEGER POWERS YIELD
;COMPLEX ANSWERS, THE MAIN ALGORITHM CALCULATES
; EXP(B*LOG(ABSF(A)))
ENTRY E.C3C3
E.C3C3: MOVE C,0(PA) ;GET EXPONENT
EC3C3A: JUMPE A,ATEST ;SPECIAL IF BASE = 0
MOVM E,C ;SET EXP. POSITIVE.
MOVEI D,0 ;CLEAR AC D TO ZERO
LSHC D,11 ;SHIFT 9 PLACES LEFT
SUBI D,200 ;TO OBTAIN SHIFTING FACTOR
JUMPLE D,EXP3GO ;IS D > 0
HRR F,D ;SET UP F AS AN INDEX REG.
SETZ D, ;CLEAR OUT AC D
LSH E,-1 ;RIGHT ADJUST EXP TO BIT 1.
ASHC D,(F) ;SHIFT LFT BY CONTENTS OF F
JFCL 1,EXP3GO ;IF OVERFLOW, GO TO EXP3GO.
JUMPN E,EXP3GO ;IS EXPONENT AN INTEGER ?
SKIPGE C ;YES, WAS IT NEGATIVE?
MOVNS D ;YES
MOVE C,D ;NO
JRST EC3C1 ;USE INTEGER ROUTINE, THEN RETURN
EXP3GO: MOVM E,A ;GET ABS(BASE) IF NE 0 OR 1.
MOVE D,C ;SAVE AC C.
PUSHJ P,ALOG. ;CALCULATE LOG OF "B"
FMPRM E,D ;CALCULATE A*LOG(B)
JRST EXPON. ;CALCULATE EXP(A*LOG(B)) AND RETURN
SUBTTL ALOG.
;FLOATING POINT SINGLE PRECISION LOGARITHM FUNCTION
;LOG(ABSF(X)) IS CALCULATED BY THE SUBROUTINE, AND AN
;ARGUMENT OF ZERO IS RETURNED AS MINUS INFINITY.
;ALOG IS THE ENTRY POINT FOR LOGE(X).
;FOR LOGE(X), THE ALGORITHM IS:
; LOGE(X) = (I + LOG2(F))*LOGE(2)
; WHERE X = (F/2)*2^(I+1), AND LOG2(F) IS GIVEN BY
; LOG2(F) = C1*Z + C3*Z^3 + C5*Z^5 - 1/2
; AND Z = (F-SQRT(2))/(F+SQRT(2))
;THE CALLING SEQUENCE FOR THE ROUTINE IS AS FOLLOWS:
; MOVE E,ARGUMENT
; PUSHJ P,ALOG.
;THE ANSWER IS RETURNED IN ACCUMULATOR E
ALOG.: JUMPE E, LZERO ;CHECK FOR ZERO ARGUMENT
CAMN E, ONE ;CHECK FOR 1.0 ARGUMENT
JRST ZERANS ;IT IS 1.0 RETURN ZERO ANS.
ASHC E, -33 ;SEPARATE FRACTION FROM EXPONENT
ADDI E, 211000 ;FLOAT THE EXPONENT AND MULT. BY 2
MOVSM E, G ;NUMBER NOW IN CORRECT FL. FORMAT
MOVSI E, 567377 ;SET UP -401.0 IN E
FADM E, G ;SUBTRACT 401 FROM EXP.*2
ASH F, -8 ;SHIFT FRACTION FOR FLOATING
TLC F, 200000 ;FLOAT THE FRACTION PART
FAD F, L1 ;F = F-SQRT(2.0)/2.0
MOVE E, F ;PUT RESULTS IN E
FAD E, L2 ;E = E+SQRT(2.0)
FDV F, E ;F = F/E
MOVE H, F ;STORE NEW VARIABLE IN H
FMP F, F ;CALCULATE Z^2
MOVE E, L3 ;PICK UP FIRST CONSTANT
FMP E, F ;MULTIPLY BY Z^2
FAD E, L4 ;ADD IN NEXT CONSTANT
FMP E, F ;MULTIPLY BY Z^2
FAD E, L5 ;ADD IN NEXT CONSTANT
FMP E, H ;MULTIPLY BY Z
FAD E, G ;ADD IN EXPONENT TO FORM LOG2(X)
FMP E, L7 ;MULTIPLY TO FORM LOGE(X)
POPJ P,
ZERANS: TDZA E, E ;MAKE ANSWER ZERO
LZERO: MOVE E,MIFI ;PICK UP MINUS INFINITY
POPJ P,
;CONSTANTS
L1: 577225754146 ;-0.707106781187
L2: 201552023632 ;1.414213562374
L3: 200462532521 ;0.5989786496
L4: 200754213604 ;0.9614706323
L5: 202561251002 ;2.8853912903
L7: 200542710300 ;0.69314718056
MIFI: 400000000001 ;LARGEST NEGATIVE FLOATING NUMBER
SUBTTL SINGLE PRECISION EXPONENTIAL FUNCTION
;FLOATING POINT SINGLE PRECISION EXPONENTIAL FUNCTION
;CALCULATE EXP(C).
;IF X<=-89.415..., THE PROGRAM RETURNS ZERO AS THE ANSWER
;IF X>=88.029..., THE PROGRAM RETURNS 377777777777 AS THE ANSWER
;THE RANGE OF THE ARGUMENT IS REDUCED AS FOLLOWS:
;EXP(X) = 2**(X*LOG(E)BASE2) = 2**(M+F)
;WHERE M IS AN INTEGER AND F IS B FRACTION
;2**M IS CALCULATED BY ALGEBRAICALLY ADDING M TO THE EXPONENT
;OF THE RESULT OF 2**F. 2**F IS CALCULATED AS
;2**F = 2(0.5+F(B+A*F^2 - F-C(F^2 +D)**-1)**-1)
;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
; MOVE PA,ARGUEMENT
; PUSHJ P,EXPON.
;THE ANSWER IS RETURNED IN ACCUMULATOR 0
EXPON.: CAMGE D,E77 ;IS EXP. < -89.41...?
JRST AZERO ;YES, GO TO EXIT.
CAMLE D,E7 ;IS EXP. > +88.029..?
JRST AINF ;YES, GET LARGEST FLOATING NUMBER
;FALL INTO STANDARD ALGORITHM
EXP1: SETZ G, ;INITIALIZE G
MULI D, 400 ;SEPARATE FRACTION AND EXPONENT
TSC D, D ;GET "B" POSITIVE EXPONENT
MUL E, E5 ;FIXED POINT MULTIPLY BY LOG2(E)
ASHC E, -242(D) ;SEPARATE FRACTION AND INTEGER
AOSG E ;ALGORITHM CALLS FOR MULT. BY 2
ADDI E,1 ;ADJUST IF FRACTION WAS NEGATIVE
HRRZM E,H ;SAVE FOR FUTURE SCALING
JUMPG F,EXP2 ;GO AHEAD IF ARG > 0.
TRNE F,377 ;ARE ALL THESE BITS 0?
ADDI F,200 ;NO. FIX UP.
EXP2: ASH F, -8 ;MAKE ROOM FOR EXPONENT
TLC F, 200000 ;PUT 200 IN EXPONENT BITS
FADB F, G ;NORMALIZE, RESULTS TO F AND G
FMP F, F ;FORM X^2
MOVE A, E2 ;GET FIRST CONSTANT
FMP A, F ;E2*X^2 IN A
FAD F, E4 ;ADD E4 TO RESULTS IN F
MOVE D, E3 ;PICK UP E3
FDV D, F ;CALCULATE E3/(F^2 + E4)
FSB A, D ;E2*F^2-E3(F^2 + E4)**-1
MOVE E, G ;GET G AGAIN
FSB A, E ;SUBTRACT FROM PARTIAL SUM
FAD A, E1 ;ADD IN E1
FDVM E, A ;DIVIDE BY F
FAD A, E6 ;ADD O.5
FSC A, (H) ;SCALE THE RESULTS
POPJ P,
E1: 204476430062 ;9.95459578
E2: 174433723400 ;0.03465735903
E3: 212464770715 ;617.97226953
E4: 207535527022 ;87.417497202
E5: 270524354513 ;LOG(E), BASE 2
E6: 0.5
E7: 207540074636 ;88.029...
E77: 570232254037 ;-89.415986
IFE BIS,<
ONE: 1.0
ENTRY E.F2D1,E.F2D2,E.F2FP,E.F2F2
E.F2D1:E.F2D2:E.F2FP:E.F2F2:
OUTSTR [ASCIZ /?EXPON is not assembled for double precision floating point.
/]
JRST STOPR.##
END
>
SUBTTL COMP-2 RAISED TO COMP POWER
;THIS ROUTINE CALCULATES A DOUBLE PRECISION NUMBER RAISED
;TO A FIXED POINT POWER. THE CALCULATION IS A**N, WHERE N
;IS AN INTEGER OF THE FORM
; N= Q(0) + Q(1)*2 + Q(2)*4 + ... WHERE Q(I) = 0 OR 1
;THE ONLY RESTRICTION ON THE BASE OR EXPONENT IS THAT
;AN EXPONENT OF 400000000000 IS NOT HANDLED CORRECTLY.
ENTRY E.F2D1 ;COMP-2 TO 1-WORD COMP
ENTRY E.F2D2 ;COMP-2 TO 2-WORD COMP
E.F2D2: SKIPN G,(PA) ;GET HIGH WORD
JRST [SKIPN G,1(PA) ;0, GET LOW WORD
JRST AB.ONE ;BOTH ZERO
JRST E.F2G] ;NO, CONTINUE WITH G SETUP
AOJE G,[SKIPGE G,1(PA) ;-1, GET LOW WORD
JRST E.F2G ;SMALL NEGATIVE NUMBER
JRST .+1] ;DON'T KNOW WHAT TO USE
SETZB A,B ;ERROR, RETURN 0 OR +INFINITY
E.F2D1: SKIPN G,(PA) ;IS EXPONENT 0?
JRST AB.ONE ;YES, A**0 GIVES 1
JUMPE A,ABTEST ;IF BASE = 0
E.F2G: JUMPGE G,[DMOVE D,A ;IS EXPONENT NEGATIVE?
JRST DEX2] ;NO, PUT ARG IN D,D+1 AND START MAIN LOOP
MOVMS G ;GET POSITIVE VALUE
DMOVE D,ONE ;GET DOUB. PRECISION 1.0
DFDV D,A ;CALCULATE (1/X)**N, SINCE N .L. 0
JOV OVR1 ;OVER/UNDERFLOW OCCURED
DEX2: DMOVE A,ONE ;GET DOUB. PREC. 1.0
JRST DEX4 ;START CALCULATING POWERS OF X (OR 1/X)
DEX3: DFMP D,D ;SQUARE X (OR 1/X) AGAIN
JOV OVR ;TEST FOR OVER/UNDERFLOW
LSH G,-1 ;LOOK AT NEXT BIT IN EXPONENT
DEX4: TRZN G,1 ;IS LO BIT ON? (I.E. NEXT POWER OF 2 WANTED)
JRST DEX5 ;NO, DON'T MULTIPLY INTO ANSWER
DFMP A,D ;MULTIPLY POWER OF X INTO ANSWER
JOV OVR ;TEST FOR OVER/UNDERFLOW
DEX5: JUMPN G,DEX3 ;IF G .N. 0, GET MORE POWERS OF X (OR 1/X)
DEX6: POPJ P,
OVR: DMOVE D,A ;PUT CURRENT RESULT IN SAFE PLACE
OVR1: JSP E,.+1 ;GET THE FLAGS
;NOTE IN EXTENDED SECTIONS USE SFM (JRST 14,)
TLNE E,(1B11) ;TEST FOR UNDERFLOW
JRST UNDERF ;YES, SET ANS TO ZERO
PUSHJ PP,ABTEST ;OVERFLOW, RETURN LARGE NUMBER
JUMPGE D,DEX6 ;IF THE ARG IS NEGATIVE
TRNN G,1 ;AND THE EXPONENT IS ODD
DMOVN A,A ;THEN RETURN -INFINITY
POPJ P,
AB.ONE: DMOVE A,ONE ;A**0 GIVES 1
POPJ P,
ABTEST: SKIPL (PA) ;IS BASE 0 WITH POSITIVE EXPONENT?
JRST ABZERO ;YES, RETURN 0
AB.INF: HRLOI A,377777 ;BASE IS 0 WITH NEG. EXP- OVERFLOW
HRLOI B,377777 ;RETURN LARGEST POSITIVE NUMBER
SETOM OVFLO. ;SET OVERFLOW ON
POPJ P,
UNDERF: SETOM OVFLO. ;UNDERFLOW
ABZERO: SETZB A,B ;OUTPUT ANS = 0
POPJ P, ;RETURN
SUBTTL COMP-2 RAISED TO FLOATING POINT POWER
;THIS PROGRAM CALCULATES A**B, WHERE A AND B ARE DOUBLE
;PRECISION FLOATING POINT NUMBERS. THE ALGORITHM USED IS
;A**B = DEXP(B*DLOG(/A/)). THE ABSOLUTE VALUE OF A IS USED
;IN THIS CALCULATION BECAUSE A NEGATIVE NUMBER TO A
;NON-INTEGER POWER PRODUCES A COMPLEX ANSWER, AND B IS
;PRESUMED TO BE NON-INTEGER.
;CHECKS TO SEE IF B IS AN INTEGER AND TREATS IT
;ACCORDINGLY BY CALLING E.F2D1.
;GIVES ERROR FOR NEGATIVE NUMBER TO NON-INTEGER POWER.
ENTRY E.F2FP ;COMP-2 RAISED TO COMP-1 POWER
ENTRY E.F2F2 ;COMP-2 RAISED TO COMP-2 POWER
E.F2FP: MOVE F,(PA) ;GET EXPONENT
TDZA G,G ;ZERO LOW WORD
E.F2F2: DMOVE F,(PA) ;GET EXPONENT
JUMPE F,AB.ONE ;IF EXPONENT = 0 RETURN ANSWER OF 1.0
JUMPE A,ABTEST ;IF BASE = 0
MOVM D,F ;MAGNITUDE OF WORD 1 OF EXP
SETZ C, ;CLEAR AC C
LSHC C,11 ;SHIFT 9 BITS
CAILE C,233 ;INTEGER MUST FIT IN 27 BITS
JRST DEXPD ;NOT AN INTEGER
SUBI C,200 ;200 COMPLEMENT ON EXP
MOVM E,G ;WORD 2 OF EXP
JUMPN E,DEXPD ;IS WORD 2 ZERO?
LSH D,-1 ;COULD STILL BE INTEGER
HRRZ E,C ;SET E AS INDEX REG
SETZ C, ;CLEAR C
ASHC C,(E) ;SHIFT LEFT BY CONTENT S OF E
JOV DEXPD ;OVERFLOW
JUMPN D,DEXPD ;IS IT INTEGER??
SKIPL (PA) ;YES-NEED SIGN OF EX
SKIPA K,C ;GET EXPONENT IN SAFE PLACE
MOVN K,C ;NEGATIVE-NEGATE INTEGER
MOVEI PA,K ;POINT TO EXPONENT
JRST E.F2D1 ;CALL INTEGER ROUTINE
DEXPD: JUMPGE A,DEXPD2 ;NEGATIVE??
DMOVN A,A ;YES-GET ABS VALUE
OUTSTR [ASCIZ /%Attempt to raise negative double precision number to non-integer power
/]
DEXPD2: PUSHJ P,DLOG. ;CALC. LOG(/A/).
DFMP A,(PA) ;CALC. B*LOG(/A/)
JOV OVUNFL
JRST DEXP. ;CALC. EXP[B*LOG(/A/)] AND
;LEAVE IT IN AC'S 0 AND 1.
OVUNFL: JUMPE A,AB.ONE ;IF EXP= 0, ANS = 1.0.
JUMPL A,UNDERF ;GO TO UNDERFLOW.
JRST AB.INF ;OVERFLOW, OUTPUT ANS = INFINITY
SUBTTL DLOG. DOUBLE PRECISION LOGARITHM FUNCTION
;THIS PROGRAM CALCULATES THE LOGARITHM OF A DOUBLE PRECISION
;ARGUMENT. THE ALGORITHM USED IS DESCRIBED ON PAGES 29-30 OF
;RALSTON AND WILF, "MATHEMATICAL METHODS FOR DIGITAL COMPUTERS".
;THE ARGUMENT X IS WRITTEN AS
; X = (2**N)*F WHERE 1/2 < F < 1
;THEN LOG(X) = (N*LOGE(2)) + LOG(F)
;F IS REDUCED BY FIXED POINT MULTIPLICATION BY NOT MORE THAN
;THREE CONSTANTS. THIS YIELDS
; 0 < T = A1*A2*A3*F - 1.0 < (2**-7)/5
;NOTE THAT NOT ALL THE A1,A2,A3 NEED BE INCLUDED IN THE PRODUCT.
;FINALLY,
; LOG(F) = LOG(1+T) - LOG(A1) - LOG(A2) - LOG(A3)
;LOG(1+T) IS CALCULATED AS A TAYLOR SERIES IN T.
;THE ROUTINE IS CALLED IN THE FOLLOWING MANNER:
; DMOVE 0,BASE
; PUSHJ 17,DLOG.
;THE RESULT IS LEFT IN ACCUMULATOR 0 AND 1.
;AVAILABLE TABLES FOR LOGARITHMS LIST 16 SIGNIFICANT DIGITS,
;ALL OF WHICH AGREE WITH VALUES PRODUCED BY THIS PROGRAM. AN
;ESTIMATE OF THE MAXIMUM THEORETICAL ERROR IS OBTAINED BY OB-
;SERVING THAT, AFTER REDUCTION, THE ARGUMENT IS LESS THAN
;(2**-7)/5 = .00625. HENCE, THE ERROR IN THE TAYLOR SERIES
;IS LESS THAN
; (.00625**9)/9 = 1.4573*10**-20/9 = .1619*10**-20
;CONSTANTS AND TEMPORARY LOCATIONS AND STUFF
LOGLST: DOUBLE 200471174064,325425031470 ;0.61180 15411 05992 8976
DOUBLE 200402252251,151350376610 ;0.50455 60107 52395 2859
DOUBLE 177637144373,057714113734 ;0.40546 51081 08164 3810
DOUBLE 177506061360,207057302360 ;0.31845 37311 18534 6147
DOUBLE 176710776761,346515041520 ;0.22314 35513 14209 7553
DOUBLE 176537746034,051711723600 ;0.17185 02569 26659 2214
DOUBLE 175557032242,271265512760 ;0.08961 21586 89687 12374
DOUBLE 173770123303,236031377700 ;0.03077 16586 66753 68689
LIST: 354000000000 ;1.11011 BINARY
324000000000 ;1.10101 BINARY
300000000000 ;1.10000 BINARY
260000000000 ;1.01100 BINARY
240000000000 ;1.01000 BINARY
230000000000 ;1.00110 BINARY
214000000000 ;1.00011 BINARY
204000000000 ;1.00001 BINARY
DLOGE2: DOUBLE 200542710277,276434757153 ;0.69314 71805 59945 30941 72321
A9: DOUBLE 175707070707,034343434344 ;1/9
A8: DOUBLE 601400000000,0 ;-1/8
A7: DOUBLE 176444444444,222222222222 ;1/7
A6: DOUBLE 601252525252,652525252526 ;-1/6
A5: DOUBLE 176631463146,146314631463 ;1/5
A4: DOUBLE 600400000000,0 ;-1/4
A3: DOUBLE 177525252525,125252525253 ;1/3
A2: DOUBLE 577400000000,0 ;-1/2
ONE:
A1: 1.0
DZERO: DOUBLE 0,0
DLOG.: JUMPG A,DLOGP ;ARGUMENT IS GREATER THAN 0
JUMPE A,LGZERO ;CHECK FOR ZERO ARGUMENT
OUTSTR [ASCIZ /%Attempt to take DLOG of Negative Arg
/]
DMOVN A,A ;GET ABSF(X)
DLOGP: CAMN A,ONE ;X PRECISELY 1.0?
JUMPE B,ABZERO ;YES, RETURN ZERO
DLOG1: LDB E,[POINT 8,A,8] ;NO,PICK UP EXPONENT FROM HI ORDER
SUBI E,200 ;GET EXPONENT EXCESS 200
FSC E,233 ;MAKE FLOATING POINT NUMBER
SETZ F, ;SET UP LO HALF
DFMP E,DLOGE2 ;CALCULATE N*LOGE(2)
DMOVE K,E ;SAVE AWAY IN SAFE PLACE
TLZ A,777000 ;MASK OUT EXPONENT
ASHC A,8 ;NORMALIZE FRACTION TO BIT 1
SETZB E,F ;INITIALIZE REDUCTION CONSTANT TO 0
DLOG2: LDB H,[POINT 3,A,4] ;GET HI 3 BITS BELOW 1/2
MOVE C,B ;GET COPY OF LOW HALF
MUL C,LIST(H) ;FIXED POINT MULTIPLY FOR LO HALF
MUL A, LIST(H) ;MULTIPLY HI ORDER, TOO
TLO B,(1B0) ;SET BIT 0, TO AVOID OVERFLOW
ADD B,C ;COMBINE RESULTS OF MULTIPLY
TLZN B,(1B0) ;CLEAR BIT 0, WAS THERE OVERFLOW?
ADDI A,1 ;YES, PROPOGATE CARRY
ASH H,1 ;TURN BITS INTO D.P. POINTER
DFAD E,LOGLST(H)
TLZE A,200000 ;IS THE PRODUCT .GE. 1.0?
JRST DLOG3 ;YES
ASHC A,1 ;NO, GET RID OF EXTRANEOUS ZERO
JRST DLOG2 ;TRY ANOTHER MULTIPLICATION
DLOG3: ASHC A,-7 ;MAKE ROOM FOR THE EXPONENT
TLC A,200000 ;INSERT EXPONENT
DFAD A,DZERO ;NORMALIZE
DMOVN E,E ;NEGATE LOG OF MAGIC NUMBERS
DFAD E,K ;AND ADD INTO FINAL SUMMATION
DMOVEM E,K
DMOVE E,A9 ;PICK UP CONSTANT TO START
MOVEI H,A8 ;INIT TABLE POINTER AT A8
DLOG4: DFMP E,A ;MULTIPLY ACCUMULATED SUM BY X
DFAD E,0(H) ;ADD NEXT CONSTANT INTO PARTIAL SUM
ADDI H,2 ;UPDATE POINTER TO NEXT CONSTANT
CAIG H,A1 ;ARE WE DONE YET?
JRST DLOG4 ;NO, LOOP BACK FOR MORE TAYLOR SERIES
DFMP A,E ;YES, ONE LAST MULTIPLICATION
DFAD A,K ;AND ADD SERIES SUM INTO FINAL ANSWER
POPJ P, ;RETURN
LGZERO: MOVSI A,400000 ;SET UP LARGE NEG. NUM.
DFAD A,A ;CAUSE NEGATIVE OVERFLOW
;TRAP ROUTINE FIXED IT UP AND
;PRINTED MESSAGE
POPJ P,
SUBTTL DEXP. DOUBLE PRECISION EXPONENTIAL FUNCTION
;ARGUMENT. AN ARGUMENT OF ZERO CAUSES AN IMMEDIATE EXIT WITH AN
;ANSWER OF 1.0 . AN ARGUMENT WHOSE MAGNITUDE EXCEEDS 88.028
;CAUSES THE ROUTINE TO EXIT WITH 0 IF THE ARGUMENT WAS NEGATIVE.
;AND 377777777777 IF THE ARGUMENT WAS POSITIVE. THIS IS
;BECAUSE DIRECT CALCULATION OF EXP(X) FOR ABSF(X)>88.028 WOULD
;CAUSE EXPONENT OVERFLOW OR UNDERFLOW.
;THE ROUTINE USES THE FOLLOWING ALGORITHM:
;EXP(X) = 2**(X*LOG2(E))
; = 2**(M+F) WHERE M IS AN INTEGER AND 0<F<1
; = 2**(M+N+R) WHERE 0<R<1/8 AND M+N+R=X*LOG2(E)
; = 2**(M+N) * EXP(R*LOG(2))
;2**M IS CALCULATED EASILY WITH THE FLOATING SCALE INSTRUCTION.
;2**N IS CALCULATED BY DETERMINING THE CORRECT INTERVAL OF N AND
;USING A TABLE OF POWERS OF TWO FROM 2**1/8 TO 2**7/8.
;FINALLY, EXP(R*LOG(2)) IS CALCULATED BY A CONTINUED FRACTION
;TAKEN FROM RALSTON AND WILF, "METHODS FOR DIGITAL COMPUTERS" :
;EXP(R*LOG(2)) = 1+AA4/((B4/R) -C4 + D4*R + H4/(R + B4/R))
;THE FOLLOWING ERRORS HAVE BEEN OBSERVED WITH DEXP:
; 1. WITH -10.0<X<10.0, ERRORS RANGED FROM 0 TO 48
; UNITS IN THE 19TH SIGNIFICANT DIGIT. THE MEAN ERROR
; FOR 20 READINGS WAS 5.4 UNITS IN THE 19TH DIGIT.
; 2. WITH 10.0<X<40.0, ERRORS RANGED FROM 0 TO 41 UNITS
; IN THE 19TH SIGNIFICANT DIGIT. THE MEAN ERROR FOR
; 7 READINGS WAS 16 UNITS IN THE 19TH SIGNIFICANT DIGIT.
; 3. WITH 40.0<X<88.0, ERRORS RANGED FROM 5 TO 57 UNITS IN
; THE 19TH SIGNIFICANT DIGIT. THE MEAN ERROR FOR 12
; READINGS WAS 44.4 UNITS IN THE 19TH SIGNIFICANT DIGIT.
;THE ERRORS REFERRED TO ABOVE ARE ABSOLUTE ERRORS. IT SHOULD
;BE NOTED THAT ADDITIONAL ERRORS ARE INTRODUCED BY ERRORS IN
;THE DOUBLE PRECISION INPUT AND OUTPUT ROUTINES.
;THE ROUTINE IS CALLED IN THE FOLLOWING MANNER:
; DMOVE 0,BASE
; PUSHJ 17,DEXP.
;CONSTANTS AND TEMPORARY LOCATIONS AND STUFF
DCON1: 88.028
DLOG2E: DOUBLE 201561250731,112701376057 ;LOG2(E) = 1.44269 50408 88963 40740
TABLE: 0 ;0
040000000000 ;1/8
100000000000 ;2/8
140000000000 ;3/8
200000000000 ;4/8
240000000000 ;5/8
300000000000 ;6/8
340000000000 ;7/8
POWERS: DOUBLE 1.0,0 ;2**0 = 1.0
DOUBLE 201427127017,037250572672 ;2**1/8 = 1.09050 77326 65257 65919
DOUBLE 201460337602,214333425134 ;2**2/8 = 1.18920 71150 02721 06671
DOUBLE 201513773265,115425047073 ;2**3/8 = 1.29683 95546 51009 66590
DOUBLE 201552023631,237635714441 ;2**4/8 = 1.41421 35623 73095 04878
DOUBLE 201612634520,212520333270 ;2**5/8 = 1.54221 08254 07940 824
DOUBLE 201656423746,126551655275 ;2**6/8 = 1.68179 28305 07429 086
DOUBLE 201725403067,076722207113 ;2**7/8 = 1.83400 80864 09342 463
AA4: DOUBLE 206744575555,062215755376 ;AA4= 60.59319 17173 36463 11080
B4: DOUBLE 207535527021,213670572221 ;B4 = 87.41749 72022 35527 474
MC4: DOUBLE 572033202222,715562022402 ;MC4=-C4 = -30.29659 58586 68231 555
D4: DOUBLE 201414631463,063146314632 ;D4 = 1.05
H4: DOUBLE 210654261010,261565402456 ;H4 = 214.17286 81454 77042 3113
DEXP.: JUMPE A,AB.ONE ;RETURN 1.0 FOR ARG OF ZERO
MOVM C,A ;GET POS VALUE OF EXPONENT
CAML C,DCON1 ;TOO BIG TO COMPUTE?
JRST [JUMPGE A,AB.INF ;YES, RETURN INFINITY
JRST ABZERO] ;OR ZERO IF NEGATIVE
DFMP A,DLOG2E
HLRE E,A ;EXTRACT EXPONENT
ASH E,-9 ;...
TSC E,E ;TAKE 1'S COMPLEMENT IF NUM .L. 0
SKIPGE A ;CHANGE EXPONENT BITS TO SIGN BITS
TLOA A,377000 ;NUMBER NEGATIVE, SET BITS
TLZ A,377000 ;NUMBER POSITIVE, CLEAR BITS
ASHC A,8 ;LEFT JUSTIFY ARG FRACTION BITS
DMOVE C,A ;GET ANOTHER COPY OF FRACTION
ASHC A,-243(E) ;SIMULATE THREE-WORD SHIFT WITH...
;TWO LONG SHIFTS. THIS LEAVES INTEGER
;IN A, FRAC IN B AND C.
LSH D,1 ;SQUEEZE OUT SIGN BIT
IFN C-B-1,<
EXCH B,C-1 ;CROCK BECAUSE B AND C ARE NOT ADJACENT
>
LSHC C,43-200(E) ;THEN DO 2ND LONG SHIFT, (THE LSHC HERE
;PREVENTS OVERFLOW GOING LEFT)
IFN C-B-1,<
EXCH B,C-1 ;PUT ACCS BACK THE WAY WE WANT
>
TLZ B,(1B0) ;CLEAR SIGN BIT. IF FRAC WAS <0,
;THIS GIVES 1-FRAC AND PROPER EXP.
HRRE H,A ;SAVE EXPONENT FOR FUTURE SCALING
MOVEI G,7 ;GET INDEX REGISTER POINTER TO TABLE
REDUCE: CAME B,TABLE(G) ;DOES ARGUMENT MATCH TABLE ENTRY?
JRST REDUC2 ; NO - KEEP LOOKING
JFFO C,.+2 ; HOW MANY LEADING 0'S
JRST REDUC1 ; NONE - C CONTAINS 0
CAIL D,^D28 ; ARE ANY LEFT OF BIT 28 ON?
REDUC1: JRST [LSH G,1 ; NO, IFF LO HALF=0!<377(OCT). CHANGE INDEX TO POINTER
DMOVE A,POWERS(G) ;PICK UP DOUBLE WORD ANSWER
JRST SCALE] ;SCALE RESULTS AND EXIT
REDUC2: CAMGE B,TABLE(G) ;IS ARGUMENT GREATER THAN ENTRY?
SOJA G,REDUCE ;NO, TRY NEXT LOWER ENTRY
SUB B,TABLE(G) ;YES, ALL DONE -REDUCE ARGUMENT
LSH G,1 ;SAVE INDEX AS A POINTER
IFN C-B-1,<
EXCH C,B+1
>
ASHC B,-8 ;MAKE ROOM FOR EXPONENT
TLO B,200000 ;INSERT EXPONENT
DFAD B,DZERO ;NORMALIZE RESULT
DMOVE A,B ;PUT RESULT AN A,B
IFN C-B-1,<
EXCH C,B+1
>
DMOVE C,B4 ;GET B4/F
DFDV C,A
DMOVEM C,K ;SAVE B4/F
DFAD C,A ;GET F+B4/F
DMOVEM C,E ;GET H4/(F+B4/F)
DMOVE C,H4
DFDV C,E
DFMP A,D4 ;GET D4*F
DFAD C,A ;GET (B4/F)-C4+D4*F+(H4/(F+B4/F))
DFAD C,MC4
DFAD C,K
DMOVE A,AA4 ;GET 1.0+AA4/REST
DFDV A,C
DFAD A,ONE
JUMPE G,SCALE ;MULTIPLY BY POWER OF TWO?
DFMP A,POWERS(G)
SCALE: FSC A,(H) ;SCALE RESULTS
POPJ P, ;RETURN
END