;XMATH 06/09/80 ;XYBASIC Interpreter Source Module ;Copyright (C) 1978, 1979, 1980 by Mark Williams Company, Chicago ;floating point extended function package if float and (not f9511) and (not fpbcd) ;8080 FLOATING POINT POLYNOMIAL EXPANDER - FPOLY ; THIS ROUTINE EVALUATES A POLYNOMIAL OF THE FORM ;F(X) = A(0) + A(1)*X^(1*(K+J)) + A(2)*X^(2*(K+J)) ; + A(2)*X^(2*(K+J)) + ... ; WHERE ; A(0) IS AN INITIAL SUM VALUE IN THE ; THE FLOATING POINT ACCUMULATOR ; A(N) IS A TERM CONSTANT FROM A LIST, THE ; ADDRESS OF WHICH IS SUPPLIED IN HL ; AND K AND J ARE INTEGER VALUES SIMULATED BY ; A SUPPLIED VALUE FOR THE INITIAL POWER TERM, AND ; A SUPPLIED VALUE FOR THE POWER TERM TO BE ; MULTIPLIED BY AFTER EACH TERM CALCULATION. ; THE EVALUATION IS TERMINATED BY ONE OF TWO EVENTS. ; 1. IF THE SUPPLIED SIGNIFICANCE STOP VALUE ; IS EVER LESS THAN THE SIGNIFICANCE INDEX ; RETURNED BY THE FLOATING POINT ADD ROUTINE ; (INDICATING THE TERM WOULD NOT AFFECT THE RESULT) ;OR ; 2. A TERM CONSTANT WITH AN EXPONENT OF ZERO ; IS ENCOUNTERED. ; ;DETAILED CALLING SEQUENCE ... ; ; LXI H,INPWR ;ADDRESS OF POWER LIST ; PUSH H ;ON STACK ; LXI H,KLIST ;ADDRESS OF CONSTANT LIST ; CALL FPOLY ; ... ;RETURN HERE WITH F(X) INFACC ; ;INPWR: DS 4 ;INITIAL POWER OF X ; DS 4 ;TERM POWER MULTIPLIER ; ; ... ; ;KLIST: DB -1 ;SIGNIFICANCE STOP VALUE ; DW K,K ;TERM CONSTANT A(1) ; DW K,K ;TERM CONSTANT A(2) ; ... ; DW K,K ;TERM CONSTANT A(N) ; DB 0 ;TERMINATOR ; ; A SIGNIFICANCE STOP VALUE OF -1 WILL CAUSE THE ;ROUTINE TO UNCONDITIONALLY CALCULATE ALL TERMS IN ;THE LIST. A VALUE OF 16 IS USED IN THE SIN/COS ;ROUTINE WITH NO EFFECT ON THE ACCURACY. ; SEE THE DESCRIPTION OF THE INDEX IN THE MATH PACKAGE WRITEUP. FPOLY: LXI D,SIG ;MOVE ARGS OVER MOV A,M ;GET SIG STAX D ;AND SAVE INX H SHLD VECT ;SAVE TERM TABLE ADDRESS POP H ;POP OFF RETURN ADDRESS XTHL ;AND EXCHANGE WITH ARGUMENT MVI C,8 FPOL0: MOV A,M INX H INX D STAX D DCR C JNZ FPOL0 CALL FTEST ;RESTORE LXI H,FPSUM ;AND SAVE PUSH H CALL FSTOR ; SUM LXI H,FPTXN CALL FLOAD ;GET CURRENT POWER FPOL1: LHLD VECT ;AND CURRENT TERM CALL FMUL ;A(N)*X^(N+K) POP H PUSH H CALL FADD ;ADD IT IN PUSH PSW LDA SIG CMP E ;CHECK FOR NO SIGNIFICANCE JC FPOL2 POP PSW POP H PUSH H CALL FSTOR LHLD VECT LXI D,4 DAD D SHLD VECT MOV A,M ANA A POP H JZ FLOAD PUSH H LXI H,FPTXN PUSH H CALL FLOAD LXI H,FPTXN+4 CALL FMUL POP H CALL FSTOR JMP FPOL1 FPOL2: POP PSW POP H RET ;8080 INTEGER/FRACTIONAL PART - FINT ; CALL FINT ;ARGUMENT IN FAC ; ... ;RETURNS WITH SIGNED INTEGER PART IN HL ; AND SIGNED FRACIONAL PART IN FAC ; ;OVERFLOW CONDITIONS SET IF ; -32767 > VALUE > 32767 FINT: CALL FTEST ;SET UP ABCD TO FAC MOV E,A ;SAVE EXP MOV A,B ;GET M0 ANI 80H ;AND ISOLATE SIGN BIT STA FINSN ;SAVE SIGN MOV A,E ;GET EXPONENT CPI 129 ;AND CHECK RANGE JNC FINT0 ;NOT TOO SMALL LXI H,0 ;VALUE .LT. 1.0 RET FINT0: CPI 144 ;SEE IF TOO LARGE JC FINT2 ;VALUE IS IN RANGE LDA FINSN ;OVERFLOW RLC ;SET HL LXI H,7FFFH ;TO +/- JNC $+4 ; FULL SCALE INX H PUSH H CALL FTEST ;RESTORE ABCD POP H ;AND STC ;SET OFLOW AND RET ;EXIT ;SET UP MANTISSA FOR SHIFTING AND CLEAR INTEGER FINT2: XCHG ;D->H PUSH B ;MOV BC POP D ;TO DE SUI 80H ;MAKE SHIFT COUNT MOV C,A ;INTO C MOV B,H ;MOVE M2 MOV A,D ;SET ORI 80H ;MANTISSA MOV D,A ; B23 LXI H,0 ;CLEAR INTEGER ;MANTISSA (M0,M1,M2) IN DEB, INTEGER (I0,I1) IN HL ;SHIFT H L D E B LEFT 'C' BITS ; I0,I1,M0,M1,M2 FINT3: DAD H ;I0,I1 LEFT XCHG DAD H ;M0,M1 LEFT INTO CY JNC $+4 INX D ;CY PROPAGATE INTO I0,I1 XCHG XRA A ;CLEAR A & CY ADD B ;GET M2 RAL ;LEFT B7->CY, 0->B0 MOV B,A ; AND RESTORE JNC $+4 INX D ;CY PROPAGATE INTO M0,M1 DCR C ;DECREMENT EXP JNZ FINT3 ;AND CONTINUE XCHG ;I0,I1 <-> M0',M1' ;CHECK NEW MANTISSA (M0',M1',M2') FOR ZERO ORA H ;CHECK ORA L ;MN' JZ FINT6 ;MN' IS ZERO ;NORMALISE NEW MANTISSA IN HLB ; DECREMENTING EXPONENT IN C FINT4: MOV A,H ;GET M0' RLC ;CHECK HIGH ORDER BIT B23 JC FINT5 ;NORMALIZATION COMPLETE DAD H ;SHIFT M0',M1' LEFT MOV A,B ;GET M2' RAL ;M2' LEFT B7->CY, 0->B0 MOV B,A ;RESTORE JNC $+4 ;NO DATA INX H ;PROPAGATE BIT DCR C ;DECREMENT EXPONENT JMP FINT4 ;AND CONTINUE ;PREPARE FRACTIONAL PART FOR STORAGE FINT5: MOV A,H ;GET M0' ANI 127 ;KILL B23 MOV H,A ;(SIGN) LDA FINSN ;GET SIGN OF INPUT ORA H ;AND MOVE TO B23 MOV H,L ;SWAP FOR MOV L,A ;SHLD LATER MOV A,C ;GET EXPONENT ADI 80H ;AND BIAS UP ;STORE FRACIONAL PART FOR RE-LOAD FINT6: STA FINFP ;STORE EXP SHLD FINFP+1 ;STORE M0',M1' MOV A,B ;GET M2' STA FINFP+3 ;AND STORE ;TRANSFER SIGN TO INTEGER LDA FINSN ;GET SIGN RLC JNC FINT7 ;POSITIVE MOV A,D ;NEGATE CMA ;INTEGER MOV D,A ; BY MOV A,E ; COMPLEMENTING CMA ; AND MOV E,A ; ADDING INX D ; ONE FINT7: PUSH D ;SAVE INTEGER LXI H,FINFP ;RE-LOAD CALL FLOAD ;FRACTION TO FAC POP H ;RESTORE INTEGER RET ;AND EXIT ;FSQR - 8080 FLOATING POINT SQUARE ROOT ;METHOD: APPROXIMATION FOLLOWED BY THREE ; NEWTON ITERATIONS ; ; SQUARE ROOT (X) ; ;LET X=2^(2B)*F WHEN .25<=F<1 ;THEN SQR(X)=2^B*SQR(F) ;WHERE SQR(F)=P(I) I=NUMBER OF ITERATION ; ; P(1)=A*F+B AS A FIRST APPROXIMATION ; WHERE A=.875, B=.27863 WHEN .25<=F<.5 ; OR ; A=.578125, B=.421875 WHEN .5<=F<1 ;AND THEN ; P(I+1)=(P(I)+F/P(I))/2 ; AND P(4) IS FINAL RESULT FOR SQR(F) ; ;CALLING SEQUENCE WITH 8008/8080 MATH PACKAGE ; CALL FSQR ;SQR(FPAC) -> FPAC ; ;NEGATIVE INPUTS WILL BE TREATED AS POSITIVE AFTER CALLING FCERN FSQR: CALL FTEST ;RESTORE RZ ;ZERO CM FCERN ;NEGATIVE ERROR CM FABS ;SET FLAG AND USE ABS VALUE LXI H,FPTXX ;AND SAVE CALL FSTOR ;X MVI E,80H ;CALCULATE RRC ;EXPONENT JNC $+4 ; FOR DCR E ; F MOV A,E ;AND LXI H,FPTF ;STORE CALL FSTOR ; F RRC ;DECIDE LXI H,FPTA1 JNC $+6 ;WHICH LXI H,FPTA2 ; A AND B PUSH H ; TO USE CALL FLOAD ; AND FETCH A LXI H,FPTF ;CALCULATE CALL FMUL ;AF POP H ; AND INX H ; THEN INX H ; GET INX H ; B INX H CALL FADD ;AF+B LXI H,FPTP ;AND PUSH H CALL FSTOR ; SAVE P1 CALL NEWTN ;THREE CALL NEWTN ;NEWTON CALL NEWTN ; ITERATIONS MOV H,A ;SAVE RESULT EXPONENT LDA FPTXX ;CALCULATE SUI 127 ;EXPONENT RAR ; FOR ADD H ; RESULT POP H MOV M,A ;AND SET IT JMP FLOAD ;LOAD RESULT NEWTN: LXI H,FPTF ;GET CALL FLOAD ;F LXI H,FPTP PUSH H CALL FDIV ;F/P POP H PUSH H CALL FADD ;F/P+P DCR A ;(F/P+P)/2 POP H ;AND JMP FSTOR ; SAVE P(I+1) FPTA1: DB 80H,60H,0,0 ;.875 DB 7FH,0EH,0A8H,98H ;.27863 FPTA2: DB 80H,14H,0,0 ;.578125 DB 7FH,58H,0,0 ;.421875 ;********* END - FSQR ********* ;8080 FLOATING POINT EXPONENTIAL - FEXP ; CALL FEXP ;E^(FAC) -> FAC ; ; **** INTERNAL OVER/UNDERFLOW WILL OCCUR FOR INPUT ; VALUES GREATER THAN LN(2^127)/(LOGBASE2(E)). ; OR APPROXIMATELY ; -61 < VALUE < 61 ; ; IF OVERFLOW OCCURS THEN 'ERROR' WILL BE CALLED AND ; FAC WILL BE SET TO FULL SCALE VALUE. ; ;METHOD: ; POLYNOMIAL APPROXIMATION ;TO FIND E^X THE FOLLOWING IDENTITY IS USED. ;TO REDUCE THE RANGE, WE LET ; X LOG2 E = N + D + Z ;WHERE ... ; N IS THE INTEGRAL PORTION OF THE REAL NUMBER, ; D IS A DISCRETE FRACTION (1/8, 3/8, 5/8, OR ; 7/8) OF THE REAL NUMBER (FRAC PART), AND ; Z IS THE REMAINDER WHICH IS IN THE RANGE ; -1/8 <= Z <= 1/8 ; THUS, ; E^X = 2^N * 2^D * 2^Z ;AND IT IS NECCESSARY TO ONLY APPROXIMATE 2^Z FOR ;-1/8 <= Z <= 1/8 BY USING THE POLYNOMIAL F(Z). ; ; F(Z) = A0 + A1*Z + A2*Z^2 + A3*Z^3 + A4*Z^4 + ; A5*Z^5 ;WHERE ... ; A0 = 1.0 ; A1 = .69314718057 ; A2 = .24022648580 ; A3 = .055504105406 ; A4 = .0096217398747 ; A5 = .0013337729375 ; FEXP: CALL FTEST ;SET ABCD TO FAC MOV A,B ;GET SIGN STA FPSGN ;AND SAVE CALL FABS ;TAKE ABS VALUE LXI H,LOG2E ;X=X*LOG2(E) CALL FMUL CALL FINT ;INTEGER TO HL, FRAC IN FAC PUSH H ;SAVE INTEGER LXI H,FPTX ;AND CALL FSTOR ; FRACTION POP H ;RESTORE INTEGER MOV A,H ;TO CHECK RANGE ANA A JZ FEXP1 ;CAN'T BE > 256 FEXP0: LDA FPSGN ;NUMBER TOO LARGE RLC ;SO RETURN JC FZRO ; 0 OR FULL SCALE JMP FOVER ;NONFATAL OV ERROR, TAKE MAX AND RETURN FEXP1: MOV A,L ;CHECK N CPI 127 ;LARGEST THAT WILL FIT JNC FEXP0 ;E^88 = 1.7 E+38 STA FPTN ;SAVE INT PART N ;SELECT VALUE FOR 'D' LDA FPTX ;GET EXPONENT OF FRAC LXI D,4 LXI H,FPTDF CPI 127 JC FEXP5 ;FRAC < .25 DAD D ;FPTDF+4 TO HL JZ FEXP5 ;.25 <= FRAC < .5 DAD D ;FPTDF+8 TO HL LDA FPTX+1 ;GET M0 RLC ;SHIFT B6 RLC ;TO CARRY JNC FEXP5 ;.5 <= FRAC < .75 DAD D ;.75 <= FRAC < 1, FPTDF+12 TO HL FEXP5: PUSH H ;HL POINTS TO SELECTED D LXI H,FPTX ;X TO FAC CALL FLOAD POP H ;ADDRESS OF D PUSH H ;SAVED AGAIN CALL FSUB ;CALCULATE Z MOV E,A ;AND LDA FPSGN ; TRANSFER RLC ; SIGN MOV A,E ; TO Z CC FCHS LXI H,FPTZ0 ;SAVE Z PUSH H ;PUSH FOR FPOLY CALL FSTOR ; IN LXI H,FPTZ1 ; POLY CALL FSTOR ; TABLE LXI H,FPONE ;GET TERM CALL FLOAD ;A0 = 1.0 LXI H,EXPLY ; AND DO THE CALL FPOLY ; POLYNOMIAL F(Z) POP H ;SELECT LXI D,16 ;SIGNED XCHG ; VALUE LDA FPSGN ; FOR RLC ; 2^D JNC $+4 DAD H ;2^-D TABLE DAD D ;HL HAS ADDRESS PUSH H ;OF 2^D CALL FTEST ;RESTORE POP H ;AND CALCULATE CALL FMUL ;2^D * F(Z) ;NOW MULTIPLY BY 2^N MOV E,A ;SAVE EXP LDA FPSGN ;AND RLC ;TRANSFER LDA FPTN ;SIGN JNC $+5 ; TO N CMA INR A ;NEGATE N ADD E ;RESULT = 2^N * 2^D * F(Z) LXI H,FPTZ1 ;STORE TEMP PUSH H CALL FSTOR ;FOR RELOAD POP H ; TO JMP FLOAD ; EXIT WITH FAC = E^X ;D TEST VALUES AND 2^D TABLES FPTDF: DW 7EH,0 ;.125 DW 407FH,0 ;.375 DW 2080H,0 ;.625 DW 6080H,0 ;.875 ;2^D FOR PLUS D VALUES DW 0B81H,0C295H ;2^(1/8) DW 2581H,0D7FEH ;2^(3/8) DW 4581H,2A67H ;2^(5/8) DW 6A81H,0C7C0H ;2^(7/8) ;2^D FOR MINUS D VALUES DW 6A80H,0C7C0H ;2^(-1/8) DW 4580H,2A67H ;2^(-3/8) DW 2580H,0D7FEH ;2^(-5/8) DW 0B80H,0C295H ;2^(-7/8) ;MISC. CONSTANTS LOG2E: DW 3881H,3BAAH ;LOG2(E) 1.442695041 FPMAX: DW 7FFFH,-1 ;POS FULL SCALE EXPLY: DB -1 ;POLY TABLE DW 3180H,1872H ;A1 .69314718057 DW 757EH,0EFFDH ;A2 .24022648580 DW 637CH,4658H ;A3 .055504105406 DW 1D7AH,81A4H ;A4 .0096217398747 DW 2E77H,0FED1H ;A5 .0013337729375 DB 0 ;TABLE TERMINATOR ;****** END - FEXP ****** ;8080 FLOATING POINT NATURAL LOGARITHM - FLN ; THIS ROUTINE CALCULATES THE NATURAL LOGARITHM ;OF THE NUMBER IN THE FLOATING POINT ACCUMULATOR. ; ; LN(FAC) -> FAC ; ;METHOD - POLYNOMIAL APPROXIMATION ; ; GIVEN A NORMALISED REAL NUMBER ; X = 2^K * F ;WHERE THE RANGE OF F IS .5 <= F < 1, AND J AND G ;ARE FOUND SUCH THAT X = 2^J * G WHERE ;(SQR(2)/2) <= G < SQR(2). THIS IS DONE BY SETTING ;J=K-1, G=2*F IF F < SQR(2)/2 AND J=K, G=F OTHER- ;WISE. ; ; THUS: ; LN(X) = J * LN(2) + LN(G) ; ; THE APPROXIMATION FOR LN(G), WHERE ;SQR(2)/2 <= G < SQR(2), IS BASED ON THE SERIES ; LN((V+X)/(V-X) = 2((X/V) + (X^3/(3*V^3)) + ; (X^5/(5*V^5)) + ... ) ;WHICH CONVERGES FOR (-V= SQR(2)/2 DCX H ;POINT TO K DCR M ;J = K-1 INX H ;BACK TO F INR M ;G = 2 * F FLNA: CALL FLOAD ;GET G LXI H,FPONE ;GET 1.0 CALL FADD ;G+1 LXI H,FLNZI ;AND PUSH H CALL FSTOR ; SAVE LXI H,FLNM2 ;GET -2.0 CALL FADD ;G-1 POP H PUSH H CALL FDIV ;(G-1)/(G+1) LXI H,FLNF ;SAVE PUSH H CALL FSTOR ;Z POP H CALL FMUL ;Z^2 POP H PUSH H ;SAVE FOR FPOLY CALL FSTOR ;SAVE LXI H,FLNZ2 ;SAVE CALL FSTOR ;TERM POWER CALL FZRO ;CLEAR SUM LXI H,FLNPL ;AND CALL FPOLY ; CALC POLYNOMIAL LXI H,FLNM2 ;ADD CALL FSUB ;2.0 LXI H,FLNF ;AND PUSH H CALL FMUL ;Z*G(Z) POP H PUSH H CALL FSTOR ;Z*G(Z) ;FLOAT J, CALCULATE J * LN(2) LDA FLNK ;GET J SUI 80H ;UN-BIAS EXP MVI E,32 ;SET INTEGER SCALING MOV D,A ;INTEGER TO I3 LXI B,-1 ;INTEGER MVI A,-1 ;NEGATIVE JM $+5 ;GO FLOAT INX B ;INTEGER INR A ;NEGATIVE CALL FFLOT ;FLOAT J LXI H,FLNL2 CALL FMUL ;J*LN(2) POP H JMP FADD ;Z*G(Z) FLNR2: DW 3580H,0F304H ;SQR(2)/2 FLNL2: DW 3180H,1872H ;LN(2) FLNM2: DW 8082H,0 ;-2.0 FLNPL: DB 16 ;SIG INDEX STOPPER DW 2A80H,0A9AAH ;B2 .666666564181 DW 4C7FH,45CFH ;B4 .400018840613 DW 117FH,0ABAEH ;B6 .28453572660 DW 7EH,0 ;B8 .125 DB 0 ;TERMINATOR ;******* END - FLN ****** ;8080 FLOATING POINT SIN/COS ROUTINE - FSIN, FCOS ; THIS ROUTINE WILL CALCULATE THE SINE OR COSINE ;OF THE ANGLE IN THE FLOATING POINT ACCUMULATOR. ; ; THE ANGLE MUST BE IN RADIANS AND MAY HAVE ANY ;MAGNITUDE. ; ; SIN(FAC) -> FAC ; COS(FAC) -> FAC ; ;METHOD: ; POLYNOMIAL APPROXIMATION ; ; GIVEN A REAL NUMBER, X, N, AND Y ARE DEFINED ;SUCH THAT ; X/(2*PI) = N + Y ;WHERE N IS AN INTEGER AND 0 <= Y < 1. ;THUS, X = 2*PI*N + 2*PI*Y, AND THE IDENTITIES ARE ; SIN(X) = SIN(2*PI*Y) AND COS(X) = COS(2*PI*Y) ; ;THE POLYNOMIAL APPROXIMATION, F(Z), FOR THE ;FUNCTION SIN(2*PI*Z)/Z IS USED WHERE ;-.25 <= Z < .25. ; ; THE PROPERTIES OF SINES AND COSINES ARE USED TO ;COMPUTE THESE FUNCTIONS AS FOLLOWS: ; ; COS(2*PI*Y) = F(Z) ;WHERE ; Z = .25-Y IN THE RANGE 0 <= Y < .5 ; Z = Y-.75 IN THE RANGE .5 <= Y < 1 ; ; SIN(2*PI*Y) = F(Z) ;WHERE ; Z = Y IN THE RANGE 0 <=Y < .25 ; Z = .5-Y IN THE RANGE .25 <= Y < .75 ; Z = Y-1 IN THE RANGE .75 <= Y < 1 ; ; F(Z) = A1*Z + A2*Z^3 + A3*Z^5 + A4*Z^7 + A5*Z^9 ;WHERE ; A1 = 6.2831853 ; A2 = -41.341681 ; A3 = 81.602481 ; A4 = -76.581285 ; A5 = 39.760722 FSIN: DB 3EH ;(MVI A,0AFH) FN CODE FOR SIN FCOS: XRA A ;0 FN CODE FOR COS STA FSCFX ;SAVE FUNCTION CODE CALL FTEST ;RESTORE MOV A,B ;GET M0 ANI 80H ;STRIP OFF STA FSCSG ; AND SAVE CALL FABS ;GET ABS(FAC) LXI H,F2PI CALL FDIV ;X/(2*PI) CPI 129 CNC FINT ;Y >= 1, GET FRACTIONAL PART MOV E,A ;SAVE EXPONENT LDA FSCFX ;CHECK RLC ;FUNCTION CODE MOV A,E ;RESTORE EXPONENT JC FSC3 ;FUNCTION IS SIN CPI 80H JC FSC2 ;Y < .5 LXI H,FSC75 ;Y - .75 JMP FSC6 FSC2: CALL FCHS ;-Y LXI H,FSC25 ;.25-Y JMP FSC4A FSC3: CPI 127 ;Y < .25 ? JC FSC7 ;YES, GO DO POLY CPI 80H JC FSC4 ;Y < .5 MOV E,A ;SAVE EXP MOV A,B ;CHECK RLC ;B22 RLC ; FOR MOV A,E ; .75, RESTORE EXPONENT JC FSC5 ;Y >= .75 FSC4: CALL FCHS ;Y < .75 SO -Y LXI H,FSC50 FSC4A: CALL FADD ;.5-Y JMP FSC7 FSC5: LXI H,FPONE FSC6: CALL FSUB ;Y-1 FSC7: LXI H,FSCZ ;SAVE PUSH H CALL FSTOR ;Z POP H PUSH H ;SAVE FOR FPOLY CALL FMUL ;Z^2 LXI H,FSCZ2 CALL FSTOR ;SAVE Z^2 CALL FZRO ;ZERO SUM LXI H,FSCPL ;AND CALL FPOLY ; DO POLY MOV E,A ;SAVE EXP LDA FSCFX ;AND CHECK RLC ; FUNCTION MOV A,E ;RESTORE EXPONENT RNC ;AND EXIT IF COS FUNCTION LDA FSCSG ;GET INPUT SIGN FSC8: RLC ;TO CY -- FATN ENTRY POINT MOV A,E ;RESTORE EXPONENT RNC ;RETURN IF INPUT + JMP FCHS ;ELSE CHANGE RESULT SIGN FPI: DW 4982H,0DB0FH ;PI FHPI: DW 4981H,0DB0FH ;PI/2 FSC75: DW 4080H,0 ;.75 FSC25: DW 7FH,0 ;.25 FSC50: DW 80H,0 ;.5 FSCPL: DB 16 ;SIN/COS POLY SIG STOPPER F2PI: DW 4983H,0DB0FH ;A1 PI*2 DW 0A586H,0E25DH ;A2 -41.341681 DW 2387H,7834H ;A3 81.602481 DW 9987H,9E29H ;A4 -76.581285 DW 1F86H,0FB0AH ;A5 39.760722 DB 0 ;TERMINATOR ;****** END - FSIN/FCOS ****** ;8080 FLOATING POINT REAL BASE TO REAL EXPONENT ; FAX ; THIS ROUTINE RAISES A REAL NUMBER TO A REAL POWER. ; ***** NOTE ***** ; BECAUSE OF OVERFLOW CONDITIONS IN THE ROUTINE ;FLN, THE FOLLOWING LIMITS ARE PLACED ON THE ;INPUTS TO FAX: ; ; A > 0 ;AND ; ABS(X*LN(A)) < 50.0 ; ;CALLING SEQUENCE ... ; LXI H,X ;ADDRESS OF X ; CALL FAX ;FAC ^ X -> FAC ; ;METHOD ;If X is a positive integer < 256, iterated multiplication is used. ;Otherwise the identity A^X = E^(X*LN(A)) is used. FATOX: XCHG INX H ;ARG2 POINTER TO HL FAX: PUSH H call ftest ;load registers lxi h,ftmpa call fstor ;save A in FTMPA pop h push h call fload ;get X to FACC call fint ;check if X is integer push psw push h ;save integer part lxi h,ftmpa call fload ;restore A to FACC pop d ;integer part to DE pop psw pop h ;X pointer to HL ora d ;Zero set iff A positive integer < 256 jnz fax2 ;raise to real power inr e push d ;save the integer power+1 lxi h,fpone ;initialize result to 1.0 call fload fax1: pop d dcr e ;decrement multiply count jz ftest ;done, set registers and return push d lxi h,ftmpa call fmul ;multiply result by A jmp fax1 ;and continue multiplying fax2: push h ;raising to real, large or negative power CALL FLN ;LN(A) POP H CALL FMUL ;X*LN(A) JMP FEXP ;E^(X*LN(A)) if false ;8080 FLOATING POINT COMMON LOGARITHM - FLOG ; THIS ROUTINE CALCULATES THE COMMON OR BASE 10 ;LOGARITHM OF THE NUMBER IN THE FLOATING POINT ;ACCUMULATOR. ; ; LOG(FAC) -> FAC ; ;METHOD ; THE FOLLOWING IDENTITY IS USED. ; ; LOG(X) = LN(X)/LN(10) FLOG: CALL FLN ;LN(ABS(X)) LXI H,FLGTI ;1/LN(10) JMP FMUL ;LN(X)/LN(10) ;8080 FLOATING POINT ANTILOG - FALOG ; ; THIS ROUTINE CALCULATES THE BASE 10 ANTILOG OF ;THE NUMBER IN THE FLOATING POINT ACCUMULATOR. ; ; ALOG(FAC) -> FAC ; ;METHOD ; THE FOLLOWING IDENTITY IS USED ; 10^X = E^(X*LN(10)) FALOG: CALL FTEST ;RESTORE LXI H,FLOGT ;LN(10) CALL FMUL ;X*LN(10) JMP FEXP ;E^(X*LN(10)) ;8080 DEGREE>RADIAN AND RADIAN>DEGREE - FCDR, FCRD ; THESE ROUTINES CONVERT THE NUMBER IN THE FLOATING ;POINT ACCUMULATOR FROM DEGREES TO RADIANS (FCDR) ;OR FROM RADIANS TO DEGREES (FCRD). ; ; FAC;DEGREES -> FAC;RADIANS ; FAC;RADIANS -> FAC;DEGREES ; ;METHOD ; THE FOLLOWING IDENTITIES ARE USED. ; ; X;RADIANS = X;DEGREES * PI/180 ; X;DEGREES = X;RADIANS * 180/PI FCDR: CALL FTEST ;RESTORE LXI H,FXDR JMP FMUL ;X*PI/180 FCRD: CALL FTEST ;RESTORE LXI H,FXRD JMP FMUL ;X*180/PI endif ;8080 FLOATING POINT TANGENT - FTAN ; THIS ROUTINE CALCULATES THE TANGENT OF THE ANGLE ;IN THE FLOATING POINT ACCUMULATOR. ; THE ANGLE MUST BE EXPRESSED IN RADIANS. ; ; TAN(FAC) -> FAC ; ;METHOD ; THE FOLLOWING IDENTITY IS USED. ; ; TAN(X) = SIN(X)/COS(X) ; ; ***** NOTE ***** ; FOR VALUES OF X CLOSE TO 0 OR CLOSE TO A MULTIPLE ;OF 2*PI, EITHER THE SIN OR COS OF X WILL BE CLOSE ;TO UNITY. THIS CONDITION WILL AFFECT THE ACCURACY ;USING THE IDENTITY ABOVE AND IT IS SUGGESTED ;THAT OTHER METHODS BE USED TO COMPUTE THE FUNCTION ;TAN(X) IF THIS INACCURACY CAN NOT BE TOLERATED. FTAN: CALL FTEST ;RESTORE RZ LXI H,FTMPA PUSH H CALL FSTOR ;SAVE X CALL FCOS ;COS(X) LXI H,FTMPB CALL FSTOR ;SAVE COS(X) POP H CALL FLOAD ;GET X CALL FSIN ;SIN(X) LXI H,FTMPB JMP FDIV ;SIN(X)/COS(X) if false FLGTI: DW 5E7FH,0D95BH ;1/LN(10) FLOGT: DW 1382H,8E5DH ;LN(10) FXDR: DW 0E7BH,35FAH ;PI/180 FXRD: DW 6586H,0E12EH ;180/PI endif ;****** END -FAX, FLOG, FALOG, FCDR, FCRD, FTAN ***** if not camac ;8080 FLOATING POINT ARCTANGENT - FATAN ; THIS ROUTINE CALCULATES THE ARCTANGENT OF THE ;NUMBER IN THE FLOATING POINT ACCUMULATOR. ; ; ATAN(FAC) -> FAC ; ;METHOD ; POLYNOMIAL APPROXIMATION ; ; THIS ROUTINE IS BUILT AROUND A POLYNOMIAL, F(X), ;THAT APPROXIMATES ATAN(Z) IN THE RANGE ;-.25 <= Z <= .25. THE ATAN(Z) FOR VALUES OF Z ;OUTSIDE THIS RANGE IS FOUND BY USING THE FOLLOWING ;IDENTITIES: ; ATAN(-Z) = -ATAN(Z) ;AND ; ATAN(Z) = A(K) + ATAN((Z-B(K))/(Z*B(K)+1)) ;WHERE ; A(K) = K*PI/7, B(K) = TAN(A(K)) ;AND K IS DETERMINED SO THAT ;TAN((2*K-1)*PI/14) <= ABS(Z) <= TAN((2*K+1)*PI/14) ;WITH K = 1, 2, OR 3. ; ; HAVING DETERMINED THE VALUE OF K APPROPRIATE ;TO Z, THE TRANSFORMATION ; X = (Z-B(K))/(Z*B(K)+1) ;PUTS X IN THE RANGE -TAN(PI/14) < X < TAN(PI/14). ;THE POLYNOMIAL WORKS OVER A SLIGHTLY LARGER RANGE ;TO ALLOW DETERMINATION OF K USING THE VALUES ; K=0 IF ABS(Z) < .25 ; K=1 IF .25 < ABS(Z) < .75 ; K=2 IF .75 < ABS(Z) < 2 ; K=3 IF ABS(Z) > 2 ;THEN ; ATAN(Z) = A(K) + F(X) IF Z >= 0 ; ATAN(Z) = -A(K) - F(X) IF Z < 0 ;F(X) = X*(1 - Q1*X^2 + Q2*X^4 - Q3*X^6) ;WHERE ; Q1 = .333329573 ; Q2 = .199641035 ; Q3 = .131779888 FATAN: CALL FTEST ;RESTORE RZ MOV E,A ;SAVE EXPONENT MOV A,B ;GET M0 ANI 80H ;ISOLATE SIGN, B23 STA FATSG ;AND SAVE CALL FABS MOV E,A ;SAVE EXP LXI H,FATK MVI M,0 ;K=0 CPI 127 JZ FAT4 ;.25 <= Z < .5, SO K=1 JC FAT7 ;Z < .25, SO K=0 AND DO F(X), X=Z MVI M,2 CPI 129 JZ FAT6 ;.75 <= Z < 1, K=2 JNC FAT4 ;Z > 2, K=3 MOV A,B RLC RLC ;CHECK 1/4 BIT DCR M JNC FAT6 ;.5 <= Z <.75, K=1 FAT4: INR M FAT6: MOV A,M ;FETCH K PUSH PSW ;SAVE K MOV A,E ;GET EXPONENT LXI H,FATZX ;SAVE CALL FSTOR ;Z POP PSW ;GET K DCR A ;K-1 ADD A ;(K-1)*2 ADD A ;(K-1)*4 ADD A ;(K-1)*8 LXI H,FATC1 ;START OF TABLE MOV E,A ;CALCULATE MVI D,0 ;ADDRESS DAD D ; OF A(K) PUSH H ;AND SAVE IT LXI D,4 DAD D ;ADDRESS OF B(K) PUSH H ;SAVED ALSO CALL FLOAD ;GET B(K) LXI H,FATZX PUSH H CALL FMUL ;Z*B(K) LXI H,FPONE CALL FADD ;Z*B(K)+1 LXI H,FATXM CALL FSTOR ;SAVE Z*B(K)+1 POP H CALL FLOAD ;GET Z POP H ;ADDRESS OF B(K) CALL FSUB ;Z-B(K) LXI H,FATXM CALL FDIV ;(Z-B(K))/(Z*B(K)+1) FAT7: LXI H,FATZX PUSH H CALL FSTOR ;SAVE X POP H CALL FMUL ;X^2 LXI H,FATXM PUSH H ;SAVE FOR FPOLY CALL FSTOR ;SAVE INITIAL POWER LXI H,FATX2 ;SAVE TERM POWER MULTIPLIER CALL FSTOR LXI H,FPONE CALL FLOAD ;INITIALISE SUM LXI H,FATPL CALL FPOLY ;F(X) LXI H,FATZX CALL FMUL ;X*F(X) MOV E,A ;SAVE EXPONENT LDA FATK ;GET K ANA A ;SEE IF 0 MOV A,E ;RESTORE EXPONENT JZ FAT8 ;K=0, NO A(K) POP H ;ADDRESS OF A(K) CALL FADD ;A(K)+F(X) FAT8: MOV E,A ;SAVE EXP LDA FATSG ;GET INPUT SIGN JMP FSC8 ;EXIT AS FOR SIN AND COS FATC1: DW 657FH,0FAC8H ;A(1) PI/7 .4487989506 DW 767FH,0F390H ;B(1) TAN(A(1)) .4815746188 DW 6580H,0FAC8H ;A(2) 2*PI/7 .8975979011 DW 2081H,0C681H ;B(2) TAN(A(2)) 1.253960337 DW 2C81H,0BB56H ;A(3) 3*PI/7 1.346396852 DW 0C83H,7F33H ;B(3) TAN(A(3)) 4.381286272 FATPL: DB -1 ;ATAN POLY TABLE DW 0AA7FH,2DAAH ;Q1 -.333329573 DW 4C7EH,0B36EH ;Q2 .199641035 DW 867EH,4FF1H ;Q3 -.131779888 DB 0 ;TERMINATOR ;****** END - FATAN ****** endif ;end of NOT CAMAC contditional if false ;8080 FLOATING POINT INVERSE SIN/COS - FASIN, FACOS ; THESE ROUTINES CALCULATE THE INVERSE ;SINE AND COSINE OF A NUMBER IN THE FLOATING ;POINT ACCUMULATOR. ; THE RESULTANT ANGLE IS IN RADIANS. ; ; ASIN(FAC) -> FAC ; ACOS(FAC) -> FAC ; ;METHOD ; THE FOLLOWING IDENTITIES ARE USED ; ; ASIN(X) = ATAN(X/SQR(1-X^2)) ; ACOS(X) = ATAN(SQR(1-X^2)/X) ;AND ; ASIN(-X) = -ASIN(X) ; ACOS(-X) = ACOS(X) + PI FASIN: CALL FTEST ;RESTORE CPI 129 JC FAS4 ;ABS(X) < 1 JZ FAS1 ;ABS(X) MAY = 1 FAS0: CALL FTEST JMP FCERN FAS1: MOV E,A ;SAVE EXP MOV A,B ;CHECK ANI 127 ;FOR JZ FAS3 ; UNITY FAS2: MOV A,E ;RESTORE EXP JMP FAS0 ;ERROR ABS(X) > 1 FAS3: ADD C JNZ FAS2 ;> 1 ADD D ;CHECK M3 JNZ FAS2 ;> 1 MOV A,B ;X = +- 1 FSETP: PUSH PSW ;SAVE SIGN LXI H,FHPI CALL FLOAD ;GET PI/2 MOV E,A ;SAVE EXP POP PSW ;GET SIGN JMP FSC8 ;EXIT AS FOR SIN AND COS FAS4: CALL FR1X2 ;SQR(1-X^2) LXI H,FASCT CALL FSTOR ;SAVE LXI H,FASCX CALL FLOAD ;GET X LXI H,FASCT CALL FDIV ;X/SQR(1-X^2) JMP FATAN ;ATAN(X/SQR(1-X^2)) FACOS: CALL FTEST ;RESTORE ANA A ;SET FLAGS JZ FSETP ;X=0, SET PI/2 CPI 129 JC FAC1 ;ABS(X) < 1 CONTINUE MOV E,A ;SAVE EXP FAC0: MOV A,E ;RESTORE EXPONENT STC ;SET ERROR FLAG JNZ FCERN ;AND RET IF ABS(X) > 1 MOV A,B ;CHECK ANI 127 ;FOR JNZ FAC0 ; X ADD C ; = JNZ FAC0 ; 1 ADD D JNZ FAC0 ;> 1 JMP FZRO ;X=1 SO RETURN 0 FAC1: CALL FR1X2 ;SQR(1-X^2) LXI H,FASCX CALL FDIV ;SQR(1-X^2)/X CALL FATAN ;ATAN(SQR(-1X^2)) RP ;LEAVE IF + LXI H,FPI JMP FADD ;ADD PI ;ROUTINE TO CALCULATE SQR(1-X^2) FR1X2: LXI H,FASCX CALL FSTOR ;SAVE X LXI H,FASCX CALL FMUL ;X^2 CALL FCHS ;-X^2 LXI H,FPONE CALL FADD ;1-X^2 JMP FSQR ;SQR(1-X^2) ;****** END - FASIN, FACOS ****** endif endif ;end of FLOAT conditional ;****** END - XMATH PACKAGE ****** ;end of XMATH PAGE