;;; -*-MIDAS-*- ;;; ************************************************************** ;;; ***** MACLISP ****** STANDARD ARITHMETIC FUNCTIONS *********** ;;; ************************************************************** ;;; ** (C) COPYRIGHT 1981 MASSACHUSETTS INSTITUTE OF TECHNOLOGY ** ;;; ****** THIS IS A READ-ONLY FILE! (ALL WRITES RESERVED) ******* ;;; ************************************************************** PGBOT ARI ;THE ARITHMETIC PAGE - ARITHMETIC SUBROUTINES IFN BIGNUM,[ SUBTTL ARITHMETIC FUNCTIONS WITH BIGNUM==1 ZEROP: MOVEI R,2 JRST ZMP MINUSP: TDZA R,R PLUSP: MOVEI R,1 ZMP: JSP T,NVSKIP JRST .+2 JFCL XCT .+2(R) JRST FALSE JUMPL TT,TRUE ;FOR MINUSP JUMPG TT,TRUE ;FOR PLUSP JUMPE TT,TRUE ;FOR ZEROP MINUS: JSP T,NVSKIP JRST MNSBG JRST MNSFX MOVNS TT JRST FLOAT1 MNSFX: CAMN TT,[400000000000] JRST ABSOV MOVNS TT JRST FIX1 ADD1: MOVEI R,1 JRST SUB11 SUB1: MOVNI R,1 SUB11: JSP T,NVSKIP JRST A1S1BG JRST A1S1FX JUMPL R,.+3 FAD TT,[1.0] JRST FLOAT1 FSB TT,[1.0] JRST FLOAT1 A1S1FX: CAMN TT,[1_43] JUMPL R,A1S11 ADD TT,R CAMN TT,[1_43] ;DONT WANT TO GET -2E35. BY ADD1 JUMPG R,ABSOV JRST FIX1 A1S11: PUSHJ P,ABSOV ;CANT SUB1 FROM -2E35. AND HRROS (A) A1S1BG: PUSH P,B ;ADD1 AND SUB1 FOR BIGNUM PUSH P,CPOPBJ MOVEI B,IN1 JUMPL R,.DIF JRST .PLUS ABSOV: PUSH P,B ;OVERFLOW FROM ADD1, SUB1, ABS, MOVEI TT,1 ; MINUS, HAIPART, GCD, ETC. PUSHJ P,C1CONS MOVE B,A MOVEI TT,0 PUSHJ P,C1CONS HRRM B,(A) PUSHJ P,BNCONS JRST POPBJ ;;; MOBY DISPATCH TABLES FOR THE VARIOUS ARITHMETIC OPERATIONS CAIA . ;UNUSED WORD JRST GRSWF COMPR: JRST GRSWX JFCL 0 JRST GRBFX JRST GRFXB JRST GRBB SKIPE VZFUZZ 0 FSBR D,TT DIFFA: SUB D,TT JRST PLOV JRST PL2BN JRST PL1BN JRST BNDF SKIPE VZFUZZ ;-3(R) SKIP UNLESS FUZZ HACK TO BE PULLED 0 ;-2(R) OPERATION IDENTITY - VALUE WHEN NO ARGS GIVEN FADR D,TT ;-1(R) FLOATING POINT INSTRUCTION FOR OPERATION PLUSA: ADD D,TT ;0(R) FIXED POINT INSTRUCTION FOR OPERATION JRST PLOV ;1(R) ACTION ON ARITHMETIC OVERFLOW JRST PL2BN ;2(R) BIGNUMBER ACCUMULATION MEETS FIXNUM ARG JRST PL1BN ;3(R) FIXNUM ACCUMULATION MEETS BIGNUM ARG JRST BNPL ;4(R) BIGNUM ACCUMULATION, BIGNUM ARG CAIA 1 FMPR D,TT TIMESA: IMUL D,TT JRST TIMOV JRST TIM2BN JRST TIM1BN JRST BNTIM CAIA 1 FDVR D,TT QUOA: JRST QUOAK JRST QUOOV JRST DV2BN JRST DV1BN JRST BNDV QUOOV: SKIPN RWG JRST OVFLER AOS D,T JFCL 8.,PLOV JRST T14E QUOAK: CAMN D,[400000,,0] ;ORDINARY FIXED POINT DIVISION JRST QUOAK1 ;DOESN'T ALWAYS WIN ON SETZ QUOAK2: IDIVM D,TT MOVE D,TT JRST T14EX2 QUOAK1: CAMN TT,XC-1 ;SETZ/(-1) => POSITIVE SETZ JRST DIVSEZ CAIN TT,1 ;SETZ/1 => SETZ JRST T14EX2 JRST QUOAK2 ;IDIVM WORKS FOR OTHER CASES T1: JUMPE T,NMCK0 ;ONLY ONE ARG GIVEN - GIVE IT OUT MOVE TT,-2(R) ;NO ARGS GIVEN - GIVE OUT OPERATORS IDENTITY JRST FIX1 .QUO: SKIPA R,[QUOA] ;C KEEPS ADDRESS OF FUNCTION TYPE .TIMES: MOVEI R,TIMESA SETZM REMFL JRST T21 .DIF: SKIPA R,[DIFFA] .PLUS: MOVEI R,PLUSA T21: MOVNI T,1 PUSH P,A PUSH P,B JRST T20 QUOTIENT: SKIPA R,[QUOA] TIMES: MOVEI R,TIMESA SETZM REMFL JRST T22 DIFFERENCE: SKIPA R,[DIFFA] PLUS: MOVEI R,PLUSA T22: AOJGE T,T1 T20: MOVE F,T ;D - ACCUMULATED VALUE ADDI F,1(P) ;TT - NEXT VALUE IN LINE HRL F,T T24: MOVNI T,-1(T) HRLS T ;R - ADDRESS OF INSTRUCTION DISPATCH TABLE MOVEM T,PLUS8 ;F - AOBJN POINTER TO ARG VECTOR ON PDL MOVE A,-1(F) JSP T,NVSKIP ;PICK UP FIRST ARG AND DISPATCH TO APPROPRIATE LOOP JRST T2 JRST T3 MOVE D,TT JRST 2,@[.+1] T4: MOVE A,(F) ;FLOATING POINT ARITHMETIC LOOP JSP T,NVSKIP JRST T6 JRST T5 T7: XCT -1(R) ;FLOATING SUM OPERATED WITH FLOATING NEXT ARG XCT -3(R) ;SKIP UNLESS ZFUZZ HACK REQUIRED JSP A,ZFZCHK T7A: AOBJN F,T4 JFCL 8.,T7O T7X: MOVE TT,D ;EXIT ARITHMETIC LOOP WITH ACCUMULATED VALUE T7X1: SUB P,PLUS8 JRST FLOAT1 T7O: JSP T,T7O0 JRST T7X1 ZFZCHK: MOVE T,D JRST 2,@[.+1] FDVR T,TT JFCL 8,ZFZCH9 MOVM T,T CAMGE T,@VZFUZZ SETZ D, ZFZCH9: JRST 2,(A) ;DON'T LET FDVR AFFECT OVERFLOW/UNDERFLOW ;;; IFN BIGNUM ;ARITH OPS FOR BIGNUM==1 CONTINUED T5: EXCH D,AGDBT JSP T,IFLOAT ;FLOATING SUM, NEXT IS FIXED POINT EXCH D,AGDBT JRST T7 T6: CAIN R,QUOA JRST T6A PUSHJ P,FLBIG ;FLOATING SUM, NEXT WAS BIGNUM JRST T7 T6A: PUSHJ P,FLBIGQ ;SPECIAL HACK FOR JPG JRST T7 SETZ D, ;IF BIGNUM TOO LARGE, WE GET JRST T7A ; UNDERFLOW, NOT OVERFLOW T3: MOVE D,TT ;FIXED POINT ARITHMETIC LOOP JRST 2,@[.+1] T15: MOVE A,(F) JSP T,NVSKIP XCT 3(R) ;DISPATCH TO CONVERT SUM TO BIGNUM JRST T14 ;OPERATE ON TWO FIXED POINT MOVEM TT,AGDBT MOVE TT,D ;FIXED POINT SUM CONVERTED TO FLOATING JSP T,IFLOAT ;AND ENTER FLOATING LOOP MOVE D,TT MOVE TT,AGDBT JRST T7 ;IFLOAT CANNOT HAVE SET OFVLO FLG T14: MOVE T,D ;SAVE OLD SUM, JUST INCASE THERE IS OVERFLO XCT 0(R) ;OPERATE FIXED POINT T14EX2: JFCL 8,1(R) ;CHECK FOR OVERFLO, IF SO DISPATCH TO BIGNUM T14E: AOBJN F,T15 T14EX: MOVE TT,D T14EX1: SUB P,PLUS8 JRST FIX1 FXIDEN: JSP T,FXNV1 JRST PDLNKJ FLIDEN: JSP T,FLNV1 JRST PDLNKJ ABS: JSP T,NVSKIP JRST ABSBG SKIPA T,CFIX1 MOVEI T,FLOAT1 JUMPGE TT,PDLNMK CAMN TT,[1_43] ;ABS OF -2**35. IS NO LONGER FIXNUM JRST ABSOV MOVMS TT JRST (T) REMAINDER: SETZB F,PLUS8 JSP T,NVSKIP JRST REMBIG ;BIGNUM SKIPA D,TT JSP T,REMAIR ;FLONUM IS ERROR - RETURNS TO THE NVSKIP EXCH A,B ;FIRST ARG IS FIXNUM JSP T,NVSKIP JRST REMAI2 ;IF SECOND IS BIGNUM NOW, MAYBE GIVE OUT FIRST SKIPA T,D JSP T,REMAIR ;FLONUM IS ERROR JUMPE TT,BPDLNKJ MOVE D,TT SETZ TT, ;IN THE CASE OF (\ SETZ 1), TRY TO WIN IDIV T,D JRST FIX1 REMAI2: SKIPL T,(B) ;WELL, IF FIRST ARG IS SETZ, AND JRST BPDLNKJ ; SECOND ARG IS +SETZ, THEN REMAINDER CAME T,[400000,,] ; SHOULD BE 0, NOT SETZ! JRST BPDLNKJ MOVE A,(A) PUSH P,AR1 ;MUST SAVE AR1 PUSHJ P,BNTRS1 ;SKIPS 2 UNLESS BIGNUM IS POP P,AR1 ; +SETZ (OR SETZ) JRST 0POPJ POP P,AR1 JRST BPDLNKJ FLOAT: TDZA R,R MOVEI R,TRUTH JSP T,NVSKIP JRST FLBIGF JRST FLOAT4 FIX4: JUMPE R,PDLNKJ ;ARG IS ALREADY OF REQUIRED TYPE. IF "CALL"ED, THEN RETURN LISP ANSWER IN A POPJ P, ;ELSE IF "NCALL"ED, RETURN NUMERIC ANSWER IN TT FLOAT4: JSP T,IFLOAT JUMPE R,FLOAT1 POPJ P, IFXERR: WTA [ARG TOO BIG FOR FIXNUM - IFIX!] JRST $IFIX1 $IFIX: PUSH P,CFIX1 $IFIX1: JSP T,FLTSKP POPJ P, CAML TT,[244000,,] JRST IFXERR JSP T,IFIX POPJ P, $FIX: JSP T,NVSKIP POPJ P, POPJ P, MOVM T,TT CAML T,[244000,,] JRST FIXBIG JRST FIX2 .GREAT: EXCH A,B .LESS: PUSH P,A PUSH P,B MOVNI T,2 LESSP: SKIPA A,[CAML D,2] GREATERP: HRLZI A,(CAMG D,) MOVEI D,GRFAIL MOVEI R,GRSUCE GTR1: MOVE F,T AOJGE T,GTR9 HRRI A,TT ADDI F,2(P) HRLI F,(T) PUSHJ FXP,SAV5M2 HRLI D,(JRST) MOVEM D,CFAIL HRLI R,(JRST) MOVEM R,CSUCE MOVEI R,COMPR MOVEM A,GRESS0 JRST T24 GTR9: MOVEI D,QMAX+1(A) SOJA T,WNALOSS MIN: SKIPA A,[CAML D,1] MAX: HRLOI A,(CAMG D,) AOJE T,NMCK0 MOVEI D,MXF MOVEI R,MXS SOJA T,GTR1 MXF: MOVE AR1,AR2A SKIPA D,TT MXS: MOVE AR2A,AR1 AOBJN F,GRSUC1 MAXFIN: MOVEI B,(AR1) PUSHJ FXP,RST5M2 2DIF JRST @(B),MAX923,QFIXNUM MAX923: T14EX ;FIXNUM T7X ;FLONUM T13X ;BIGNUM GRSUC2: MOVE D,TT GRSUC1: 2DIF JRST @(AR2A),GRS923,QFIXNUM GRS923: T15 ;FIXNUM T4 ;FLONUM T12 ;BIGNUM GRSUCE: AOBJN F,GRSUC2 GRSFIN: MOVEI A,TRUTH GRSF1: PUSHJ FXP,RST5M2 SUB P,PLUS8 POPJ P, GRFAIL: MOVEI A,NIL JRST GRSF1 GRSWF: SKIPA AR1,[QFLONUM] GRSWX: MOVEI AR1,QFIXNUM MOVE AR2A,AR1 JRST GRESS0 ] ;END OF ARITH OPS WITH BIGNUM==1 IFE BIGNUM,[ SUBTTL ARITHMETIC FUNCTIONS WITH BIGNUM==0 ADD1: JSP T,FLTSKP AOJA TT,FIX1 FAD TT,[1.0] JRST FLOAT1 SUB1: JSP T,FLTSKP SOJA TT,FIX1 FSB TT,[1.0] JRST FLOAT1 REMAINDER: JSP T,FXNV1 JSP T,FXNV2 IDIV TT,TT+1 MOVE TT,TT+1 JRST FIX1 MINUS: JSP T,FLTSKP SKIPA T,CFIX1 MOVEI T,FLOAT1 MOVNS TT JRST (T) ABS: JSP T,FLTSKP SKIPA T,CFIX1 MOVEI T,FLOAT1 MOVMS TT JRST (T) MINUSP: SKIPA R,[JUMPGE TT,FALSE] PLUSP: MOVE R,[JUMPLE TT,FALSE] JSP T,FLTSKP JFCL XCT R JRST TRUE ZEROP: JSP T,FLTSKP JFCL JUMPE TT,TRUE JRST FALSE $IFIX: $FIX: TDZA R,R MOVEI R,TRUTH JSP T,FIXFLO TLNN T,FL ;FIXFLO LEFT TYPE BITS IN T JRST FIX4 JSP T,IFIX JUMPE R,FIX1 POPJ P, FIX4: JUMPE R,PDLNKJ POPJ P, FLOAT: TDZA R,R MOVEI R,TRUTH JSP T,FIXFLO TLNN T,FX ;FIXFLO LEFT TYPE BITS IN T JRST FIX4 JSP T,IFLOAT JUMPE R,FLOAT1 POPJ P, FIXFLO: PUSH P,A LSH A,-SEGLOG HLL T,ST(A) ;LEAVES TYPE BITS IN T TLNN T,FX+FL JRST FLOAT3 POP P,A MOVE TT,(A) JRST (T) FLOAT3: POP P,A %WTA NMV3 JRST FIXFLO MIN: SKIPA A,[CAMLE F,1] MAX: HRLOI A,(CAMGE F,) AOJE T,NMCK0 MOVEI D,MINMAX SOJA T,MNMX1 MINMAX: XCT MNMX0 ;CAMG F,TT OR CAML F,TT MOVE F,TT JRST PLUS4 .GREAT: EXCH A,B .LESS: PUSH P,A PUSH P,B MOVNI T,2 LESSP: SKIPA A,[CAML F,2] GREATERP: HRLZI A,(CAMG F,) MOVEI D,GRESS MNMX1: HRLI D,(JRST) MOVEM D,PLUS3 MOVNM T,PLUS8 MOVE R,T AOJGE T,MNMX9 HRRI A,TT MOVEM A,GRESS0 ;THIS IS ALSO MNMX0 ADD R,P MOVE A,1(R) SETOM PLUS0 JSP T,FLTSKP SETZM PLUS0 MOVE F,TT AOJA R,PLUS7 MNMX9: MOVEI D,QMAX+1(A) SOJA T,WNALOSS GRESS: XCT GRESS0 JRST GRUSE MOVE F,TT CAME P,R JRST PLUS9 SUB P,PLUS8 JRST TRUE GRUSE: SUB P,PLUS8 JRST FALSE .DIF: PUSH P,A PUSH P,B MOVNI T,2 DIFFERENCE: MOVE R,[JRST DIF2] MOVE D,R SOJA D,DIF1 SKIPA D,[FSBR F,TT] DIF2: MOVE D,[SUB F,TT] MOVEM D,PLUS3 MOVE D,[FSBR F,TT] MOVEM D,PLUS6 MOVE F,TT JRST PLUS4 .QUO: PUSH P,A PUSH P,B MOVNI T,2 QUOTIENT: MOVE R,[JRST QUO2] MOVE D,R SOJA D,QUO1 SKIPA D,[FDVR F,TT] QUO2: MOVE D,[JRST QUO3] MOVEM D,PLUS3 MOVE D,[FDVR F,TT] MOVEM D,PLUS6 MOVE F,TT JRST PLUS4 QUO3: CAIN TT,1 CAME F,[400000,,0] CAIA SKIPA TT,F IDIVM F,TT EXCH F,TT ;ALL THIS LOSSAGE SO THAT F+1 WONT BE DISTURBED JFCL 8.,.+2 JRST PLUS4 SKIPN RWG JRST OVFLER SKIPGE TT SOSA F,TT AOS F,TT JFCL 8.,OVFLER JRST PLUS4 .TIMES: PUSH P,A PUSH P,B MOVNI T,2 TIMES: MOVE R,[IMUL F,TT] MOVE D,[FMPR F,TT] QUO1: MOVEI F,1 JRST PLUS1 .PLUS: PUSH P,A PUSH P,B MOVNI T,2 PLUS: MOVE R,[ADD F,TT] MOVE D,[FADR F,TT] DIF1: MOVEI F,0 PLUS1: MOVNM T,PLUS8 JUMPE T,PLUS2 ADD T,P MOVEM R,PLUS3 SETZM PLUS0 MOVE R,T PLUS7: MOVEM D,PLUS6 HRLS PLUS8 JRST 2,@[PLUS4] PLUS5: MOVE D,PLUS6 ;FAD F,TT OR FMP F,TT OR ETC. MOVEM D,PLUS3 SETOM PLUS0 EXCH F,TT JSP T,IFLOAT EXCH F,TT PLUS3A: XCT PLUS3 PLUS4: CAMN P,R JRST PLUS2 PLUS9: MOVE A,1(R) JSP T,FLTSKP JRST .+4 SKIPE PLUS0 AOJA R,PLUS3A AOJA R,PLUS5 SKIPE PLUS0 JSP T,IFLOAT AOJA R,PLUS3A PLUS2: MOVE TT,F JFCL 8.,PLUS2V PLUS2A: SUB P,PLUS8 ;FALL THRU TO MAKNUM SKIPN PLUS0 JRST FIX1 JRST FLOAT1 PLUS2V: JSP T,T7O0 JRST PLUS2A ] ;END OF ARITH OPS WITH BIGNUM=0 T7O0: SKIPE VZUNDERFLOW ;NON-NIL => FLOATING UNDERFLOW TLNN T,100 .SEE %PCFXU ; YIELDS ZERO RESULT INSTEAD OF ERROR JRST UNOVER MOVEI TT,0 JRST (T) SUBTTL GENERAL EXPONENTIATION ROUTINE EXPT: JRST 2,@[.+1] ;SUBR 2 - COMPUTE A^B EXCH A,B ;FIND TYPE OF EXPONENT FIRST IFN BIGNUM,[ JSP T,NVSKIP ;EXPONENT IS . . . JRST XPT.B ;IT'S A BIGNUM JRST XPT.X ;IT'S A FIXNUM EXCH A,B ;IT'S A FLONUM JSP T,NVSKIP ;BASE IS . . . JRST XPTBL ;BIGNUM BASE JSP T,IFLOAT ;FIXNUM BASE - FLOAT IT ] ;END OF IFN BIGNUM IFE BIGNUM,[ JSP T,FLTSKP ;EXPONENT IS . . . JRST XPT.X ;IT'S A FIXNUM EXCH A,B ;IT'S A FLONUM JSP T,FLTSKP ;BASE IS . . . JSP T,IFLOAT ;FIXNUM BASE - FLOAT IT ] ;END OF IFE BIGNUM XPTLL: PUSH P,CFLOAT1 ;FLONUM^FLONUM SKIPN (B) ; X^0.0 => 1.0 JRST 1.0PJ JUMPE TT,CPOPJ ; 0.0^X => 0.0 PUSHJ P,LOG.. ;SO COMPUTE FLONUM^FLONUM BY USING THE FORMULA: FMPR TT,(B) ; B (B LOG A) JRST EXP.. ; A = E XPT.X: EXCH A,B ;FIXNUM EXPONENT FOUND MOVE D,TT BG$ JSP T,NVSKIP ;CHECK BASE FOR FIXNUM EXPONENET BG$ JRST XPTBX ;BIGNUM BASE BG% JSP T,FLTSKP JRST XPTXX0 ;FIXNUM BASE PUSH P,CFLOAT1 ;FLONUM BASE => FLONUM RESULT XPTLX: JSP R,XPTZL ;CHECK EASY CASES SKIPA R,TT ;NORMAL CASE - USE THE MULTIPLY XPTLX1: FMPR R,R ; AND SQUARE HACK TRNE D,1 FMPR T,R JFCL 8,XPTOV ;CHECK FOR OVERFLOW LSH D,-1 JUMPN D,XPTLX1 XPTLX2: MOVE TT,T ;ANSWER GOES IN TT POPJ P, XPTOV: JSP T,T7O0 POPJ P, XPTXX0: PUSHJ P,XPTXX JRST FIX1 POPJ P, ;;; SKIPS IF ANSWER IS A BIGNUM XPTXX: JSP R,XPTZX ;FIXNUM^FIXNUM - CHECK EASY CASES JUMPL D,ZPOPJ IFE BIGNUM,[ SKIPA R,TT XPTXX5: IMUL R,R TRNE D,1 IMUL T,R LSH D,-1 JUMPN D,XPTXX5 MOVE TT,T JFCL 8,XPTOV POPJ P, ] ;END OF IFE BIGNUM IFN BIGNUM,[ SKIPGE R,TT JRST XPTXX3 JFFO R,.+1 LSH R,1(F) JUMPE R,2XPT ;XPTZX HAS CHECKED BASE, SO IT'S NOT 0/1/-1 MOVE R,TT XPTXX3: MOVE TT,T ;HERE YOU GO FANS, YOU BASIC MULTIPLY BY SQUARING LOOP. MOVEM D,NORMF TRNE D,1 IMUL T,R JFCL 8.,EXPT6C LSH D,-1 JUMPN D,XPTXX4 MOVE TT,T POPJ P, XPTXX4: MOVE F,R IMUL R,R JFCL 8.,EXPT6B JRST XPTXX3 2XPT: MOVNI F,(F) IMULI D,36.-1(F) MOVEI TT,1 CAIL D,35. JRST 2BGXPT ASH TT,(D) POPJ P, 2BGXPT: IDIVI D,35. ASH TT,(R) JSP T,FIX1A PUSHJ P,NCONS 2BGXP1: MOVE B,CIN0 PUSHJ P,XCONS SOJG D,2BGXP1 PUSHJ P,BGNMAK JRST POPJ1 ] ;END OF IFN BIGNUM IFN BIGNUM,[ XPTBL: PUSH P,A ;BIGNUM^FLONUM PUSHJ P,FLBIG ;SO FLOAT THE BIGNUM, THEN USE SUB P,R70+1 ; FLONUM^FLONUM JRST XPTLL XPT.B: EXCH A,B ;BIGNUM FOUND AS EXPONENT HLRZ D,(TT) HRRZ D,(D) TLNE TT,400000 TLO D,400000 ;D GETS SIGN-BIT IN 4.9, RANDOM-NON-ZERO-BIT IN 3.1 TLO D,1 ;AND ODDP-BIT IN 1.1 JSP T,NVSKIP JRST OVFLER JRST XPTZX0 PUSH P,CFLOAT1 JSP R,XPTZL ;FLONUM^BIGNUM -- CHECK EASY CASES MOVMS TT CAML TT,T ;T SUPPOSED TO HAVE 1.0 JRST OVFLER SKIPN VZUNDERFLOW JRST UNFLER JRST ZPOPJ ;PUTS A RANDOM ZERO IN TT, AND POPJS XPTZX0: PUSH P,CFIX1 JSP R,XPTZX ;FIXNUM^BIGNUM -- CHECK EASY CASES JUMPL D,ZPOPJ ;N^- ==> 0 JRST OVFLER ;;; MUST SKIP 1 AS POPJ SINCE ONLY COME HERE FROM XPTXX EXPT6B: MOVE R,F ;RESTORE R, AND LEAVE OLD D IN NORMF EXPT6C: PUSHJ FXP,SAV5 ;EXPECTS RUNNING SQUARER IN R, ACCUMULATION IN TT PUSHJ P,BNCV ;NOTE THAT D CANT BE ZERO WHEN WE COME HERE MOVE B,A ;ACCUMULATION AS BIGNUM IN B MOVE TT,R PUSHJ P,BNCVTM MOVE A,TT ;RUNNING SQUARER IN A EXPT1A: MOVEM A,-4(P) MOVE D,NORMF EXPT1: TRNN D,1 ;-4(P) AND A HAVE RUNNING SQUARER, B HAS ACCUMULATION JRST EXPT2 MOVEM D,NORMF PUSHJ P,BNMUL MOVE D,NORMF EXCH A,-4(P) EXPT3: LSH D,-1 ;-4(P) NOW HAS ACCUMULATION, A HAS RUNNING SQUARER JUMPE D,EXPT4 MOVE B,A MOVEM D,NORMF PUSHJ P,BNMUL MOVE B,-4(P) JRST EXPT1A EXPT2: MOVEM B,-4(P) JRST EXPT3 EXPT4: JSP R,RSTR5 PUSHJ P,BNCONS JRST POPJ1 XPTBX: SOJG D,XPTBX1 ;BIGNUM^FIXNUM AOJG D,CPOPJ ; X^1 => X MOVEI A,IN0 JUMPL D,CPOPJ ; X^-N => 0 AOJA A,CPOPJ ; X^0 => 1 ;HACK HACK - IN0 => IN1 XPTBX1: MOVE A,TT ;EXPONENT > 1 SOS (P) ;COUNTERACT POPJ1 IN EXPT1 PUSHJ FXP,SAV5 MOVE B,BN.1 ;1, STORED AS A BIGNUM AOJA D,EXPT1 ;RESTORE VALUE OF D ] ;END OF IFN BIGNUM XPTII: PUSH P,CFIX1 ;SUBR 2 NCALLABLE (REAL NAME: ^) JSP T,FXNV1 JSP T,FXNV2 JRST 2,@[.+1] PUSHJ P,XPTXX POPJ P, LERR [SIXBIT \RESULT LARGER THAN FIXNUM - #^!\] XPTI$: PUSH P,CFLOAT1 ;SUBR 2, NCALLABLE (REAL NAME: ^$) JSP T,FLNV1 JSP T,FXNV2 JRST 2,@[XPTLX] ;OVERFLOW MUST BE CLEAR ON ENTRY TO XPTLX XPTZL: JUMPN TT,XPTZL1 ;FLONUM BASE (CFLOAT1 ON PDL) SKIPN D ; 0.0^X => 0.0, 1.0PJ: MOVSI TT,(1.0) ; EXCEPT 0.0^0.0 => 1.0 POPJ P, XPTZL1: JUMPGE D,XPTZL2 ; -Y 1 Y MOVSI T,(1.0) ; X = (---) FDVR T,TT ; X MOVE TT,T MOVMS D XPTZL2: CAMN TT,[-1.0] JRST XPTM1 ;BASE IS -1.0 CAMN TT,[1.0] POPJ P, ;BASE IS 1.0 MOVSI T,(1.0) ;T GETS 1.0 IN ANY CASE JRST (R) XPTZX: JUMPN TT,XPTZX1 ;FIXNUM BASE - PDL HAS CFIX1 JUMPN D,CPOPJ ; 0^X => 0, AOJA TT,CPOPJ ; EXCEPT 0^0 => 1 XPTZX1: CAMN TT,XC-1 ;BASE = -1 JRST XPTM1 CAIN TT,1 ;FOR BASE = 1, ALSO EASY POPJ P, MOVEI T,1 ;T GETS 1 IN ANY CASE JRST (R) XPTM1: TRNN D,1 ;FOR BASE = -1 OR -1.0, SIMPLY MOVMS TT ; ASCERTAIN PARITY OF EXPONENT POPJ P, SUBTTL RANDOM RANDOM: SKIPA F,CFIX1 MOVEI F,CPOPJ AOJG T,RNDM0 AOJLE T,RAND9 POP P,A JUMPE A,IRAND ;ONE ARG OF NIL CAUSES INITIALIZATION PUSH P,F JSP F,RNDM0 MOVE D,TT ;ANY OTHER ARGUMENT SHOULD BE A JSP T,FXNV1 ; FIXNUM N, AND WE GENERATE A JUMPLE TT,RAND1 ; FIXNUM IN THE RANGE 0 TO N-1 TLZ D,400000 IDIV D,TT SKIPA TT,R RAND1: SETZ TT, ;RETURN 0 FOR NON-POSITIVE ARGUMENTS POPJ P, IRAND: MOVE TT,[171622221402] ;A GOOD STARTING NUMBER IRAND0: MOVEI T,LRBLOCK-1 ;INITIALIZE THE RANDOMNESS IRAND3: MOVE D,TT MULI D,3125. DIV D,[377777777741] MOVEM R,TT TLCE T,400000 JRST IRAND5 HRLM R,RBLOCK(T) JRST IRAND3 IRAND5: HRRM R,RBLOCK(T) SOJGE T,IRAND3 MOVEI D,ROFSET MOVEM D,RNOWS RNDM1: MOVEI T,LRBLOCK-1 MOVEM T,RBACK JRST RNDM1A RNDM2: MOVEI D,LRBLOCK-1 MOVEM D,RNOWS JRST RNDM2A RNDM0: SOSGE T,RBACK ;BASIC COMBINATION FOR RANDOMNESS JRST RNDM1 RNDM1A: SOSGE D,RNOWS JRST RNDM2 RNDM2A: MOVE TT,RBLOCK(T) ADDB TT,RBLOCK(D) JRST (F) SUBTTL HAULONG FUNCTION HAULONG: PUSH P,CFIX1 .HAU: BG$ JSP T,NVSKIP BG$ JRST 1HAU BG% JSP T,FLTSKP JRST 4HAU %WTA FXNMER JRST .HAU 4HAU: MOVM D,TT MOVEI TT,35.+1 3HAU1: JFFO D,.+2 TDZA TT,TT SUBI TT,(R) POPJ P, IFN BIGNUM,[ 1HAU: MOVEI F,(TT) ;RECEIVES BN HEADER IN TT HRRZ R,(F) ;LEAVES HAULONG IN TT, PTR TO NEXT TO LAST MOVEI TT,35.+1 ;IN F, CNT OF # OF ZEROS FOR LAST WD IN R JUMPE R,3HAU 2HAU: ADDI TT,35. HRRZ D,(R) JUMPE D,3HAU MOVEI F,(R) MOVEI R,(D) JRST 2HAU 3HAU: HLRZ T,(R) MOVE D,(T) JRST 3HAU1 ] ;END OF IFN BIGNUM SUBTTL HAIPART FUNCTION HAIPART: IFN BIGNUM,[ JSP T,NVSKIP JRST 1HAI ] IFE BIGNUM, JSP T,FLTSKP JRST 0HAI %WTA FXNMER JRST HAIPART 0HAI: MOVM TT,TT JFFO TT,.+2 JRST 0POPJ ;FOR ZERO ARG, JUST RETURN ARG! HRREI F,-36.(D) ;-<# OF BITS IN ARG> NO IN AC F JSP T,FXNV2 JUMPLE D,0HAI1 ADD D,F JUMPG D,0HAI4 ;MORE DIGITS REQUESTED THAN ARE AVAILABLE? LSH TT,(D) ;GETTING HAI PART INTO AC TT JUMPGE TT,FIX1 IFN BIGNUM, JRST ABSOV IFE BIGNUM, JRST OVFLER 0HAI4: SKIPL (A) JRST PDLNKJ JRST MNSFX ;NEGATE THE FIXNUM, TO GET "ABS" 0HAI1: JUMPE D,0POPJ ;RETURNS A FIXNUM ZERO CAMGE D,F JRST 0HAI4 MOVNS D 0HAI2: SETO F, ;REQUESTING LOW PART BY NEG COUNT LSH F,(D) ;CREATE MASK TO LET PROPER BITS THRU ANDCM TT,F JRST FIX1 IFN BIGNUM*USELESS,[ 3HAI: MOVNS D ;ACTUALLY ASKING FOR LOW PART CAILE D,35. JRST 3HAI1 JUMPE D,0POPJ HLRZ TT,(TT) MOVE TT,(TT) JRST 0HAI2 3HAI1: PUSH FXP,D PUSHJ P,1HAU POP FXP,D CAIL D,(TT) JRST PDLNKJ IDIVI D,35. PUSH P,C MOVEI F,C ;F WILL BE POINTER TO LAST OF FORMNG LIST MOVE C,(A) ;C HOLDS POINTER TO FNAL RESULT MOVEI B,(C) ;B GOES CDR'ING DOW INPUT ARG 3HAI2: HLRZ TT,(B) MOVE TT,(TT) PUSHJ P,C1CONS HRRM A,(F) MOVEI F,(A) HRRZ B,(B) SOJG D,3HAI2 ;D HOLDS HOW MANY WORDS TO USE JUMPE R,3HAI3 ;R HOLDS HOW MANY LEFT OVER BITS FROM D WORDS HLRZ TT,(B) MOVE TT,(TT) MOVNI D,1 LSH D,(R) ANDCM TT,D JUMPE TT,3HAI3 PUSHJ P,C1CONS HRRM A,(F) 3HAI3: MOVEI A,(C) PUSH P,AR1 PUSHJ P,BNTRUN ;IN LOPART CASE, MAY NEED TO GET POP P,AR1 ; RID OF LEADING ZEROS POP P,C HRRZ B,(A) ;MAYBE WHAT WE HAVE IS SHORT ENOUGH JUMPN B,BGNMAK ; TO FIT IN A FIXNUM; IF SO, WE CAN JRST CAR ; USE ONE WE JUST CONSED FOR BIGNUM! ] ;END OF IFN BIGNUM*USELESS SUBTTL LENGTH AND BIGP FUNCTIONS LNGTER: WTA [NON-LIST - LENGTH!] JRST LNGTH0 LENGTH: SKIPA T,CFIX1 MOVEI T,CPOPJ LNGTH0: SKIPE V.RSET JRST LNGTH5 ;FOR *RSET MODE, USE SLOW ERROR-CHECKING LOOP LNG1A: MOVEI TT,777777 .SEE $LISTEN ;SAVES R LNGTH1: JUMPE A,LNGTH2 HRRZ A,(A) SOJG TT,LNGTH1 LNGTE1: MOVEI TT,(A) ;MAKNUM JSP T,FXCONS WTA [CIRCULAR POINTER - LENGTH!] JRST LNGTH0 LNGTH2: XORI TT,777777 ;ONE'S COMPLEMENT! JRST (T) LNGTH5: MOVEI TT,777777 LNGTH6: SKIPN D,A ;DONE IF NIL SEEN JRST LNGTH2 LSH D,-SEGLOG SKIPL ST(D) .SEE LS JRST LNGTER HRRZ A,(A) SOJG TT,LNGTH6 JRST LNGTE1 IFE BIGNUM, BIGP==:FALSE IFN BIGNUM,[ BIGP: PUSHJ P,TYPEP ;SUBR 1 - IS IT A BIGNUM? CAIE A,QBIGNUM SETZ A, ;RETURNS T OR NIL JRST NOTNOT ] ;END OF IFN BIGNUM SUBTTL BOOLE AND ODDP FUNCTIONS BOOLE: SKIPA F,CFIX1 MOVEI F,CPOPJ MOVE R,T ADDI R,2(P) HRLI T,-1(T) MOVEM T,PLUS8 MOVE A,-1(R) JSP T,FXNV1 DPB TT,[350400,,BOOLI] PUSHJ P,BOOLG MOVE D,TT BOOLL: PUSHJ P,BOOLG XCT BOOLI JRST BOOLL BOOLG: CAIL R,(P) JRST BOOL1 MOVE A,(R) JSP T,FXNV1 AOJA R,CPOPJ BOOL1: ADD P,PLUS8 POP P,B JRST (F) ODDP1: %WTA FXNMER ODDP: SKOTT A,FX IFN BIGNUM, JRST ODDP4 IFE BIGNUM, JRST ODDP1 ODDP2: MOVE TT,(A) ODDP21: TRNN TT,1 JRST FALSE JRST TRUE IFN BIGNUM,[ ODDP4: TLNN TT,BN JRST ODDP1 MOVE TT,(A) ODDP3: HLRZ TT,(TT) MOVE TT,(TT) JRST ODDP21 ] ;END OF IFN BIGNUM SUBTTL FSC, ROT, LSH, AND GCD FUNCTIONS $FSC: JSP T,FLTSKP ;SUBR 2 JFCL JSP T,FXNV2 CAIG D,-1 FSC TT,(D) JRST FLOAT1 $ASH: HRLZI R,(ASH TT,(D)) JRST SHIFTX $ROT: SKIPA R,[ROT TT,(D)] ;SUBR 2 $LSH: HRLZI R,(LSH TT,(D)) ;SUBR 2 SHIFTX: PUSH P,CFIX1 SHIFTY: JSP T,FLTSKP JFCL JSP T,FXNV2 XCT R POPJ P, $LOADB: PUSH P,CFIX1 JSP T,FXNV1 JSP T,FXNV2 JSP T,FXNV3 JRST .LODB1 %LOADB: PUSH P,CFIX1 .LODB1: MOVE D,(C) ROT D,-6 HRR D,(B) ROT D,-6 MOVE TT,(A) JRST .LDB2 $LDB: PUSH P,CFIX1 JSP T,FXNV1 JSP T,FXNV2 EXCH TT,D LSH D,30 JRST .LDB2 %LDB: PUSH P,CFIX1 MOVE D,(A) MOVE TT,(B) .LDB2: HRRI D,TT LDB TT,D POPJ P, $DEPOB: PUSH P,CFIX1 JSP T,FXNV1 JSP T,FXNV2 JSP T,FXNV3 JSP T,FXNV4 JRST .DPOB1 %DEPOB: PUSH P,CFIX1 .DPOB1: MOVE R,(AR1) MOVE D,(C) ROT D,-6 HRR D,(B) ROT D,-6 MOVE TT,(A) JRST .DPB2 $DPB: PUSH P,CFIX1 JSP T,FXNV1 JSP T,FXNV2 JSP T,FXNV3 LSH D,30 JRST .DPB1 %DPB: PUSH P,CFIX1 ;Args = ( newbyte position_30.+size_24. word ) MOVE D,(B) .DPB1: MOVE TT,(C) MOVE R,(A) .DPB2: HRRI D,TT DPB R,D ;puts result in TT POPJ P, IFN USELESS,[ IFE BIGNUM, GCD: .GCD: PUSH P,CFIX1 ;SUBR 2 - NCALLABLE JSP T,FXNV1 ;GCD OF FIXNUM ARGS ONLY JSP T,FXNV2 MOVM TT,TT ;GCD(-X,Y) = GCD(X,Y) MOVM D,D ;GCD(X,-Y) = GCD(X,Y) .GCD0: JUMPE TT,.GCD2 ;GCD(0,Y) = ABS(Y) JUMPE D,CPOPJ ;GCD(X,0) = ABS(X) CAMGE D,TT EXCH D,TT JRST .GCD1 .GCD3: MOVE D,TT MOVE TT,R .GCD1: IDIV D,TT ;GOOD OLD EUCLIDEAN ALGORITHM JUMPN R,.GCD3 POPJ P, .GCD2: MOVE TT,D POPJ P, IFN BIGNUM,[ GCD0: %WTA FXNMER ;NON-FIXNUM VALUE GCD: SETZ R, ;SUBR 2 - GCD, EVEN OF BIGNUM ARGS JSP T,NVSKIP TRO R,1 ;TURN ON BIT IF BIGNUM JRST .+2 ;FIXNUMS ARE OK TOO JRST GCD0 ;DON'T LIKE FLONUMS EXCH A,B MOVE D,TT JSP T,NVSKIP ;NOW CHECK OTHER ARG TRO R,2 JRST .+2 JRST GCD0 ;I TOLD YOU, I DON'T LIKE FLONUMS! JRST .+1(R) ;SO FIGURE OUT THIS MESS JRST GCDXX ;FIXNUM AND FIXNUM EXCH A,B ;FIXNUM AND BIGNUM JRST GCDBX ;BIGNUM AND FIXNUM JRST GCDBG ;BIGNUM AND BIGNUM GCDXX: MOVM TT,TT ;GCD OF TWO FIXNUMS JUMPL TT,GCDOV1 ;CHECK OUT -400000000000 CASES MOVM D,D JUMPL D,GCDOV PUSH P,CFIX1 ;EVERYTHING OKAY - CAN USE .GCD0 JRST .GCD0 ] ;END OF IFN BIGNUM ] ;END OF IFN USELESS SUBTTL FUNCTIONS: = < > 1+ 1+$ 1- 1-$ $EQUAL: JSP T,FLTSKP ;NUMERIC EQUAL = JRST IEQUAL EXCH A,B MOVE D,TT $EQL1: JSP T,FLTSKP JRST 2EQNF $IEQ: CAME D,TT JRST FALSE JRST TRUE IEQUAL: EXCH A,B MOVE D,TT JSP T,FLTSKP JRST $IEQ JRST 1EQNF $LESS: EXCH A,B $GREAT: JSP T,FLTSKP ;NUMERIC GREATERP AND LESSP <,> JRST IGRT MOVE D,TT EXCH A,B $IGL1: JSP T,FLTSKP JRST 2GPNF $IGL: CAMG D,TT JRST FALSE JRST TRUE IGRT: MOVE D,TT MOVE A,B JSP T,FLTSKP JRST $IGL JRST 1GPNF IADD1: JSP T,FLTSKP ;FIXNUM ADD1 1+ AOJA TT,FIX1 %WTA IARERR JRST IADD1 %WTA $ARERR $ADD1: JSP T,FLTSKP ;FLONUM ADD1 1+$ JRST $ADD1-1 FADRI TT,(1.0) JRST FLOAT1 ISUB1: JSP T,FLTSKP ;FIXNUM SUB1 1- SOJA TT,FIX1 %WTA IARERR JRST ISUB1 %WTA $ARERR $SUB1: JSP T,FLTSKP ;FLONUM SUB1 1-$ JRST $SUB1-1 FSBRI TT,(1.0) JRST FLOAT1 SUBTTL FUNCTIONS: + +$ - -$ * *$ // //$ $ARITH: SETOM PLUS0 SKIPA IARITH: SETZM PLUS0 ;SET UP FOR FIXNUM ARITHMETIC AOJGE T,ARIT0 I$B: JRST 2,@[.+1] SKIPA B,T I$ART2: XCT R POP P,A ;MAIN LOOP FOR FIXNUM AND FLONUM ARITHMETIC ARITH: JSP T,FLTSKP ;MAKE SURE NO MIXED MODES, RETURN MACHINE NUMBER IN TT TDZA T,T MOVNI T,1 CAME T,PLUS0 JRST ARTHER AOJLE B,I$ART2 CAIN B,69.+1 ;SIGNAL FOR CASE WITH ONE ARG EXCH TT,D XCT F IARDS: SKIPE PLUS0 ;DISPATCH TO CONS UP FINAL ANSWER JRST FLOAT1 JRST FIX1 ARIT0: MOVE TT,D JUMPN T,IARDS MOVEI T,69. JRST I$B IDIFFERENCE: SKIPA F,[SUB TT,D] ;- IPLUS: MOVE F,[ADD TT,D] ;+ MOVE R,[ADD D,TT] MOVEI D,0 JRST IARITH IQUOTIENT: SKIPA F,[IDIV TT,D] ;/ ITIMES: MOVE F,[IMUL TT,D] ;* MOVE R,[IMUL D,TT] MOVEI D,1 JRST IARITH $DIFFERENCE: SKIPA F,[FSBR TT,D] ;-$ $PLUS: MOVE F,[FADR TT,D] ;+$ MOVE R,[FADR D,TT] MOVEI D,0 JRST $ARITH $QUOTIENT: SKIPA F,[FDVR TT,D] ;/$ $TIMES: MOVE F,[FMPR TT,D] ;*$ MOVE R,[FMPR D,TT] MOVSI D,(1.0) JRST $ARITH IARZAR: MOVE TT,D JRST FIX1 ;;; ********** NUMBER SUBRS FOR LISP ********** SUBTTL SIN AND COS FUNCTIONS ;;; SIN IS A TOPS-10/TENEX JSYS, SO MUST CALL THIS $SIN. FOO! - GLS $SIN: PUSH P,CFLOAT1 SIN.: JSP T,FLTSKP JSP T,IFLOAT MOVM T,TT ;SIN(-X)=-SIN(X) CAMLE T,C1.0E5 ;ARG SHOULD BE <= 1.0E5 (ELSE RESULT JRST SIN.ER ; WOULD BE GROSSLY INACCURATE) CAMG T,[.001] ;THE RELATIVE ERROR OF APPROXIMATION [BY THIS RATIONAL ; ; FUNCTION] IS BOUNDED BY ABOUT 2.0E-7, BUT OCCASIONALLY ; ; COMES CLOSE TO THIS. SINCE THE ERROR OF TRUNCATION ; ; INHERENT IN TAKING X-(1/6)*X**3 FOR THE TAYLOR SERIES ; ; OF SIN(X) IS MUCH LESS THAN 2.0E-7, IT WILL BE SUFFICIENT ; ; TO TAKE X FOR SIN(X) WHENEVER THE RELATIVE ERROR TERM ; ; [(1/6)*X**3] IS LESS THAN 2.0E-7. SOLVING, WE FIND JRST SIN.XT ; X=.001 WILL DO. EXCH T,TT SIN.0: FDVR TT,PI%2 ;DIVIDE ARG BY PI/2 (ARG IS NOW IN QUADRANTS) MULI TT,400 ;TT GETS CHARACTERISTIC, R GETS MANTISSA SETZB R,F ASHC D,-243(TT) ;D GETS INTEGER PART, R GETS FRACTION (OF ARG) ASHC R,-8. ;R GETS HIGH 27. BITS OF FRACTION, F GETS REST TLO R,200000 ;FLOAT R LSH F,-8. TLO F,145000 ;FLOAT F (NOTE: 145=200-33; R,F NOW FORM 2-WORD FLOATING NUMBER) FADR R,F ;ADD F TO R (THIS WHOLE MESS PRESERVES PRECISION AND NORMALIZES) TRCN D,3 ;R IS NOW A QUADRANT 1 ANGLE - WHAT WAS ORIGINAL QUADRANT? JRST SIN.1 ;QUADRANT 1 - ALL IS WELL TRCE D,3 MOVN T,T ;QUADRANT 2 OR 3 - MUST REVERSE SIGN: SIN(X)=-SIN(X-PI) TRNE D,1 FSBR R,FPWUN ;QUADRANT 2 OR 4 - SUBTRACT 1 TO PUT IN RANGE -1.0 TO 0 SIN.1: SKIPGE T ;TEST SINE SIGN FLAG MOVN R,R ;IF NEGATIVE, RESULT MUST BE NEGATIVE MOVE D,R FMPR D,D ;D <- R*R IS ALWAYS NON-NEGATIVE MOVE TT,SIN.CF+4 ;MOBY APPROXIMATION MOVEI T,3 SIN.2: FMPR TT,D FADR TT,SIN.CF(T) SOJGE T,SIN.2 FMPR TT,R SIN.XT: CAMLE TT,[1.0] ;THIS IS A CROCK TO MAKE SURE ABS(RESULT) NOT >1 MOVSI TT,(1.0) CAMGE TT,[-1.0] MOVSI TT,(-1.0) POPJ P, ;RETURN - RESULT IS IN TT PI%2: 1.570796326 ;A PIECE OF PI (ABOUT 50%) SIN.CF: 1.5707963185 ;COEFFICIENTS FOR SIN APPROXIMATION -0.6459637111 0.07968967928 -0.00467376557 0.00015148419 COS: PUSH P,CFLOAT1 COS.: JSP T,FLTSKP JSP T,IFLOAT SKIPLE T,TT MOVN T,TT FADR T,PI%2 ;PI/2-X IN T, SINCE COS(X) = SIN(PI/2-X) MOVM TT,T ;|PI/2-X| IN TT CAMLE TT,C1.0E5 JRST COS.ER JRST SIN.0 SUBTTL SQRT FUNCTION COMMENT | OLD SQRT ALGORITHM SQRT: PUSH P,CFLOAT1 SQRT.: JSP T,FLNV1 JUMPL TT,SQR$ER ;NEGATIVE ARG IS AN ERROR SQRT..: MOVE D,TT ;D GETS ARG LDB T,[341000,,TT] ;FOR FIRST APPROXIMATION, TRY ADDI T,100 ; HALVING CHARACTERISTIC OF ARGUMENT, DPB T,[331100,,TT] ; AND USE SAME MANTISSA MOVEI T,5 ;NOW DO MOBY ITERATION SQRT.1: MOVE R,TT ; R <- TT MOVE TT,D FDVR TT,R ; R + D/R FADR TT,R ; TT <- --------- FSC TT,-1 ; 2 SOJN T,SQRT.1 POPJ P, | ;END OF OLD SQRT ALGORITHM COMMENT | ANOTHER OLD SQRT ALGORITHM ;;; THIS SQRT ALGORITHM IS BASED ON ONE BY KAHAN, ORIGINALLY ;;; DESIGNED FOR THE IBM 7094. THAT VENERABLE MACHINE LOOKED ;;; LIKE THE PDP-10 (27.-BIT MANTISSA AND 8-BIT EXPONENT). ;;; (THANKS TO RJF FOR HELP IN CODING THIS.) ;;; ;;; THE IDEA IS TO DECOMPOSE THE ARGUMENT X INTO: ;;; F * 2.0^(2*I - J) ;;; WHERE THE FRACTION F IS BETWEEN 0.5 (INCLUSIVE) AND 1.0 ;;; (EXCLUSIVE), AND I AND J ARE INTEGERS, J BEING 0 OR 1. ;;; ONE THEN COMPUTES THE INITIAL APPROXIMATION AS: ;;; A0 = (C + F/2.0 - J/4.0) * 2.0^I ;;; WHERE C IS THE MAGIC CONSTANT 0.4826004, CHOSEN FOR THE ;;; BEST POSSIBLE FIT TO A CURVE. ONE THEN PERFORMS AN ;;; ITERATION CALCULATING: ;;; A = (A + X/A)/2.0 ;;; ALL ARITHMETIC IS DONE WITHOUT ROUNDING EXCEPT LAST ADD. ;;; THREE ITERATIONS SHOULD SUFFICE; A3 IS THE RESULT. ;;; THE INITIAL APPROXIMATION CAN BE CALCULATED QUICKLY BY ;;; MEANS OF THE FOLLOWING TRICK. LET THE EXPONENT BE ;;; E = 2*I - J = 2*N + M ;;; SUCH THAT M IS 0 OR 1; THEN J=M AND I=N+M. MOREOVER, ;;; NOTE THAT THE PDP-10 EXPONENT X=E+200 (OCTAL), BECAUSE ;;; OF EXCESS-200 NOTATION. HENCE X=2*(N+100)+M. ;;; WE FIRST PICK OFF THE M BIT AS A SEPARATE WORD AND ;;; SHIFT IT RIGHT. THANKS TO THE PARTICULAR REPRESENTATION ;;; OF EXPONENT AND FRACTION, THIS PRODUCES A WORD WITH ;;; A FRACTION OF M/2. NOW WE WILL ADD TOGETHER THIS WORD, ;;; THE ORIGINAL ARGUMENT, AND A MAGIC CONSTANT, AND SHIFT ;;; THE SUM RIGHT BY 1. SHIFTING AFTERWARDS GIVES GREATER ;;; ACCURACY AND TAKES FEWER INSTRUCTIONS, BUT FOR PURPOSES ;;; OF EXPOSITION LET US ASSUME THE THREE SUMMANDS TO HAVE ;;; BEEN PRE-SHIFTED. ;;; SHIFTING THE ORIGINAL ARGUMENT RIGHT PRODUCES A WORD WITH ;;; FRACTION F/2+M/2 AND MACHINE EXPONENT N+100. SHIFTING ;;; THE M/2 PRODUCES M/4. THE MAGIC CONSTANT IS CHOSEN SUCH ;;; THAT, WHEN SHIFTED, ITS FRACTION IS C (0.4826004) AND ;;; ITS MACHINE EXPONENT IS 100. ADDING THESE TOGETHER ;;; PRODUCES FRACTION F/2 + 3*M/4 + C AND MACHINE EXPONENT ;;; N+200. HOWEVER, SINCE F IS NORMALIZED, THE ADDITION ;;; OF 3*M/4 IS GUARANTEED TO OVERFLOW INTO THE EXPONENT FIELD; ;;; THIS RESULTS IN SUBTRACTING M/4 FROM THE FRACTION, AND ;;; ADDING M INTO THE MACHINE EXPONENT. THE RESULT IS THUS: ;;; (C + F/2 - M/4) * 2.0^(N+M) ;;; WHICH IS THE DESIRED VALUE. SQRT: PUSH P,CFLOAT1 SQRT.: JSP T,FLNV1 JUMPG TT,SQRT.. JUMPL TT,SQR$ER ;NEGATIVE ARG IS AN ERROR POPJ P, ;ZERO ARGUMENT => ZERO ;;; POSITIVE ARGUMENT IS IN TT NOW SQRT..: MOVE R,TT ;SAVE ARGUMENT IN R FOR LATER MOVS D,TT ANDI D,1000 LSH D,22-1 ;D HAS M/2 AS A SINGLE BIT ADD TT,D ;ADD INTO ORIGINAL ARGUMENT ADD TT,[200756135462] ;EXPONENT 200, FRACTION 2*0.4826004 LSH TT,-1 ;NOW WE HAVE INITIAL APPROXIMATION IRPC ROUND,,[ R]AC,,[DDR] IFSN AC,R, MOVE D,R ; TT + R/TT FDV AC,TT ;COMPUTE TT <- --------- FAD!ROUND TT,AC ; 2 FSC TT,-1 ;LAST TIME ONLY, ADD ROUNDED TERMIN POPJ P, | ;END OF ANOTHER OLD SQRT ALGORITHM ;;; I HAVE NO IDEA HOW THIS WORKS! - GLS ;;; THANKS TO RJF AND KAHAN. ;;; KAHAN CLAIMS THE ERROR LIES BETWEEN -.5 AND +.516 LSB'S SQRT: PUSH P,CFLOAT1 SQRT.: PUSHJ P,NUMFLT JUMPG TT,SQRT.. JUMPL TT,SQR$ER ;NEGATIVE ARG IS AN ERROR POPJ P, ;ZERO ARGUMENT => ZERO ;;; POSITIVE ARGUMENT IS IN TT NOW SQRT..: MOVE R,TT ;SAVE ARG FOR LATER ASH TT,-1 ADD TT,[265116421] ;THAT'S 265116421 (KAHAN BLACK MAGIC) TLON TT,400 JRST SQRT.2 FMPRI TT,301456 ;(301456=(FSC 1.1796875 100) ;; BKPH suggested using 301456 instead of 301461 in mid 1980 (JONL 10/16/80) JRST SQRT.3 SQRT.2: FMPRI TT,300653 ;(300653)=(FSC 0.833984375 100) ;NOW TWO NEWTON ITERATIONS, MODIFIED SQRT.3: MOVE D,R FDV D,TT ;UNROUNDED DIVIDE FAD TT,D ;UNROUNDED ADD ; FSC TT,-1 SUB TT,[1000002645] ;KAHAN SEZ: INSTEAD OF DIVISION BY 2, SUBTRACT 1000002645 FDV R,TT ;UNROUNDED DIVIDE FADR TT,R ;ROUNDED ADD! FSC TT,-1 POPJ P, ;;; A FEW HINTS, PAINFULLY WORKED OUT BY GLS AND RZ: ;;; THE ASH BY -1 DIVIDES THE EXPONENT BY 2, AND MUNCHES ;;; THE MANTISSA IN A BIZARRE WAY. ;;; THE ADDITION OF 265116421 IS GUARANTEED TO CARRY ;;; INTO THE 3.9 BIT, ASSUMING A NORMALIZED INPUT. THIS ;;; WILL COMPLEMENT THE ORIGINAL LOW EXPONENT BIT. ;;; THIS IS THEN TESTED BY THE TLON, WHICH ALSO FORCES ;;; THE 3.9 BIT ON, MAKING THE NEW NUMBER NORMALIZED. ;;; THE SUBTRACTION OF 1000002645 INDEED DIVIDES BY 2, ;;; BY SUBTRACTING 1 FROM THE EXPONENT; AND THE REST DOES ;;; A WEIRD LITTLE PERTURBATION WHICH, HOWEVER, CANNOT ;;; BORROW FROM THE EXPONENT. SUBTTL LOG FUNCTION LOG: PUSH P,CFLOAT1 LOG.: PUSHJ P,NUMFLT LOG..: JUMPLE TT,LOG.ER ;NON-POSITIVE ARG IS AN ERROR MULI TT,400 HRREI TT,-201(TT) ;SAVE CHARACTERISTIC IN TT LSH D,-8. ;REDUCE ARG TO VALUE X BETWEEN 1.0 AND 2.0 TLO D,201000 MOVEI R,0 CAMN D,FPWUN ;LOG(1.0)=0.0 (ALSO FOR WHOLE POWERS OF 2 THIS SAVES TIME) JRST LOG.2 MOVE T,D ; X - SQRT(2) FSBR T,ROOT2 ; T <- ------------- FADR D,ROOT2 ; X + SQRT(2) FDVRB T,D FMPR D,D ; D <- T*T MOVEI F,3 ;MOBY APPROXIMATION TO LOG BASE 2 LOG.1: FMPR R,D FADR R,LOG.CF(F) SOJGE F,LOG.1 FMPR R,T FADR R,[0.5] LOG.2: JSP T,IFLOAT ;FLOAT CHARACTERISTIC FADR TT,R ;ADD TO LOG OF MANTISSA FMPR TT,[0.6931471806] ;MULTIPLY BY LN 2 TO GET LOG BASE E POPJ P, ROOT2: 1.4142135625 ;SQRT(2) LOG.CF: 2.885390073 ;COEFFICIENTS FOR LOG APPROXIMATION 0.9618007623 0.5765843421 0.4342597513 NUMFLT: IFE BIGNUM, JSP T,FLTSKP IFN BIGNUM, JSP T,NVSKIP IFN BIGNUM, JRST NUMFL3 JSP T,IFLOAT POPJ P, IFN BIGNUM,[ NUMFL3: PUSH P,A PUSHJ P,FLBIG JRST POPAJ ] ;END OF IFN BIGNUM SUBTTL ATAN FUNCTION ATAN: PUSH P,CFLOAT1 ATAN.: EXCH A,B PUSHJ P,NUMFLT PUSH FXP,TT MOVEI A,(B) PUSHJ P,NUMFLT POP FXP,D MOVM R,TT ;GET ABSOLUTE VALUE OF Y MOVM F,D ;GET ABSOLUTE VALUE OF X MOVEM R,ATAN.Y ;SAVE ABS(Y) MOVEM F,ATAN.X ;SAVE ABS(X) HLR D,TT ;D HAS ,, MOVEM D,ATAN.S ;SAVE THAT MESS (HAS SIGNS OF X AND Y) MOVE T,R JFCL 8,.+1 FSBR T,F ; ABS(Y)-ABS(X) FADR R,F ; T <- ----------------- FDVRB T,R ; ABS(Y)+ABS(X) FMPR R,R ; R <- T*T MOVE D,ATAN.C+7 ;MOBY APPROXIMATION MOVEI F,6 ATAN.1: FMPR D,R FADR D,ATAN.C(F) SOJGE F,ATAN.1 FMPR D,T MOVM TT,D CAMGE TT,[.7855] CAMGE TT,[.7853] JRST ATAN.3 JUMPGE D,ATAN.2 ;PATCH UP FOR WHEN RATIONAL APPROXIMATION NOT VERY GOOD MOVE D,ATAN.Y ;WE CAN USE Y/X FOR ATAN (Y/X) FDVR D,ATAN.X JRST ATAN.4 ATAN.2: MOVN D,ATAN.X FDVR D,ATAN.Y FADR D,PI%2 JRST ATAN.4 ATAN.3: FADR D,[0.7853981634] ;PI/4 ATAN.4: MOVN TT,D ;NOW WE HAVE A QUADRANT 1 RESULT (CALL IT Q) FADR TT,PI% ;PATCH-UP STUFF TO GET RIGHT QUADRANT SKIPL F,ATAN.S ; X>0 I X<0 EXCH D,TT ;-------------------------I------------------------- FSC D,1 ; D <- PI-Q I D <- Q TRNE F,400000 ; TT <- Q I TT <- PI-Q FADR TT,D ; Y>0 I Y<0 I Y>0 I Y<0 JFCL 8,ATAN.7 ;------------I------------I------------I------------ POPJ P, ; TT<-Q I TT<-2*PI-Q I TT<-PI-Q I TT<-PI+Q PI%: 3.1415926536 ;A WELL-KNOWN NUMBER ATAN.C: 0.9999993329 ;COEFFICIENTS FOR ATAN APPROXIMATION -0.3332985605 0.1994653599 -0.139085335 0.0964200441 -0.0559098861 0.0218612288 -0.004054058 SUBTTL EXP FUNCTION EXP: PUSH P,CFLOAT1 EXP.: JSP T,FLTSKP JSP T,IFLOAT EXP..: SETZ R, MOVEM TT,EXP.S ;SAVE SIGN OF ARG ON PDL MOVM TT,TT ;GET ABSOLUTE VALUE OF ARG CAMLE TT,[88.0] ;WAS REQUESTED POWER > 88.0? JRST EXP.A ;YES, CAN'T REPRESENT SOMETHING THIS BIG FMPR TT,[0.4342944819] ;LOG BASE 10. OF E ;FROM NOW ON WE DO 10.^X, NOT E^X MOVE F,FPWUN ;F HOLDS 10.^ CAMG TT,FPWUN ;IF ARG <=1.0 GO DO RATIONAL APPROXIMATION JRST EXP.RX MULI TT,400 ASHC D,-243(TT) ;D GETS INTEGER PART OF ARG ; CAIG D,43 ;THIS IS OLD CHECK, JONL SAYS OK TO ALLOW JRST EXP.1 ; LARGER RANGE EXP.A: SKIPGE TT,EXP.S ;TOO LARGE - RESULT CAN'T BE REPRESENTED TDZA TT,TT JRST EXP.ER POPJ P, ;NEGATIVE ARG PRODUCES ZERO (UNDERFLOW) EXP.1: CAIG D,7 ;SKIP IF INTEGER PART OF ARG > 7 JRST EXP.2 LDB T,[030300,,D] ;GET TOP 3 BITS OF 6 BIT INTEGER PART ANDI D,7 ;AND THEM OUT OF D MOVE F,INTLG(T) ;F GETS (10.^T)^8. = 10.^(T*8.) FMPR F,F FMPR F,F FMPR F,F EXP.2: FMPR F,INTLG(D) ;MULTIPLY F BY APPROPRIATE 10.^D (0<=D<=7) LDB TT,[103300,,R] ;NOW GET FRACTION PART OF ARG TLO TT,177000 ;THIS STRANGENESS FLOATS FADR TT,TT ; AND NORMALIZES THE FRACTION EXP.RX: MOVEI T,6 ;MOBY APPROXIMATION SKIPA R,EXP.CF+6 EXP.3: FADR R,EXP.CF(T) FMPR R,TT SOJGE T,EXP.3 FADR R,FPWUN FMPR R,R FMPR F,R ;MULTIPLY FRACTION APPROXIMATION BY 10.^ MOVE TT,FPWUN SKIPL EXP.S SKIPA TT,F ;IF ARG>0, RETURN RESULT FDVR TT,F ;IF ARG<0, RETURN 1.0/RESULT POPJ P, EXP.CF: 1.151292776 ;COEFFICIENTS FOR EXP APPROXIMATION 0.6627308843 0.2543935748 0.07295173666 0.01742111988 2.55491796^-3 9.3264267^-4 FPWUN: ;FLOATING POINT 1.0 INTLG: 1.0 ;TABLE OF 10.^X FOR INTEGRAL 0<=X<=7 REPEAT 7, 1.0^<.RPCNT+1> C1.0E5=FPWUN+5 PGTOP ARI,[ARITHMETIC SUBROUTINES]