\ Horner          Evaluation of a polynomial by the Horner method

\ Forth Scientific Library Algorithm #3

\ This routine evaluates an Nth order polynomial Y(X) at point X
\ Y(X) = \sum_i=0^N a[i] x^i                  (NOTE: N+1 COEFFICIENTS)
\ by the Horner scheme.  This algorithm minimizes the number of multiplications
\ required to evaluate the polynomial.
\ The implementation demonstrates the use of array aliasing.

\ This code conforms with ANS requiring:
\      1. The Floating-Point word set
\      2. The immediate word '%' which takes the next token
\         and converts it to a floating-point literal
\      3. Uses words 'Private:', 'Public:' and 'Reset_Search_Order'
\         to control the visibility of internal code.
\      4. The immediate word '&' to get the address of an array
\         function at either compile or run time.
\      5. Uses the words 'DArray' and '&!' to alias arrays.
\      6. Test code uses 'ExpInt' for real exponential integrals

\  (code for the dependencies , 3, 4, and 5 above are in the file 'fsl_util' )

\ This algorithm is described in many places, e.g.,
\ Conte, S.D. and C. deBoor, 1972; Elementary Numerical Analysis, an algorithmic
\ approach, McGraw-Hill, New York, 396 pages

\ (c) Copyright 1994 Everett F. Carter.  Permission is granted by the
\ author to use this software for any application provided this
\ copyright notice is preserved.

\ ( HORNER            V1.5           31 October 1994   EFC )
\ include lib/ansfloat.4th
\ include lib/fsl-util.4th
\ include lib/fexpint.4th
\ include lib/fexp.4th
\ include lib/flnflog.4th

[UNDEFINED] }horner [IF]
[UNDEFINED] fsl-array [IF] [ABORT] [THEN]

: }Horner ( &a n -- , f: x -- y[x] )
  SWAP VALUE ha{
        0 s>f
        -1 SWAP DO
                     FOVER F*
                     ha{ I } F@ F+
        -1 +LOOP
         FSWAP FDROP
;

false [IF]     \ test code =============================================
6 FLOAT [*] 1 [+] ARRAY ArrayZ{
5 FLOAT [*] 1 [+] ARRAY ArrayY{
5 FLOAT [*] 1 [+] ARRAY ArrayW{
FLOAT ArrayZ{ FSL-ARRAY
FLOAT ArrayY{ FSL-ARRAY
FLOAT ArrayW{ FSL-ARRAY
:THIS ArrayZ{ DOES> (FSL-ARRAY) ;
:THIS ArrayY{ DOES> (FSL-ARRAY) ;
:THIS ArrayW{ DOES> (FSL-ARRAY) ;

\ initialize with data for real exponential integral
: horner_init ( -- )

      S" 0.00107857" s>float ArrayZ{ 5 } F!
     S" -0.00976004" s>float ArrayZ{ 4 } F!
      S" 0.05519968" s>float ArrayZ{ 3 } F!
     S" -0.24991055" s>float ArrayZ{ 2 } F!
      S" 0.99999193" s>float ArrayZ{ 1 } F!
     S" -0.57721566" s>float ArrayZ{ 0 } F!

                   1  s>f     ArrayY{ 4 } F!
     S" 8.5733287401" s>float ArrayY{ 3 } F!
    S" 18.059016973"  s>float ArrayY{ 2 } F!
     S" 8.6347608925" s>float ArrayY{ 1 } F!
     S" 0.2677737343" s>float ArrayY{ 0 } F!

                    1  s>f     ArrayW{ 4 } F!
      S" 9.5733223454" s>float ArrayW{ 3 } F!
     S" 25.6329561486" s>float ArrayW{ 2 } F!
     S" 21.0996530827" s>float ArrayW{ 1 } F!
      S" 3.9584969228" s>float ArrayW{ 0 } F!

;

: local_exp ( -- , f: x -- expint[x] )

        FDUP
        f1.0 F< IF
                    FDUP ArrayZ{  5 }Horner
                    FSWAP FLN F-
                ELSE
                    FDUP  ArrayY{ 4 }Horner
                    FOVER ArrayW{ 4 }Horner
                    F/
                    FOVER F/
                    FSWAP -1 s>f F* FEXP F*
                THEN
;

: do_both ( -- , f: x -- )
   FDUP FDUP
   ." X = " F.
   ." Horner: " local_exp F.
   ." ExpInt: " fexpint   F.
   CR
;

\ compare ExpInt as coded in V1.0 against the general purpose
\ Horner routine

: horner_test ( -- )
    horner_init

    CR
      s" 0.2" s>float do_both
      s" 0.5" s>float do_both
      1 s>f do_both
      2 s>f do_both
      5 s>f do_both
     10 s>f do_both
;
fclear
10 set-precision
horner_test
[THEN]
[THEN]