5

In 1980's ROM (Apple IIe, Commodore 64, VIC-20, ...) which algorithm is used to compute exp(x), and where do the coefficients below come from? (Chebyshev, Remez, Pade, ...)

.byte   $71,$34,$58,$3E,$56 ; 2.14987637E-5
.byte   $74,$16,$7E,$B3,$1B ; 1.43523140E-4
.byte   $77,$2F,$EE,$E3,$85 ; 1.34226348E-3
.byte   $7A,$1D,$84,$1C,$2A ; 9.61401701E-3
.byte   $7C,$63,$59,$58,$0A ; 5.55051269E-2
.byte   $7E,$75,$FD,$E7,$C6 ; 2.40226385E-1
.byte   $80,$31,$72,$18,$10 ; 6.93147186E-1
.byte   $81,$00,$00,$00,$00 ; 1.00000000

PS: See https://math.stackexchange.com/questions/2858662/expx-approximation-in-old-1980s-computer-rom in Mathematics.SE

  • 3
    Such as it may help others better at reading this stuff than I, those coefficients are from Microsoft BASIC so I think the whole question is about Microsoft BASIC, and commented, disassembled source can be found at pagetable.com/?p=46 (serch for msbasic.zip). Since I know none of the approximate exponential algorithms, I'm not much help beyond that. – Tommy Jul 22 '18 at 21:14
  • (oh, and check out float.s inside that zip file; the quoted table is from line 1758 and the implementation of EXP is immediately below) – Tommy Jul 23 '18 at 2:06
  • @Tommy : thank you, I will look at these informations. – jpcohet Jul 23 '18 at 18:39
7

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)
  • 1
    Incidentally, if you want to see a Chebyshev generator (that does not take this corrected approach), Richard Russell posted his code that he's used with generations of BBC BASIC ports here. You'll need a BBC BASIC interpreter to run it, though, as it uses the built in function evaluator EVAL() – scruss Jul 23 '18 at 3:08
  • thank you very much. I will look at these informations. My main expectation is to know how these Chebyshev coefficients are compute to be use in this EXP(x) routine, because all my tries with traditional method have failed. – jpcohet Jul 23 '18 at 18:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.