\ Win32Forth version 6.14 has a bug in FEXPM1 that is fixed below.
\ This code was written by George Hubert and should only be used with Win32Forth 6.14.

\ Include this file prior to including novice.4th so the hyperbolic trig functions will work properly.
\ This presumably won't be necessary in future versions of Win32Forth.


code FEXPM1     ( fs: r1 -- r2 )          \ ANSI      Floating ext
\ *G Raise e to the power r1 and subtract one, giving r2.
\ *P This function allows accurate computation when its arguments are close to zero, and
\ **  provides a useful base for the standard exponential functions. Hyperbolic functions
\ **  such as cosh(x) can be efficiently and accurately implemented by using FEXPM1;
\ **  accuracy is lost in this function for small values of x if the word FEXP is used.
                fstack-check_1
                >fpu
                fldl2e                                             \ 2
                fmulp st(1), st(0)                                 \ 1
                fld1                                               \ 2
                fcom    st(1)                                      \ 2
                fstsw  ax                                          \ 2
                sahf                                               \ 2
                jbe     short L$1       \ arg > 1                  \ 2
                fchs                                               \ 2
                fcomp   st(1)                                      \ 1
                fstsw  ax                                          \ 1
                sahf                                               \ 1
                jnc     short L$2       \ arg <= -1                \ 1
                je      short L$4       \ is NAN                   \ 1
                f2xm1                                              \ 1
                jmp     short L$4
L$1:            fstp    st(0)                                      \ 1
L$2:            fxam
                fstsw  ax
                and     ax, # FPU_STATUS_CCF_MASK
                cmp     ax, # FPU_STATUS_CCF_INFINITY
                je      short L$4
                cmp     ax, #  0x700 \ FPU_STATUS_CCF_INFINITY FPU_STATUS_CCF_MASK_1 or
                jne     short L$3
                fstp    st(0)
                fld1
                fchs
                jmp     short L$4
L$3:            fld     st(0)           \ duplicate exponent       \ 2
                frndint                 \ take integer part        \ 2
\               jmp     short L$4                                           \ these three lines replaced
\               fld     st(0)           \ duplicate exponent       \ 2      \ by the previous three lines
\ L$3:          frndint                 \ take integer part        \ 2      \ per George Hubert
                fsub    st(1), st       \ get fractional part      \ 2
                fld1                                               \ 3
                fscale                  \ 2**int                   \ 3
                fstp    st(1)           \ remove unneeded part     \ 2
                fxch    st(1)           \ frac                     \ 2
                f2xm1                   \ (2**frac) - 1            \ 2
                fld1                    \ 1.0                      \ 3
                faddp   st(1), st       \ 2**frac                  \ 2
                fmulp   st(1), st       \ 2**(int + frac)          \ 1
                fld1
                fsubp   st(1), st
L$4:            fpu>
                float;