SKP
*               ******************************
*               SUBROUTINE TO CALCULATE SIN(X)
*               ******************************
*
*               CALLED BY A JMP ESIN WITH THE ARGUMENT
*               IN FLOATING RADIANS IN THE REGISTERS.
*               THE FLOATING RESULT IS RETURNED IN A 8 B
*
ECOS  JSB .FADA,I
      DEF PIBY2
ESIN  JSB .FMPA,I
      DEF TOPI
      STA XTEMP
      STB XTEMP+1   X=2*X/PI
      JSB .FADA,I
      DEF K1
      JSB .PWR2
      DEC -2
      JSB .IENT
      JMP TRGER-1   ERROR IF EXPONENT >= 15
      JSB FLOAT
      JSB .FMPA,I
      DEF MM4
      JSB .FADA,I
      DEF XTEMP
      STA XTEMP
      STB XTEMP+1   X=X-4*ENTIER((X+1)/4)
      JSB .FSBA,I
      DEF K1
      SSA           X<1 ?
      JMP PAST      YES
      LDA K2
      LDB K2+1
      JSB .FSBA,I
      DEF XTEMP
      STA XTEMP
      STB XTEMP+1   X=2-X
PAST  LDA XTEMP
      LDB XTEMP+1
      JSB .FMPA,I
      DEF XTEMP
      JSB .PWR2
      DEC 1
      JSB .FSBA,I
      DEF K1
      JSB .CHEB
      DEF COEF1
      JSB .FMPA,I
      DEF XTEMP
      JMP FR12A,I   ANS=X+CHEBYI2*X**2-1)
*
TOPI  DEC .636619772     2/PI
MM4   DEC -4.
COEF1 DEC 1.18496E-6
      DEC -1.365875E-4
      DEC 9.118016E-3
      DEC -.2852615692
      DEC 2.5525579248
      OCT 0
      SPC 5
*               *****************************
*               SUSROUTINE TO COMPUTE RND(X)
*               *****************************
*
*               THE ARGUMENT OF RND IS A DUMMY ONE
*               THE ROUTINE GENERATES A RANDOM NUMBER
*               IN THE A & B REGISTERS
*
*  R=X/M, X=C*X MOD M, M=2^31, C=2^15 + 3
*
ERND  CLA           X IS INITIALLY 1
      STA EXP       INITIALIZE EXPONENT
      LDA XH        COMPUTE
      ALS             HIGH
      ADA XH            PART
      LDB XL              2*XH
      CLE,ERB               + XH +
      ADA 1                   XL*2^15
      LDB XL
      RBL,CLE,SLB,ERB    ADD XL[15] TO
      INA                  (A) (FROM 2*XL)
      CLE,ELB            2*XL
      ADB XL               + XL
      ELA,CLE,SLA,ERA    ADD OVERFLOW
      CLE,INA              TO (A)
      ADB FLGBT     ADD IN TRAILING BIT OF XL*2^15
      SEZ           ADD OVERFLOW
      INA             TO (A)
      ELA,CLE,ERA   ERASE A[15]
      STA XH        STORE
      STB XL          INTEGER
      JSB .PACK     NORMALIZE AND PACK
      JMP FR12A,I
      SKP
*               ******************************
*               SUBROUTINE TO CALCULATE SQR(X)
*               ******************************
*
*               CALLED BY A JMP ESQR WITH THE ARGUMENT
*               IN FLOATING POINT FORM IN THE REGISTERS.
*               THE FLOATING RESULT IS RETURNED IN A & B
*
ESQR  SZA,RSS       X=0 ?
      JMP FR12A,I   YES, ANS = 0
      SSA           X<O ?
      JSB ERROR     YES ERROR
SQRER STA XTEMP
      JSB .FLUN
      SLA,ARS       EXP(X) ODD?
      JMP ODD
      ADA M1
      STA SBOX      SBOX=EXPO(X)/2-1
      STB XTEMP+1   LOW MANTISSA/2
      LDA XTEMP
      JSB .FMPA,I
      DEF SA2
      JSB .FADA,I
      DEF SB2       Y=SB2+SA2*X
