Category : Forth Source Code
Archive   : VPSFP102.ZIP
Filename : FPTRANS.SEQ

\ FILE: FPTRANS.SEQ TRANSCENDENTAL FUNCTIONS
\ VP Planner floating point for F-PC.

/*
Kent Brothers and Stephenson Software Ltd. are making these floating
point routines available to the Forth language community to further the
growth of Forth as an application development language. It is licensed
for incorporation into programs (with or without modifications) with no
fee as long as the following notice is displayed on screen and in the
documentation:

VP-Planner floating point routines copyright 1985-1989
Kent Brothers, Stephenson Software Ltd.

The source code may be distributed free of charge provided this entire
message (including the notice) is left at the top.
*/

CR .( Loading FPTRANS.SEQ ... )


\ KMB 850818 ATAN2 corrected to clean up stack.
\ KMB 850819 All functions modified to handle NANs.
\ KMB 850820 MAXEXPARG and MINEXPARG revised for extended precision.
\ KMB 850902 Changed EXP and EXPC to return 0 and -1 when underflow occurs.
\ JWB 860225 Vectored kernel transendental functions and constants.
\ KMB 861109 EXPBC vectored and EXP and ^^ de-vectored.
\ KMB 861113 ATAN2 changed to call ATAN instead of (ATAN).
\ DJL 870206 ?8087 changed to load HFP.OVL from path.
\ DJL 870506 (EXP) fix for XF' < 0.
\ DJL 870506 |(LN) split from |LN to match with HFPPRIM2.
\ DJL 870706 Changed HFP.OVL to VP87.OVR
\ DJL 870722 Changed VP87.OVR to VPP87.OVR
\ JWB 890331 Converted to F-PC
\ JWB 890404 Made non user words headerless
\ JWB 890405 Renamed trig functions to confrom to FVG fp standard.
\ JWB 890405 Set up some aliases for LN EXP LOG10 to have FVG names
\ available for these functions.

FORTH DEFINITIONS
FLOATING DECIMAL

HEADERLESS
\ Leave r2 on floating point stack, the square root of r1.
: |SQRT ( F: r1 -- r2 ) ( KMB 85 08 19 )
?NONAN1
IF FDUP F0<
IF FDROP FERR
ELSE FDUP F0= NOT
IF EXP/FRAC 2 /MOD SWAP IF -1 2SCALE 1+ THEN
1. 5 0 DO FOVER FOVER F/ F+ 2./ LOOP FSWAP FDROP 2SCALE
THEN
THEN
THEN ;

' |SQRT IS FSQRT


2. FSQRT FCONSTANT |SQRT2 ' |SQRT2 IS SQRT2

: FP# ( NE -- ) ( F: r1 r2 r3 r4 -- r5 ) ( KMB 85 06 10 )
2 0 DO 1.E-5 F* F+ LOOP 1.E-5 F* FSWAP FSIGN FABS F+ ( NE FS | r )
0< IF FNEGATE THEN ?DUP
IF DUP 0<
IF 10. ABS 1- ?DUP IF 0 DO 10.* LOOP THEN F/
ELSE 10. 1- ?DUP IF 0 DO 10.* LOOP THEN F*
THEN
THEN ;

\ : FLIT] ] COMPILE FLIT XF, ; IMMEDIATE ( KMB 85 06 11 )

\ : FP#] ( NE -- ) ( F: r1 r2 r3 r4 r5 -- ) ( KMB 85 06 11 )
\ FP# [COMPILE] FLIT] ; IMMEDIATE

\ The following polynomial and rational approximations are taken from
\ "Computer Approximations" by J.F. Hart et al. (1968).

