TITLE LIBFUN - PPL LIBRARY FUNCTIONS SUBTTL E.A.TAFT/ 9-SEP-71 V.2 HISEG SEARCH PPL ;THE SYSTEM FUNCTIONS IN THIS SECTION ARE MATHEMATICAL FUNCTIONS ;WHICH TAKE ARGS OF TYPE REAL AND RETURN REALS. AN INT OR DBL IS CONVERTED ;TO A REAL FOR CONVENIENCE ;ROUTINE TO PREPARE A SINGLE ARGUMENT FOR USE BY LIBRARY FUNCTIONS LIBPRP: CALL GRAB1 ;PICK UP AND CHECK FOR SINGLE ATOMIC ARG CAIN T,U.INT/2 ;AN INTEGER? CALL INTRL1 ;YES, CONVERT TO A REAL CAIN T,U.DBL/2 ;A DBL? CALL DBLRL1 ;YES. CONV. TO REAL CAIE T,U.REAL/2 ;A REAL? ILLTYPE ;NO, IMPROPER CALLING ARG RETURN ;OK, RETURN TO CALLER ;THE ALGORITHMS FOR THE SQRT,EXP,LN,SIN,COS, AND ATAN SYSTEM ;FUNCTIONS ARE TAKEN FROM THE DEC-SUPPLIED FORTRAN LIBRARY (LIB40). ;EXPLANATIONS FOR THE ALGORITHMS ARE CONTAINED IN THE MANUAL ;'SCIENCE LIBRARY AND FORTRAN UTILITY SUBPROGRAMS'. ;SQRT(A) - RETURN THE SQUARE ROOT OF A. ;A MUST BE NON-NEGATIVE SSQRT: EXP 1 ;TAKES ONE ARG CALL LIBPRP ;MAKE SURE ARG IS A REAL JUMPGE R,.+2 ;LEGAL ARG? SFNERR MSG(SQRNG) ;SQUARE ROOT OF NEGATIVE NUMBER JUMPE R,SRETRN ;SQRT(0)=0 IDIV R,[XWD 1000,0] ;SEPARATE EXP AND FRACTION HRREI AC1,-201(R) ;SAVE TRUE EXPONENT -1 ROT AC1,-1 ;DIVIDE EXP BY 2 SAVING REMAINDER JUMPL AC1,SQRT3 ;JUMP IF FRACTION > .5 ;HERE FOR RESULTING FRACTION < .5 FSC R2,177 ;FIX UP EXPONENT SO NUMBER IN RANGE .25 TO .5 MOVE R,R2 ;SAVE RESULT FMPRI R2,200640 ;TAKE FIRST LINEAR APPROXIMATION FADRI R2,177465 JRST SQRT1 ;GO PERFORM NEWTON'S METHOD ;HERE FOR FRACTION > .5 SQRT3: FSC R2,200 ;FIX UP EXPONENT SO NUMBER IN RANGE .5<= R2 < 1 MOVE R,R2 ;SAVE RESULT FMPRI R2,200450 ;COMPUTE LINEAR APPROXIMATION FADRI R2,177660 ;HERE TO PERFORM TWO ITERATIONS OF NEWTON'S METHOD SQRT1: SAVE R ;REMEMBER FRACTION FDV R,R2 ;COMPUTE .5*(APPROX+F/APPROX) FAD R2,R FSC R2,-1 RESTORE R ;GET BACK FRACTION FDV R,R2 ;COMPUTE APPROX+F/APPROX FADR R,R2 FSC R,(AC1) ;HALVE AND SCALE RESULT JRST SRETRN ;RETURN REAL RESULT ;EXP(A) - RETURNS EXPONENTIAL OF A SEXP: EXP 1 ;TAKES ONE ARG CALL LIBPRP ;ENSURE A REAL ARG EXP1: CAML R,[570232254037];CHECK FOR WITHIN RANGE -89.410 TRNE AC3,377 ;MAKE CORRECTION TO SCALED ARG ADDI AC3,200 ASHH: ASH AC3,-^D8 ;MAKE ROOM FOR EXPONENT IN SCALED ARG TLC AC3,200000 ;PUT 200 IN EXPONENT BITS FSC AC3,0 ;NORMALIZE SCALED ARG, CALL IT X MOVE AC1,AC3 ;SAVE X FMP AC3,AC3 ;TAKE X^2 MOVE R,AC3 ;MOVE X^2 TO R FMP R,[174433723400];PERFORM POWER SERIES EXPANSION FAD AC3,[207535527022] MOVE AC4,[212464770715] FDV AC4,AC3 FSB R,AC4 FSB R,AC1 FAD R,[204476430062] FDVM AC1,R FADRI R,(0.5) FSC R,(AC2) ;SCALE RESULTS JRST SRETRN ;RETURN REAL RESULT ;LN(A) - RETURN THE NATURAL LOGARITHM OF A SLN: EXP 1 ;TAKES ONE ARG CALL LIBPRP ;ENSURE REAL ARG CALL LNSUB ;CALL LN AS SUBROUTINE SO POWER CAN USE JRST SRETRN ;RETURN REAL RESULT LNSUB: JUMPG R,.+2 ;ARG MUST BE POSITIVE AND NONZERO SFNERR MSG(ARGNG) ;ARG MUST BE POSITIVE AND NONZERO CAMN R,[1.0] ;LN(1)? JRST RETZER ;YES. RETURN ZERO RESULT IDIV R,[XWD 1000,0] ;SEPARATE EXPONENT AND FRACTION ADDI R,211000 ;FLOAT THE EXPONENT AND MULT BY 2 MOVS AC1,R ;MOVE FLOATED EXPONENT TO AC1 FSBRI AC1,210401 ;SUBTRACT FLOATING 401 OCTAL FROM IT TLC R2,200000 ;FLOAT THE FRACTION PART FAD R2,[577225754146] ;SUBTRACT SQRT(2)/2 MOVE AC2,R2 ;SAVE RESULTS FAD AC2,[201552023632] ;ADD SQRT(2) FDVB R2,AC2 ;SAVE REDUCED ARG, CALL IT Z FMP R2,R2 ;CALCULATE Z^2 MOVE R,[200462532521] ;PICK UP FIRST CONSTANT FMP R,R2 ;MULT BY Z^2 FAD R,[200754213604] ;ADD IN NEXT CONSTANT FMP R,R2 ;MULT BY Z^2 FAD R,[202561251002] ;ADD IN NEXT CONSTANT FMP R,AC2 ;MULT BY Z FAD R,AC1 ;ADD IN EXPONENT TO FORM LOG2(X) FMP R,[200542710300] ;CONVERT TO LOGE(X) RETURN ;HERE TO RETURN ZERO FOR LN(1.0) RETZER: SETZ R, RETURN ;COS(A) - RETURNS COSINE OF A SCOS: EXP 1 ;TAKES ONE ARG CALL LIBPRP ;ENSURE REAL ARG FADR R,[201622077325] ;COS(X)=SIN(X+PI/2) JRST SIN1 ;ENTER SINE ROUTINE ;SIN(A) - RETURNS SINE OF A SSIN: EXP 1 ;TAKES ONE ARG CALL LIBPRP ;ENSURE REAL ARG SIN1: MOVM AC1,R ;TAKE MAGNITUDE OF ARG CAMG AC1,[170000000000] ;COMPARE TO 2^-9 JRST SRETRN ;SIN(X)=X IF X IS VERY SMALL FDV AC1,[201622077325] ;DIVIDE X BY PI/2 CAMG AC1,[1.0] ;IS X/(PI/2) < 1.0? JRST SIN2 ;YES, ALREADY IN FIRST QUADRANT MULI AC1,400 ;NO. SEPARATE EXP FROM FRACTION ASH AC2,-202(AC1) ;GET X MOD 2 PI MOVEI AC1,200 ;PREPARE FLOATING EXPONENT ROT AC2,3 ;SAVE THREE BITS TO DETERMINE QUADRANT LSHC AC1,^D27 ;ARG NOW IN RANGE (-1,1) FSC AC1,0 ;NORMALIZE THE ARG JUMPE AC2,SIN2 ;ALREADY IN FIRST QUADRANT IF BITS 000 TLCE AC2,1000 ;SUBTRACT 1.0 IF BITS ARE 001 OR 011 FSBRI AC1,(1.0) TLCE AC2,3000 ;CHECK FOR FIRST QUADRANT, 001 TLNN AC2,3000 ;CHECK FOR THIRD QUADRANT, 010 MOVN AC1,AC1 SIN2: JUMPGE R,.+2 ;CHECK SIGN OF ORIGINAL ARG MOVN AC1,AC1 ;SIN(-X)=-SIN(X) MOVE R,AC1 ;SAVE REDUCED ARG FMPR AC1,AC1 ;COMPTE X^2 MOVE AC2,[164475536722] ;GET 1ST CONSTANT FMP AC2,AC1 ;MULTIPLY BY X^2 FAD AC2,[606315546346] ;ADD NEXT CONSTANT FMP AC2,AC1 ;MULTIPLY BY X^2 FAD AC2,[175506321276] ;ADD NEXT CONSTANT FMP AC2,AC1 ;MULTIPLY BY X^2 FAD AC2,[577265210372] ;ADD PI/2 FMP AC2,AC1 ;MULTIPLY BY X^2 FAD AC2,[201622077325] ;ADD NEXT CONSTANT FMPR R,AC2 ;MULTIPLY BY X JRST SRETRN ;RETURN REAL RESULT ;ATAN(A) - RETURNS ARCTANGENT OF A SATAN: EXP 1 ;TAKES ONE ARG CALL LIBPRP ;PREPARE SINGLE REAL ARG MOVM R2,R ;GET MAGNITUDE CAMG R2,[145000000000] ;COMPARE TO 2^-33 JRST SRETRN ;ATAN(X)=X IF X IS VERY SMALL HLLO AC2,R ;SAVE SIGN OF ARG. RH USED AS FLAG CAML R2,[233000000000] ;IF ARG>2^33 THEN RETURN WITH JRST AT4 ;ATAN(X)=PI/2 MOVSI AC1,(1.0) ;FORM 1.0 CAMG R2,AC1 ;IS ABS(X)>1? TRZA AC2,-1 ;NO. SWAP FLAG FDVM AC1,R2 ;YES. REPLACE X BY 1/X TLC AC2,(AC2) ;SWAP SIGN IF NECESSARY MOVE AC3,R2 ;SAVE ARG FMP R2,R2 ;TAKE X^2 MOVE AC1,[201562663021] ;GET FIRST CONSTANT FAD AC1,R2 ;MULT BY X^2 MOVE R,[600360700773] ;GET NEXT CONSTANT FDVM R,AC1 ;FORM -A3/(X^2+B3) FAD AC1,R2 ;ADD X^2 TO PARTIAL SUM FAD AC1,[202650373270] ;ADD NEXT CONSTANT MOVE R,[574071125540] ;GET NEXT CONSTANT FDVM R,AC1 ;DIVIDE BY PARTIAL SUM FAD AC1,R2 ;ADD X^2 TO PARTIAL SUM FAD AC1,[203660615617] ;ADD NEXT CONSTANT MOVE R,[202732621643] ;GET NEXT CONSTANT FDV R,AC1 ;DIVIDE BY PARTIAL SUM FAD R,[176545543401] ;ADD NEXT CONSTANT FMP R,AC3 ;ADD SAVED X TRNE AC2,-1 ;CHECK '>1.0' INDICATOR FSB R,[201622077325] ;ATAN(A)=-ATAN(1/A)-PI/2) JRST AT5 ;HERE TO RETURN PI/2 AS ANSWER AT4: MOVE R,[201622077325] AT5: JUMPGE AC2,.+2 ;CHECK SIGN FLAG MOVN R,R JRST SRETRN ;RETURN RESULT ;RANDOM(A) - RETURNS A RANDOM REAL BETWEEN 0 AND 1. ; IF A=0, THE NEXT RANDOM NUMBER IN SEQUENCE IS RETURNED ; IF A<0, THE RANDOM NUMBER SEQUENCE IS INITIALIZED WITH THE ; TIME OF DAY. ; IF A>0, THE RANDOM NUMBER SEQUENCE IS INITIALIZED WITH A. SRANDO: EXP 1 ;TAKES ONE ARG CALL LIBPRP ;CONVERT ARG TO A REAL JUMPE R,NXTRAN ;IF ZERO, JUST GET NEXT RANDOM NUMBER JUMPG R,.+2 ;USER SUPPLYING OWN ARG? MSTIME R, ;NO, INITIALIZE WITH TIME OF DAY TLZ R,760000 ;MASK 5 BITS FOR SAFETY MOVEM R,RANDOM ;STORE PRESENT RANDOM NUMBER ;HERE TO GENERATE NEXT RANDOM NUMBER NXTRAN: MOVE R,[^D630360016] ;SETUP MAGIC CONSTANT MUL R,RANDOM ;MULTIPLY BY PREVIOUS RANDOM NUMBER ASHC R,4 ;SEPARATE RESULT INTO 2 31-BIT WORDS LSH R2,-4 ADD R,R2 ;ADD RESULTS TOGETHER TLZE R,760000 ;PERFORM END-AROUND CARRY MOD 31 BITS ADDI R,1 MOVEM R,RANDOM ;STORE NEW RANDOM INTEGER HLRZ R2,R ;CONVERT INTEGER TO FLOATING POINT HLLI R, ; BETWEEN ZERO AND ONE. FSC R2,216 FSC R,174 FADR R,R2 ;COMBINE HIGH AND LOW COMPONENTS JRST SRETRN ;RETURN RANDOM REAL END