BTH2  STA YTEMP
      STB YTEMP+1
      LDA XTEMP
      LDB XTEMP+1
      JSB .FDVA,I
      DEF YTEMP
      JSB .FADA,I
      DEF YTEMP
      JSB .PWR2
      DEC -1
      STA YTEMP
      STB YTEMP+1   Y=(Y+X/Y)/2
      LDA XTEMP
      LDB XTEMP+1
      JSB .FDVA,I
      DEF YTEMP
      JSB .FADA,I
      DEF YTEMP
      STA YTEMP
      STB YTEMP+1
      JSB .PWR2
SBOX  OCT 0
      JMP FR12A,I   ANS=(P+F/P)*2**SBOX
ODD   STA SBOX      SBOX = EXP(X)/2
      ADB MSK0
      STB XTEMP+1   LOW MANTISSA/2
      LDA XTEMP
      JSB .FMPA,I
      DEF SA1
      JSB .FADA,I
      DEF SB1
      JMP BTH2      Y=SB1+SA1*X
*
SA1   DEC .875
SA2   DEC .578125
SB1   DEC .27863
SB2   DEC .421875
      SPC 5
*               ******************************
*               SUBROUTINE TO CALCULATE INT(X)
*               ******************************
*
*               CALLED BY A JMP EINT WITH THE ARGUMENT
*               IN FLOATING POINT FORM IN THE REGISTERS.
*               ENTIER(X) IS FLOATED AND RETURNED IN A & B
*
*               CALCULATES THE FLOATING POINT EQUIVALENT
*               OF ENTIER(X)
*
EINT  STB SBOX      SAVE LOW WORD
      LDB .31       PRESET BIAS FOR
      STB EXP         BINARY POINT
      LDB SBOX      24-BIT
      JSB IFIX        INTEGER?
      JMP EINT1     NO
      JSB .PACK     YES, PACK IT
      JMP FR12A,I
EINT1 LDA GETCR     RETRIEVE
      LDB SBOX        NUMBER
      JMP FR12A,I
      SKP
*               ******************************
*               SUBROUTINE TO CALCULATE LOG(X)
*               ******************************
*
*               CALLED BY A JMP ELOG WITH THE ARGUMENT
*               IN FLOATING POINT FORM IN THE REGISTERS.
*               THE FLOATING RESULT IS RETURNED IN A & B
*
ELOG  JSB .LOG
      JMP FR12A,I
.LOG  NOP
      SZA,RSS       NON-ZERO ARGUMENT?
      JMP .LOG1     NO
      SSA           YES, POSITIVE ARGUMENT?
      JSB ERROR     NO
LOGER STA XTEMP     YES
      JSB .FLUN
      STB XTEMP+1
      JSB FLOAT
      STA YTEMP
      STB YTEMP+1   Y=EXPO(X)
      LDA XTEMP
      LDB XTEMP+1
      JSB .FADA,I
      DEF R22
      STA UTEMP
      STB UTEMP+1   U=X+SQR(0.5)
      LDA XTEMP
      LDB XTEMP+1
      JSB .FSBA,I
      DEF R22
      JSB .FDVA,I
      DEF UTEMP
      STA UTEMP
      STB UTEMP+1   U=(X-SQR(0.5))/U
      JSB .FMPA,I
      DEF UTEMP
      JSB .FSBA,I
      DEF CCC
      STA XTEMP
      STB XTEMP+1
      LDA MB
      LDB MB+1
      JSB .FDVA,I
      DEF XTEMP
      JSB .FADA,I
      DEF AAA
      JSB .FMPA,I
      DEF UTEMP
      JSB .FSBA,I
      DEF HALF
      JSB .FADA,I
      DEF YTEMP
      JSB .FMPA,I
      DEF LE2
      JMP .LOG,I    ANS=LOG(2)*(EXPO(X)-0.5+U*
*                                   (A-B/X))
.LOG1 JSB ERROR     LOG OF ZERO
LNZR  LDA MNEG      RETURN
      LDB B776        NEGATIVE
      JMP .LOG,I        INFINITY
