home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / forth / compiler / fpc / source / tfloat.seq < prev    next >
Text File  |  1990-04-27  |  50KB  |  1,655 lines
  1. \\ High Level FLoating Point .  IEEE Format.     11:30 07Nov88RS)
  2.  
  3.           Copyright 1988 by Robert L. Smith
  4.       2300 St. Francis Drive, Palo Alto, CA 94303
  5.                  (415) 856-9321
  6.  
  7.  These routines may be freely used, provided only that the
  8.  copyright notice be displayed and preserved.
  9.  
  10. \                                               09:02 23Oct88RS)
  11.  
  12. Please: If you wish to distribute any of your changes to this
  13.  package, please give your name, address and telephone number
  14.  so that any other users can contact you if they find problems.
  15.  
  16.       ( Approximate times on 4.77 MHz 8088 are: )
  17.       ( F+    10 milliseconds )
  18.       ( F/    20 milliseconds )
  19.       ( FLN   30 milliseconds )
  20.  
  21. \ Load Block                                    14:33 08Nov88RS)
  22.  
  23. : COPYRIGHT
  24.      CR ." Floating Forth  Version 1.2   Jan. 25, 1990 " CR
  25.      CR ." written by Robert L. Smith "
  26.      CR ." 2300 St. Francis Drive, Palo Alto, CA 94303 " CR ;
  27.  
  28. {
  29.  
  30. DECIMAL
  31. >FORTH
  32.  
  33. ONLY FORTH ALSO DEFINITIONS
  34.  
  35. FLOAD FLOAT4TH.SEQ
  36.  
  37. ONLY FORTH ALSO COMPILER ALSO TARGET ALSO DEFINITIONS ASSEMBLER ALSO
  38.  
  39. FORTH DECIMAL TARGET >LIBRARY
  40.  
  41. ' FCON ALIAS FCONSTANT          \ make compilers FCONSTANT be the version
  42.                                 \ of FCONSTANT used in the compiler
  43. FORTH
  44. ' FPOP IS FLOAT_POP             \ link to compilers FCONSTANT
  45. TARGET
  46.  
  47. CODE DSHFT8     ( lo hi -- lowest middle high )
  48.                 LODSW           \ ax=low, bx=hi
  49.                 SUB CL, CL
  50.                 MOV CH, AL      MOV AL, AH
  51.                 MOV AH, BL      MOV BL, BH
  52.                 SUB BH, BH
  53.                 DEC SI          DEC SI          MOV 0 [SI], CX
  54.                 DEC SI          DEC SI          MOV 0 [SI], AX
  55.                 RET             END-CODE
  56.  
  57. }
  58. \ 8-bit shift functions
  59.  
  60. DSHFT8   Perform an 8 bit shift on a double precision number,
  61.    returning a triple word result (no loss of bits).  The
  62.    results can be considered as either right or left shifted.
  63. {
  64.  
  65. IMACRO D1+      ( d1 -- d1+ )
  66.                 ADD 0 [SI], # 1 WORD
  67.                 ADC BX, # 0     END-IMACRO
  68.  
  69. \  MA  MB  MC  XMD  Y*  FRACT*                 07:07 15Oct88RS)
  70.  
  71. 8 ARRAY MA
  72.  
  73. : Y*   ( d1 d2 -- t )  ( High order 8 bits of d1 and d2 = 0 )
  74.         MA 2!  MA 4 + 2!  MA 2+ @ MA 6 + @ UM* 0
  75.         MA @ MA 6 + @ UM* D+  MA 2+ @ MA 4 + @ UM* D+
  76.         MA @ MA 4 + @ * + ;
  77.  
  78. : FRACT*   ( n1 n2 -- n3 )
  79.         UM* NIP ;
  80.  
  81. }
  82. \  Partial and fractional multiplication.       20:51 23Oct88RS)
  83.  
  84. MA  MB  MC  XMD  Variables used for multiplication.
  85.  
  86. Y*   Take the product of two 24-bit numbers, returning a triple
  87.    precision result.  The 24-bit multiplicands are double
  88.    numbers with the most significant 8 bits cleared to 0.
  89.  
  90. FRACT*  Take the product of two single numbers, treated as
  91.    fractions with the binary point at the left.
  92.  
  93. {
  94. \  F#BYTES  FPSTAT  FPSTAT  FPS0  FDEPTH        16:17 15Nov88RS)
  95.  
  96. 4 CONSTANT F#BYTES           ( Size of floating point number )
  97.  
  98. TABLE FPSTAT    0 , 0 ,          ( Holds error information )
  99.                 END-TABLE
  100.  
  101. 20 F#BYTES * CONSTANT FPSSIZE           \ Size of floating pt. stack
  102.  
  103. FPSSIZE 3 F#BYTES * + ARRAY FPSTACK     \ Floating Point Stack
  104.                                         \ with room for underflow
  105.  
  106. : FSP0          ( -- a1 )               \ Point to base of stack.
  107.                 FPSTACK FPSSIZE + ;
  108.  
  109. VARIABLE FSP                            \ Points to top of stack
  110.                                         \ initialized by FCLEAR
  111.  
  112. : FDEPTH   ( -- n )
  113.         FSP0 FSP @ - F#BYTES / ;
  114.  
  115. }
  116. \ Definition of the Floating Point stack.       21:36 23Oct88RS)
  117.  
  118. F#BYTES   Size of a floating point number in bytes.
  119. FPSTAT    A variable to hold the Floating Point Status.
  120.  
  121. FPSSIZE   The size of the floating point stack, in bytes.
  122.  
  123. FPSTACK   The floating point stack.
  124.  
  125. FSP0      The base of the floating point stack.
  126. FSP       Contains the pointer to the F.P. Stack top.
  127.  
  128. FDEPTH    The depth of the floating point stack, in FPSSIZE
  129.           units.
  130.  
  131. {
  132. \  XBIAS  FCLEAR  FDROP  FPUSH                  14:40 08Nov88RS)
  133.  
  134. $3F80 CONSTANT XBIAS             ( Exponent bias )
  135.  
  136. : FPUSH   ( F: -- r ; d -- )
  137.         -4 FSP +!  FSP @ 2! ;           EXECUTES> FPUSH
  138.  
  139. ' FPUSH ALIAS FPUSHER
  140.  
  141. FORTH   >FORTH
  142.  
  143. : %FPUSH        ( | string" -- )
  144.                 [FORTH]
  145.                 ['] FPUSHER COMP_CALL ; IMMEDIATE
  146.  
  147. ' %FPUSH FORTH IS COMP_FPUSH    \ link into defered word
  148.  
  149. : [F#]          ( | floating_number -- )
  150.                 [FORTH]
  151.                 F# FPOP SWAP COMP_SINGLE COMP_SINGLE
  152.                 COMP_FPUSH ; IMMEDIATE
  153.  
  154. FORTH DECIMAL TARGET >LIBRARY       \ Library continues
  155.  
  156. : FCLEAR   ( -- )
  157.         FSP0 FSP ! ;
  158.  
  159. : FDROP   ( F: r -- )
  160.         F#BYTES FSP +! ;
  161.  
  162. : F2DROP  ( -- )
  163.         FDROP FDROP ;
  164.  
  165. }
  166. \  XBIAS  FCLEAR  FDROP                         14:17 08Nov88RS)
  167.  
  168. XBIAS   The exponent bias.
  169.  
  170. FCLEAR  Empties the floating point stack.
  171.  
  172. FDROP   Drop one floating point element from the F.P. stack.
  173.  
  174. FNSWAP  Exchange the n-th item on the F.P. stack with the
  175.         zeroth item.
  176.  
  177. {
  178. \ FPICK  FDUP  FOVER  FSWAP                     14:32 08Nov88RS)
  179.  
  180. : FPICK   ( F: rn rn-1 ... r0 -- fn rn-1 ... r0 rn ; n -- )
  181.         F#BYTES * FSP @ + 2@ FPUSH ;
  182.  
  183. : FDUP   ( F: r -- r r )
  184.         FSP @ 2@  F#BYTES NEGATE FSP +!
  185.         FSP @ 2! ;
  186.  
  187. : FOVER   ( F: r1 r2 -- r1 r2 r1 )
  188.         F#BYTES NEGATE FSP +!
  189.         FSP @ 8 + 2@ FSP @ 2! ;
  190.  
  191. : F2DUP   ( ? -- ? )
  192.         FOVER FOVER ;
  193.  
  194. : FSWAP   ( F: r1 r2 -- r2 r1 )
  195.         FSP @ 2@  FSP @ 4 + 2@
  196.         FSP @ 2!  FSP @ 4 + 2! ;
  197. }
  198. \ FDUP  FOVER  FSWAP                            14:35 08Nov88RS)
  199.  
  200. FPICK   Push a copy of the nth element of the F.P. stack
  201.         onto the F.P. stack.
  202.  
  203. FDUP    Duplicate the top element on the F.P. stack.
  204.  
  205. FOVER   Push a copy of the second element on the F.P. stack
  206.         onto the F.P. stack.
  207.  
  208. FSWAP   Interchange the top two elements on the F.P. stack.
  209.  
  210. {
  211. \  FNEGATE  FNIP  FROT  F-ROT                   21:40 23Oct88RS)
  212.  
  213. : FNEGATE   ( F: r1 -- r2 )
  214.         FSP @ DUP @ $8000 XOR SWAP ! ;
  215.  
  216. : FNIP   ( F: r1 r2 -- r2 )
  217.         FSP @ 2@  4 FSP +!  FSP @ 2! ;
  218.  
  219. : FROT   ( F: r1 r2 r3 -- r2 r3 r1 )
  220.         FSP @ 2@  FSP @ 4 + 2@  FSP @ 8 + 2@
  221.         FSP @ 2!  FSP @ 8 + 2!  FSP @ 4 + 2! ;
  222.  
  223. : F-ROT   ( F: r1 r2 r3 -- r3 r1 r2 )
  224.         FSP @ 2@  FSP @ 4 + 2@  FSP @ 8 + 2@
  225.         FSP @ 4 + 2!  FSP @ 2!  FSP @ 8 + 2! ;
  226.  
  227. }
  228. \  FNEGATE  FNIP  FROT  F-ROT                   21:44 23Oct88RS)
  229.  
  230. FNEGATE   Reverse the sign of the element at the top of the
  231.           F.P. stack.
  232.  
  233. FNIP      Drop the second item from the F.P. stack.
  234.  
  235. FROT      Rotate the top three items on the F.P. stack,
  236.           bring the third element to the top.
  237.  
  238. F-ROT     Rotate the top three items on the F.P. stack,
  239.           bringing the second item to the top, and moving the
  240.           former top item to the third position.
  241.  
  242. {
  243. \ FPOP  FPCOPY  F0=  FNSWAP              14:39 014:45 08Nov88RS)
  244.  
  245. : FPOP   ( F: r -- ; -- d )
  246.         FSP @ 2@ FDROP ;                EXECUTES> FPOP
  247.  
  248. : FPCOPY   ( F: r -- r ; -- d )
  249.         FSP @ 2@ ;
  250.  
  251. : F0=   ( F: r -- ; -- flag )
  252.         FPOP $7FFF AND OR 0= ;
  253.  
  254. : FPOP0=   ( F: r -- ; -- d flag )
  255.         FPOP 2DUP $7FFF AND OR 0= ;
  256.  
  257. }
  258. \ FPOP  FPUSH  FPCOPY  F0=                      23:02 23Oct88RS)
  259.  
  260. FPOP    Pop the top number from the F.P. stack, and push it
  261.         onto the parameter stack as a double number.
  262.  
  263. FPUSH   Pop the top double number and
  264.  
  265. FPCOPY  Get a copy of the top number on th F.P. stack and push
  266.         it on the parameter stack.
  267.  
  268. F0=     Pop the top member of the F.P. stack and test its value.
  269.         If the value is zero, push true, else push false onto
  270.         the parameter stack.
  271.  
  272. {
  273. \  F0<  F=                                      14:45 08Nov88RS)
  274.  
  275. : F0<   ( F: r -- ; -- flag )
  276.         FPOP DUP 0<
  277.         IF      $7FFF AND OR 0= 0=
  278.         ELSE    2DROP 0   THEN ;
  279.  
  280. : F=   ( F: r1 r2 -- ; -- flag )
  281.         FPOP0=
  282.         IF      2DROP FPOP0= NIP NIP
  283.         ELSE    FPOP D=   THEN ;
  284.  
  285. : FNSWAP   ( F: rn rn-1 ... r0 -- r0 rn-1 ... rn ; n -- )
  286.         F#BYTES * FSP @ + DUP >R 2@  FSP @ 2@
  287.         R> 2!  FSP @ 2! ;
  288.  
  289. }
  290. \  F0<  FPOP0=  F=                              23:13 23Oct88RS)
  291.  
  292. F0<     Pop and test the number at the top of the F.P. stack.
  293.         If the sign is negative and the value is non-zero, push
  294.         a true flag.  Otherwise push a false flag.
  295.  
  296. FPOP0=  Pop the top member of the F.P. stack, test, and push it
  297.         onto the parameter stack.  If the number has a value of
  298.         zero, also push a true flag; otherwise push a false
  299.         flag.
  300. F=      Pop the top two elements from the F.P. stack.  If the
  301.         two numbers are equal, push a true flag.  Otherwise,
  302.         push a false flag.
  303.  
  304. {
  305. \  F<                                           16:07 14Oct88RLS
  306.  
  307. : F<    ( F: r1 r2 -- ; -- flag )
  308.         FPOP FPOP DUP 0<
  309.         IF      DU<
  310.         ELSE    2SWAP D<
  311.         THEN ;
  312.  
  313. }
  314. \  F<                                           23:16 23Oct88RS)
  315.  
  316. F<      Pop the top two F.P. numbers from the F.P. stack and
  317.         compare them.  If the second number is arithmetically
  318.         less than the first number, push a true flag.
  319.         Otherwise push a false flag.
  320.  
  321. {
  322. \ F>                                            06:44 18Jul89RS)
  323.  
  324. : F>    ( F: r1 r2 -- ; -- flag )
  325.         FPOP FPOP DUP 0<
  326.         IF      2SWAP DU<
  327.         ELSE    D<
  328.         THEN ;
  329.  
  330. }
  331. \ F>                                            23:18 23Oct88RS)
  332.  
  333. F>      Pop the top two numbers from the F.P. stack and
  334.         compare them.  If the second number is arithmetically
  335.         greater than the top, push a true flag onto the
  336.         parameter stack.  Otherwise push a false flag.
  337.  
  338. {
  339. \  FABS  FMIN  FMAX                             06:31 18Jul89RS)
  340.  
  341.  
  342. : FABS  ( F: r1 -- r2 )
  343.         FSP @ DUP @ $7FFF AND SWAP ! ;
  344.  
  345. : FMIN   ( F: r1 r2 -- r3 )
  346.         FPOP FPOP 2OVER 2OVER DUP 0<
  347.         IF   DU<
  348.         ELSE 2SWAP D<
  349.         THEN
  350.         IF   2SWAP  THEN
  351.         2DROP FPUSH ;
  352.  
  353. : FMAX   ( F: r1 r2 -- r3 )
  354.         FPOP FPOP 2OVER 2OVER DUP 0<
  355.         IF      2SWAP DU<
  356.         ELSE    D<
  357.         THEN
  358.         IF  2SWAP  THEN
  359.         2DROP FPUSH ;
  360.  
  361. }
  362. \  FABS  FMIN  FMAX                             23:26 23Oct88RS)
  363.  
  364. FABS    Set the sign of the top of the F.P. stack to 0.
  365.  
  366. FMIN    Pop the top two members from the F.P. stack.  Push the
  367.         arithmetically smaller back onto the F.P. stack.
  368.  
  369. FMAX    Pop the top two members from the F.P. stack.  Push the
  370.         arithmetically larger back onto the F.P. stack.
  371.  
  372. {
  373. \ F@  F!  FCONSTANT  FVARIABLE                  08:03 20Oct88RS)
  374.  
  375. : F@   ( F: -- r ; addr -- )
  376.         2@ FPUSH ;
  377.  
  378. : F!   ( F: r -- ; addr -- )
  379.         >R FPOP R> 2! ;
  380.  
  381. }
  382. \ F@  F!  FCONSTANT  FVARIABLE                  00:03 23Oct88RS)
  383.  
  384. F@      Fetch the F.P. variable at the address specified by the
  385.         top of the parameter stack.  Push the contents of the
  386.         variable onto the F.P. stack.  Pop and discard the
  387.         address at the top of the parameter stack.
  388. F!      Store the number at the top of the F.P. stack into
  389.         memory at the address at the top of the parameter stack.
  390.         Pop the number from the F.P. stack, and pop the address
  391.         from the parameter stack.
  392. FCONSTANT  Create a F.P. constant with a value equal to the
  393.         number poped off the F.P. stack.
  394. FVARIABLE  Create a F.P. variable.
  395.  
  396. {
  397. \ Various FCONSTANTs: F1.0 PI F0.0 FLOG10E      07:20 15Oct88RS)
  398.  
  399. F# 2.0E0 FCONSTANT F2.0
  400.  
  401. $0000 $3F80 FPUSH FCONSTANT F1.0
  402. $0FDB $4049 FPUSH FCONSTANT PI
  403. $0000 $0000 FPUSH FCONSTANT F0.0
  404. $5BD9 $3EDE FPUSH FCONSTANT FLOG10E
  405. $5D8E $4013 FPUSH FCONSTANT FLN10.0
  406. $0000 $4120 FPUSH FCONSTANT F10.0
  407. $0000 $3F00 FPUSH FCONSTANT F0.5
  408.  
  409. }
  410. \ Various FCONSTANTs: F1.0 PI F0.0 FLOG10E      00:07 23Oct88RS)
  411.  
  412. F1.0      Floating point 1.
  413. PI        Floating point pi
  414. F0.0      Floating point 0.
  415. FLOG10E   Floating point log base 10 of e
  416. FLN10.0   Floating point natural log of 10
  417. F10.0     Floating point 10.
  418. F0.5      Floating point 0.5
  419.  
  420. {
  421. \ T2/  T2*  T>SHIFT                             07:21 15Oct88RS)
  422.  
  423. FALSE #IF       \ FALSE = load the CODE equivelant words
  424.                 \ TRUE  = load the hi-level versions
  425.  
  426. : T2/   ( t1 -- t2 )
  427.         >R D2/ R@ 1 AND
  428.         IF    $8000 OR   ELSE   $7FFF AND   THEN    R> 2/ ;
  429.  
  430. : T2*   ( t1 -- t2 )
  431.         >R DUP 0<
  432.         IF   D2* R> 2* 1 OR   ELSE   D2* R> 2*  THEN ;
  433.  
  434. : T>SHIFT   ( t1 n -- t2 )   ( n is multiple of 80 hex )
  435.         ?DUP
  436.         IF   0  DO  T2/  $80 +LOOP  THEN ;
  437.  
  438. #ELSE
  439.  
  440. MACRO  T2/
  441.         SAR   BX, # 1
  442.         RCR   0 [SI], # 1 WORD
  443.         RCR   2 [SI], # 1 WORD
  444.         END-MACRO
  445.  
  446. MACRO  T2*
  447.         SHL   2 [SI], # 1 WORD
  448.         RCL   0 [SI], # 1 WORD
  449.         RCL   BX, # 1
  450.         END-MACRO
  451.  
  452. CODE T>SHIFT
  453.         MOV CX, BX
  454.         JCXZ 0 $
  455.         MOV AX, CX
  456.         SUB DX, DX
  457.         MOV BX, # $80
  458.         IDIV BX                 \ divide divisor by $80
  459.         MOV CX, AX
  460.         LODSW
  461.         MOV BX, AX
  462.    1 $: SAR BX, # 1
  463.         RCR 0 [SI], # 1 WORD
  464.         RCR 2 [SI], # 1 WORD
  465.         LOOP 1 $
  466.         RET
  467.    0 $: LODSW
  468.         MOV BX, AX
  469.         RET             END-CODE
  470.  
  471. #ENDIF
  472.  
  473. }
  474. \ T2/  T2*  T>SHIFT                             00:10 23Oct88RS)
  475.  
  476. T2/     Shift right a triple precision number.
  477.  
  478. T2*     Shift left a triple precision number.
  479.  
  480. T>SHIFT  Shift right a triple precision number by a count equal
  481.          to the argument divided by 128 (80 hex).
  482.  
  483. {
  484. \ D>SHIFT  D<SHIFT  D-CY                        07:21 15Oct88RS)
  485.  
  486. \ : D>SHIFT   ( d1 n -- d2 )   ( n a multiple of 80 hex )
  487. \         ?DUP  IF  0  DO D2/  $80 +LOOP  THEN ;
  488. \
  489. \ : D<SHIFT   ( d1 n -- d2 )  ( n is multiple of 80 hex )
  490. \         ?DUP
  491. \         IF      0 DO  D2*  $80 +LOOP
  492. \         THEN ;
  493.  
  494. CODE D>SHIFT
  495.         MOV   CX, BX
  496.         LODSW
  497.         MOV   BX, AX
  498.         JCXZ  0 $
  499.    1 $: SAR   BX, # 1
  500.         RCR   0 [SI], # 1 WORD
  501.         LOOP 1 $
  502.    0 $: RET             END-CODE
  503.  
  504. CODE D<SHIFT
  505.         MOV   CX, BX
  506.         LODSW
  507.         MOV   BX, AX
  508.         JCXZ  3 $
  509.    4 $: SHL   0 [SI], # 1 WORD
  510.         RCL   BX, # 1
  511.         LOOP 4 $
  512.    3 $: RET             END-CODE
  513.  
  514. : D-CY   ( d1 d2 -- d3 n )   ( n is borrow or carry )
  515.         2OVER 2OVER DU<
  516.         IF      D- -1
  517.         ELSE    D- 0
  518.         THEN ;
  519.  
  520. }
  521. \ D>SHIFT  D<SHIFT  D-CY                        00:14 23Oct88RS)
  522.  
  523. D>SHIFT  Shift a double precision number right by a count equal
  524.          to the top parameter divided by 128 (80 hex).
  525.  
  526. D<SHIFT  Shift a double precision number left by a count equal
  527.          to the top parameter divided by 128 (80 hex).
  528.  
  529. D-CY     Subtract the top double number from the second.
  530.          Return the difference and a borrow flag.
  531.  
  532. {
  533. \  D+CY  ZSIGN  ZEXP  FIXGRS                    07:21 15Oct88RS)
  534.  
  535. : D+CY   ( d1 d2 -- d3 n )
  536.         DNEGATE D-CY ;
  537.  
  538. TABLE ZSIGN     0 ,     END-TABLE
  539. TABLE ZEXP      0 ,     END-TABLE
  540.  
  541. : FIXGRS   ( n1 -- n2 )
  542.         DUP $3FFF AND
  543.         IF      $C000 AND $0F OR
  544.         THEN ;
  545.  
  546. }
  547. \  D+CY  ZSIGN  ZEXP  FIXGRS                    00:18 23Oct88RS)
  548.  
  549. D+CY    Add the top two double numbers and return a carry flag.
  550.  
  551. ZSIGN   A variable to carry the sign of a result.
  552. ZEXP    A vaiable to carry the resultant exponent.
  553.  
  554. FIXGRS  Set "sticky" bit flags.
  555.  
  556. {
  557. \  UNNORMALIZE                                  07:21 15Oct88RS)
  558.  
  559. : UNNORMALIZE   ( d1 n -- d2 grs )  ( n is multiple of 80X )
  560.         NO_INLINE
  561.         >R R@ $480 <
  562.         IF      0 -ROT R> T>SHIFT ROT FIXGRS EXIT   THEN
  563.         R@ $880 <
  564.         IF      $800 R> - D<SHIFT SWAP 0 SWAP FIXGRS EXIT  THEN
  565.         R@ $D80 <
  566.         IF      SWAP DUP $3FF AND
  567.                 IF      $4000 OR   THEN  SWAP
  568.                 R> $800 - D>SHIFT 0 ROT FIXGRS EXIT
  569.         THEN
  570.         R> DROP OR
  571.         IF      0 0 $0F  ELSE   0 0 0  THEN ;
  572.  
  573. }
  574. \  UNNORMALIZE                                  00:21 23Oct88RS)
  575.  
  576. UNNORMALIZE   Unnormalize the double number by a count specified
  577.               at the top of the stack.  The number at the top of
  578.               the stack is a signed count multiplied by 128
  579.               (128 hex).
  580.  
  581. {
  582. \  >NORMALIZE                                   07:22 15Oct88RS)
  583.  
  584. : >NORMALIZE   ( d1 -- grs d2 n )
  585.         0 -ROT 0 $480 0
  586.         DO      DROP DUP $100 U<
  587.                 IF      I LEAVE
  588.                 ELSE    T2/ I
  589.                 THEN
  590.        $80 +LOOP ;
  591.  
  592. }
  593. \  >NORMALIZE                                   00:26 23Oct88RS)
  594.  
  595. >NORMALIZE   Normalize the double number at the top of the
  596.              parameter stack.  Return GRS (Guard, Round, and
  597.              Sticky) bits and a normalizing count multiplied
  598.              by 128.
  599.  
  600. {
  601. \ MROUND  EVROUND  DENORMALIZE1                 07:57 31Oct88RS)
  602.  
  603. : MROUND   ( d1 evenflg -- d2 )
  604.         IF      D1+ SWAP $FFFE AND SWAP
  605.         ELSE    D1+
  606.         THEN ;
  607.  
  608. : EVROUND   ( d1 grs -- d2 )
  609.         DUP 0<
  610.         IF      $8000 =   MROUND
  611.         ELSE    DROP
  612.         THEN ;
  613.  
  614. : DENORMALIZE1   ( d1 n -- d2 )
  615.                 NEGATE $80 + UNNORMALIZE EVROUND ;
  616.  
  617. }
  618. \ MROUND  EVROUND  DENORMALIZE1                 07:57 31Oct88RS)
  619.  
  620. MROUND  Round the double number on the stack.  If the flag at
  621.         the top of the stack is true, then round to even.
  622.  
  623. EVROUND  If the guard bit is set, round the double number on the
  624.          stack.
  625.  
  626. DENORMALIZE1  Denormalize the double number according to the
  627.         shifted count at the top of the stack.  Round toward
  628.         even.
  629. {
  630. \ DENORMALIZE2  1NORMALIZE                      07:56 31Oct88RS)
  631.  
  632. : DENORMALIZE2   ( grs d1 n -- d2 )
  633.         NEGATE $80 + UNNORMALIZE >R ROT R> OR EVROUND ;
  634.  
  635. : 1NORMALIZE    ( d1 -- d2 n )
  636.         $8000 $C00 0
  637.         DO      DROP DUP $7F >
  638.                 IF      $7F AND I NEGATE LEAVE
  639.                 ELSE    D2* $8000
  640.                 THEN
  641.        $80 +LOOP ;
  642.  
  643. }
  644.                                                 08:04 31Oct88RS)
  645.  
  646. DENORMALIZE2  Denormalize the combination of d1 and GRS
  647.               according to the shifted count n .
  648.  
  649. 1NORMALIZE    Shift the double number left until it is in a
  650.               normalized form.  A shifted form of the count is
  651.               left on the stack.  The shifted form is used to
  652.               speed up the conversion process.
  653.  
  654. {
  655. \ NORMALIZE  4NORMALIZE                         08:16 31Oct88RS)
  656.  
  657. : NORMALIZE   ( d1 -- d2 n )
  658.         NO_INLINE
  659.         2DUP D0=
  660.         IF      0 EXIT  THEN
  661.         1NORMALIZE ;
  662.  
  663. : 4NORMALIZE   ( d1 -- d2 )
  664.         NO_INLINE
  665.         2DUP OR 0=  IF   EXIT   THEN
  666.         $1000 0
  667.         DO      D2* DUP $7F >
  668.                 IF      $7F AND $3F00 I - OR LEAVE
  669.                 THEN
  670.        $80 +LOOP ;
  671.  
  672. }
  673.                                                 08:26 31Oct88RS)
  674.  
  675. NORMALIZE  Perform a normalization process on the double number.
  676.            Push the shifted count on the stack.
  677.  
  678. 4NORMALIZE  Normalize a double number by shifting left.  This is
  679.             used only by 2NORMALIZE (and ultimately by FLN).
  680.             The result is really a floating point number on the
  681.             parameter stack.
  682.  
  683. {
  684. \ 3NORMALIZE                                    08:18 31Oct88RS)
  685.  
  686. : 3NORMALIZE   ( d1 -- d2 )
  687.         $4480 $4000
  688.         DO      DUP $200 <
  689.                 IF      D1+ D2/ DUP $0FF >
  690.                         IF      D2/
  691.                         ELSE    $7F AND
  692.                         THEN
  693.                         I + LEAVE
  694.                 ELSE    D2/
  695.                 THEN
  696.        $80 +LOOP ;
  697.  
  698. }
  699.                                                 08:25 31Oct88RS)
  700.  
  701. 3NORMALIZE  Normalize a double number by shifting right.  This
  702.             routine is used only by 2NORMALIZE (and ultimately
  703.             by FLN).  The result is really a floating point
  704.             number on the parameter stack.
  705.  
  706. {
  707. \ 2NORMALIZE                                    08:20 31Oct88RS)
  708.  
  709. : 2NORMALIZE   ( d1 -- d2 )
  710.         $7FFF AND DUP $0FF >
  711.         IF      3NORMALIZE
  712.         ELSE    DUP $80 <
  713.                 IF      4NORMALIZE
  714.                 ELSE
  715.                         $3F80 OR
  716.                 THEN
  717.         THEN ;
  718.  
  719. }
  720.                                                 08:27 31Oct88RS)
  721.  
  722. 2NORMALIZE  Normalize the double number, performing either a
  723.             left or right shift, as required.
  724.  
  725. {
  726. \ FLOAT                                         07:15 19Oct88RS)
  727.  
  728. : FLOAT         ( F: -- r ; d -- )
  729.                 NO_INLINE
  730.                 2DUP OR 0= IF  FPUSH EXIT THEN
  731.                 DUP 0< $8000 AND >R DABS DUP $80 U<
  732.                 IF      NORMALIZE >R $7F AND $4B00 R> + OR
  733.                         R> OR FPUSH EXIT
  734.                 THEN    >NORMALIZE $4B00 + >R ROT DUP 0<
  735.                 IF      $8000 =  MROUND
  736.                         DUP $100 U< 0=
  737.                         IF      D2/ $7F AND R> $80 + OR
  738.                                 R> OR FPUSH EXIT
  739.                         THEN
  740.                 ELSE    DROP
  741.                 THEN    $7F AND R> OR
  742.                 R> OR FPUSH ;
  743.  
  744. }
  745.                                                 08:28 31Oct88RS)
  746.  
  747. FLOAT   Convert the double number on the parameter stack to a
  748.         floating point number on the floating point stack.
  749.  
  750. {
  751. \ DINTABS                                       05:49 01Nov88RS)
  752.  
  753. : DINTABS   ( F: r -- ; -- d flag )
  754.         NO_INLINE
  755.         FPOP DUP $7F80 AND DUP $3F80 <
  756.         IF      DROP 2DROP 0 0 0 EXIT  THEN
  757.         SWAP $7F AND $80 OR SWAP $4B00 - DUP 0<
  758.         IF      0 SWAP DO   D2/  $80 +LOOP  0
  759.         ELSE    $0400 MIN DUP
  760.                 IF      0 DO   D2*   $80 +LOOP
  761.                 ELSE    DROP
  762.                 THEN    DUP 0<
  763.         THEN  ;
  764.  
  765. }
  766.                                                 05:45 01Nov88RS)
  767.  
  768. DINTABS  Pop the top number from the floating point stack.
  769.          Take the absolute value and convert it to a double
  770.          number on the parameter stack.  If the resulting
  771.          number is positive, push a 0 onto the stack.
  772.          Otherwise, push a true flag (-1) onto the stack.
  773.  
  774. {
  775. \ INT   BMASK                                   07:23 15Oct88RS)
  776.  
  777. : INT   ( F: r -- ; -- d )
  778.         FDUP F0<
  779.         IF      DINTABS >R DNEGATE R>
  780.         ELSE    DINTABS
  781.         THEN    ABORT"  Out of range in INT " ;
  782.  
  783. TABLE BMASK     $0000 , $8000 , $C000 , $E000 , $F000 , $F800 ,
  784.                 $FC00 , $FE00 , $FF00 , $FF80 , $FFC0 , $FFE0 ,
  785.                 $FFF0 , $FFF8 , $FFFC , $FFFE , $FFFF ,
  786. END-TABLE
  787.  
  788. }
  789.                                                 05:54 01Nov88RS)
  790.  
  791. INT     Pop a floating point number from the f.p. stack and
  792.         convert it to a double number.
  793.  
  794. BMASK   An array of bit masks used to obtain the integer part
  795.         of a floating point number.
  796.  
  797. {
  798. \ FINT                                          07:24 15Oct88RS)
  799.  
  800. : FINT   ( F: r1 -- r2 )
  801.         NO_INLINE
  802.         FPOP DUP $7F80 AND DUP $3F80 <
  803.         IF      DROP $8000 AND NIP 0 SWAP FPUSH EXIT  THEN
  804.         DUP $4B00 <
  805.         IF      $3B00 - 2* $100 / 2* DUP $20 <
  806.                 IF      BMASK + @ AND
  807.                 ELSE    $1F AND BMASK + @ ROT AND SWAP
  808.                 THEN
  809.         ELSE    DROP
  810.         THEN
  811.         FPUSH ;
  812.  
  813. }
  814.                                                 05:54 01Nov88RS)
  815.  
  816. FINT    Convert the number at the top of the f.p. stack to its
  817.         integer part represented as a floating point number.
  818.  
  819. {
  820. \ ROUND1                                        07:24 15Oct88RS)
  821.  
  822. : ROUND1   ( grs d1 n1 -- d2 n2 )
  823.         >R ROT DUP 0<
  824.         IF      $8000 = MROUND
  825.                 DUP $0FF >
  826.                 IF      D2/ R> $80 +
  827.                 ELSE    R>
  828.                 THEN
  829.         ELSE    DROP R>
  830.         THEN ;
  831.  
  832. }
  833.                                                 05:58 01Nov88RS)
  834.  
  835. ROUND1  The input parameters represent a floating point number
  836.         broken into an exponent at the top, a double number,
  837.         and the GRS (guard, round, sticky) bits.  Round the
  838.         number according to GRS.
  839.  
  840. {
  841. \ ROUND2                                        07:24 15Oct88RS)
  842.  
  843. : ROUND2   ( d1 grs -- d2 n )
  844.         $8000 =
  845.         IF      SWAP $FFFE AND SWAP
  846.         THEN
  847.         DUP $0FF >
  848.         IF      D2/ $80
  849.         ELSE    0
  850.         THEN ;
  851.  
  852. }
  853.                                                 05:58 01Nov88RS)
  854.  
  855. ROUND2  Round the double number according to the GRS bits at
  856.         the top of the stack.
  857.  
  858. {
  859. \ Aux for F+                                    07:24 15Oct88RS)
  860.  
  861. : (F-X1=X2)   ( F: -- r ; d1 d2 -- )  ( F- : Equal exponents )
  862.         NO_INLINE
  863.         D- DUP 0<
  864.         IF      ( mantissa1 < mantissa2 )
  865.                 DNEGATE ZSIGN @ $8000 XOR ZSIGN !
  866.         ELSE    2DUP D0=  IF   EXIT   THEN
  867.         THEN   ZEXP @ $80 >
  868.         IF      ( Equal exponents, normal r1' )
  869.                 NORMALIZE ZEXP @ + DUP $80 <
  870.                 IF      ( denormalize )
  871.                         DENORMALIZE1
  872.                 ELSE    SWAP $7F  AND OR
  873.                 THEN
  874.          THEN ( ZSIGN @ OR ) ;
  875.  
  876. }
  877.                                                 06:02 01Nov88RS)
  878.  
  879. (F-X1=X2)  Auxilliary floating point subtraction of magnitudes
  880.            function for the case of equal exponents for the two
  881.            operands.
  882.  
  883. {
  884. \ Auxilliary for F+                             08:33 20Oct88RS)
  885.  
  886. : (F+X1=X2)   ( F: -- r ; d1 d2 -- )
  887.         D+ DUP $0FF >
  888.         IF      ( normal case )  OVER 1 AND
  889.                 IF      ( Round to even case )
  890.                         D1+ D2/ SWAP $FFFE AND SWAP
  891.                 ELSE    D2/
  892.                 THEN   ZEXP @ + DUP 0< ABORT"  Overflow in F+ "
  893.         ELSE    DUP $7F >
  894.                 IF      $7F AND ZEXP @ DUP 0=
  895.                         IF      DROP $80   THEN  OR
  896.                 THEN
  897.         THEN  ( ZSIGN @ OR ) ;
  898.  
  899. }
  900.                                                 06:01 01Nov88RS)
  901.  
  902. (F+X1=X2)  Auxillary floating point addition of magnitudes
  903.            function for the case of equal exponents of the
  904.            operands.
  905.  
  906. {
  907. \ Auxilliary for F+                             07:25 15Oct88RS)
  908.  
  909. : (1F-)   ( d1 d2 x1-x2 -- grs d3 n )
  910.         NO_INLINE
  911.         UNNORMALIZE NEGATE DUP >R
  912.         IF    D1+  THEN   D- R> -ROT        ( grs d3 )
  913.         DUP $7F >  IF   0 EXIT  THEN
  914.         T2* DUP $7F >
  915.         IF      $FF80 EXIT  THEN
  916.         T2* 0 $C800 $0100
  917.         DO      DROP DUP $7F >
  918.                 IF      I NEGATE LEAVE THEN
  919.                 D2* I NEGATE $80
  920.         +LOOP ;
  921.  
  922. }
  923.                                                 06:03 01Nov88RS)
  924.  
  925. (1F-)   An auxillary function for floating point subtraction of
  926.         magnitudes.
  927.  
  928. {
  929. \ Aux for F+                                    07:25 15Oct88RS)
  930.  
  931. : (F+-AUX)   ( d2 sx2 -- d3 sx2 x1-x2 flag )
  932.         DUP $7F80 AND DUP 0=
  933.         IF      DROP >R $7F AND R> $80 OR $80   THEN
  934.         ZEXP @ SWAP - DUP 0< ;
  935.  
  936. }
  937.                                                 06:05 01Nov88RS)
  938.  
  939. (F+-AUX)  An auxillary function for floating point addition.
  940.  
  941. {
  942. \ Aux for F+                                    07:26 15Oct88RS)
  943.  
  944. : (F-)   ( d1 d2 sx2 -- d3 )        ( Subtract magnitudes )
  945.         NO_INLINE
  946.         (F+-AUX)
  947.         IF      ( x2 > x1 )
  948.                 NEGATE SWAP DUP $7F80 AND ZEXP !
  949.                 $8000 AND ZSIGN !  >R 2SWAP R>
  950.         ELSE    ( x2 <= x1 )  NIP DUP 0=
  951.                 IF      ( Exponents are equal )
  952.                         DROP (F-X1=X2) EXIT
  953.                 THEN
  954.         THEN    ( d1 d2 x1-x2 )
  955.         (1F-) ZEXP @ + DUP $80 <
  956.         IF   DENORMALIZE2  ELSE   ROUND1 SWAP $7F AND OR  THEN ;
  957.  
  958. }
  959.                                                 06:06 01Nov88RS)
  960.  
  961. (F-)    Auxilliary function for the subtraction of magnitudes
  962.         of floating point numbers.
  963.  
  964. {
  965. \ Auxilliary for F+                             08:46 20Oct88RS)
  966.  
  967. : (F+)   ( d1 d2 sx2 -- d3 )
  968.         NO_INLINE
  969.         (F+-AUX)
  970.         IF      ( x2 > x1 )  NEGATE SWAP DUP $7F80 AND ZEXP !
  971.                 $8000 AND ZSIGN !  >R 2SWAP R>
  972.         ELSE    ( x2 <= x1 )   NIP DUP 0=
  973.                 IF      ( x1 = x2 ) DROP (F+X1=X2) EXIT  THEN
  974.         THEN    ( d1 d2 x1-x2 ) UNNORMALIZE
  975.         >R D+ DUP $0FF >   IF  D1+ R> 0=
  976.                 IF      SWAP $FFFC AND SWAP  THEN  D2/ $80
  977.         ELSE    R@ 0<
  978.                 IF   D1+ R> ROUND2  ELSE   R> DROP 0  THEN
  979.         THEN  ZEXP @ + DUP 0< ABORT"  Overflow in F+ " DUP
  980.         IF      SWAP $7F AND SWAP  THEN  OR ;
  981. }
  982. (F+)                                            06:11 01Nov88RS)
  983.  
  984. (F+)    Auxillary function for addition of floating point
  985.         numbers having the same sign.
  986.  
  987. {
  988. \  F+                                           07:26 15Oct88RS)
  989.  
  990. : F+   ( F: r1 r2 -- r3 )
  991.         NO_INLINE
  992.         FSP @ 4 + 2@ DUP $8000 AND ZSIGN !
  993.         DUP $7F80 AND DUP ZEXP !
  994.         IF      $7F AND $80 OR
  995.         ELSE    $7F AND 2DUP OR 0=
  996.                 IF   2DROP FNIP EXIT   THEN
  997.         THEN
  998.         FPOP0=  IF    2DROP 2DROP EXIT   THEN
  999.         FDROP DUP $7F AND $80 OR SWAP $FF80 AND DUP
  1000.         ZSIGN @ XOR 0<
  1001.         IF      (F-)  ELSE    (F+)   THEN
  1002.         ZSIGN @ XOR FPUSH ;
  1003.  
  1004. }
  1005. F+                                              06:11 01Nov88RS)
  1006.  
  1007. F+      Floating point addition.
  1008.  
  1009. {
  1010. \ F-  FIX  ZDEN  ZQUOT                          08:46 20Oct88RS)
  1011.  
  1012. : F-   ( F: r1 r2 -- r3 )
  1013.         FNEGATE F+ ;
  1014.  
  1015. : FIX   ( F: r -- ; -- d )
  1016.         FDUP F0<
  1017.         IF      FABS F0.5 F+ DINTABS >R DNEGATE R>
  1018.         ELSE    F0.5 F+ DINTABS
  1019.         THEN    ABORT"  Out of range in FIX " ;
  1020.  
  1021. 2VARIABLE ZDEN   2VARIABLE ZQUOT
  1022.  
  1023. }
  1024. F-  FIX  ZDEN  ZQUOT                            06:10 01Nov88RS)
  1025.  
  1026. F-      Floating point subtraction.
  1027.  
  1028. FIX     Pop a number from the floating point stack, convert it
  1029.         to a double number and push the result on the parameter
  1030.         stack.  Issue an error message if the number cannot be
  1031.         properly converted.
  1032.  
  1033. ZDEN    A variable for temporary results in f.p. division.
  1034. ZQUOT   A variable for temporary results in f.p. division.
  1035.  
  1036. {
  1037. \ Aux for F/                                    07:27 15Oct88RS)
  1038.  
  1039. : (1F/)   ( 0 d1 d2 -- t )
  1040.         ( Set quotient exponent, possible num adjust )
  1041.         FDROP DUP $7F AND $80 OR SWAP $7F80 AND DUP
  1042.         IF      XBIAS - NEGATE
  1043.         ELSE    DROP NORMALIZE
  1044.         THEN
  1045.         2/ ZEXP +!
  1046.         2OVER 2OVER DU< 0=
  1047.         IF      ZDEN 2!  T2/
  1048.         ELSE    ZDEN 2!  $-40 ZEXP +!
  1049.         THEN
  1050.         ZDEN 2@ DSHFT8 DROP ZDEN 2! ;  ( Shift den left by 8 )
  1051.  
  1052. }
  1053. \ Aux for F/                                    06:13 01Nov88RS)
  1054.  
  1055. (1F/)   An auxillary function for floating point division.
  1056.         Set the quotient exponent, possibly adjust the
  1057.         numerator.
  1058.  
  1059. {
  1060. \ Aux for F/                                    06:15 01Nov88RS)
  1061.  
  1062. : (2F/)   ( t -- d )
  1063.         ZDEN @ UM/MOD DUP ZQUOT !
  1064.         ZDEN 2+ @ UM* D-CY
  1065.         IF      ZDEN 2@ D+CY
  1066.                 IF      ZDEN 2@ D+ -2
  1067.                 ELSE    -1
  1068.                 THEN
  1069.                 ZQUOT +!
  1070.         THEN ;
  1071.  
  1072. }
  1073. \ Aux for F/                                    06:14 01Nov88RS)
  1074.  
  1075. (2F/)   An auxilliary function for floating point division.
  1076.         In this routine we obtain the most significant part
  1077.         of the quotient.
  1078.  
  1079. {
  1080. \ Aux for F/                                    15:24 14Oct88RLS
  1081.  
  1082. : (3F/)   ( 0 d1 n -- d2 )  ( Usual generate low quotient )
  1083.         UM/MOD DUP ZQUOT 2+ !
  1084.         ZDEN 2+ @ UM* D-CY
  1085.         IF      ZDEN 2@ D+CY
  1086.                 IF      ZDEN 2@ D+ -2
  1087.                 ELSE    -1
  1088.                 THEN
  1089.                 ZQUOT 2+ +!
  1090.         THEN ;
  1091.  
  1092. }
  1093.                                                 06:16 01Nov88RS)
  1094.  
  1095. (3F/)   Auxillary function for floating point division.  This
  1096.         routine is normally called to obtain the low order part
  1097.         of the quotient.
  1098.  
  1099. {
  1100. \ Aux for F/                                    15:23 14Oct88RLS
  1101.  
  1102. : (4F/)   ( 0 d1 n -- d2 )
  1103.         ( Gen low quotient for hi num = hi den )
  1104.         DROP NIP ZDEN 2+ @ SWAP 0 ZDEN @ 0 D+
  1105.         ZDEN 2+ @ 0 D- NIP 0<           \ TJZ 01/25/90 17:37:40.03 ADDED NIP
  1106.         IF      ZDEN 2@ D+CY
  1107.                 IF      ZDEN 2@ D+ -3
  1108.                 ELSE    -2
  1109.                 THEN
  1110.         ELSE    -1
  1111.         THEN
  1112.         ZQUOT 2+ ! ;
  1113.  
  1114. }
  1115. \ Aux for F/                                    06:18 01Nov88RS)
  1116.  
  1117. (4F/)   Auxilliary function for floating point division.
  1118.         Generate the low part of the quotient for the case of
  1119.         the high part of the numerator equal to the high part
  1120.         of the denominator.
  1121.  
  1122. {
  1123. \ SROUND                                        07:27 15Oct88RS)
  1124.  
  1125. : SROUND        ( d1 -- d2 )
  1126.                 NO_INLINE
  1127.                 DUP 0<
  1128.                 IF      ( Round up )  2DROP ZQUOT 2@ D1+ EXIT
  1129.                 THEN    D2* 2DUP ZDEN 2@ DU< 0=
  1130.                 IF      ( Do we round to even? )
  1131.                         ZDEN 2@ D=
  1132.                         IF      ( Yes )
  1133.                                 ZQUOT 2@ D1+
  1134.                                 SWAP $FFFE AND SWAP
  1135.                         ELSE    ZQUOT 2@ D1+
  1136.                         THEN
  1137.                 ELSE    2DROP ZQUOT 2@
  1138.                 THEN    ;
  1139.  
  1140. }
  1141.                                                 06:20 01Nov88RS)
  1142.  
  1143. SROUND  A rounding function used in floating point division.
  1144.  
  1145. {
  1146. \ Auxilliary for F/                             07:27 15Oct88RS)
  1147.  
  1148. : (UF/)   ( d1 -- d2 )   ( Generate unnormalized quotient )
  1149.         ZEXP @ 2* NEGATE $80 + UNNORMALIZE DUP 0<
  1150.         IF      $8000 = MROUND
  1151.         ELSE    DROP
  1152.         THEN ;
  1153.  
  1154. VARIABLE ZTEMP
  1155.  
  1156. }
  1157.                                                 06:22 01Nov88RS)
  1158.  
  1159. (UF/)   An auxilliary function used in floating point division
  1160.         to generate an unnormalized quotient.
  1161.  
  1162. ZTEMP   Another variable for temporary storage.
  1163.  
  1164. {
  1165. \  F/                                           08:47 20Oct88RS)
  1166.  
  1167. : F/   ( F: r1 r2 -- r3 )
  1168.         NO_INLINE
  1169.       0 0 FSP @ 4 + 2@ DUP FSP @ @ XOR $8000 AND ZSIGN !
  1170.       DUP $7F80 AND DUP 2/ ZEXP !
  1171.       IF      $7F AND $80 OR
  1172.       ELSE    $7F AND 2DUP OR 0=
  1173.               IF  2DROP FDROP FDROP 0 ZSIGN @ FPUSH EXIT
  1174.               THEN  NORMALIZE $80 + 2/ ZEXP !
  1175.       THEN  FPOP0= ABORT"  Floating Division by 0" (1F/) (2F/)
  1176.       ZDEN @ 2DUP U<  IF   (3F/)  ELSE  (4F/)  THEN
  1177.       ZEXP @ $40 <    IF  (UF/)  ELSE  SROUND
  1178.         ZEXP @ $3FC0 > ABORT"  Overflow in F/ "
  1179.               $7F AND ZEXP @ 2* OR
  1180.       THEN  ZSIGN @ OR FPUSH ;
  1181. }
  1182.                                                 06:31 01Nov88RS)
  1183.  
  1184. F/      The floating point division routine.  The rather sneaky
  1185.         technique used here is attributable to Roedy Green, the
  1186.         author of BBL/Abundance.  To do a division by a double
  1187.         number, shift the denominator until the m.s. bit is
  1188.         set.  Use the high order part to obtain the first
  1189.         approximation, along with its remainder.  If neccessary,
  1190.         make one or two stages of correction to obtain the
  1191.         exact high order part and temporary remainder.  Repeat
  1192.         the process to obtain the low order part.
  1193.  
  1194. {
  1195. \ (F*CLEANUP)    Aux for F*                     08:47 20Oct88RS)
  1196.  
  1197. : (F*CLEANUP)   ( t1 flag -- d2 )
  1198.                 NO_INLINE
  1199.                 IF    ZEXP @ 2* NEGATE $80 + UNNORMALIZE DUP 0<
  1200.                         IF    $8000 = ZTEMP @ 0= AND MROUND
  1201.                         ELSE  DROP
  1202.                         THEN  EXIT            \ leave here
  1203.                 THEN    ROT DUP 0<
  1204.                 IF      $8000 = ZTEMP @ 0= AND MROUND DUP $0FF >
  1205.                         IF      D2/ $40 ZEXP +!  ZEXP @ $3FC0 >
  1206.                                 ABORT"  Overflow in F* "
  1207.                         THEN
  1208.                 ELSE  DROP
  1209.                 THEN  $7F AND ZEXP @ 2* OR ;
  1210.  
  1211. }
  1212. \ (F*CLEANUP)    Aux for F*                     06:32 01Nov88RS)
  1213.  
  1214. (F*CLEANUP)  Auxillary function for F*
  1215.  
  1216. {
  1217. \ 1F*  Aux for F*                               07:30 15Oct88RS)
  1218.  
  1219. : 1F*   ( n1 -- n1 )
  1220.      DUP $7F80 AND DUP 2/ $40 + ZEXP !
  1221.      IF      $7F AND $80 OR
  1222.      ELSE    $7F AND NORMALIZE $80 + 2/ ZEXP !
  1223.      THEN ;
  1224.  
  1225. }
  1226.                                                 06:32 01Nov88RS)
  1227.  
  1228. 1F*     Auxillary function for F*
  1229.  
  1230. {
  1231. \  F*  Floating point multiply                  08:47 20Oct88RS)
  1232.  
  1233. : F*   ( F: r1 r2 -- r3 )
  1234.         NO_INLINE
  1235.      FPOP FSP @ @ OVER XOR $8000 AND ZSIGN !  $7FFF AND 2DUP D0=
  1236.      IF      2DROP FDROP 0 ZSIGN @ FPUSH EXIT  THEN
  1237.      1F* FPOP $7FFF AND 2DUP D0=
  1238.      IF      2DROP 2DROP 0 ZSIGN @ FPUSH EXIT  THEN
  1239.      DUP $7F80 AND DUP
  1240.      IF      XBIAS - 2/ ZEXP +! $7F AND $80 OR
  1241.      ELSE    $7F AND NORMALIZE $80 + XBIAS - 2/ ZEXP +!
  1242.      THEN   Y* DUP 0< 0=
  1243.      IF      T2* $-40 ZEXP +!  THEN
  1244.      ROT ZTEMP !  DSHFT8 ZEXP @ DUP $3FC0 > ABORT"   Overflow in F* "
  1245.      $40 < (F*CLEANUP)
  1246.      ZSIGN @ OR FPUSH ;
  1247.  
  1248. }
  1249. F*                                              06:32 01Nov88RS)
  1250.  
  1251. F*      Floating point multiplication function.
  1252.  
  1253. {
  1254. \ F**+N   Raise fp to positive integer power.   07:30 15Oct88RS)
  1255.  
  1256. : F**+N     ( F: r1 -- r2 ; n -- )
  1257.         $7FFF AND >R F1.0
  1258.         BEGIN R@ 1 AND
  1259.               IF   FOVER F*  THEN
  1260.               R> 2/ DUP
  1261.         WHILE >R FSWAP FDUP F* FSWAP
  1262.         REPEAT
  1263.         DROP FNIP ;
  1264.  
  1265. }
  1266.                                                 06:34 01Nov88RS)
  1267.  
  1268. F**+N   Raise the floating point number at the top of the f.p.
  1269.         stack to the positive integer power at the top of the
  1270.         parameter stack.
  1271.  
  1272. {
  1273. \ F**N  F**N*  MF**2                            07:31 15Oct88RS)
  1274.  
  1275. : F**N     ( F: r1 -- r2 ; n -- )  ( r1^n )
  1276.         DUP 0<
  1277.         IF    ABS F**+N F1.0 FSWAP F/
  1278.         ELSE  F**+N
  1279.         THEN ;
  1280.  
  1281. : F**N*    ( F: r1 r2 -- r3 ; n -- ) ( r1 * [r2^n] )
  1282.         DUP 0< IF   ABS F**+N F/ ELSE    F**+N F*  THEN ;
  1283.  
  1284. : MF**2   ( xlo lhi -- x^2lo x^2hi )
  1285.         DUP >R UM* NIP 0 D2*
  1286.         R> DUP UM* D+ ;
  1287.  
  1288. }
  1289. \ F**N  F**N*  MF**2                            06:39 01Nov88RS)
  1290.  
  1291. F**N    Raise the number at the top of the f.p. stack to the
  1292.         power specified at the top of the parameter stack.
  1293.  
  1294. F**N    Raise the number at the top of the f.p. stack to the
  1295.         power specified at the top of the parameter stack,
  1296.         then multiply by the number second on the f.p. stack.
  1297. MF**2   Square the double number mantissa on the parameter
  1298.         stack.
  1299.  
  1300. {
  1301. \  D2**N                                        07:31 15Oct88RS)
  1302.  
  1303. TABLE 2**NTAB   $0001 , $0002 , $0004 , $0008 , $0010 , $0020 ,
  1304.                 $0040 , $0080 , $0100 , $0200 , $0400 , $0800 ,
  1305.                 $1000 , $2000 , $4000 , $8000 ,
  1306. END-TABLE
  1307.  
  1308. : D2**N    ( n -- d )
  1309.         $1F AND DUP $10 <
  1310.         IF      2* 2**NTAB + @ 0
  1311.         ELSE    $0F AND 2* 2**NTAB + @ 0 SWAP
  1312.         THEN ;
  1313.  
  1314. 2VARIABLE DROOT
  1315.  
  1316. }
  1317.                                                 06:47 01Nov88RS)
  1318.  
  1319. 2**NTAB  A table of 2 raised to various powers.
  1320.  
  1321. D2**N    Return a double number representing 2 raised to the
  1322.          power specified.
  1323.  
  1324. DROOT    A double number variable for temporary use with
  1325.          square roots.
  1326.  
  1327. {
  1328. \ Auxilliaries for Square Root.                 14:18 14Oct88RLS
  1329.  
  1330. : SQRTSTEP   ( d1 n1 -- d2 n2 )
  1331.         2* >R D2* D2* DUP R@ >
  1332.         IF      R@ 1 OR -  R> 2 OR
  1333.         ELSE    R> THEN ;
  1334.  
  1335. : DSQRTSTEP   ( d1 d2 -- d3 d4 )
  1336.         D2* DROOT 2!  D2* D2* 2DUP DROOT 2@ 2SWAP D<
  1337.         IF      DROOT 2@ D1+ D-
  1338.                 DROOT 2@ D1+ D1+
  1339.         ELSE
  1340.                 DROOT 2@
  1341.         THEN ;
  1342.  
  1343. }
  1344.                                                 06:48 01Nov88RS)
  1345.  
  1346. SQRTSTEP  Auxillary function for FSQRT .
  1347.  
  1348. DSQRTSTEP  Another auxillary function for FSQRT .
  1349.  
  1350. {
  1351. \ Auxilliary Square Root function.              08:45 20Oct88RS)
  1352.  
  1353. : FSQRT1   ( F: r1 -- ; -- n1 n2 n3 )
  1354.         NO_INLINE
  1355.         FPOP0=  IF  FPUSH EXIT  THEN  DUP 0<
  1356.         ABORT"  Negative argument for FSQRT "
  1357.         DUP $7F80 AND DUP ZEXP !
  1358.         IF      $7F AND $80 OR
  1359.         ELSE    $7F AND NORMALIZE $80 + ZEXP +!
  1360.         THEN   DSHFT8 DROP  ZEXP @ XBIAS - DUP $80 AND
  1361.         IF      $80 - 2/ XBIAS + ZEXP !  0 D2* D2*
  1362.         ELSE    2/ XBIAS + ZEXP !  0 D2*
  1363.         THEN
  1364.         1- 2  7 0
  1365.         DO      SQRTSTEP  LOOP ;
  1366.  
  1367. }
  1368.                                                 06:48 01Nov88RS)
  1369.  
  1370. FSQRT1  Yet another auxillary function for FSQRT .
  1371.  
  1372. {
  1373. \ FSQRT                                         07:32 15Oct88RS)
  1374.  
  1375. : FSQRT   ( F: r1 -- r2 )
  1376.         FSQRT1 ROT DROP  5 0  DO  SQRTSTEP  LOOP
  1377.         ROT DROP 0 SWAP 0  $0C 0
  1378.         DO      DSQRTSTEP  LOOP
  1379.         D2/ OVER 1 AND
  1380.         IF      D2/ 2SWAP D0=
  1381.                 IF      D1+ SWAP $FFFE AND SWAP
  1382.                 ELSE    D1+
  1383.                 THEN
  1384.         ELSE    D2/ 2SWAP 2DROP
  1385.         THEN
  1386.         DUP $0FF >
  1387.         IF      D2/ -1 ZEXP +!   THEN
  1388.         $7F AND ZEXP @ OR FPUSH ;
  1389. }
  1390.                                                 06:49 01Nov88RS)
  1391.  
  1392. FSQRT   Replace the number at the top of the f.p. stack with
  1393.         its square root.
  1394.  
  1395. {
  1396. \ Tables for logarithms.                        07:32 15Oct88RS)
  1397.  
  1398. TABLE LOGTAB1
  1399.         $0000 , $0000 ,   $00FC , $14D8 ,   $01F0 , $A30C ,
  1400.         $02DE , $1A51 ,   $03C4 , $E0EE ,   $04A5 , $54BE ,
  1401.         $057F , $CC1C ,   $0654 , $96A7 ,   $0723 , $FDF2 ,
  1402.         $07EE , $461B ,   $08B3 , $AE56 ,   $0974 , $715D ,
  1403.         $0A30 , $C5E1 ,   $0AE8 , $DEE0 ,   $0B9C , $EBFB ,
  1404.         $0C4D , $19C3 ,   $0CF9 , $91F6 ,   $0D22 , $7BBE ,
  1405.         $0E47 , $FBE4 ,
  1406. END-TABLE
  1407.  
  1408. TABLE LOGTAB2
  1409.         $0000 , $0000 ,   $0208 , $2BB1 ,   $0421 , $662D ,
  1410.         $064C , $D797 ,   $088B , $C741 ,   $0ADF , $A036 ,
  1411.         $0D49 , $F69E ,   $0FCC , $8E36 ,
  1412. END-TABLE
  1413.  
  1414. }
  1415.                                                 06:50 01Nov88RS)
  1416.  
  1417. LOGTAB1  Auxillary table used for logarithm functions.
  1418.  
  1419. LOGTAB2  Auxillary table used for logarithm functions.
  1420.  
  1421. {
  1422. \ Auxilliary functions                          08:09 31Oct88RS)
  1423.  
  1424. : YLN2*   ( n -- d )
  1425.         2* DUP >R $B172 UM* R> $17F8 UM*
  1426.         $8000 0 D+ NIP 0 D+ 2NORMALIZE
  1427.         DUP IF  $80 -  THEN ;
  1428.  
  1429. : XLN2*   ( n -- d )  ( Multiply a shifted exponent by ln 2 )
  1430.         DUP 0<
  1431.         IF      ABS YLN2* $8000 OR
  1432.         ELSE    YLN2*
  1433.         THEN ;
  1434.  
  1435. }
  1436.                                                 08:12 31Oct88RS)
  1437.  
  1438. YLN2*   Multiply the positive shifted exponent by the natural
  1439.         logarithm of 2.
  1440.  
  1441. XLN2*   Multiply the shifted exponent by the natural log of 2.
  1442.  
  1443. {
  1444. \ Auxilliary function for FLN                   07:32 15Oct88RS)
  1445.  
  1446. : FLN+   ( F: -- r ; d1 n -- )
  1447.         $3F80 - XLN2* FPUSH
  1448.         DUP $FC AND DUP $7C AND >R  $100 * >R      ( R: div 4J )
  1449.         3 AND >R 0 SWAP D2/ D2/ D2/ $1FFF AND
  1450.         0 R> D2/ D2/ D2/ DROP OR
  1451.         R@ UM/MOD SWAP 0 SWAP R> UM/MOD NIP SWAP
  1452.         DUP >R $20 R@ 0<  IF  1-  THEN
  1453.         R@ FRACT* $0555 SWAP - DUP 0=
  1454.         IF      DROP R@
  1455.         ELSE    NEGATE R@ FRACT*
  1456.         THEN    R> FRACT* 0 SWAP DSHFT8 ROT DROP D2* D2*
  1457.         D- DSHFT8 ROT DROP R> LOGTAB1 + 2@ D+
  1458.         2NORMALIZE DUP IF  $300 -  THEN
  1459.         FPUSH F+ ;
  1460. }
  1461.                                                 06:56 01Nov88RS)
  1462.  
  1463. FLN+    Auxillary function for the natural logarithm function
  1464.         used when the mantissa is less than or equal to 1.5625 .
  1465.  
  1466. {
  1467. \ Auxilliary function for Logarithm             07:33 15Oct88RS)
  1468.  
  1469. : FLN-   ( F: -- r ; d n -- )
  1470.         $3F00 - XLN2* FPUSH DUP $F8 AND 8 + $100 * >R
  1471.         0 $100 2SWAP D- DUP $F8 AND 2/ R> SWAP >R >R
  1472.         7 AND >R 0 SWAP D2/ D2/ D2/ D2/ $FFF AND
  1473.         0 R> D2/ D2/ D2/ D2/ DROP OR R@ DUP
  1474.         IF  UM/MOD SWAP 0 SWAP R> UM/MOD NIP SWAP
  1475.         ELSE    R> 2DROP
  1476.         THEN   DUP >R 2DUP MF**2 DUP R@ $033 FRACT* $400 +
  1477.         R@ FRACT* $5555 + R> FRACT* FRACT* 0 SWAP D2/ D2/ D2/
  1478.         D+ D2/ D2/ D2/ D2/ D2/ D+  D2/ $7FFF AND
  1479.         D2/ D2/ D2/ D2/ D2/ R> LOGTAB2 + 2@ D+ 2DUP D0= 0=
  1480.         IF      2NORMALIZE $380 - $8000 OR  THEN
  1481.         FPUSH F+ ;
  1482.  
  1483. }
  1484.                                                 06:56 01Nov88RS)
  1485.  
  1486. FLN-    Auxillary function for natural logarithm, used when
  1487.         the mantissa is greater than 1.5625 and less than 2.0 .
  1488.  
  1489. {
  1490. \ FLN   Natural Logarithm                       08:47 20Oct88RS)
  1491.  
  1492. : FLN   ( F: r1 -- r2 )
  1493.         NO_INLINE
  1494.         FPOP0=
  1495.         IF      CR ."  Zero argument for FLN "
  1496.                 2DROP -1 -1 FPUSH EXIT
  1497.         THEN  DUP 0< ABORT"  Negative argument for FLN "
  1498.         DUP $7F80 AND DUP 0=
  1499.         IF      1NORMALIZE
  1500.         ELSE    SWAP $7F AND $80 OR SWAP
  1501.         THEN
  1502.         OVER $C8 >
  1503.         IF      FLN-   ELSE    FLN+   THEN ;
  1504.  
  1505. }
  1506.                                                 06:57 01Nov88RS)
  1507.  
  1508. FLN     Replace the number at the top of the f.p. stack with
  1509.         its natural logarithm.
  1510.  
  1511. {
  1512. \ FLOG ( Common Logarithm ) and FPARTS          08:05 20Oct88RS)
  1513.  
  1514. : FLOG   ( F: r1 -- r2 )
  1515.         FLN FLOG10E F* ;
  1516.  
  1517. : FPARTS   ( F: r1 -- ; n -- d exp sign )   ( Aux for E. )
  1518.         NO_INLINE
  1519.         8 MIN 1 MAX FDUP F0=
  1520.         IF      FDROP DROP 0 0 0 0 EXIT  THEN
  1521.         FSP @ @ 0< >R
  1522.         FABS FDUP FLOG INT DROP
  1523.         2DUP - F10.0 F**N* F0.5 F+ FINT
  1524.         SWAP FDUP F10.0 F**N F< 0=
  1525.         IF      F10.0 F/ 1+  THEN
  1526.         >R INT R> R> ;
  1527.  
  1528. }
  1529.                                                 06:59 01Nov88RS)
  1530.  
  1531. FLOG    Replace the number at the top of the f.p. stack with
  1532.         its common (base 10) logarithm.
  1533.  
  1534. FPARTS  An auxillary function for  (E.) .
  1535.  
  1536. {
  1537. \ Numeric output E.                             14:28 08Nov88RS)
  1538.  
  1539. VARIABLE F#PLACES
  1540.  
  1541. : (E.)          ( F: r -- ; n -- addr cnt )
  1542.                 F#PLACES @
  1543.                 FPARTS BASE @ >R >R DECIMAL <# DUP ABS 0 # # 2DROP 0<
  1544.                 IF     $2D   ELSE $2B  THEN     ( Send "-" or "+" )
  1545.                 HOLD $45 HOLD                   ( Send the "E" )
  1546.                 F#PLACES @ 0 DO  #  LOOP        ( Send the fraction )
  1547.                 $2E HOLD R> 0<                  ( Send "." Check sign )
  1548.                 IF      $2D   ELSE  $20  THEN
  1549.                 HOLD                            ( Send "-" or space )
  1550.                 R> BASE ! #> ;                  ( Restore BASE )
  1551.  
  1552. : E.            ( F: r -- )
  1553.                 8 F#PLACES ! (E.) TYPE SPACE ;
  1554.  
  1555. : E.R           ( F: r -- ; places width -- )
  1556.                 >R 8 min 1 max F#PLACES ! (E.)   R> OVER - SPACES TYPE   ;
  1557.  
  1558. : .FS           ( -- )
  1559.                 FDEPTH  IF   CR FSP @ FDEPTH F#BYTES * 0 DO  DUP I + 2@
  1560.                 FPUSH E. F#BYTES  +LOOP  DROP  ELSE ."  Empty"  THEN ;
  1561.  
  1562. : (F.)          ( F: r -- ; n -- addr cnt )
  1563.                 F#PLACES @ FPARTS BASE @ >R >R DECIMAL <#
  1564.                 F#PLACES @ SWAP - 0max 0 ?DO # LOOP
  1565.                 $2E HOLD #S R> 0<               ( Send "." Check sign )
  1566.                 IF      $2D   ELSE  $20  THEN
  1567.                 HOLD                            ( Send "-" or space )
  1568.                 R> BASE ! #> ;                  ( Restore BASE )
  1569.  
  1570. : F.            ( F: r -- )
  1571.                 8 F#PLACES ! (F.) TYPE SPACE ;
  1572.  
  1573. : F.R           ( F: r -- ; places width -- )
  1574.                 >R 8 min 1 max F#PLACES ! (F.)   R> OVER - SPACES TYPE   ;
  1575.  
  1576. }
  1577.  
  1578. (E.)  E.                                        14:29 08Nov88RS)
  1579.  
  1580. (E.)    Auxillary function for E.
  1581.  
  1582. E.      The floating point output routine.
  1583.  
  1584. .FS     A utility for checking the contents of the F.P. stack.
  1585.  
  1586. {
  1587. \ Auxilliary finctions for numeric input.       14:05 14Oct88RLS
  1588.  
  1589. : Ee(?   ( n -- flag )
  1590.         DUP 40 =                     ( Check for left paren )
  1591.         IF    DROP -1
  1592.         ELSE  DUP 69 = SWAP 101 = OR  ( Check for "e" or "E" )
  1593.         THEN ;
  1594.  
  1595. : -?    ( addr1 -- addr2 flag )
  1596.         DUP 1+ C@ 45 =
  1597.         IF      1+ -1
  1598.         ELSE    0
  1599.         THEN ;
  1600.  
  1601. }
  1602.                                                 07:03 01Nov88RS)
  1603.  
  1604. E.(?    Auxillary function to test for exponential indicator
  1605.         in floating point input.  The indicator should be one
  1606.         of the following:   E  e  (
  1607.  
  1608. -?      Check the character at the specified address for
  1609.         ASCII - sign.  If found, increment the address pointer
  1610.         and return a true flag.  Otherwise, return a false flag.
  1611.  
  1612. {
  1613. \ Numeric Input Conversion                      11:27 07Nov88RS)
  1614.  
  1615. : (FNUMBER?)    ( a1 -- f1 ; F: -- r )  \ convert string a1 to floating point #
  1616.         0 0     ROT -? >R DUP 1+ C@ Ee(?
  1617.         IF      1+ ROT DROP 1 -ROT 0 >R   ( 1 to mantissa )
  1618.         ELSE    CONVERT DUP C@ 46 =          ( Check for "." )
  1619.                 IF      DUP >R CONVERT DUP R> - 1-
  1620.                 ELSE    0  THEN    >R
  1621.         THEN
  1622.         DUP C@ Ee(?
  1623.         IF      -? >R >R 0 0 R> CONVERT NIP
  1624.         ELSE    0 SWAP 0 >R THEN
  1625.         C@ DUP 0= SWAP DUP 32 = SWAP 41 = OR OR 0=
  1626.         IF      DROP 2DROP R>DROP R>DROP R>DROP F0.0 FALSE
  1627.         ELSE    R> IF   NEGATE  THEN  R> - BASE @ 0 FLOAT F**N
  1628.                 FLOAT F* R>  IF  FNEGATE  THEN
  1629.                 ( TRUE ) 1      \ return 1 to be compatible with SFLOAT
  1630.         THEN    ;
  1631.  
  1632. : $F#   ( a1 -- F: -- r )       \ convert string a1 to floating point #
  1633.         (FNUMBER?) 0= ABORT"  Bad Floating Point Input" ;
  1634.  
  1635. : F#    ( F: -- r )             \ convert string from input stream to
  1636.                                 \ a floating point nubmer
  1637.         $20 WORD $F# ;          EXECUTES> F#
  1638.  
  1639. FORTH TARGET >TARGET
  1640.  
  1641. }
  1642.                                                 07:07 01Nov88RS)
  1643.  
  1644. $F#     This function converts a counted string into a floating point
  1645.         number.
  1646.  
  1647. F#      This is the function used to get a floating point number
  1648.         from the input stream.  Usage examples follow:
  1649.  
  1650.                 F# -2.34
  1651.                 F# 34.5e6
  1652.                 F# -.1E-2
  1653.                 F# 2.34(5)
  1654.  
  1655.