\ Aitken             Aitken Interpolation        ACM Algorithm #70

\ Forth Scientific Library Algorithm #9

\ Evaluates the (N-1)th degree Lagrange polynomial given N data coordinates
\ and the value where interpolation is desired.  The polynomial is
\ generated by Aitken's iterative scheme.
\
\ This is an ANS Forth program 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. Uses the words 'DArray' and '&!' to alias arrays.
\      5. The immediate word '&' to get the address of an array
\         at either compile or run time.
\      6. Uses '}malloc' and '}free' to allocate and release memory
\         for dynamic arrays ( 'DArray' ).
\      7. The compilation of the test code is controlled by the VALUE TEST-CODE?
\         and the conditional compilation words in the Programming-Tools wordset
\      8. The second test uses 'Logistic' for the logistic function.


\ 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 .( AITKEN            V1.1.2           20 August 1994   EFC )

Private:

FLOAT DARRAY x{                 \ array pointer
FLOAT DARRAY fx{                \ scratch array

: }fcopy ( n &src &dest -- )

    ROT 0 DO
                OVER I } F@
                DUP  I } F!
          LOOP

          DROP DROP
;

Public:

: Aitken ( &x &f n -- ) ( f: u -- ans )
        ROT & x{  &!                       \ point to x{} data

        & fx{ OVER }malloc

        \ copy passed F data into the local F array
        \ because it gets modified by the subsequent calculation.

        DUP ROT fx{ }fcopy


        DUP 1- 0 DO
                    DUP I 1+ DO
                               FDUP x{ J } F@ F- fx{ i } F@ F*
                               FOVER x{ I } F@ F-
                               fx{ J } F@ F* F-

                               x{ I } F@ x{ J } F@ F- F/

                               fx{ I } F!

                              LOOP
             LOOP

        FDROP

        fx{ SWAP 1- } F@

        & fx{ }free
;


Reset_Search_Order


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


9 FLOAT ARRAY x{
9 FLOAT ARRAY y{

: A_coords1 ( -- )

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

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

           A_coords1

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

           FSIN FSWAP             \ get exact value for later

           x{  y{ 9 Aitken

          ."      interpolated value: " F.
          ."   exact value: " F. CR


;

: A_coords2 ( -- )

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

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

           A_coords2

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

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

           x{ y{ 5 Aitken

          ."      interpolated value: " F.
          ."   exact value: " F. CR


;


[THEN]