*
R22   DEC .707106781     SQR(0.5)
LE2   DEC .6931471806    LOG BASE E OF 2
AAA   DEC 1.2920070987
MB    DEC -2.6398577035
CCC   DEC 1.6567626301
      SPC 5
*               ******************************
*               SUBROUTINE TO COMPUTE SGN(X)
*               ******************************
*
*               CALLED BY A JMP ESGN WITH THE ARGUMENT
*               IN FLOATING POINT FORM IN THE REGISTERS.
*               ON RETURN A & B CONTAIN THE FOLLOWING:
*
*               IF X>0 THEN +1.
*               IF X=0 THEN  0
*               IF X<0 THEN -1.
*
ESGN  CLB
      SZA,RSS       ZERO?
      JMP FR12A,I   YES
      SSA,RSS       NO, POSITIVE?
      LDB .2        YES, SET EXPONENT
      LDA FLGBT     LOAD MANTISSA
      SZB           POSITIVE?
      RAR           YES, CORRECT MANTISSA
      JMP FR12A,I
      SKP
*               ******************************
*               SUBROUTINE TO CALCULATE EXP(X)
*               ******************************
*
*               CALLED BY A JMP EEXP WITH THE ARGUMENT
*               IN FLOATING POINT FORM IN THE REGISTERS.
*               THE FLOATING RESULT IS RETURNED IN A & B
*
EEXP  JSB .EXP
      JMP FR12A,I
.EXP  NOP
      JSB .FMPA,I
      DEF L2E
      STA XTEMP
      STB XTEMP+1   X=Z*LOG2(E)
      JSB .IENT
      JMP .EXP1
      STA INTE      INTE = ENTIER(X)
      JSB FLOAT
      STA YTEMP
      STB YTEMP+1   Y=ENTIER(X)
      LDA INTE
      ADA M124
      SSA,RSS       INTE >=124 ?
      JMP EXPER-1   YES,ERROR
      ADA .244      INTE <-120 ?
      SSA
      JMP ZERE      YES,ANS=0
      LDA XTEMP
      LDB XTEMP+1
      JSB .FSBA,I
      DEF YTEMP
      STA XTEMP
      STB XTEMP+1   X=X-ENTIER(X)
      JSB .FMPA,I
      DEF XTEMP
      STA UTEMP
      STB UTEMP+1   U=X**2
      JSB .FADA,I
      DEF AAAA
      STA YTEMP
      STB YTEMP+1   Y=X**2+AAAA
      LDA BBBB
      LDB BBBB+1
      JSB .FDVA,I
      DEF YTEMP
      STA YTEMP
      STB YTEMP+1   Y=BBBB/Y
      LDA CCCC
      LDB CCCC+1
      JSB .FMPA,I
      DEF UTEMP
      JSB .FADA,I
      DEF DDDD
      JSB .FSBA,I
      DEF XTEMP
      JSB .FSBA,I
      DEF YTEMP
      STA YTEMP
      STB YTEMP+1   Y=-X+DDDD+CCCC*X**2-Y
      LDA XTEMP
      LDB XTEMP+1
      JSB .FDVA,I
      DEF YTEMP
      JSB .FADA,I
      DEF HALF
      ISZ INTE
      NOP
      JSB .PWR2
INTE  OCT 0
      JMP .EXP,I    ANS=(0.5+X/Y)*2**INTE
ZERE  CLA
      CLB
      JMP .EXP,I    ANS=0
.EXP1 LDA X2TMP
      SSA
      JMP ZERE
      JSB ERROR
EXPER LDA INF
      LDB M2
      JMP .EXP,I
*
M124  DEC -124
.244  DEC 244
AAAA  DEC 87.417497202
BBBB  DEC 617.9722695
CCCC  DEC .03465735903
DDDD  DEC 9.9545957821
L2E   DEC 1.4426950409
      SKP
*               ******************************
*               SUBROUTINE TO COMPUTE CHEBY(X)
*               ******************************
*
*               CALLING SEQUENCE:
*
*               LDA X         (FLOATING)
*               LDB X+1
*               JSB .CHEB     (RESULT FLOATING)
*               DEF C         (TABLE OF CHEBY.COEFFS.,FLOATING)
*
.CHEB NOP
      JSB .FMPA,I
      DEF K2
      STA X2TMP
      STB X2TMP+1   X2 =X*2
      LDB .CHEB,I
      STB CTMP      C POINTS TO COEFFICIENT TABLE
      LDA 1,I
      INB
      LDB 1,I       GET FIRST COEFF
      STA DTMP
      STB DTMP+1    D=C(N)
      CLA
      STA BTMP
      STA BTMP+1    B=0
