PCOEF()

Fortran version: 

    SUBROUTINE PCOEF( N, XN, YN, CN )
      INTEGER, INTENT(IN   ) :: N            ! length of input vector
      REAL   , INTENT(IN   ) :: XN( N )      ! input vector of X-values
      REAL   , INTENT(IN   ) :: YN( N )      ! input vector of Y-values
      REAL   , INTENT(  OUT) :: CN( N )      ! output vector of polynomial coefficients

C version: none

Summary:

PCOEF() finds the array CN of coefficients for the polynomial P going through the sets of points <XN(k),YN(k)>, k = 1,...,N. Must have N l< 16; in practice, N should not exceed 8, and the points XN(k) should be well-spaced-out, or else numerical instabilities may arise. To evaluate P at X, evaluate the sum
SUM( CN( K ) * X**(K-1) ; K = 1,...,N )
Note that the following code ( the "Horner trick") is an ostensibly-efficient way to evaluate this in Fortran; however, it does create long dependency-chains that greatly slow down modern deeply-pipelined superscalar processors: 
        ...
        Y = CN( N )
        DO  K = N-1, 1, -1
            Y = X * Y  +  CN( K )
        END DO
        ...
For modern (deeply-pipelined superscalar) processors, it may be significantly more efficient to "unroll" this loop in order to obtain a form which pipelines better, e.g.: 
        ...
        K = MOD( N, 2 )
        IF ( K .EQ. 0 ) THEN
            Y = CN( N )
            M = N
        ELSE
            Y = 0.0
            M = N - 1
        END IF
        X2 = X*X       
        DO  K = M, 1, -2
            Y = X2 * Y  +  X2 * CN( K ) + X * CN( K-1 ) + CN( K-1 )
        END DO
        ...
(If N is fixed and known in advance, there are (processor specific) forms that are better yet, if the operation of polynomial-evaluation is to be performed many, many times.)

See also POLY().

Fortran Usage:

!! under construction !!
Previous: NAMEVAL

Next: PERMUTI

Up: Utility Routines

To: Models-3/EDSS I/O API: The Help Pages