.ref-Research-V6/usr/source/s3/exp.s

.globl exp, _exp
/
ldfps = 170100^tst
stfps = 170200^tst
ldexp = 176400^movif
stexp = 175000^movfi
/
/       exp accepts its argument and returns its result
/       in fr0.  The carry bit is set if the result overflows.
/       The coefficients are #1067 from Hart & Cheney.
/
/       movf    arg,fr0
/       jsr     pc,exp
/       movf    fr0,result
/
_exp:
        mov     r5,-(sp)
        mov     sp,r5
        mov     4(r5),fr0
        jsr     pc,exp
        mov     (sp)+,r5
        rts     pc

exp:
        stfps   -(sp)
        ldfps   $200            /di mode
        movf    fr2,-(sp)
        movf    fr1,-(sp)
        tstf    fr0
        cfcc
        bne     1f
        movf    $one,fr0        /exp(0)
        clc
        br      out
1:
        modf    log2e,fr0       /exp(x) = 2^(x*log2(e))
        cfcc
        bmi     2f
        movfi   fr1,-(sp)       /save integer part
        subf    $half,fr0
        br      3f
2:
        movfi   fr1,-(sp)
        dec     (sp)
        addf    $half,fr0
3:
        movf    fr0,fr1         / -.5 < x < +.5
        mulf    fr1,fr1         /arg**2
/
        movf    P2,fr2
        mulf    fr1,fr2
        addf    P1,fr2
        mulf    fr1,fr2
        addf    P0,fr2
        mulf    fr2,fr0         /xP(x**2)
/
        movf    fr1,fr2
        addf    Q1,fr2
        mulf    fr1,fr2
        addf    Q0,fr2          /Q(x**2)
/
        movf    fr2,fr1
        subf    fr0,fr1
        addf    fr2,fr0
        divf    fr1,fr0         /(Q+xP)/(Q-xP)
        mulf    sqrt2,fr0
/
        stexp   fr0,-(sp)
        add     (sp)+,(sp)
/
        cmp     (sp),$177
        ble     2f
        tst     (sp)+
        movf    big,fr0         /overflow
        sec
        br      1f
2:
        cmp     (sp),$-177
        bge     2f
        tst     (sp)+
        clrf    fr0             /underflow
        clc
        br      1f
2:
        ldexp   (sp)+,fr0
        clc
1:
out:
        movf    (sp)+,fr1
        movf    (sp)+,fr2
        ldfps   (sp)+
        rts     pc
/
/
        .data
P0:      42675;  36404;  77563;  46675
P1:      41241; 116724; 114237;  60333
P2:      36675;  27102; 125560; 136652
Q0:      43210; 100661;  76072;  62453
Q1:      42151;  27450;  75350; 112503
log2e:   40270; 125073;  24534;  13761
sqrt2:   40265;  02363;  31771; 157144
half    = 40000
one     = 40200
/
big:     77777; 177777; 177777; 177777
/
/ P0 = .15139 06799 05433 89158 94328 d4
/ P1 = .20202 06565 12869 27227 886   d2
/ P2 = .23093 34775 37502 33624       d-1
/
/ Q0 = .43682 11662 72755 84984 96814 d4
/ Q1 = .23318 42114 27481 62379 0295  d3
/ Q2 = .1                             d1
/
/ log2e = 1.44269 50408 88963 40735 99246
/ sqrt2 = 1.41421 35623 73095 04880 16887