LOPC  ISZ CTMP
      ISZ CTMP      N=N-1
      LDA CTMP,I
      SZA,RSS       C(N)=0 ?
      JMP COUT      ZERO FLAGS END OF TABLE
      LDA BTMP      NO
      LDB BTMP+1
      STA ATMP
      STB ATMP+1    A=B
      LDA DTMP
      LDB DTMP+1
      STA BTMP
      STB BTMP+1    B=D
      JSB .FMPA,I
      DEF X2TMP
      JSB .FSBA,I
      DEF ATMP
      JSB .FADA,I
      DEF CTMP,I
      STA DTMP
      STB DTMP+1    D=C(N) -A+B*X2
      JMP LOPC
COUT  LDA DTMP
      LDB DTMP+1
      JSB .FSBA,I
      DEF ATMP
      JSB .FMPA,I
      DEF HALF
      ISZ .CHEB
      JMP .CHEB,I   ANS=(D-A)/2
*
X2TMP BSS 2
ATMP  BSS 2
BTMP  BSS 2
CTMP  BSS 2
DTMP  BSS 2
      SPC 5
*               ********************************************
*               SUBROUTINE TO COMPUTE THE ENTIER OF A NUMBER
*               WHOSE EXPONENT IS LESS THAN 15
*               ********************************************
*
*               CALLING SEQUENCE:
*
*               LDA X         (FLOATING)
*               LDA X+1
*               JSB .IENT.    (RESULT INTERGER)
*               JSB ERROR     (EXIT IF EXPO(X)>14)
*
*
.IENT NOP
      STA X2TMP     STORE HIGH PART
      LDA 1         MOVE LOW PART TO A
      AND MSK0      ISOLATE EXPONENT
      SLA,RAR
      JMP *+4       IF NEGATIVE OK
      ADA M15
      SSA,RSS       EXPO(X) > 14
      JMP .IENT,I   YES, ERROR RETURN
      ISZ .IENT     NO BUMP RETURN POINT
      LDA X2TMP     RESTORE HIGH PART
      JSB IFIX      CALL ENTIER
      NOP
      LDA 1         PUT RESULT INTO (A)
      JMP .IENT,I
      SKP
*               ******************************
*               SUBROUTINE TO FLOAT AN INTEGER
*               ******************************
*
*               CALLED BY JSB FLOAT WITH INTEGER IN A
*               THE FLUATING POINT EQUIVALENT IS RETURNED
*               IN A & B
*
FLOAT NOP
      LDB .15
      STB EXP
      CLB
      JSB .PACK
      JMP FLOAT,I
      SPC 5
*               ****************************************
*               SUBROUTINE TO MULTIPLY BY A POWER OF TWO
*               ****************************************
*
*               CALLING SEQUENCE
*
*               LDA X         (FLOATING)
*               LDB X+1
*               JSB .RWR2     (RESULT FLOATING)
*               DEC N         (INTEGER POWER)
*
*               RETURNS WITH X*2^N IN A&B
*               NO CHECK IS MADE FOR EXPONENT
*               OVERFLOW OR UNDERFLOW
*
.PWR2 NOP
      SZA,RSS       X=0 ?
      JMP .RET      YES, ANS=0
      STA X2TMP
      JSB .FLUN
      STB X2TMP+1
      ADA .PWR2,I
      RAL
      AND MSK0      NEW EXPO = (OLD EXPO) +N
      STA 1
      ADB X2TMP+1   KEEP OLD MANTISSA
      LDA X2TMP
.RET  ISZ .PWR2
      JMP .PWR2,I
      SPC 2
TT1   EQU .FDV
TT2   EQU IDIV
TT3   EQU TEMPS+4
TT4   EQU TEMPS+5
FFLAG EQU SBOX