\ BNewt FNewt     Newton Interpolation           ACM Algorithm #168 and #169

\ Forth Scientific Library Algorithm #13

\ BNewt   Newton Interpolation with backward differences   ACM Algorithm #168
\ FNewt   Newton Interpolation with forward  differences   ACM Algorithm #169

\ These routines use the Newton interpolation scheme to interpolate
\ within N data points (X).  The interpolated data (dif) has to be passed
\ as either backward divided differences or forward divided differences
\ depending upon which routine is used.  The interpolation point (z) is expected
\ on the float stack.

\ The routine returns 3 floating point values, the interpolated value (p),
\ the estimate of the first derivative (d), and an estimated error bound (e).
\ In practice the error bound estimate does not seem very accurate.

\ This is an ANS Forth program requiring:
\      1. The Floating-Point word set
\      2. Uses words 'Private:', 'Public:' and 'Reset_Search_Order'
\         to control the visibility of internal code.
\      3. Uses the words 'DARRAY' and '&!' to alias arrays.
\      4. The immediate word '&' to get the address of an array
\         at either compile or run time.
\      5. Uses Divided Differences words, 'forward_div_difs' and
\         'backward_div_difs'
\      6. The compilation of the test code is controlled by the VALUE TEST-CODE?
\         and the conditional compilation words in the Programming-Tools wordset.
\      7. The test code uses the immediate word '%' which takes the next token
\         and converts it to a floating-point literal
\      8. The second test uses 'logistic' for the logistic equation.
\

\ Collected Algorithms from ACM, Volume 1 Algorithms 1-220,
\ 1980; Association for Computing Machinery Inc., New York,
\ ISBN 0-89791-017-6

\ (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.

CR .( Newton            V1.4           31 October 1994   EFC )

Private:

FLOAT DARRAY xn{
FLOAT DARRAY dif{

FVARIABLE z

: z1 ( i -- ) ( f: -- z1 )
      z F@
      xn{ SWAP } F@ F-
;

: Newton ( i -- ) ( f: e d p -- e d p )

             \ calculate new D
             DUP z1
             FROT F* FOVER F+

             \ calculate new P
             FSWAP
             DUP z1
             F* dif{ OVER } F@ F+

             \ calculate new E
             FROT
             z1 FABS
             F*
             FOVER FABS F+

             \ restore stack order
             FROT FROT

;

: scale_error ( -- ) ( f: e d p -- e d p )
          FROT
          1.5e0 F*
          FOVER FABS F-
          1.0e-5 F*
          FROT FROT
;


Public:

: FNewt ( &xn &dif n -- ) ( f: z -- e d p)
        >R
        & dif{ &!
        & xn{  &!
        R>

        z F!        
        0.0e0  0.0e0 0.0e0

        \ this loop direction is opposite what is in the original
        \ publication, but apparently the original pub is wrong here.
        0 DO
             I Newton
          LOOP

          \ scale error estimate
          scale_error
;

: BNewt ( &xn &dif n -- ) ( f: z -- e d p )
        >R
        & dif{ &!
        & xn{  &!
        R>

        z F!        
        0.0e0  0.0e0 0.0e0
        
        \ this loop direction is opposite what is in the original
        \ publication, but apparently the original pub is wrong here.
        1- 0 SWAP DO
                    I Newton
                   -1 +LOOP

          \ scale error estimate
          scale_error
;


Reset_Search_Order


TEST-CODE? [IF]     \ test code =============================================


9 FLOAT ARRAY x{
9 FLOAT ARRAY y{
9 FLOAT ARRAY dif{

: test_coords1 ( n -- )

     0 DO I S>F % 0.25 F*
          FDUP x{ I } F!
          FSIN y{ I } F!
     LOOP

;

: test_coords2 ( n -- )

      0 DO I 2* S>F % -4.0 F+
           FDUP x{ I } F!
           % 1.0 % 1.0 logistic  y{ I } F!
      LOOP
;

\ forward differences tests

: fnewt-test1 ( -- ) ( f: u -- )   \ u can be in the range 0..2 for this test

           9 test_coords1

           \ convert y{} to forward divided differences
           x{ y{ dif{ 9 forward_div_difs

           FDUP FDUP CR ." Interpolation point: " F. CR

           FSIN FSWAP             \ get exact value for later
           FDUP FCOS FSWAP

           x{ dif{ 9 FNewt

          ."      interpolated value, deriv, err: " F. F. F. CR
          ."      exact value: " FSWAP F. F. CR


;


: fnewt-test2 ( -- ) ( f: u -- )   \ u is in the range -4..4 for this test

           5 test_coords2

           \ convert y{} to forward divided differences
           x{  y{ dif{ 5 forward_div_difs

           FDUP FDUP CR ." Interpolation point: " F. CR

           % 1.0 % 1.0 logistic FSWAP             \ get exact value for later
           FDUP % 1.0 % 1.0 d_logistic FSWAP

           x{ dif{ 5 FNewt

          ."      interpolated value, deriv, err: " F. F. F. CR
          ."      exact value: " FSWAP F. F. CR


;


\ backwards differences tests

: bnewt-test1 ( -- , f: u -- )   \ u can be in the range 0..2 for this test

           9 test_coords1

           \ convert y{} to backward divided differences
           x{  y{ dif{ 9 backward_div_difs

           FDUP FDUP CR ." Interpolation point: " F. CR

           FSIN FSWAP             \ get exact value for later
           FDUP FCOS FSWAP

           x{  dif{ 9 BNewt

          ."      interpolated value, deriv, err: " F. F. F. CR
          ."      exact value: " FSWAP F. F. CR


;

: bnewt-test2 ( -- , f: u -- )   \ u is in the range -4..4 for this test

           5 test_coords2

           \ convert y{} to backward divided differences
           x{  y{ dif{ 5 backward_div_difs

           FDUP FDUP CR ." Interpolation point: " F. CR

           % 1.0 % 1.0 logistic FSWAP             \ get exact value for later
           FDUP % 1.0 % 1.0 d_logistic FSWAP

           x{  dif{ 5 BNewt

          ."      interpolated value, deriv, err: " F. F. F. CR
          ."      exact value: " FSWAP F. F. CR


;


[THEN]