Monte Davidoff's floating point routines for early Microsoft BASIC used Chebyshev Modified Taylor series for EXP(x)
. There's a very helpful disassembly of the TRS-80 MC-10 ROM here: http://www.roust-it.dk/coco/mc10/romlist.txt. It's 6800 assembly, and the whole commented routine (using the same constants) is:
TBLF59B FCB $81,$38,$AA,$3B,$29 ;1.44269504 (CF) correction factor for EXP function
TBLF5A0 FCB $07 ;eight coeff's... tchebyshev modified taylor series coeffs for exp(x)
FCB $71,$34,$58,$3E,$56 ;0.00002150 1/(7! * CF^7)
FCB $74,$16,$7E,$B3,$1B ;0.00014352 1/(6! * CF^6)
FCB $77,$2F,$EE,$E3,$85 ;0.00134226 1/(5! * CF^5)
FCB $7A,$1D,$84,$1C,$2A ;0.00961402 1/(4! * CF^4)
FCB $7C,$63,$59,$58,$0A ;0.05550513 1/(3! * CF^3)
FCB $7E,$75,$FD,$E7,$C6 ;0.24022638 1/(2! * CF^2)
FCB $80,$31,$72,$18,$10 ;0.69314719 1/(1! * CF^1)
FCB $81,$00,$00,$00,$00 ; 1.0
; --- EXP function ---
LBLF5C9:
FNC_EXP LDX #TBLF59B ;Get correction factor
BSR LBLF604 ;Multiply FPA0 by X
JSR LBLF26C ;pack fpa0 and store in fpa3
LDAA ramC9 ;get exponent of fpa0 and compare to max value
CMPA #$88 ; (128)
BLO LBLF5DA ;br if fpa0 < 128
LBLF5D7 JMP LBLF190 ;set fpa0 = 0 or ?OV ERROR
LBLF5DA JSR FNC_INT ;convert fpa0 to integer
LDAA ram80 ;get least significant byte of integer
ADDA #$81 ; =127?
BEQ LBLF5D7 ; ?OV ERROR
DECA ; adds bias of 80 (since 81 used above)
PSHA ;save exponent on stack
LDX #TBL00BA ;point (x) to FPa3
JSR LBLEF72 ;subtract fpa0 from (x)
LDX #TBLF5A0 ;point x to coeffs
BSR LBLF607 ;eval polynomial for frac part
CLR ramDC ;force mantissa to be positive
PULA
JSR LBLF179 ;calc exp of new fpa0 by adding exps of integer and frac'l parts.
RTS
LBLF5F8 STX ramDE
JSR LBLF26C
BSR LBLF604
BSR LBLF609
LDX #TBL00BA
LBLF604 JMP LBLF0EF
LBLF607 STX ramDE
LBLF609 JSR LBLF267
LDX ramDE
LDAB ,X
STAB ramCF
INX
STX ramDE
LBLF615 BSR LBLF604
LDX ramDE
LDAB #$05
ABX
STX ramDE
JSR LBLEF7D
LDX #ramBF
DEC ramCF
BNE LBLF615
RTS
Your listing missed off the all-important correction factor, 1.44269504 or 1/LOG(2)
in BASIC. The coefficients take the form
1 / ( n! * (1/LOG(2))^n)
msbasic.zip
). Since I know none of the approximate exponential algorithms, I'm not much help beyond that. – Tommy Jul 22 '18 at 21:14float.s
inside that zip file; the quoted table is from line 1758 and the implementation ofEXP
is immediately below) – Tommy Jul 23 '18 at 2:06