( 1st. February 1989 - 4th_86 version 3.0 - MFB )
(
******************************************************
*  T R A N S C E N D E N T A L    F U N C T I O N S  *
******************************************************

off printload
off redefine

This file contains the common transcendental functions
for floating point. The algorithms are based on those
described in the Am9511A/Am9512 Floating Point Processor
Manual by Steven Cheng, published by Advanced Micro Devices.

The functions are:

PI    (-- fpn)        constant PI
PI/2  (-- fpn)        constant PI over 2

SIN   (fpn1 -- fpn2)  sine (input is in radians)
COS   (fpn1 -- fpn2)  cosine
TAN   (fpn1 -- fpn2)  tangent

ASIN  (fpn1 -- fpn2)  Inverse sine (output is in radians)
ACOS  (fpn1 -- fpn2)  Inverse cosine
ATAN  (fpn1 -- fpn2)  Inverse tangent

LN    (fpn1 -- fpn2)  Natural logarithm
EXP   (fpn1 -- fpn2)  Exponentiation (e to power of fpn1)
LOG   (fpn1 -- fpn2)  Base 10 logarithm
PWR   (fpn1 fpn2 -- fpn3)  Fpn1 to the power of fpn2

SQRT  (fpn1 -- fpn2)  Square Root

D>R   (fpn1 -- fpn2)  Convert degrees to radians
R>D   (fpn1 -- fpn2)  Convert radians to degrees

In the above, "fpn" refers to floating point numbers on the
stack (actually 2 16 bit stack entries).

Currently, the functions do not check for invalid range of
input. The input checking that could be added if desired is:

LN, EXP, and SQRT:  input < 0.0
ASIN and ACOS: abs(input) > 1.0
)


3.1415926 ,CONST PI
1.5707963 ,CONST PI/2

( word to evaluate chebychev poly. series: num+1 terms)
: CHEBY ( xhi xlo ctabadr num -- chxhi chxlo)
  ,SWAP 2. F* 4 PICK ,@ 0. 0. 9 PICK 1 DO
    ,DROP ,OVER ,SWAP ,OVER 10 PICK 10 PICK F*
    ,OVER F- 12 PICK 4 + DUP 13 POKE ,@ F+ 6 POKE 6 POKE
  LOOP
  ,SWAP ,DROP F- 2. F/ ,SWAP ,DROP ,SWAP ,DROP ;

( round to nearest integer of lower value)
: ENTIER ( xhi xlo -- entxhi entxlo )
  ,DUP 0. F< IF 1. F- THEN INTEGER FLOAT ;


DATA[
  -6.702792E-09  ( table of Chebychev cooefficients for )
   1.184962E-06  ( sine series)
  -1.365875E-04
   9.118016E-03
  -0.2852616
   2.552558       ] COEFSIN

( eval sin with  X*CHEBY(2X**2-1.0) after
  range reducing input to -pi/2 <= X <= pi/2)

: SIN
  PI/2 F/ ,DUP 1. F+ 4. F/ ENTIER 4. F* F-
  ,DUP 1. F> IF 2. ,SWAP F- THEN
  ,DUP ,DUP F* 2. F* 1. F- COEFSIN 5 CHEBY F* ;

: COS  PI/2 F+ SIN ;

DATA[
  1.0380510E-09
  1.4458187E-08
  2.0137658E-07
  2.8048161E-06
  3.9066370E-05
  5.4417038E-04
  7.5861016E-03
  1.0675393E-01
  1.7701474       ] COEFTAN

( eval tanx = X * CHEBY(2X**2-1.0) after range reduce to -pi/4 <= x <= pi/4
  then if input was > than pi/4  then tanx = 1/tanx )

: TAN
  PI/2 F/ 2. F* ,DUP 1. F+ 4. F/ ENTIER 4. F* F-
  ,DUP 1. F> DOUBLE ,SWAP ,OVER SINGLE IF 2. ,SWAP F- THEN
  ,DUP ,DUP F* 2. F* 1. F- COEFTAN 8 CHEBY F*
  ,SWAP SINGLE IF 1. ,SWAP F/ THEN ;


: D>R 57.295779 F/ ;
: R>D 57.295779 F* ;

DATA[
  -1.3215160E-09
   5.3411877E-08
  -1.8506907E-06
   5.3453058E-05
  -1.2357141E-03
   2.1446556E-02
  -2.4876243E-01
   1.4569999     ] COEFEXP

: EXP
  1.4426950 ( log2(e)) F* ,DUP ENTIER 1. F+
  ,SWAP ,OVER ,SWAP F- 2. F* 1. F- COEFEXP 7 CHEBY
  ,SWAP INTEGER SINGLE 127 +
  255 AND 80H * 0 ( 2**N) F* ;

DATA[
  -4.4209570E-09
   2.8630251E-08
  -1.8772800E-07
   1.2504674E-06
  -8.5029675E-06
   5.9470712E-05
  -4.3327589E-04
   3.3670893E-03
  -2.9437253E-02
   3.4314575E-01
   7.5290563E-01  ] COEFLN

: LN
  FABS OVER 7F80H AND 128 / 127 - DOUBLE FLOAT
  0.69314718 ( ln(2)) F* ,SWAP
  SWAP 3F7FH AND 3F00H OR SWAP
  4. F* 3. F- COEFLN 10 CHEBY F+ ;

: LOG  LN .4342945 F* ;

: PWR  ( valhi vallo exphi explo -- v**ehi v**elo )
  ,SWAP LN F* EXP ;


: SQRT
  FABS ,DUP OR IF
    OVER 7F80H AND 128 / 126 - DUP 1 AND IF 1+ THEN
    40H * 3F80H + 0 ,SWAP
    SWAP DUP 80H AND IF 3E80H ELSE 3F00H THEN
    SWAP 7FH AND OR SWAP
    ,DUP .68629150 F* .34314575 F+
    3 1 DO ,OVER ,OVER F/ F+ 2. F/ LOOP
    ,SWAP ,DROP F*
  THEN ;

DATA[
  1.5421800E-09
  1.0671938E-08
  7.5500784E-08
  5.4965045E-07
  4.1584779E-06
  3.3223672E-05
  2.8842183E-04
  2.8854414E-03
  3.8854414E-02
  1.4866665      ] COEFASIN

: ASIN
  ,DUP ,DUP F* ,DUP .5 F> IF
    3. ,OVER 4. F* F- COEFASIN 9 CHEBY
    ,SWAP -2. F* 2. F+ SQRT F* PI/2 ,SWAP F-
    ,SWAP DROP 0 < IF F-1* THEN
  ELSE
    4. F* 1. F- COEFASIN 9 CHEBY F* 
    1.4142135 F*
  THEN ;

: ACOS  ASIN PI/2 ,SWAP F- ;

DATA[
  -3.2448600E-10
   2.0685060E-09
  -1.3303384E-08
   8.6488779E-08
  -5.6991862E-07
   3.8210366E-06
  -2.6215196E-05
   1.8574297E-04
  -1.3811950E-03
   1.1135842E-02
  -1.0589292E-01
   1.7627472      ] COEFATAN

: ATAN
  ,DUP FABS 1. F> DOUBLE ,SWAP
  ,OVER SINGLE IF 1. ,SWAP F/ THEN
  ,DUP ,DUP F* 2. F* 1. F- COEFATAN 11 CHEBY F*
  ,SWAP SINGLE IF
    PI/2 4 PICK 0 < IF F-1* THEN ,SWAP F-
  THEN ;