: (LN(1+X)) ( F: r1 -- r2 ) ( KMB 85 06 11 )
\ Evaluates natural logarithm of 1+X1 to 19+ significant digits for
\ 1/2^.5 <= 1+X1 <= 2^.5. Approximation 2705.
FDUP 2. F+ F/ FDUP FDUP F*
FDUP
[ .42108 .73712 .17979 .7145 0 FP# ] FLITERAL FOVER F*
[ -.96376 .90933 .68686 .59324 1 FP# ] FLITERAL F+ FOVER F*
[ .30957 .29282 .15376 .50062 2 FP# ] FLITERAL F+ F*
[ -.24013 .91795 .59210 .50987 2 FP# ] FLITERAL F+
FSWAP FDUP
[ -.89111 .09027 .93783 .12337 1 FP# ] FLITERAL F+ FOVER F*
[ .19480 .96607 .00889 .73052 2 FP# ] FLITERAL F+ F*
[ -.12006 .95897 .79605 .25472 2 FP# ] FLITERAL F+
F/ F* ;

.69314 .71805 .59945 .30942 0 FP# FCONSTANT |LN2 ' |LN2 IS LN2


: LN(1+X) ( F: r1 -- r2 ) ( KMB 85 08 19 )
?NONAN1
IF FDUP -1. F>
IF FSIGN 0< DUP IF FDUP 1. F+ F/ FNEGATE THEN ( FS | r3 )
FDUP 1. F+ EXP/FRAC FDUP SQRT2 F> IF 1+ 2./ THEN ( FS N | r3 r4 )
DUP IF ( r3 not small ) 1. F- FSWAP THEN FDROP
(LN(1+X)) ( FS N | r5 )
S>F LN2 F* F+ IF FNEGATE THEN
ELSE FDROP FERR
THEN
THEN ;

: |(LN)
EXP/FRAC FDUP SQRT2 F> IF 2./ 1+ THEN ( N F | r )
1. F- (LN(1+X)) S>F LN2 F* F+ ;

' |(LN) IS (LN)

HEADERS

\ Leave r2 = ln(r1) on floating point stack.
: LN ( F: r1 -- r2 )( KMB 850819 )( DJL 870506 )
?NONAN1
IF FSIGN ?DUP
IF 0>
IF (LN)
ELSE FDROP FERR ( should be replaced by -infinity )
THEN
ELSE FDROP FERR
THEN
THEN ;


10. LN FCONSTANT LN10

: LOG10 LN LN10 F/ ; ( F: r1 -- r2 )

LN2 16384. F* FCONSTANT MAXEXPARG
LN2 -16382. F* FCONSTANT MINEXPARG

HEADERLESS
: |EXPBC ( F: r1 -- r2 ) ( KMB 85 08 19 )
\ Evaluates 2^X1-1 for 19+ significant digits for 0 <= r1 < .5.
\ Approximation 1324.
?NONAN1
IF FDUP FDUP F*
FDUP
[ .60613 .30790 .74800 .42575 2 FP# ] FLITERAL FOVER F*
[ .30285 .61978 .21164 .59206 5 FP# ] FLITERAL F+ F*
[ .20802 .83036 .50596 .27129 7 FP# ] FLITERAL F+
FSWAP FDUP
[ .17492 .20769 .51057 .14559 4 FP# ] FLITERAL F+ FOVER F*
[ .32770 .95471 .93281 .18053 6 FP# ] FLITERAL F+ F*
[ .60024 .28040 .82517 .36653 7 FP# ] FLITERAL F+
F/ F* FDUP FNEGATE 1. F+ F/ 2.*
THEN ;

' |EXPBC IS EXPBC ( KMB 861109 )

: (EXP) ( -- N )( F: r1 -- r2 )( KMB 850902 )( DJL 870506 )
LN2 F/ INTFRAC FSWAP F>D DROP ( N' | rF )
FDUP FSIGN 0<
IF -.5 F< IF 1. F+ 1- 0 ELSE FABS -1 THEN
ELSE .5 F> ( ANDIF) DUP IF DROP 1. F- FABS 1+ -1 THEN
THEN ( N F | |XF'| )
EXPBC IF FDUP 1. F+ F/ FNEGATE THEN ; ( N | 2^XF'-1 )

HEADERS

: EXPC ( F: r1 -- r2 ) ( KMB 85 09 02 )
?NONAN1
IF FDUP MINEXPARG F>
IF FDUP MAXEXPARG F<
IF FSIGN
IF (EXP) ?DUP IF 1. F+ 2SCALE 1. F- THEN
THEN
ELSE FDROP FERR
THEN
ELSE FDROP -1.
THEN
THEN ;

\ Leave r2 = exp(r1), the anti natural logarithm on the floating
\ point stack
: EXP ( F: r1 -- r2 ) ( KMB 85 09 02 )
?NONAN1
IF FDUP MINEXPARG F>
IF FDUP MAXEXPARG F<
IF FSIGN
IF (EXP) 1. F+ 2SCALE
ELSE FDROP 1.
THEN
ELSE FDROP FERR
THEN
ELSE FDROP 0.
THEN
THEN ;


/*
0. FVAR EXP

: EXP ( F: r1 -- r2 )
FDUP MINEXPARG F> ( ANDIF) DUP IF DROP FDUP MAXEXPARG F< THEN
IF FDUP F0=
IF FDROP 1.
ELSE FACTOR-LOG EF@ F/ INTFRAC FSWAP F>D DROP ( N ) ( F: r )
FDUP F0> IF 1. F- 1+ THEN 1. 2SCALE =: EXP ( F: r' )
FNEGATE FACTOR-LOG EF@ F* 2
BEGIN DUP 21 <
WHILE DUP 1- FPSIZE * FACTOR-LOG + EF@ FOVER FOVER F<
IF FDROP 1+
ELSE F- DUP NEGATE EXP 2SCALE FNEGATE +: EXP
THEN
REPEAT DROP FNEGATE 1. F+ EXP F*
THEN
ELSE FDROP FERR
THEN ;
*/

\ Leave r1 raised to the power r2 on the floating point stack.
: ^^ ( F: r1 r2 -- r1^r2 ) ( KMB 85 08 19 )
?NONAN2
IF FOVER FSIGN FDROP ?DUP
IF 0< ( F1 | r1 r2 )
IF ( r1 < 0 ) INTFRAC F0=
IF ( r2 is an integer ) FDUP F>D DROP 1 AND
FSWAP FABS LN F* EXP IF FNEGATE THEN
ELSE ( r2 is not an integer ) FDROP FDROP FERR
THEN
ELSE ( r1 > 0 ) FSWAP LN F* EXP
THEN
ELSE ( r1 = 0 ) FSIGN 0= ( F2 | 0 r2 )
IF ( r2 = 0 ) FDROP FDROP FERR
ELSE ( r2 <> 0 ) FDROP
THEN
THEN
THEN ;

\ Leave r2, the value of 10 raised to the power of r1
\ or the anti base 10 logarithm.
: FALOG ( r1 -- r2 )
10. FSWAP ^^ ;


HEADERLESS
3.14159 .26535 .89793 .23846 0 FP# FCONSTANT |PI

' |PI IS PI
HEADERS

PI 2./ FCONSTANT PI/2
PI 2.* FCONSTANT 2PI
PI 2./ 2./ FCONSTANT PI/4
HEADERLESS

\ Evaluates the sine of r1 to 17+ significant digits for |r1| <= PI/4.
\ Approximation 3043.
: |(SIN) ( F: r1 -- r2 ) ( KMB 85 06 10 )
FDUP FDUP F*
[ .68771 .00349 .0 .0 -11 FP# ] FLITERAL FOVER F*
[ -.17571 .49292 .755 .0 -8 FP# ] FLITERAL F+ FOVER F*
[ .31336 .16216 .61904 .0 -6 FP# ] FLITERAL F+ FOVER F*
[ -.36576 .20415 .84556 .95 -4 FP# ] FLITERAL F+ FOVER F*
[ .24903 .94570 .18873 .6117 -2 FP# ] FLITERAL F+ FOVER F*
[ -.80745 .51218 .82805 .30192 -1 FP# ] FLITERAL F+ F*
[ .78539 .81633 .97448 .30701 0 FP# ] FLITERAL F+ F* ;

' |(SIN) IS (SIN)


\ Evaluates the cosine of r1 to 17+ significant digits for |r1| <= PI/4.
\ Approximation 3824.
: |(COS) ( F: r1 -- r2 ) ( KMB 85 06 10 )
FDUP F*
[ -.38577 .62037 .2 .0 -12 FP# ] FLITERAL FOVER F*
[ .11500 .49702 .4263 .0 -9 FP# ] FLITERAL F+ FOVER F*
[ -.24611 .36382 .63700 .5 -7 FP# ] FLITERAL F+ FOVER F*
[ .35908 .60445 .88581 .953 -5 FP# ] FLITERAL F+ FOVER F*
[ -.32599 .18869 .26687 .55044 -3 FP# ] FLITERAL F+ FOVER F*
[ .15854 .34424 .38154 .10898 -1 FP# ] FLITERAL F+ FOVER F*
[ -.30842 .51375 .34042 .45242 0 FP# ] FLITERAL F+ F*
[ .99999 .99999 .99999 .99996 0 FP# ] FLITERAL F+ ;

' |(COS) IS (COS)

HEADERS
\ Leave r2 = sin(r1) where r1 is any radian angle.
: FSIN ( F: r1 -- r2 ) ( KMB 85 08 19 )
?NONAN1
IF FSIGN 0< FABS PI 2.* F/ INTFRAC FSWAP FDROP
[ .5 2./ 2./ ] FLITERAL F/ INTFRAC FSWAP F>D DROP
CASE
0 OF (SIN) ENDOF
1 OF FNEGATE 1. F+ (COS) ENDOF
2 OF (COS) ENDOF
3 OF FNEGATE 1. F+ (SIN) ENDOF
4 OF (SIN) NOT ENDOF
5 OF FNEGATE 1. F+ (COS) NOT ENDOF
6 OF (COS) NOT ENDOF
7 OF FNEGATE 1. F+ (SIN) NOT ENDOF
DROP
ENDCASE IF FNEGATE THEN
THEN ;

\ Leave r2 = cos(r1) where r1 is any radian angle.
: FCOS PI/2 F+ FSIN ; ( F: r1 -- r2 ) ( KMB 85 06 10 )

HEADERLESS
\ Evaluates the tangent of r1 to 19+ significant digits for |r1| <= PI/4.
\ Approximation 4285.
: |(TAN) ( F: r1 -- r2 ) ( KMB 85 06 11 )
FDUP FDUP F*
FDUP
[ .33866 .38642 .67717 .20961 -4 FP# ] FLITERAL FOVER F*
[ .34225 .54387 .24100 .34353 -1 FP# ] FLITERAL F+ FOVER F*
[ -.15506 .85653 .48326 .63769 2 FP# ] FLITERAL F+ FOVER F*
[ .10559 .70901 .71495 .31936 4 FP# ] FLITERAL F+ F*
[ -.13068 .20264 .75482 .56683 5 FP# ] FLITERAL F+
FSWAP FDUP
[ -.15550 .33164 .03170 .99669 3 FP# ] FLITERAL F+ FOVER F*
[ .47657 .51362 .91648 .36989 4 FP# ] FLITERAL F+ F*
[ -.16638 .95238 .94711 .90019 5 FP# ] FLITERAL F+
F/ F* ;

' |(TAN) IS (TAN)

HEADERS

\ Leave r2 = tan(r1) where r1 is any radian angle.
: FTAN ( F: r1 -- r2 ) ( KMB 85 08 19 )
?NONAN1
IF FSIGN 0< FABS PI F/ INTFRAC FSWAP FDROP
[ .5 2./ ] FLITERAL F/ INTFRAC FSWAP F>D DROP
CASE
0 OF (TAN) ENDOF
1 OF FNEGATE 1. F+ (TAN) 1/X ENDOF
2 OF FSIGN 0=
IF FDROP FERR
ELSE (TAN) 1/X NOT
THEN ENDOF
3 OF FNEGATE 1. F+ (TAN) NOT ENDOF
DROP
ENDCASE IF FNEGATE THEN
THEN ;

HEADERLESS
\ Evaluates the arcsine of r1 to 17+ significant digits for |X1| <= .5.
\ Approximation 4698.
: |(ASIN) ( F: r1 -- r2 ) ( KMB 85 06 11 )
FDUP FDUP F*
FDUP
[ -.36148 .64568 .03475 .23002 2 FP# ] FLITERAL FOVER F*
[ .49908 .74735 .18143 .34756 3 FP# ] FLITERAL F+ FOVER F*
[ -.19037 .55915 .75077 .92670 4 FP# ] FLITERAL F+ FOVER F*
[ .27058 .67326 .43406 .43538 4 FP# ] FLITERAL F+ F*
[ -.12828 .25499 .97869 .27732 4 FP# ] FLITERAL F+
FSWAP FDUP
[ -.75411 .43644 .19617 .07887 2 FP# ] FLITERAL F+ FOVER F*
[ .71974 .04229 .53630 .34267 3 FP# ] FLITERAL F+ FOVER F*
[ -.22941 .55932 .65797 .84211 4 FP# ] FLITERAL F+ FOVER F*
[ .29196 .71576 .43051 .75556 4 FP# ] FLITERAL F+ F*
[ -.12828 .25499 .97869 .27795 4 FP# ] FLITERAL F+
F/ F* ;

' |(ASIN) IS (ASIN)
HEADERS

\ Leave r2 = Arcsin(r1) where r2 is a radian angle.
: FASIN ( F: r1 -- r2 ) ( KMB 85 08 19 )
?NONAN1
IF FSIGN 0< FABS FDUP 1. F>
IF FDROP DROP FERR
ELSE FDUP .5 F>
IF FNEGATE 1. F+ 2./ FSQRT (ASIN) 2.* FNEGATE PI/2 F+
ELSE (ASIN)
THEN IF FNEGATE THEN
THEN
THEN ;

\ Leave r2 = Arccos(r1) where r2 is a radian angle
: FACOS FASIN FNEGATE PI/2 F+ ; ( F: r1 -- r2 ) ( KMB 85 06 11 )

HEADERLESS
\ Evaluates the arctangent of r1 to 17+ significant digits for |r1| <= 1.
\ Approximation 5100.
: |(ATAN) ( F: r1 -- r2 ) ( KMB 85 06 11 )
FDUP FDUP F*
FDUP
[ .97627 .21591 .71763 .30370 -1 FP# ] FLITERAL FOVER F*
[ .11322 .15941 .16764 .65524 2 FP# ] FLITERAL F+ FOVER F*
[ .19257 .92014 .48155 .96135 3 FP# ] FLITERAL F+ FOVER F*
[ .11141 .29072 .84551 .83546 4 FP# ] FLITERAL F+ FOVER F*
[ .27617 .19824 .61388 .34959 4 FP# ] FLITERAL F+ FOVER F*
[ .30310 .74595 .61150 .83044 4 FP# ] FLITERAL F+ F*
[ .12097 .47001 .75809 .07217 4 FP# ] FLITERAL F+
FSWAP FDUP
[ .39917 .88424 .86537 .98150 2 FP# ] FLITERAL F+ FOVER F*
[ .42307 .16464 .80904 .78045 3 FP# ] FLITERAL F+ FOVER F*
[ .18216 .00339 .29184 .64942 4 FP# ] FLITERAL F+ FOVER F*
[ .36645 .44956 .32837 .49894 4 FP# ] FLITERAL F+ FOVER F*
[ .34343 .23596 .19753 .51717 4 FP# ] FLITERAL F+ F*
[ .12097 .47001 .75809 .07287 4 FP# ] FLITERAL F+
F/ F* ;

' |(ATAN) IS (ATAN)

HEADERS
\ Leave r2 = Arctan(r1) where r2 is a radian angle.
: FATAN ( F: r1 -- r2 ) ( KMB 85 06 11 )
?NONAN1
IF FSIGN 0< FABS FDUP 1. F>
IF 1/X (ATAN) FNEGATE PI/2 F+
ELSE (ATAN)
THEN IF FNEGATE THEN
THEN ;

\ Leave angle, in radians, where tan(angle) = y/x
: FATAN2 ( F: r y -- angle ) ( KMB 861113 )
?NONAN2
IF FSIGN ?DUP ( NSY NSY or 0 | r Y )
IF ( Y<>0 ) FSWAP FSIGN ?DUP ( NSY NSX NSX or NSY 0 | y r )
IF ( r<>0 ) F/ FATAN 0< ( NSY F | ANGLE' )
IF ( r<0 ) PI 0< IF ( y<0 ) F- ELSE ( y>0 ) F+ THEN ( | ANGLE )
ELSE DROP ( | ANGLE )
THEN
ELSE ( r=0 ) FDROP FDROP PI/2 0< IF FNEGATE THEN ( | ANGLE )
THEN
ELSE ( y=0 ) FDROP FSIGN FDROP ?DUP ( NSX NSX or 0 )
IF ( r<>0 ) 0< IF PI ELSE 0. THEN ( | ANGLE )
ELSE ( r=0 ) FERR ( | ERR )
THEN
THEN
THEN ;

\ Convert a floating-point number to the nearest signed double integer
\ equivalent, removing the real number from the floating-point stack and
\ leaving the double number result on the forth parameter stack.
' F>D ALIAS FIX ( F: r -- ) ( P: -- dn )


\ Truncate (round towards zero) to a signed double integer removing the
\ real number from the floating point stack and leaving the result on
\ the parameter stack.
: INT ( F: r -- ) ( P: -- dn )
INTFRAC FDROP F>D ;

: F-ROT ( F: r1 r2 r3 -- r3 r1 r2 )
FROT FROT ;

: FNIP ( F: r1 r2 -- r2 )
FSWAP FDROP ;

: FTUCK ( F: r1 r2 -- r2 r1 r2 )
FSWAP FOVER ;

: F>S ( F: r -- )( P: -- n )
F>D DROP ;

: F2DUP ( F: r -- r r )
FOVER FOVER ;

: RADIANS ( -- )
NOOP ;

\ Used by Z* in COMPLEX.SEQ
: F\- ( F: r1 r2 -- r3=r2-r1 )
F- FNEGATE ;


\ Leaves the address of the base of the floating point stack.
: FSP0 ( -- adr )
F0 @ ;

\ Leaves the address of the top of the floating point stack.
: FSP@ ( -- adr )
FP @ ;

\ Check floating point stack for overflow and underflow.
: ?FSTACK ( -- )
FDEPTH 0<
IF FP! TRUE ABORT" Floating point stack underflow."
ELSE FDEPTH 100 >
IF FP! TRUE ABORT" Floating point stack overflow."
THEN
THEN ;



\ ALIASes to match FVG fp standard names and F-PC SFLOAT names.

' FLOATING ALIAS FLOATS
' S>F ALIAS IFLOAT
' S>F ALIAS N>R
' D>F ALIAS D>R
' EF@ ALIAS F@
' EF! ALIAS F!
' FP ALIAS FSP
' FP! ALIAS FCLEAR
' FPSIZE ALIAS F#BYTES
' EXP ALIAS FALN
' EXP ALIAS FEXP
' 2SCALE ALIAS F2**N*
' LN ALIAS FLN
' LOG10 ALIAS FLOG
' 1/X ALIAS 1/F
' ^^ ALIAS F**
' 2.* ALIAS F2*
' 2./ ALIAS F2/
' 0. ALIAS F0.0
' 1. ALIAS F1.0
' 10. ALIAS F10.0


/*
VARIABLE CTL-WRD

HEX

\ Leave a true flag if the 8087 is present.
CODE ?8087 ( -- flag )
BX PUSH \ Make room for flag
BX, BX XOR \ Set flag to false.
FNINIT
AX, AX XOR \ Set AX to zero
WORD CS: CTL-WRD , AX MOV \ Clear CTL-WRD
CS: CTL-WRD FNSTCW \ Store control word.
WORD CS: AX, CTL-WRD MOV
AH, # 03 CMP \ Will be 03 if 8087 is present.
1$ JNE \ Jump if no 8087 present.
BX, # FFFF MOV \ Set flag true, 8087 found.
1$: NEXT
END-CODE

: OVL8087 ( -- ) ( DJL 870722 )
?8087 IF " VPP87.OVR" CS ' FPSIZE "PATH-LOAD THEN ;
*/
DECIMAL