# include # include # include # include # include # include "divdif.h" /******************************************************************************/ double *cheby_t_zero ( int n ) /******************************************************************************/ /* Purpose: CHEBY_T_ZERO returns zeroes of the Chebyshev polynomial T(N)(X). Discussion: The I-th zero of T(N)(X) is cos((2*I-1)*PI/(2*N)), I = 1 to N Licensing: This code is distributed under the GNU LGPL license. Modified: 24 May 2011 Author: John Burkardt Parameters: Input, int N, the order of the polynomial. Output, double CHEBY_T_ZERO[N], the zeroes of T(N)(X). */ { double angle; int i; static double pi = 3.141592653589793; double *z; z = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { angle = ( double ) ( 2 * i + 1 ) * pi / ( double ) ( 2 * n ); z[i] = cos ( angle ); } return z; } /******************************************************************************/ double *cheby_u_zero ( int n ) /******************************************************************************/ /* Purpose: CHEBY_U_ZERO returns zeroes of the Chebyshev polynomial U(N)(X). Discussion: The I-th zero of U(N)(X) is cos((I-1)*PI/(N-1)), I = 1 to N Licensing: This code is distributed under the GNU LGPL license. Modified: 24 May 2011 Author: John Burkardt Parameters: Input, int N, the order of the polynomial. Output, double CHEBY_U_ZERO[N], the zeroes of U(N)(X). */ { double angle; int i; static double pi = 3.141592653589793; double *z; z = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { angle = ( double ) ( i + 1 ) * pi / ( double ) ( n + 1 ); z[i] = cos ( angle ); } return z; } /******************************************************************************/ void data_to_dif ( int ntab, double xtab[], double ytab[], double diftab[] ) /******************************************************************************/ /* Purpose: DATA_TO_DIF sets up a divided difference table from raw data. Discussion: Space can be saved by using a single array for both the DIFTAB and YTAB dummy parameters. In that case, the difference table will overwrite the Y data without interfering with the computation. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Reference: Carl deBoor, A Practical Guide to Splines, Springer, 2001, ISBN: 0387953663, LC: QA1.A647.v27. Parameters: Input, int NTAB, the number of pairs of points (XTAB[I],YTAB[I]) which are to be used as data. Input, double XTAB[NTAB], the X values at which data was taken. These values must be distinct. Input, double YTAB[NTAB], the corresponding Y values. Output, double DIFTAB[NTAB], the divided difference coefficients corresponding to the input (XTAB,YTAB). */ { int i; int j; /* Copy the data values into DIFTAB. */ for ( i = 0; i < ntab; i++ ) { diftab[i] = ytab[i]; } /* Make sure the abscissas are distinct. */ for ( i = 0; i < ntab; i++ ) { for ( j = i + 1; j < ntab; j++ ) { if ( xtab[i] - xtab[j] == 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "DATA_TO_DIF - Fatal error!\n" ); fprintf ( stderr, " Two entries of XTAB are equal!\n" ); fprintf ( stderr, " XTAB[%d] = %f\n", i, xtab[i] ); fprintf ( stderr, " XTAB[%d] = %f\n", j, xtab[j] ); exit ( 1 ); } } } /* Compute the divided differences. */ for ( i = 1; i <= ntab - 1; i++ ) { for ( j = ntab - 1; i <= j; j-- ) { diftab[j] = ( diftab[j] - diftab[j-1] ) / ( xtab[j] - xtab[j-i] ); } } return; } /******************************************************************************/ void data_to_dif_display ( int ntab, double xtab[], double ytab[], double diftab[] ) /******************************************************************************/ /* Purpose: DATA_TO_DIF_DISPLAY sets up a divided difference table and prints out intermediate data. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Parameters: Input, int NTAB, the number of pairs of points (XTAB[I],YTAB[I]) which are to be used as data. Input, double XTAB[NTAB], the X values at which data was taken. These values must be distinct. Input, double YTAB[NTAB], the corresponding Y values. Output, double DIFTAB[NTAB], the divided difference coefficients corresponding to the input (XTAB,YTAB). */ { int i; int j; if ( !r8vec_distinct ( ntab, xtab ) ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "DATA_TO_DIF_DISPLAY - Fatal error!\n" ); fprintf ( stderr, " Two entries of XTAB are equal!\n" ); exit ( 1 ); } printf ( "\n" ); printf ( " The divided difference table:\n" ); printf ( "\n" ); printf ( " " ); for ( i = 0; i < ntab; i++ ) { printf ( "%10f ", xtab[i] ); } printf ( "\n" ); printf ( "\n" ); printf ( "%6d ", 0 ); for ( i = 0; i < ntab; i++ ) { printf ( "%10f ", ytab[i] ); } printf ( "\n" ); /* Copy the data values into DIFTAB. */ for ( i = 0; i < ntab; i++ ) { diftab[i] = ytab[i]; } /* Compute the divided differences. */ for ( i = 1; i <= ntab-1; i++ ) { printf ( "%6d ", i ); for ( j = ntab-1; i <= j; j-- ) { diftab[j] = ( diftab[j] - diftab[j-1] ) / ( xtab[j] - xtab[j-i] ); } for ( j = i; j < ntab; j++ ) { printf ( "%10f ", diftab[j] ); } printf ( "\n" ); } return; } /******************************************************************************/ void data_to_r8poly ( int ntab, double xtab[], double ytab[], double c[] ) /******************************************************************************/ /* Purpose: DATA_TO_R8POLY computes the coefficients of a polynomial interpolating data. Discussion: Space can be saved by using a single array for both the C and YTAB parameters. In that case, the coefficients will overwrite the Y data without interfering with the computation. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Reference: Carl deBoor, A Practical Guide to Splines, Springer, 2001, ISBN: 0387953663, LC: QA1.A647.v27. Parameters: Input, int NTAB, the number of data points. Input, double XTAB[NTAB], YTAB[NTAB], the data values. Output, double C[NTAB], the coefficients of the polynomial that passes through the data (XTAB,YTAB). C(0) is the constant term. */ { if ( !r8vec_distinct ( ntab, xtab ) ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "DATA_TO_R8POLY - Fatal error!\n" ); fprintf ( stderr, " Two entries of XTAB are equal.\n" ); exit ( 1 ); } data_to_dif ( ntab, xtab, ytab, c ); dif_to_r8poly ( ntab, xtab, c, c ); return; } /******************************************************************************/ void dif_antideriv ( int ntab, double xtab[], double diftab[], int *ntab2, double xtab2[], double diftab2[] ) /******************************************************************************/ /* Purpose: DIF_ANTIDERIV integrates a polynomial in divided difference form. Discussion: This routine uses the divided difference representation (XTAB, DIFTAB) of a polynomial to compute the divided difference representation (XTAB, ANTTAB) of the antiderivative of the polynomial. The antiderivative of a polynomial P(X) is any polynomial Q(X) with the property that d/dX Q(X) = P(X). This routine chooses the antiderivative whose constant term is zero. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Parameters: Input, int NTAB, the size of the input table. Input, double XTAB[NTAB], the abscissas for the divided difference table. Input, double DIFTAB[NTAB], the divided difference table. Output, int *NTAB2, the size of the output table, which is NTAB+1. Input, double XTAB2[NTAB2], the abscissas for the divided difference table for the antiderivative. Output, double DIFTAB2[NTAB2], the divided difference table for the antiderivative. */ { int i; double *xtab1; double *diftab1; /* Using a temporary copy of the difference table, shift the abscissas to zero. */ xtab1 = ( double * ) malloc ( ntab * sizeof ( double ) ); diftab1 = ( double * ) malloc ( ntab * sizeof ( double ) ); for ( i = 0; i < ntab; i++ ) { xtab1[i] = xtab[i]; } for ( i = 0; i < ntab; i++ ) { diftab1[i] = diftab[i]; } dif_shift_zero ( ntab, xtab1, diftab1 ); /* Append a final zero to XTAB. */ *ntab2 = ntab + 1; for ( i = 0; i < *ntab2; i++ ) { xtab2[i] = 0.0; } /* Get the antiderivative of the standard form polynomial. */ r8poly_ant_cof ( ntab, diftab1, diftab2 ); free ( xtab1 ); free ( diftab1 ); return; } /******************************************************************************/ void dif_append ( int ntab, double xtab[], double diftab[], double xval, double yval, int *ntab2, double xtab2[], double diftab2[] ) /******************************************************************************/ /* Purpose: DIF_APPEND adds a pair of data values to a divided difference table. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Parameters: Input, int NTAB, the size of the difference table. Input, double XTAB[NTAB], the abscissas of the table. Input, double DIFTAB[NTAB], the difference table. Input, double XVAL, the data abscissa to be added to the table. Input, double YVAL, the data value to be added to the table. Output, int *NTAB2, the updated size of the difference table. Output, double XTAB2[*NTAB2], the updated abscissas. Output, double DIFTAB2[*NTAB2], the updated difference table. */ { int i; *ntab2 = ntab + 1; /* Move the original data up one index. */ for ( i = *ntab2-1; 1 <= i; i-- ) { diftab2[i] = diftab[i-1]; xtab2[i] = xtab[i-1]; } /* Recompute the data. */ xtab2[0] = xval; diftab2[0] = yval; for ( i = 1; i < *ntab2; i++ ) { diftab2[i] = ( diftab2[i] - diftab2[i-1] ) / ( xtab2[i] - xtab2[0] ); } return; } /******************************************************************************/ void dif_basis ( int ntab, double xtab[], double *diftab ) /******************************************************************************/ /* Purpose: DIF_BASIS: all Lagrange basis polynomials in divided difference form. Discussion: The I-th Lagrange basis polynomial for a set of NTAB X values XTAB, L(I,NTAB,XTAB)(X) is a polynomial of order NTAB-1 which is zero at XTAB(J) for J not equal to I, and 1 when J is equal to I. The Lagrange basis polynomials have the property that the interpolating polynomial through a set of NTAB data points (XTAB,YTAB) may be represented as P(X) = Sum ( 1 <= I <= N ) YTAB(I) * L(I,NTAB,XTAB)(X) Higher order interpolation at selected points may be accomplished using repeated X values, and scaled derivative values. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Reference: Carl deBoor, A Practical Guide to Splines, Springer, 2001, ISBN: 0387953663, LC: QA1.A647.v27. Parameters: Input, int NTAB, the number of X data points XTAB, and the number of basis polynomials to compute. Input, double XTAB[NTAB], the X values upon which the Lagrange basis polynomials are to be based. Output, double DIFTAB[NTAB*NTAB], points to a list of NTAB * NTAB values, the set of divided difference tables, stored as consecutive rows. Logical row I of DIFTAB contains the table for the I-th Lagrange basis polynomial. */ { int i; int j; double *pointer1; double *pointer2; pointer1 = diftab; for ( i = 0; i < ntab; i++ ) { pointer2 = pointer1; for ( j = 0; j < ntab; j++ ) { if ( j == i ) { *pointer1 = 1.0; } else { *pointer1 = 0.0; } pointer1++; } data_to_dif ( ntab, xtab, pointer2, pointer2 ); } return; } /******************************************************************************/ void dif_basis_deriv ( int nd, double xd[], double xdp[], double ddp[] ) /******************************************************************************/ /* Purpose: DIF_BASIS_DERIV: Lagrange basis derivative difference tables. Discussion: Given ND points XD, a Lagrange basis polynomial L(J)(X) is associated with each point XD(J). This function computes a table DDP(*,*) whose J-th column contains the difference table for the first derivative of L(J)(X). Licensing: This code is distributed under the GNU LGPL license. Modified: 01 June 2013 Author: John Burkardt Reference: Carl deBoor, A Practical Guide to Splines, Springer, 2001, ISBN: 0387953663, LC: QA1.A647.v27. Parameters: Input, int ND, the number of data points. Input, double XD[ND], the X values upon which the Lagrange basis polynomials are to be based. Output, double XDP[ND-1], the X values upon with the derivative difference table is based. In fact, these are all 0. Output, double DDP[(ND-1)*ND], the divided difference tables for all the Lagrange basis polynomials. Column J of DDP contains the table for basis polynomial associated with XD(J). */ { double *dd; int i; int j; double *yd; /* Process the vectors one column at a time. */ dd = ( double * ) malloc ( nd * sizeof ( double ) ); yd = ( double * ) malloc ( nd * sizeof ( double ) ); for ( j = 0; j < nd; j++ ) { /* Set the data. */ for ( i = 0; i < nd; i++ ) { yd[i] = 0.0; } yd[j] = 1.0; /* Compute the divided difference table. */ data_to_dif ( nd, xd, yd, dd ); /* Compute the divided difference table for the derivative. */ dif_deriv_table ( nd, xd, dd, xdp, ddp + j * ( nd - 1 ) ); } free ( dd ); free ( yd ); return; } /******************************************************************************/ void dif_basis_derivk ( int nd, double xd[], int k, double xdp[], double ddp[] ) /******************************************************************************/ /* Purpose: DIF_BASIS_DERIVK: Lagrange basis K-th derivative difference tables. Discussion: Given ND points XD, a Lagrange basis polynomial L(J)(X) is associated with each point XD(J). This function computes a table DDP(*,*) whose J-th column contains the difference table for the K-th derivative of L(J)(X). Licensing: This code is distributed under the GNU LGPL license. Modified: 03 June 2013 Author: John Burkardt Reference: Carl deBoor, A Practical Guide to Splines, Springer, 2001, ISBN: 0387953663, LC: QA1.A647.v27. Parameters: Input, int ND, the number of data points. Input, double XD[ND], the X values upon which the Lagrange basis polynomials are to be based. Input, int K, the index of the derivative. Output, double XDP[ND-1], the X values upon with the derivative difference table is based. In fact, these are all 0. Output, double DDP[(ND-1)*ND], the divided difference tables for all the Lagrange basis polynomials. Column J of DDP contains the table for basis polynomial associated with XD(J). */ { double *dd; int i; int j; double *yd; /* Process the vectors one column at a time. */ dd = ( double * ) malloc ( nd * sizeof ( double ) ); yd = ( double * ) malloc ( nd * sizeof ( double ) ); for ( j = 0; j < nd; j++ ) { /* Set the data. */ for ( i = 0; i < nd; i++ ) { yd[i] = 0.0; } yd[j] = 1.0; /* Compute the divided difference table. */ data_to_dif ( nd, xd, yd, dd ); /* Compute the divided difference table for the derivative. */ dif_derivk_table ( nd, xd, dd, k, xdp, ddp + j * ( nd - k ) ); } free ( dd ); free ( yd ); return; } /******************************************************************************/ void dif_basis_i ( int ival, int ntab, double xtab[], double diftab[] ) /******************************************************************************/ /* Purpose: DIF_BASIS_I: I-th Lagrange basis polynomial in divided difference form. Discussion: The I-th Lagrange basis polynomial for a set of NTAB X values XTAB, L(I,NTAB,XTAB)(X) is a polynomial of order NTAB-1 which is zero at XTAB(J) for J not equal to I, and 1 when J is equal to I. The Lagrange basis polynomials have the property that the interpolating polynomial through a set of NTAB data points (XTAB,YTAB) may be represented as P(X) = Sum ( 1 <= I <= N ) YTAB(I) * L(I,NTAB,XTAB)(X) Higher order interpolation at selected points may be accomplished using repeated X values, and scaled derivative values. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Reference: Carl deBoor, A Practical Guide to Splines, Springer, 2001, ISBN: 0387953663, LC: QA1.A647.v27. Parameters: Input, int IVAL, the index of the desired Lagrange basis polynomial. IVAL should be between 1 and NTAB. Input, int NTAB, the number of data points XTAB. Input, double XTAB[NTAB], the X values upon which the Lagrange basis polynomial is to be based. Output, double DIFTAB[NTAB], the divided difference table for the IVAL-th Lagrange basis polynomial. */ { int i; /* Check IVAL. */ if ( ival < 1 || ntab < ival ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "DIF_BASIS_I - Fatal error!\n" ); fprintf ( stderr, " IVAL must be between 1 and %d.\n", ntab ); fprintf ( stderr, " but your value is %d\n", ival ); exit ( 1 ); } /* Initialize DIFTAB to Delta(I,J). */ for ( i = 0; i <= ntab-1; i++ ) { diftab[i] = 0.0; } diftab[ival] = 1.0; /* Compute the IVAL-th Lagrange basis polynomial. */ data_to_dif ( ntab, xtab, diftab, diftab ); return; } /******************************************************************************/ void dif_deriv_table ( int nd, double xd[], double yd[], double xdp[], double ydp[] ) /******************************************************************************/ /* Purpose: DIF_DERIV_TABLE computes the divided difference table for a derivative. Licensing: This code is distributed under the GNU LGPL license. Modified: 01 June 2013 Author: John Burkardt Reference: Carl deBoor, A Practical Guide to Splines, Springer, 2001, ISBN: 0387953663, LC: QA1.A647.v27. Parameters: Input, int ND, the size of the input table. Input, double XD[ND], the abscissas for the divided difference table. Input, double YD[ND], the divided difference table. Output, double XDP[ND-1], the abscissas for the divided difference table for the derivative. Output, double YDP[ND-1], the divided difference table for the derivative. */ { int i; double *xd_temp; double *yd_temp; /* Using a temporary copy of the difference table, shift the abscissas to zero. */ xd_temp = ( double * ) malloc ( nd * sizeof ( double ) ); yd_temp = ( double * ) malloc ( nd * sizeof ( double ) ); for ( i = 0; i < nd; i++ ) { xd_temp[i] = xd[i]; } for ( i = 0; i < nd; i++ ) { yd_temp[i] = yd[i]; } dif_shift_zero ( nd, xd_temp, yd_temp ); /* Construct the derivative. */ for ( i = 0; i < nd - 1; i++ ) { xdp[i] = 0.0; } for ( i = 0; i < nd - 1; i++ ) { ydp[i] = ( double ) ( i + 1 ) * yd_temp[i+1]; } free ( xd_temp ); free ( yd_temp ); return; } /******************************************************************************/ void dif_derivk_table ( int nd, double xd[], double dd[], int k, double xdk[], double ddk[] ) /******************************************************************************/ /* Purpose: DIF_DERIVK_TABLE computes the divided difference table for K-th derivative. Licensing: This code is distributed under the GNU LGPL license. Modified: 01 June 2013 Author: John Burkardt Reference: Carl deBoor, A Practical Guide to Splines, Springer, 2001, ISBN: 0387953663, LC: QA1.A647.v27. Parameters: Input, int ND, the size of the input table. Input, double XD[ND], the abscissas for the divided difference table. Input, double DD[ND], the divided difference table. Input, int K, the index of the derivative. 0 <= K. Input, double XDK[ND-K], the abscissas for the divided difference table for the derivative. Output, double DDK[NDP], the divided difference table for the derivative. */ { double *dd_temp; int i; int j; int ndk; double *xd_temp; if ( k < 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "DIF_DERIVK_TABLE - Fatal error!\n" ); fprintf ( stderr, " K < 0.\n" ); exit ( 1 ); } if ( nd <= k ) { return; } /* Shift the abscissas to zero. */ ndk = nd; xd_temp = ( double * ) malloc ( ndk * sizeof ( double ) ); dd_temp = ( double * ) malloc ( ndk * sizeof ( double ) ); for ( i = 0; i < ndk; i++ ) { xd_temp[i] = xd[i]; } for ( i = 0; i < ndk; i++ ) { dd_temp[i] = dd[i]; } dif_shift_zero ( ndk, xd_temp, dd_temp ); /* Repeatedly differentiate. */ for ( j = 1; j <= k; j++ ) { ndk = ndk - 1; for ( i = 0; i < ndk; i++ ) { dd_temp[i] = ( double ) ( i + 1 ) * dd_temp[i+1]; } } for ( i = 0; i < ndk; i++ ) { ddk[i] = dd_temp[i]; xdk[i] = 0.0; } free ( xd_temp ); free ( dd_temp ); return; } /******************************************************************************/ void dif_print ( int ntab, double xtab[], double diftab[], char *title ) /******************************************************************************/ /* Purpose: DIF_PRINT prints the polynomial represented by a divided difference table. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Parameters: Input, int NTAB, the dimension of the arrays DIFTAB and XTAB. Input, double XTAB[NTAB], the X values for the polynomial. Input, double DIFTAB[NTAB], the divided difference table for the polynomial. Input, char *TITLE, a title. */ { int i; printf ( "\n" ); printf ( "%s\n", title ); printf ( "\n" ); printf ( " p(x) = %14f\n", diftab[0] ); for ( i = 1; i < ntab; i++ ) { printf ( " + ( x - %10f ) * ( %14f\n", xtab[i-1], diftab[i] ); } printf ( " " ); for ( i = 1; i < ntab; i++ ) { printf ( ")" ); } printf ( "\n" ); return; } /******************************************************************************/ void dif_shift_x ( int nd, double xd[], double yd[], double xv ) /******************************************************************************/ /* Purpose: DIF_SHIFT_X replaces one abscissa of a divided difference table with a new one. Discussion: This routine shifts the representation of a divided difference polynomial by dropping the last X value in XD, and adding a new X value to the beginning of the XD array, suitably modifying the coefficients stored in YD. The representation of the polynomial is changed, but the polynomial itself should be identical. Licensing: This code is distributed under the GNU LGPL license. Modified: 23 June 2011 Author: John Burkardt Reference: Carl deBoor, A Practical Guide to Splines, Springer, 2001, ISBN: 0387953663, LC: QA1.A647.v27. Parameters: Input, int ND, the number of divided difference coefficients, and the number of entries in XD. Input/output, double XD[ND], the X values used in the representation of the divided difference polynomial. After a call to this routine, the last entry of XD has been dropped, the other entries have shifted up one index, and XV has been inserted at the beginning of the array. Input/output, double YD[ND], the divided difference coefficients corresponding to the XD array. On output, this array has been adjusted. Input, double XV, a new X value which is to be used in the representation of the polynomial. On output, XD[0] equals XV and the representation of the polynomial has been suitably changed. Note that XV does not have to be distinct from any of the original XD values. */ { int i; /* Recompute the divided difference coefficients. */ for ( i = nd - 2; 0 <= i; i-- ) { yd[i] = yd[i] + ( xv - xd[i] ) * yd[i+1]; } /* Shift the XD values up one position and insert XV. */ for ( i = nd - 1; 0 < i; i-- ) { xd[i] = xd[i-1]; } xd[0] = xv; return; } /******************************************************************************/ void dif_shift_zero ( int nd, double xd[], double yd[] ) /******************************************************************************/ /* Purpose: DIF_SHIFT_ZERO shifts a divided difference table so that all abscissas are zero. Discussion: When the abscissas are changed, the coefficients naturally must also be changed. The resulting pair (XD, YD) still represents the same polynomial, but the entries in YD are now the standard polynomial coefficients. Licensing: This code is distributed under the GNU LGPL license. Modified: 01 November 2011 Author: John Burkardt Reference: Carl deBoor, A Practical Guide to Splines, Springer, 2001, ISBN: 0387953663, LC: QA1.A647.v27. Parameters: Input, int ND, the length of the XD and YD arrays. Input/output, double XD[ND], the X values that correspond to the divided difference table. On output, XD contains only zeroes. Input/output, double YD[ND], the divided difference table for the polynomial. On output, YD is also the coefficient array for the standard representation of the polynomial. */ { int i; int j; for ( j = 1; j <= nd; j++ ) { /* Recompute the divided difference coefficients. */ for ( i = nd - 2; 0 <= i; i-- ) { yd[i] = yd[i] - xd[i] * yd[i+1]; } /* Shift the XD values up one position and insert XV. */ for ( i = nd - 1; 0 < i; i-- ) { xd[i] = xd[i-1]; } xd[0] = 0.0; } return; } /******************************************************************************/ void dif_to_r8poly ( int ntab, double xtab[], double diftab[], double c[] ) /******************************************************************************/ /* Purpose: DIF_TO_R8POLY converts a divided difference table to a standard polynomial. Discussion: The vector DIFTAB, containing the divided difference polynomial coefficients is overwritten with the standard form polynomial coefficients, but the abscissa vector XTAB is unchanged. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Reference: Carl deBoor, A Practical Guide to Splines, Springer, 2001, ISBN: 0387953663, LC: QA1.A647.v27. Parameters: Input, int NTAB, the number of coefficients, and abscissas. Input, double XTAB[NTAB], the X values used in the divided difference representation of the polynomial. Input, double DIFTAB[NTAB], the divided difference table. Output, double C[NTAB], the standard form polyomial coefficients. C[0] is the constant term, and C[NTAB-1] is the coefficient of X**(NTAB-1). */ { int i; int j; for ( i = 0; i < ntab; i++ ) { c[i] = diftab[i]; } /* Recompute the divided difference coefficients. */ for ( j = 1; j <= ntab-1; j++ ) { for ( i = 1; i <= ntab-j; i++ ) { c[ntab-i-1] = c[ntab-i-1] - xtab[ntab-i-j] * c[ntab-i]; } } return; } /******************************************************************************/ double dif_val ( int ntab, double xtab[], double diftab[], double xv ) /******************************************************************************/ /* Purpose: DIF_VAL evaluates a divided difference polynomial at a point. Discussion: DATA_TO_DIF must be called first to set up the divided difference table. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Reference: Carl deBoor, A Practical Guide to Splines, Springer, 2001, ISBN: 0387953663, LC: QA1.A647.v27. Parameters: Input, integer NTAB, the number of divided difference coefficients in DIFTAB, and the number of points XTAB. Input, double XTAB[NTAB], the X values upon which the divided difference polynomial is based. Input, double DIFTAB[NTAB], the divided difference table. Input, double XV, a value of X at which the polynomial is to be evaluated. Output, double DIF_VAL, the value of the polynomial at XV. */ { int i; double yv; yv = diftab[ntab-1]; for ( i = 2; i <= ntab; i++ ) { yv = diftab[ntab-i] + ( xv - xtab[ntab-i] ) * yv; } return yv; } /******************************************************************************/ double *dif_vals ( int nd, double xd[], double yd[], int nv, double xv[] ) /******************************************************************************/ /* Purpose: DIF_VALS evaluates a divided difference polynomial at a set of points. Licensing: This code is distributed under the GNU LGPL license. Modified: 23 May 2011 Author: John Burkardt Reference: Carl deBoor, A Practical Guide to Splines, Springer, 2001, ISBN: 0387953663, LC: QA1.A647.v27. Parameters: Input, int ND, the order of the difference table. Input, double XD[ND], the X values of the difference table. Input, double YD[ND], the divided differences. Input, int NV, the number of evaluation points. Input, double XV[NV], the evaluation points. Output, double DIF_VALS[NV], the value of the divided difference polynomial at the evaluation points. */ { int i; int j; double *yv; yv = ( double * ) malloc ( nv * sizeof ( double ) ); for ( j = 0; j < nv; j++ ) { yv[j] = yd[nd-1]; for ( i = 2; i <= nd; i++ ) { yv[j] = yd[nd-i] + ( xv[j] - xd[nd-i] ) * yv[j]; } } return yv; } /******************************************************************************/ int i4_max ( int i1, int i2 ) /******************************************************************************/ /* Purpose: I4_MAX returns the maximum of two I4's. Licensing: This code is distributed under the GNU LGPL license. Modified: 29 August 2006 Author: John Burkardt Parameters: Input, int I1, I2, are two integers to be compared. Output, int I4_MAX, the larger of I1 and I2. */ { int value; if ( i2 < i1 ) { value = i1; } else { value = i2; } return value; } /******************************************************************************/ int i4_min ( int i1, int i2 ) /******************************************************************************/ /* Purpose: I4_MIN returns the smaller of two I4's. Licensing: This code is distributed under the GNU LGPL license. Modified: 29 August 2006 Author: John Burkardt Parameters: Input, int I1, I2, two integers to be compared. Output, int I4_MIN, the smaller of I1 and I2. */ { int value; if ( i1 < i2 ) { value = i1; } else { value = i2; } return value; } /******************************************************************************/ double *lagrange_rule ( int n, double x[] ) /******************************************************************************/ /* Purpose: LAGRANGE_RULE computes the weights of a Lagrange interpolation rule. Discussion: Given N abscissas X, an arbitrary function F(X) can be interpolated by a polynomial P(X) of order N (and degree N-1) using weights that depend only on X. Standard Lagrange interpolation can be rewritten into this form, which is more economical than evaluating the individual Lagrange basis polynomials. If we define W(I) = 1 / product ( 1 <= J <= N, J /= I ) ( X(J) - X(I) ) then P(XV) = sum ( 1 <= I <= N ) W(I) * F( X(I) ) / ( XV - X(I) ) / sum ( 1 <= I <= N ) W(I) / ( XV - X(I) ) except when XV = X(J), for some J, when we set: P(X(J)) = F(X(J)) Modified: 30 December 2004 Author: John Burkardt Reference: Jean-Paul Berrut, Lloyd Trefethen, Barycentric Lagrange Interpolation, SIAM Review, Volume 46, Number 3, September 2004, pages 501-517. Parameters: Input, int N, the order of the rule. Input, double X[N], the abscissas of the rule. Output, double LAGRANGE_RULE[N], the weights of the rule. */ { int i; int j; double *w; w = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { w[i] = 1.0; } for ( i = 0; i < n; i++ ) { for ( j = 0; j < n; j++ ) { if ( i != j ) { w[j] = w[j] / ( x[i] - x[j] ); } } } return w; } /******************************************************************************/ double lagrange_sum ( int n, double x[], double w[], double y[], double xv ) /******************************************************************************/ /* Purpose: LAGRANGE_SUM carries out a Lagrange interpolation rule. Discussion: It is assumed that LAGRANGE_RULE has already been called to compute the appropriate weights for the given set of abscissas. Licensing: This code is distributed under the GNU LGPL license. Modified: 24 May 2001 Author: John Burkardt Reference: Jean-Paul Berrut, Lloyd Trefethen, Barycentric Lagrange Interpolation, SIAM Review, Volume 46, Number 3, September 2004, pages 501-517. Parameters: Input, int N, the order of the rule. Input, double X[N], the abscissas of the rule. Input, double W[N], the weights of the rule. Input, double Y[N], the function values at the abscissas. Input, double XV, a point where an interpolated value is needed. Output, double LAGRANGE_SUM, the interpolated function value. */ { double bot; int i; double top; double yv; for ( i = 0; i < n; i++ ) { if ( xv == x[i] ) { yv = y[i]; return yv; } } top = 0.0; bot = 0.0; for ( i = 0; i < n; i++ ) { top = top + w[i] * y[i] / ( xv - x[i] ); bot = bot + w[i] / ( xv - x[i] ); } yv = top / bot; return yv; } /******************************************************************************/ double lagrange_val ( int n, double x[], double y[], double xv ) /******************************************************************************/ /* Purpose: LAGRANGE_VAL applies a naive form of Lagrange interpolation. Discussion: Given N abscissas X, an arbitrary function Y(X) can be interpolated by a polynomial P(X) of order N (and degree N-1) using Lagrange basis polynomials of degree N-1. Standard Lagrange interpolation can be rewritten into this form, which is more economical than evaluating the individual Lagrange basis polynomials. If we define L(I)(XV) = product ( 1 <= J <= N, J /= I ) ( XV - X(J) ) / ( X(I) - X(J) ) then P(XV) = sum ( 1 <= I <= N ) Y( X(I) ) * L(I)(XV) Applying this form of the interpolation rule directly involves about N^2 work. There are more efficient forms of the rule. Licensing: This code is distributed under the GNU LGPL license. Modified: 24 May 2011 Author: John Burkardt Parameters: Input, int N, the number of data points. Input, double X[N], the abscissas. Input, double Y[N], the function values at the abscissas. Input, double XV, a point where an interpolated value is needed. Output, double LAGRANGE_VAL, the interpolated function value. */ { int i; int j; double poly; double yv; yv = 0.0; for ( i = 0; i < n; i++ ) { poly = 1.0; for ( j = 0; j < n; j++ ) { if ( j != i ) { poly = poly * ( xv - x[j] ) / ( x[i] - x[j] ); } } yv = yv + y[i] * poly; } return yv; } /******************************************************************************/ void nc_rule ( int norder, double a, double b, double xtab[], double weight[] ) /******************************************************************************/ /* Purpose: NC_RULE computes the weights of a Newton-Cotes quadrature rule. Discussion: For the interval [A,B], the Newton-Cotes quadrature rule estimates Integral ( A <= X <= B ) F(X) dX using NORDER equally spaced abscissas XTAB(I) and a weight vector WEIGHT(I): Sum ( 1 <= I <= N ) WEIGHT(I) * F ( XTAB(I) ). For the CLOSED rule, the abscissas include the points A and B. For the OPEN rule, the abscissas do not include A and B. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Parameters: Input, int NORDER, the order of the rule. Input, double A, B, the left and right endpoints of the interval over which the quadrature rule is to be applied. Input, double XTAB[NORDER], the abscissas of the rule. Output, double WEIGHT[NORDER], the weights of the rule. */ { int i; double *poly_cof; double yvala; double yvalb; /* Allocate temporary space for POLY_COF. */ poly_cof = ( double * ) malloc ( norder * sizeof ( double ) ); for ( i = 1; i <= norder; i++ ) { /* Compute the Lagrange basis polynomial which is 1 at XTAB(I), and zero at the other nodes. */ r8poly_basis_1 ( i, norder, xtab, poly_cof ); /* Evaluate the antiderivative of the polynomial at the left and right endpoints. */ yvala = r8poly_ant_val ( norder-1, poly_cof, a ); yvalb = r8poly_ant_val ( norder-1, poly_cof, b ); weight[i-1] = yvalb - yvala; } /* Free up POLY_COF space. */ free ( poly_cof ); return; } /******************************************************************************/ void ncc_rule ( int norder, double xtab[], double weight[] ) /******************************************************************************/ /* Purpose: NCC_RULE computes the coefficients of a Newton-Cotes closed quadrature rule. Discussion: For the interval [-1,1], the Newton-Cotes quadrature rule estimates Integral ( -1 <= X <= 1 ) F(X) dX using NORDER equally spaced abscissas XTAB(I) and a weight vector WEIGHT(I): Sum ( 1 <= I <= N ) WEIGHT(I) * F ( XTAB(I) ). For the CLOSED rule, the abscissas include A and B. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Parameters: Input, int NORDER, the order of the rule. Output, double XTAB[NORDER], the abscissas of the rule. Output, double WEIGHT[NORDER], the weights of the rule. */ { double a; double b; int i; /* Compute a closed quadrature rule. */ a = -1.0; b = 1.0; for ( i = 1; i <= norder; i++ ) { xtab[i-1] = ( ( double ) ( norder - i ) * a + ( double ) ( i - 1 ) * b ) / ( double ) ( norder - 1 ); } nc_rule ( norder, a, b, xtab, weight ); return; } /******************************************************************************/ void nco_rule ( int norder, double xtab[], double weight[] ) /******************************************************************************/ /* Purpose: NCO_RULE computes the coefficients of a Newton-Cotes open quadrature rule. Discussion: For the interval [-1,1], the Newton-Cotes quadrature rule estimates Integral ( -1 <= X <= 1 ) F(X) dX using NORDER equally spaced abscissas XTAB(I) and a weight vector WEIGHT(I): Sum ( 1 <= I <= N ) WEIGHT(I) * F ( XTAB(I) ). For the OPEN rule, the abscissas do not include A and B. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Parameters: Input, int NORDER, the order of the rule. Output, double XTAB[NORDER], the abscissas of the rule. Output, double WEIGHT[NORDER], the weights of the rule. */ { double a; double b; int i; a = -1.0; b = 1.0; for ( i = 1; i <= norder; i++ ) { xtab[i-1] = ( ( double ) ( norder + 1 - i ) * a + ( double ) ( i ) * b ) / ( double ) ( norder + 1 ); } nc_rule ( norder, a, b, xtab, weight ); return; } /******************************************************************************/ void r8_swap ( double *x, double *y ) /******************************************************************************/ /* Purpose: R8_SWAP swaps two R8's. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Parameters: Input/output, double *X, *Y. On output, the values of X and Y have been interchanged. */ { double z; z = *x; *x = *y; *y = z; return; } /******************************************************************************/ void r8mat_transpose_print ( int m, int n, double a[], char *title ) /******************************************************************************/ /* Purpose: R8MAT_TRANSPOSE_PRINT prints an R8MAT, transposed. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. Licensing: This code is distributed under the GNU LGPL license. Modified: 28 May 2008 Author: John Burkardt Parameters: Input, int M, N, the number of rows and columns. Input, double A[M*N], an M by N matrix to be printed. Input, char *TITLE, a title. */ { r8mat_transpose_print_some ( m, n, a, 1, 1, m, n, title ); return; } /******************************************************************************/ void r8mat_transpose_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi, int jhi, char *title ) /******************************************************************************/ /* Purpose: R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT, transposed. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. Licensing: This code is distributed under the GNU LGPL license. Modified: 20 August 2010 Author: John Burkardt Parameters: Input, int M, N, the number of rows and columns. Input, double A[M*N], an M by N matrix to be printed. Input, int ILO, JLO, the first row and column to print. Input, int IHI, JHI, the last row and column to print. Input, char *TITLE, a title. */ { # define INCX 5 int i; int i2; int i2hi; int i2lo; int inc; int j; int j2hi; int j2lo; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); if ( m <= 0 || n <= 0 ) { fprintf ( stdout, "\n" ); fprintf ( stdout, " (None)\n" ); return; } for ( i2lo = i4_max ( ilo, 1 ); i2lo <= i4_min ( ihi, m ); i2lo = i2lo + INCX ) { i2hi = i2lo + INCX - 1; i2hi = i4_min ( i2hi, m ); i2hi = i4_min ( i2hi, ihi ); inc = i2hi + 1 - i2lo; fprintf ( stdout, "\n" ); fprintf ( stdout, " Row:" ); for ( i = i2lo; i <= i2hi; i++ ) { fprintf ( stdout, " %7d ", i - 1 ); } fprintf ( stdout, "\n" ); fprintf ( stdout, " Col\n" ); fprintf ( stdout, "\n" ); j2lo = i4_max ( jlo, 1 ); j2hi = i4_min ( jhi, n ); for ( j = j2lo; j <= j2hi; j++ ) { fprintf ( stdout, "%5d:", j - 1 ); for ( i2 = 1; i2 <= inc; i2++ ) { i = i2lo - 1 + i2; fprintf ( stdout, " %14f", a[(i-1)+(j-1)*m] ); } fprintf ( stdout, "\n" ); } } return; # undef INCX } /******************************************************************************/ void r8poly_ant_cof ( int n, double poly_cof[], double poly_cof2[] ) /******************************************************************************/ /* Purpose: R8POLY_ANT_COF integrates an R8POLY in standard form. Discussion: The antiderivative of a polynomial P(X) is any polynomial Q(X) with the property that d/dX Q(X) = P(X). This routine chooses the antiderivative whose constant term is zero. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Parameters: Input, int N, the order of the polynomial. Input, double POLY_COF[N], the polynomial coefficients. POLY_COF[0] is the constant term, and POLY_COF[N-1] is the coefficient of X**(N-1). Output, double POLY_COF2[N+1], the coefficients of the antiderivative polynomial, in standard form. The constant term is set to zero. */ { int i; /* Set the constant term. */ poly_cof2[0] = 0.0; /* Integrate the polynomial. */ for ( i = 1; i <= n; i++ ) { poly_cof2[i] = poly_cof[i-1] / ( double ) i; } return; } /******************************************************************************/ double r8poly_ant_val ( int n, double poly_cof[], double xval ) /******************************************************************************/ /* Purpose: R8POLY_ANT_VAL evaluates the antiderivative of an R8POLY in standard form. Discussion: The constant term of the antiderivative is taken to be zero. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Parameters: Input, int N, the order of the polynomial. Input, double POLY_COF[N], the polynomial coefficients. POLY_COF[0] is the constant term, and POLY_COF[N-1] is the coefficient of X**(N-1). Input, double XVAL, the point where the antiderivative is to be evaluated. Output, double R8POLY_ANT_VAL, the value of the antiderivative of the polynomial at XVAL. */ { int i; double value; value = 0.0; for ( i = n-1; 0 <= i; i-- ) { value = ( value + poly_cof[i] / ( double ) ( i + 1 ) ) * xval; } return value; } /******************************************************************************/ void r8poly_basis ( int ntab, double xtab[], double *poly_cof ) /******************************************************************************/ /* Purpose: R8POLY_BASIS computes all Lagrange basis polynomials in standard form. Discussion: The I-th Lagrange basis polynomial for a set of NTAB X values XTAB, L(I,NTAB,XTAB)(X) is a polynomial of order NTAB-1 which is zero at XTAB(J) for J not equal to I, and 1 when J is equal to I. The Lagrange basis polynomials have the property that the interpolating polynomial through a set of NTAB data points (XTAB,YTAB) may be represented as P(X) = Sum ( 1 <= I <= N ) YTAB(I) * L(I,NTAB,XTAB)(X) Higher order interpolation at selected points may be accomplished using repeated X values, and scaled derivative values. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Parameters: Input, int NTAB, the number of data points XTAB. Input, double XTAB[NTAB], the X values upon which the Lagrange basis polynomial is to be based. Output, double *POLY_COF, points to NTAB * NTAB values. The polynomial coefficients for the I-th Lagrange basis polynomial are stored in (logical) row I. POLY_COF[0,*] is the constant term, and POLY_COF[NTAB-1,*] is the coefficient of X**(NTAB-1). */ { int i; int j; double *pointer1; double *pointer2; pointer1 = poly_cof; for ( i = 0; i < ntab; i++ ) { pointer2 = pointer1; for ( j = 0; j < ntab; j++ ) { if ( j == i ) { *pointer1 = 1.0; } else { *pointer1 = 0.0; } pointer1++; } /* Compute the divided difference table for the IVAL-th Lagrange basis polynomial. */ data_to_dif ( ntab, xtab, pointer2, pointer2 ); /* Convert the divided difference table coefficients to standard polynomial coefficients. */ dif_to_r8poly ( ntab, xtab, pointer2, pointer2 ); } return; } /******************************************************************************/ void r8poly_basis_1 ( int ival, int ntab, double xtab[], double poly_cof[] ) /******************************************************************************/ /* Purpose: R8POLY_BASIS_1 computes the I-th Lagrange basis polynomial in standard form. Discussion: The I-th Lagrange basis polynomial for a set of NTAB X values XTAB, L(I,NTAB,XTAB)(X) is a polynomial of order NTAB-1 which is zero at XTAB(J) for J not equal to I, and 1 when J is equal to I. The Lagrange basis polynomials have the property that the interpolating polynomial through a set of NTAB data points (XTAB,YTAB) may be represented as P(X) = Sum ( 1 <= I <= N ) YTAB(I) * L(I,NTAB,XTAB)(X) Higher order interpolation at selected points may be accomplished using repeated X values, and scaled derivative values. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Parameters: Input, int IVAL, the index of the desired Lagrange basis polynomial. IVAL should be between 1 and NTAB. Input, int NTAB, the number of data points XTAB. Input, double XTAB[NTAB], the X values upon which the Lagrange basis polynomial is to be based. Output, double POLY_COF[NTAB], the polynomial coefficients for the IVAL-th Lagrange basis polynomial. */ { int i; /* Check IVAL. */ if ( ival < 1 || ntab < ival ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8POLY_BASIS_1 - Fatal error!\n" ); fprintf ( stderr, " IVAL must be between 1 and %d.\n", ntab ); fprintf ( stderr, " but your value is %d.\n", ival ); exit ( 1 ); } /* Initialize POLY_COF to the IVAL-th column of the identity matrix. */ for ( i = 0; i <= ntab-1; i++ ) { poly_cof[i] = 0.0; } poly_cof[ival-1] = 1.0; /* Compute the divided difference table for the IVAL-th Lagrange basis polynomial. */ data_to_dif ( ntab, xtab, poly_cof, poly_cof ); /* Convert the divided difference table coefficients to standard polynomial coefficients. */ dif_to_r8poly ( ntab, xtab, poly_cof, poly_cof ); return; } /******************************************************************************/ void r8poly_der_cof ( int n, double poly_cof[], double poly_cof2[] ) /******************************************************************************/ /* Purpose: R8POLY_DER_COF computes the coefficients of the derivative of a real polynomial. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Parameters: Input, int N, the order of the polynomial. Input, double POLY_COF[N], the coefficients of the polynomial to be differentiated. POLY_COF[0] is the constant term, and POLY_COF[N-1] is the coefficient of X**(N-1). Output, double POLY_COF2[N-1], the coefficients of the derivative of the polynomial. */ { int i; for ( i = 0; i < n; i++ ) { poly_cof2[i] = ( double ) ( i + 1 ) * poly_cof[i+1]; } return; } /******************************************************************************/ double r8poly_der_val ( int n, double poly_cof[], double xval ) /******************************************************************************/ /* Purpose: R8POLY_DER_VAL evaluates the derivative of a real polynomial in standard form. Discussion: A polynomial in standard form, with coefficients POLY_COF(*), may be written: P(X) = POLY_COF[0] + POLY_COF[1] * X ... + POLY_COF[N-1] * X**(N-1) Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Parameters: Input, int N, the order of the polynomial. Input, double POLY_COF[N], the polynomial coefficients. POLY_COF[0] is the constant term, and POLY_COF[N-1] is the coefficient of X**(N-1). Input, double XVAL, a value where the derivative of the polynomial is to be evaluated. Output, double R8POLY_DER_VAL, the value of the derivative of the polynomial at XVAL. */ { int i; double value; value = ( double ) ( n - 1 ) * poly_cof[n-1]; for ( i = n-2; 1 <= i; i-- ) { value = value * xval + ( double ) i * poly_cof[i]; } return value; } /******************************************************************************/ int r8poly_order ( int na, double a[] ) /******************************************************************************/ /* Purpose: R8POLY_ORDER returns the order of a polynomial. Discussion: The order of a polynomial is the degree plus 1. The order of a constant polynomial is 1. The order of the zero polynomial is debatable, but this routine returns the order as 1. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Parameters: Input, int NA, the order of the polynomial (ignoring zero coefficients). Input, double A[NA], the coefficients of the polynomial. Output, int R8POLY_ORDER, the degree of the polynomial. */ { int value; value = na; while ( 1 < value ) { if ( a[value-1] != 0.0 ) { return value; } value = value - 1; } return value; } /******************************************************************************/ void r8poly_print ( int n, double poly_cof[], char *title ) /******************************************************************************/ /* Purpose: R8POLY_PRINT prints out a real polynomial in standard form. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Parameters: Input, int N, the order of the polynomial. Input, double POLY_COF[N], the polynomial coefficients. POLY_COF[0] is the constant term and POLY_COF[N-1] is the coefficient of X**(N-1). Input, char *TITLE, a title to be printed first. TITLE may be blank. */ { int i; int n2; /* Print the title. */ printf ( "\n" ); printf ( "%s\n", title ); /* Determine the "real" degree of the polynomial. N2 is the index of the highest order nonzero coefficient. */ n2 = r8poly_order ( n, poly_cof ); printf ( "\n" ); printf ( " p(x) = " ); if ( poly_cof[n2-1] < 0.0 ) { printf ( "- " ); } else { printf ( " " ); } printf ( "%8f", fabs ( poly_cof[n2-1] ) ); if ( 2 <= n2 ) { printf ( " * x" ); } if ( 3 <= n2 ) { printf ( " ^ %d", n2 - 1 ); } printf ( "\n" ); for ( i = n2 - 1; 1 <= i; i-- ) { if ( poly_cof[i-1] < 0.0 ) { printf ( " - " ); } else if ( poly_cof[i-1] == 0.0 ) { break; } else if ( 0.0 <= poly_cof[i-1] ) { printf ( " + " ); } printf ( "%10f", fabs ( poly_cof[i-1] ) ); if ( 2 <= i ) { printf (" * x" ); } if ( 3 <= i ) { printf ( " ^ %d", i - 1 ); } printf ( "\n" ); } return; } /******************************************************************************/ void r8poly_shift ( double scale, double shift, int n, double poly_cof[] ) /******************************************************************************/ /* Purpose: R8POLY_SHIFT adjusts the coefficients of a polynomial for a new argument. Discussion: Assuming P(X) is a polynomial in the argument X, of the form: P(X) = C(N-1) * X**(N-1) + ... + C(1) * X + C(0), and that Z is related to X by the formula: Z = SCALE * X + SHIFT then this routine computes coefficients C for the polynomial Q(Z): Q(Z) = C(N-1) * Z**(N-1) + ... + C(1) * Z + C(0) so that: Q(Z(X)) = P(X) Example: P(X) = 2 * X**2 - X + 6 Z = 2.0 * X + 3.0 Q(Z) = 0.5 * Z**2 - 3.5 * Z + 12 Q(Z(X)) = 0.5 * ( 4.0 * X**2 + 12.0 * X + 9 ) - 3.5 * ( 2.0 * X + 3 ) + 12 = 2.0 * X**2 - 1.0 * X + 6 = P(X) Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Parameters: Input, double SHIFT, SCALE, the shift and scale applied to X, so that Z = SCALE * X + SHIFT. Input, int N, the order of the polynomial. Input/output, double POLY_COF[N]. On input, the coefficient array in terms of the X variable. On output, the coefficient array in terms of the Z variable. */ { int i; int j; for ( i = 1; i <= n; i++ ) { for ( j = i+1; j <= n; j++ ) { poly_cof[j-1] = poly_cof[j-1] / scale; } } for ( i = 1; i <= n; i++ ) { for ( j = n-1; i <= j; j-- ) { poly_cof[j-1] = poly_cof[j-1] - shift * poly_cof[j]; } } return; } /******************************************************************************/ double r8poly_val_horner ( int n, double poly_cof[], double xval ) /******************************************************************************/ /* Purpose: R8POLY_VAL_HORNER evaluates a real polynomial in standard form. Discussion: A polynomial in standard form, with coefficients POLY_COF(*), may be written: P(X) = POLY_COF[0] + POLY_COF[1] * X ... + POLY_COF[N-1] * X**(N-1) Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Parameters: Input, int N, the order of the polynomial. Input, double POLY_COF[N], the polynomial coefficients. POLY_COF[0] is the constant term, and POLY_COF[N-1] is the coefficient of X**(N-1). Input, double XVAL, a value where the polynomial is to be evaluated. Output, double R8POLY_VAL_HORNER, the value of the polynomial at XVAL. */ { int i; double value; value = poly_cof[n-1]; for ( i = n-2; 0 <= i; i-- ) { value = value * xval + poly_cof[i]; } return value; } /******************************************************************************/ int r8vec_distinct ( int n, double x[] ) /******************************************************************************/ /* Purpose: R8VEC_DISTINCT is true if the entries in an R8VEC are distinct. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Input, double X[N], the vector to be checked. Output, int R8VEC_DISTINCT is 1 if all N elements of X are distinct. */ { int i; int j; for ( i = 1; i <= n-1; i++ ) { for ( j = 1; j <= i - 1; j++ ) { if ( x[i] == x[j] ) { return 0; } } } return 1; } /******************************************************************************/ void r8vec_indicator ( int n, double a[] ) /******************************************************************************/ /* Purpose: R8VEC_INDICATOR sets an R8VEC to the indicator vector {1,2,3...}. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Parameters: Input, int N, the number of elements of A. Output, double A[N], the array to be initialized. */ { int i; for ( i = 0; i <= n-1; i++ ) { a[i] = ( double ) ( i + 1 ); } return; } /******************************************************************************/ void r8vec_print ( int n, double a[], char *title ) /******************************************************************************/ /* Purpose: R8VEC_PRINT prints an R8VEC. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Parameters: Input, int N, the number of components of the vector. Input, double A[N], the vector to be printed. Input, char *TITLE, a title to be printed first. TITLE may be blank. */ { int i; printf ( "\n" ); printf ( "%s\n", title ); printf ( "\n" ); for ( i = 0; i <= n - 1; i++ ) { printf ( "%6d %10f\n", i + 1, a[i] ); } return; } /******************************************************************************/ void roots_to_dif ( int nroots, double roots[], int *ntab, double xtab[], double diftab[] ) /******************************************************************************/ /* Purpose: ROOTS_TO_DIF sets a divided difference table for a polynomial from its roots. Discussion: This turns out to be a simple task, because of two facts: * The divided difference polynomial of one smaller degree which passes through the values ( ROOT(I), 0 ) is the zero polynomial, and hence has a zero divided difference table. * We want a polynomial of one degree higher, but we don't want it to pass through an addditional point. Instead, we specify that the polynomial is MONIC. This means that the divided difference table is almost the same as for the zero polynomial, except that there is one more pair of entries, an arbitrary X value, and a Y value of 1. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Parameters: Input, int NROOTS, is the number of roots. Input, double ROOTS[NROOTS], the roots of the polynomial. Output, int *NTAB, the size of the divided difference table. This is NROOTS+1 Output, double XTAB[NTAB], the abscissas of the table. Output, double DIFTAB[NTAB], the divided difference table. */ { int i; *ntab = nroots + 1; /* Build the appropriate difference table for the polynomial through ( ROOTS(I), 0 ) of degree NTAB-2. */ for ( i = 0; i < *ntab-1; i++ ) { diftab[i] = 0.0; } for ( i = 0; i < *ntab-1; i++ ) { xtab[i] = roots[i]; } /* Append the extra data to make a monic polynomial of degree NTAB-1 which is zero at the NTAB-1 roots. */ xtab[*ntab-1] = 0.0; diftab[*ntab-1] = 1.0; return; } /******************************************************************************/ void roots_to_r8poly ( int nroots, double roots[], int *nc, double c[] ) /******************************************************************************/ /* Purpose: ROOTS_TO_R8POLY converts polynomial roots to polynomial coefficients. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Parameters: Input, int NROOTS, the number of roots specified. Input, double ROOTS[NROOTS], the roots. Output, integer *NC, the order of the polynomial, which will be NROOTS+1. Output, double C[*NC], the coefficients of the polynomial. */ { int i; double *xtab; *nc = nroots + 1; /* Initialize C to (0, 0, ..., 0, 1). Essentially, we are setting up a divided difference table. */ xtab = ( double * ) malloc ( ( nroots + 1 ) * sizeof ( double ) ); for ( i = 0; i < *nc-1; i++ ) { xtab[i] = roots[i]; } xtab[*nc-1] = 0.0; for ( i = 0; i < *nc-1; i++ ) { c[i] = 0.0; } c[*nc-1] = 1.0; /* Convert to standard polynomial form by shifting the abscissas of the divided difference table to 0. */ dif_shift_zero ( *nc, xtab, c ); free ( xtab ); return; } /******************************************************************************/ int s_len_trim ( char *s ) /******************************************************************************/ /* Purpose: S_LEN_TRIM returns the length of a string to the last nonblank. Discussion: It turns out that I also want to ignore the '\n' character! Licensing: This code is distributed under the GNU LGPL license. Modified: 05 October 2014 Author: John Burkardt Parameters: Input, char *S, a pointer to a string. Output, int S_LEN_TRIM, the length of the string to the last nonblank. If S_LEN_TRIM is 0, then the string is entirely blank. */ { int n; char *t; n = strlen ( s ); t = s + strlen ( s ) - 1; while ( 0 < n ) { if ( *t != ' ' && *t != '\n' ) { return n; } t--; n--; } return n; } /******************************************************************************/ void timestamp ( ) /******************************************************************************/ /* Purpose: TIMESTAMP prints the current YMDHMS date as a time stamp. Example: 31 May 2001 09:45:54 AM Licensing: This code is distributed under the GNU LGPL license. Modified: 24 September 2003 Author: John Burkardt Parameters: None */ { # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct tm *tm; time_t now; now = time ( NULL ); tm = localtime ( &now ); strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm ); fprintf ( stdout, "%s\n", time_buffer ); return; # undef TIME_SIZE }