# include # include # include # include # include # include # include "r8poly.h" /******************************************************************************/ int i4_max ( int i1, int i2 ) /******************************************************************************/ /* Purpose: i4_max() returns the maximum of two I4's. Licensing: This code is distributed under the MIT license. Modified: 29 August 2006 Author: John Burkardt 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 MIT license. Modified: 29 August 2006 Author: John Burkardt 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 r8_sign ( double x ) /******************************************************************************/ /* Purpose: r8_sign() returns the sign of an R8. Licensing: This code is distributed under the MIT license. Modified: 08 May 2006 Author: John Burkardt Input: double X, the number whose sign is desired. Output: double R8_SIGN, the sign of X. */ { double value; if ( x < 0.0 ) { value = - 1.0; } else { value = + 1.0; } return value; } /******************************************************************************/ void r82poly2_print ( double a, double b, double c, double d, double e, double f ) /******************************************************************************/ /* Purpose: r82poly2_print() prints a second order polynomial in two variables. Licensing: This code is distributed under the MIT license. Modified: 21 May 2011 Author: John Burkardt Input: double A, B, C, D, E, F, the coefficients. */ { printf ( " %g * x^2 + %g * y^2 + %g * xy +\n", a, b, c ); printf ( " %g * x + %g * y + %g\n", d, e, f ); return; } /******************************************************************************/ int r82poly2_type ( double a, double b, double c, double d, double e, double f ) /******************************************************************************/ /* Purpose: r82poly2_type() analyzes a second order polynomial in two variables. Discussion: The polynomial has the form A x^2 + B y^2 + C xy + Dx + Ey + F = 0 The possible types of the solution set are: 1: a hyperbola; 9x^2 - 4y^2 -36x - 24y - 36 = 0 2: a parabola; 4x^2 + 1y^2 - 4xy + 3x - 4y + 1 = 0; 3: an ellipse; 9x^2 + 16y^2 +36x - 32y - 92 = 0; 4: an imaginary ellipse (no real solutions); x^2 + y^2 - 6x - 10y + 115 = 0; 5: a pair of intersecting lines; xy + 3x - y - 3 = 0 6: one point; x^2 + 2y^2 - 2x + 16y + 33 = 0; 7: a pair of distinct parallel lines; y^2 - 6y + 8 = 0 8: a pair of imaginary parallel lines (no real solutions); y^2 - 6y + 10 = 0 9: a pair of coincident lines. y^2 - 2y + 1 = 0 10: a single line; 2x - y + 1 = 0; 11; all space; 0 = 0; 12; no solutions; 1 = 0; Licensing: This code is distributed under the MIT license. Modified: 21 May 2011 Author: John Burkardt Reference: Daniel Zwillinger, editor, CRC Standard Mathematical Tables and Formulae, CRC Press, 30th Edition, 1996, pages 282-284. Input: double A, B, C, D, E, F, the coefficients. Output: int TYPE, indicates the type of the solution set. */ { double delta; double j; double k; int type; /* Handle the degenerate case. */ if ( a == 0.0 && b == 0.0 && c == 0.0 ) { if ( d == 0.0 && e == 0.0 ) { if ( f == 0.0 ) { type = 11; } else { type = 12; } } else { type = 10; } return type; } delta = 8.0 * a * b * f + 2.0 * c * e * d - 2.0 * a * e * e - 2.0 * b * d * d - 2.0 * f * c * c; j = 4.0 * a * b - c * c; if ( delta != 0.0 ) { if ( j < 0.0 ) { type = 1; } else if ( j == 0.0 ) { type = 2; } else if ( 0.0 < j ) { if ( r8_sign ( delta ) != r8_sign ( a + b ) ) { type = 3; } else if ( r8_sign ( delta ) == r8_sign ( a + b ) ) { type = 4; } } } else if ( delta == 0.0 ) { if ( j < 0.0 ) { type = 5; } else if ( 0.0 < j ) { type = 6; } else if ( j == 0.0 ) { k = 4.0 * ( a + b ) * f - d * d - e * e; if ( k < 0.0 ) { type = 7; } else if ( 0.0 < k ) { type = 8; } else if ( k == 0.0 ) { type = 9; } } } return type; } /******************************************************************************/ void r82poly2_type_print ( int type ) /******************************************************************************/ /* Purpose: r82poly2_type_print() prints the meaning of the output from R82POLY2_TYPE. Licensing: This code is distributed under the MIT license. Modified: 21 May 2011 Author: John Burkardt Input: int TYPE, the type index returned by R82POLY2_TYPE. */ { if ( type == 1 ) { printf ( " The set of solutions forms a hyperbola.\n" ); } else if ( type == 2 ) { printf ( " The set of solutions forms a parabola.\n" ); } else if ( type == 3 ) { printf ( " The set of solutions forms an ellipse.\n" ); } else if ( type == 4 ) { printf ( " The set of solutions forms an imaginary ellipse.\n" ); printf ( " (There are no real solutions).\n" ); } else if ( type == 5 ) { printf ( " The set of solutions forms a pair of intersecting lines.\n" ); } else if ( type == 6 ) { printf ( " The set of solutions is a single point.\n" ); } else if ( type == 7 ) { printf ( " The set of solutions form a pair of distinct parallel lines.\n" ); } else if ( type == 8 ) { printf ( " The set of solutions forms a pair of imaginary parallel lines.\n" ); printf ( " (There are no real solutions).\n" ); } else if ( type == 9 ) { printf ( " The set of solutions forms a pair of coincident lines.\n" ); } else if ( type == 10 ) { printf ( " The set of solutions forms a single line.\n" ); } else if ( type == 11 ) { printf ( " The set of solutions is all space.\n" ); } else if ( type == 12 ) { printf ( " The set of solutions is empty.\n" ); } else { printf ( " This type index is unknown.\n" ); } return; } /******************************************************************************/ double r8mat_det_3d ( double a[] ) /******************************************************************************/ /* Purpose: r8mat_det_3d() computes the determinant of a 3 by 3 R8MAT. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. The determinant of a 3 by 3 matrix is a11 * a22 * a33 - a11 * a23 * a32 + a12 * a23 * a31 - a12 * a21 * a33 + a13 * a21 * a32 - a13 * a22 * a31 Licensing: This code is distributed under the MIT license. Modified: 15 May 2010 Author: John Burkardt Input: double A[3*3], the matrix whose determinant is desired. Output: double R8MAT_DET_3D, the determinant of the matrix. */ { double det; det = a[0+0*3] * ( a[1+1*3] * a[2+2*3] - a[1+2*3] * a[2+1*3] ) + a[0+1*3] * ( a[1+2*3] * a[2+0*3] - a[1+0*3] * a[2+2*3] ) + a[0+2*3] * ( a[1+0*3] * a[2+1*3] - a[1+1*3] * a[2+0*3] ); return det; } /******************************************************************************/ double *r8mat_inverse_3d ( double a[] ) /******************************************************************************/ /* Purpose: r8mat_inverse_3d() inverts a 3 by 3 R8MAT using Cramer's rule. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. If the determinant is zero, A is singular, and does not have an inverse. In that case, the output is set to NULL. If the determinant is nonzero, its value is an estimate of how nonsingular the matrix A is. Licensing: This code is distributed under the MIT license. Modified: 14 May 2010 Author: John Burkardt Input: double A[3*3], the matrix to be inverted. Output: double R8MAT3_INVERSE[3*3], the inverse of the matrix A. */ { double *b; double det; /* Compute the determinant of A. */ det = a[0+0*3] * ( a[1+1*3] * a[2+2*3] - a[1+2*3] * a[2+1*3] ) + a[0+1*3] * ( a[1+2*3] * a[2+0*3] - a[1+0*3] * a[2+2*3] ) + a[0+2*3] * ( a[1+0*3] * a[2+1*3] - a[1+1*3] * a[2+0*3] ); if ( det == 0.0 ) { return NULL; } b = ( double * ) malloc ( 3 * 3 * sizeof ( double ) ); b[0+0*3] = ( a[1+1*3] * a[2+2*3] - a[1+2*3] * a[2+1*3] ) / det; b[0+1*3] = - ( a[0+1*3] * a[2+2*3] - a[0+2*3] * a[2+1*3] ) / det; b[0+2*3] = ( a[0+1*3] * a[1+2*3] - a[0+2*3] * a[1+1*3] ) / det; b[1+0*3] = - ( a[1+0*3] * a[2+2*3] - a[1+2*3] * a[2+0*3] ) / det; b[1+1*3] = ( a[0+0*3] * a[2+2*3] - a[0+2*3] * a[2+0*3] ) / det; b[1+2*3] = - ( a[0+0*3] * a[1+2*3] - a[0+2*3] * a[1+0*3] ) / det; b[2+0*3] = ( a[1+0*3] * a[2+1*3] - a[1+1*3] * a[2+0*3] ) / det; b[2+1*3] = - ( a[0+0*3] * a[2+1*3] - a[0+1*3] * a[2+0*3] ) / det; b[2+2*3] = ( a[0+0*3] * a[1+1*3] - a[0+1*3] * a[1+0*3] ) / det; return b; } /******************************************************************************/ void r8mat_print ( int m, int n, double a[], char *title ) /******************************************************************************/ /* Purpose: r8mat_print() prints an R8MAT. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. Entry A(I,J) is stored as A[I+J*M] Licensing: This code is distributed under the MIT license. Modified: 28 May 2008 Author: John Burkardt Input: int M, the number of rows in A. int N, the number of columns in A. double A[M*N], the M by N matrix. char *TITLE, a title. */ { r8mat_print_some ( m, n, a, 1, 1, m, n, title ); return; } /******************************************************************************/ void r8mat_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi, int jhi, char *title ) /******************************************************************************/ /* Purpose: r8mat_print_some() prints some of an R8MAT. 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 MIT license. Modified: 26 June 2013 Author: John Burkardt Input: int M, the number of rows of the matrix. M must be positive. int N, the number of columns of the matrix. N must be positive. double A[M*N], the matrix. int ILO, JLO, IHI, JHI, designate the first row and column, and the last row and column to be printed. char *TITLE, a title. */ { # define INCX 5 int i; int i2hi; int i2lo; 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; } /* Print the columns of the matrix, in strips of 5. */ for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX ) { j2hi = j2lo + INCX - 1; if ( n < j2hi ) { j2hi = n; } if ( jhi < j2hi ) { j2hi = jhi; } fprintf ( stdout, "\n" ); /* For each column J in the current range... Write the header. */ fprintf ( stdout, " Col: "); for ( j = j2lo; j <= j2hi; j++ ) { fprintf ( stdout, " %7d ", j - 1 ); } fprintf ( stdout, "\n" ); fprintf ( stdout, " Row\n" ); fprintf ( stdout, "\n" ); /* Determine the range of the rows in this strip. */ if ( 1 < ilo ) { i2lo = ilo; } else { i2lo = 1; } if ( m < ihi ) { i2hi = m; } else { i2hi = ihi; } for ( i = i2lo; i <= i2hi; i++ ) { /* Print out (up to) 5 entries in row I, that lie in the current strip. */ fprintf ( stdout, "%5d:", i - 1 ); for ( j = j2lo; j <= j2hi; j++ ) { fprintf ( stdout, " %14g", a[i-1+(j-1)*m] ); } fprintf ( stdout, "\n" ); } } return; # undef INCX } /******************************************************************************/ double *r8poly_add ( int na, double a[], int nb, double b[] ) /******************************************************************************/ /* Purpose: r8poly_add() adds two R8POLY's. Discussion: The polynomials are in power sum form. The power sum form is: p(x) = a(0) + a(1)*x + ... + a(n-1)*x^(n-1) + a(n)*x^(n) Licensing: This code is distributed under the MIT license. Modified: 21 November 2013 Author: John Burkardt Input: int NA, the degree of polynomial A. double A[NA+1], the coefficients of the first polynomial factor. int NB, the degree of polynomial B. double B[NB+1], the coefficients of the second polynomial factor. Output: double C[max(NA,NB)+1], the coefficients of A + B. */ { int i; double *c; c = ( double * ) malloc ( ( i4_max ( na, nb ) + 1 ) * sizeof ( double ) ); if ( nb == na ) { for ( i = 0; i <= na; i++ ) { c[i] = a[i] + b[i]; } } else if ( nb < na ) { for ( i = 0; i <= nb; i++ ) { c[i] = a[i] + b[i]; } for ( i = nb + 1; i <= na; i++ ) { c[i] = a[i]; } } else if ( na < nb ) { for ( i = 0; i <= na; i++ ) { c[i] = a[i] + b[i]; } for ( i = na + 1; i <= nb; i++ ) { c[i] = b[i]; } } return c; } /******************************************************************************/ double *r8poly_ant_coef ( int n, double poly_cof[] ) /******************************************************************************/ /* Purpose: r8poly_ant_coef() integrates a polynomial 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 MIT license. Modified: 15 October 2019 Author: John Burkardt Input: int N, the order of the polynomial. double POLY_COF[0:N], the polynomial coefficients. POLY_COF[0] is the constant term, and POLY_COF[N] is the coefficient of X^(N). Output: double POLY_COF2[0:N+1], the coefficients of the antiderivative polynomial, in standard form. The constant term is set to zero. */ { int i; double *poly_cof2; poly_cof2 = ( double * ) malloc ( ( n + 2 ) * sizeof ( double ) ); /* Set the constant term. */ poly_cof2[0] = 0.0; /* Integrate the polynomial. */ for ( i = 1; i < n + 2; i++ ) { poly_cof2[i] = poly_cof[i-1] / ( double ) ( i ); } return poly_cof2; } /******************************************************************************/ double r8poly_ant_value ( int n, double poly_cof[], double xval ) /******************************************************************************/ /* Purpose: r8poly_ant_value() evaluates the antiderivative of a polynomial. Discussion: The constant term of the antiderivative is taken to be zero. Licensing: This code is distributed under the MIT license. Modified: 16 October 2019 Author: John Burkardt Input: int N, the order of the polynomial. double POLY_COF[0:N], the polynomial coefficients. POLY_COF[0] is the constant term, and POLY_COF[N] is the coefficient of X^(N). real XVAL, the evaluation point. Output: real r8poly_ant_value, the value of the antiderivative at XVAL. */ { int i; double yval; yval = 0.0; for ( i = n; 0 <= i; i-- ) { yval = ( yval + poly_cof[i] / ( double ) ( i + 1 ) ) * xval; } return yval; } /******************************************************************************/ int r8poly_degree ( int na, double a[] ) /******************************************************************************/ /* Purpose: r8poly_degree() returns the degree of a polynomial. Discussion: The degree of a polynomial is the index of the highest power of X with a nonzero coefficient. The degree of a constant polynomial is 0. The degree of the zero polynomial is debatable, but this routine returns the degree as 0. Licensing: This code is distributed under the MIT license. Modified: 15 May 2010 Author: John Burkardt Input: int NA, the dimension of A. double A[NA+1], the coefficients of the polynomials. Output: int R8POLY_DEGREE, the degree of A. */ { int degree; degree = na; while ( 0 < degree ) { if ( a[degree] != 0.0 ) { return degree; } degree = degree - 1; } return degree; } /******************************************************************************/ double *r8poly_deriv ( int n, double c[], int p ) /******************************************************************************/ /* Purpose: r8poly_deriv() returns the derivative of a polynomial. Licensing: This code is distributed under the MIT license. Modified: 16 May 2010 Author: John Burkardt Input: int N, the degree of the polynomial. double C[N+1], the polynomial coefficients. C[I] is the coefficient of X^I. int P, the order of the derivative. 0 means no derivative is taken. 1 means first derivative, 2 means second derivative and so on. Values of P less than 0 are meaningless. Values of P greater than N are meaningful, but the code will behave as though the value of P was N+1. Output: double R8POLY_DERIV CP[N-P+1], the polynomial coefficients of the derivative. */ { double *cp; double *cp_temp; int d; int i; if ( n < p ) { return NULL; } cp_temp = r8vec_copy_new ( n + 1, c ); for ( d = 1; d <= p; d++ ) { for ( i = 0; i <= n-d; i++ ) { cp_temp[i] = ( double ) ( i + 1 ) * cp_temp[i+1]; } cp_temp[n-d+1] = 0.0; } cp = r8vec_copy_new ( n - p + 1, cp_temp ); free ( cp_temp ); return cp; } /******************************************************************************/ void r8poly_division ( int na, double a[], int nb, double b[], int *nq, double q[], int *nr, double r[] ) /******************************************************************************/ /* Purpose: r8poly_division() computes the quotient and remainder of two R8POLY's. Discussion: The polynomials are assumed to be stored in power sum form. The power sum form is: p(x) = a(0) + a(1)*x + ... + a(n-1)*x^(n-1) + a(n)*x^(n) Licensing: This code is distributed under the MIT license. Modified: 08 October 2019 Author: John Burkardt Input: int NA, the dimension of A. double A[NA+1], the coefficients of the polynomial to be divided. int NB, the dimension of B. double B[NB+1], the coefficients of the divisor polynomial. Output: int *NQ, the degree of Q. If the divisor polynomial is zero, NQ is returned as -1. double Q[NA+1-NB], contains the quotient of A/B. If A and B have full degree, Q should be dimensioned Q(0:NA-NB). In any case, Q(0:NA) should be enough. int *NR, the degree of R. If the divisor polynomial is zero, NR is returned as -1. double R[NB], contains the remainder of A/B. If B has full degree, R should be dimensioned R(0:NB-1). Otherwise, R will actually require less space. */ { double *a2; int i; int j; int na2; int nb2; na2 = r8poly_degree ( na, a ); nb2 = r8poly_degree ( nb, b ); if ( b[nb2] == 0.0 ) { *nq = -1; *nr = -1; return; } a2 = ( double * ) malloc ( ( na + 1 ) * sizeof ( double ) ); for ( i = 0; i <= na; i++ ) { a2[i] = a[i]; } *nq = na2 - nb2; *nr = nb2 - 1; for ( i = *nq; 0 <= i; i-- ) { q[i] = a2[i+nb2] / b[nb2]; a2[i+nb2] = 0.0; for ( j = 0; j <= nb2-1; j++ ) { a2[i+j] = a2[i+j] - q[i] * b[j]; } } for ( i = 0; i <= *nr; i++ ) { r[i] = a2[i]; } free ( a2 ); return; } /******************************************************************************/ double r8poly_lagrange_0 ( int npol, double xpol[], double xval ) /******************************************************************************/ /* Purpose: r8poly_lagrange_0() evaluates the Lagrange factor at a point. Discussion: W(X) = Product ( 1 <= I <= NPOL ) ( X - XPOL(I) ) Discussion: For a set of points XPOL(I), 1 <= I <= NPOL, the IPOL-th Lagrange basis polynomial L(IPOL)(X), has the property: L(IPOL)( XPOL(J) ) = delta ( IPOL, J ) and may be expressed as: L(IPOL)(X) = W(X) / ( ( X - XPOL(IPOL) ) * W'(XPOL(IPOL)) ) Licensing: This code is distributed under the MIT license. Modified: 25 August 2011 Author: John Burkardt Input: int NPOL, the number of abscissas. NPOL must be at least 1. double XPOL[NPOL], the abscissas, which should be distinct. double XVAL, the point at which the Lagrange factor is to be evaluated. Output: double R8POLY_LAGRANGE_0, the value of the Lagrange factor at XVAL. */ { int i; double wval; wval = 1.0; for ( i = 0; i < npol; i++ ) { wval = wval * ( xval - xpol[i] ); } return wval; } /******************************************************************************/ double r8poly_lagrange_1 ( int npol, double xpol[], double xval ) /******************************************************************************/ /* Purpose: r8poly_lagrange_1() evaluates the first derivative of the Lagrange factor. Discussion: W(XPOL(1:NPOL))(X) = Product ( 1 <= I <= NPOL ) ( X - XPOL(I) ) W'(XPOL(1:NPOL))(X) = Sum ( 1 <= J <= NPOL ) Product ( I /= J ) ( X - XPOL(I) ) We also have the recursion: W'(XPOL(1:NPOL))(X) = d/dX ( ( X - XPOL(NPOL) ) * W(XPOL(1:NPOL-1))(X) ) = W(XPOL(1:NPOL-1))(X) + ( X - XPOL(NPOL) ) * W'(XPOL(1:NPOL-1))(X) Licensing: This code is distributed under the MIT license. Modified: 25 August 2011 Author: John Burkardt Input: int NPOL, the number of abscissas. double XPOL[NPOL], the abscissas, which should be distinct. double XVAL, the point at which the Lagrange factor is to be evaluated. Output: double R8POLY_LAGRANGE_1, the derivative of W with respect to XVAL. */ { double dwdx; int i; double w; dwdx = 0.0; w = 1.0; for ( i = 0; i < npol; i++ ) { dwdx = w + ( xval - xpol[i] ) * dwdx; w = w * ( xval - xpol[i] ); } return dwdx; } /******************************************************************************/ double r8poly_lagrange_2 ( int npol, double xpol[], double xval ) /******************************************************************************/ /* Purpose: r8poly_lagrange_2() evaluates the second derivative of the Lagrange factor. Discussion: W(X) = Product ( 1 <= I <= NPOL ) ( X - XPOL(I) ) W'(X) = Sum ( 1 <= J <= NPOL ) Product ( I /= J ) ( X - XPOL(I) ) W"(X) = Sum ( 1 <= K <= NPOL ) Sum ( J =/ K ) Product ( I /= K, J ) ( X - XPOL(I) ) For a set of points XPOL(I), 1 <= I <= NPOL, the IPOL-th Lagrange basis polynomial L(IPOL)(X), has the property: L(IPOL)( XPOL(J) ) = delta ( IPOL, J ) and may be expressed as: L(IPOL)(X) = W(X) / ( ( X - XPOL(IPOL) ) * W'(XPOL(IPOL)) ) Licensing: This code is distributed under the MIT license. Modified: 24 January 2004 Author: John Burkardt Input: int NPOL, the number of abscissas. NPOL must be at least 1. double XPOL[NPOL], the abscissas, which should be distinct. double XVAL, the point at which the Lagrange factor is to be evaluated. Output: double R8POLY_LAGRANGE_2, the second derivative of W with respect to XVAL. */ { double dw2dx2; int i; int j; int k; double term; dw2dx2 = 0.0; for ( k = 0; k < npol; k++ ) { for ( j = 0; j < npol; j++ ) { if ( j != k ) { term = 1.0; for ( i = 0; i < npol; i++ ) { if ( i != j && i != k ) { term = term * ( xval - xpol[i] ); } } dw2dx2 = dw2dx2 + term; } } } return dw2dx2; } /******************************************************************************/ double *r8poly_lagrange_coef ( int npol, int ipol, double xpol[] ) /******************************************************************************/ /* Purpose: r8poly_lagrange_coef() returns the coefficients of a Lagrange polynomial. Discussion: Given NPOL distinct abscissas, XPOL(*), the IPOL-th Lagrange polynomial P(IPOL)(X) is defined as the polynomial of degree NPOL - 1 which is 1 at XPOL(IPOL) and 0 at the NPOL - 1 other abscissas. Licensing: This code is distributed under the MIT license. Modified: 25 August 2011 Author: John Burkardt Input: int NPOL, the number of abscissas. NPOL must be at least 1. int IPOL, the index of the polynomial to evaluate. IPOL must be between 1 and NPOL. double XPOL[NPOL], the abscissas of the Lagrange polynomials. The entries in XPOL must be distinct. Output: double R8POLY_LAGRANGE_COEF[NPOL], the polynomial coefficients of the IPOL-th Lagrange polynomial. */ { int i; int index; int j; double *pcof; /* Make sure IPOL is legal. */ if ( ipol < 1 || npol < ipol ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8POLY_LAGRANGE_COEF(): Fatal error!\n" ); fprintf ( stderr, " 1 <= IPOL <= NPOL is required.\n" ); fprintf ( stderr, " but IPOL = %d\n", ipol ); fprintf ( stderr, " and NPOL = %d\n", npol ); exit ( 1 ); } /* Check that the abscissas are distinct. */ if ( ! r8vec_is_distinct ( npol, xpol ) ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8POLY_LAGRANGE_COEF(): Fatal error!\n" ); fprintf ( stderr, " Two entries of XPOL are equal:\n" ); exit ( 1 ); } pcof = ( double * ) malloc ( npol * sizeof ( double ) ); pcof[0] = 1.0; for ( i = 1; i < npol; i++ ) { pcof[i] = 0.0; } index = 0; for ( i = 1; i <= npol; i++ ) { if ( i != ipol ) { index = index + 1; for ( j = index; 0 <= j; j-- ) { pcof[j] = - xpol[i-1] * pcof[j] / ( xpol[ipol-1] - xpol[i-1] ); if ( 0 < j ) { pcof[j] = pcof[j] + pcof[j-1] / ( xpol[ipol-1] - xpol[i-1] ); } } } } return pcof; } /******************************************************************************/ void r8poly_lagrange_factor ( int npol, double xpol[], double xval, double *wval, double *dwdx ) /******************************************************************************/ /* Purpose: r8poly_lagrange_factor() evaluates the polynomial Lagrange factor at a point. Discussion: Suppose F(X) is at least N times continuously differentiable in the interval [A,B]. Pick NPOL distinct points XPOL(I) in [A,B] and compute the interpolating polynomial P(X) of order NPOL ( and degree NPOL-1) which passes through all the points ( XPOL(I), F(XPOL(I)) ). Then in the interval [A,B], the maximum error abs ( F(X) - P(X) ) is bounded by: C * FNMAX * W(X) where C is a constant, FNMAX is the maximum value of the NPOL-th derivative of F in [A,B], W(X) is the Lagrange factor. Thus, the value of W(X) is useful as part of an estimated bound for the interpolation error. The formula is: W(X) = Product ( 1 <= I <= NPOL ) ( X - XPOL(I) ) Note that the Chebyshev abscissas have the property that they minimize the value of W(X) over the interval [A,B]. Hence, if the abscissas may be chosen arbitrarily, the Chebyshev abscissas have this advantage over other choices. For a set of points XPOL[I], 0 <= I <= NPOL-1, the IPOL-th Lagrange basis polynomial L(IPOL)(X), has the property: L(IPOL)( XPOL(J) ) = delta ( IPOL, J ) and may be expressed as: L(IPOL)(X) = W(X) / ( ( X - XPOL[IPOL] ) * W'(XPOL[IPOL]) ) Licensing: This code is distributed under the MIT license. Modified: 25 August 2011 Author: John Burkardt Input: int NPOL, the number of abscissas. NPOL must be at least 1. double XPOL[NPOL], the abscissas, which should be distinct. double XVAL, the point at which the Lagrange factor is to be evaluated. Output: double *WVAL, the value of the Lagrange factor at XVAL. double *DWDX, the derivative of W with respect to XVAL. */ { int i; int j; double term; *wval = 1.0; for ( i = 0; i < npol; i++ ) { *wval = *wval * ( xval - xpol[i] ); } *dwdx = 0.0; for ( i = 0; i < npol; i++ ) { term = 1.0; for ( j = 0; j < npol; j++ ) { if ( i != j ) { term = term * ( xval - xpol[j] ); } } *dwdx = *dwdx + term; } return; } /******************************************************************************/ int r8poly_lagrange_value ( int npol, int ipol, double xpol[], double xval, double *pval, double *dpdx ) /******************************************************************************/ /* Purpose: r8poly_lagrange_value() evaluates the IPOL-th Lagrange polynomial. Discussion: Given NPOL distinct abscissas, XPOL[*], the IPOL-th Lagrange polynomial P(IPOL)(X) is defined as the polynomial of degree NPOL - 1 which is 1 at XPOL[IPOL] and 0 at the NPOL - 1 other abscissas. Licensing: This code is distributed under the MIT license. Modified: 25 August 2011 Author: John Burkardt Input: int NPOL, the number of abscissas. NPOL must be at least 1. int IPOL, the index of the polynomial to evaluate. IPOL must be between 0 and NPOL-1. double XPOL[NPOL], the abscissas of the Lagrange polynomials. The entries in XPOL must be distinct. double XVAL, the point at which the IPOL-th Lagrange polynomial is to be evaluated. Output: double *PVAL, the value of the IPOL-th Lagrange polynomial at XVAL. double *DPDX, the derivative of the IPOL-th Lagrange polynomial at XVAL. int r8poly_lagrange_value, 0 if no error. */ { int i; int j; double p2; /* Make sure IPOL is legal. */ if ( ipol < 0 || npol-1 < ipol ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "r8poly_lagrange_value(): Fatal error!\n" ); fprintf ( stderr, " 0 <= IPOL <= NPOL-1 is required.\n" ); exit ( 1 ); } /* Check that the abscissas are distinct. */ for ( i = 1; i < npol; i++ ) { for ( j = 0; j < i; j++ ) { if ( xpol[i] == xpol[j] ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "r8poly_lagrange_value(): Fatal error!\n" ); fprintf ( stderr, " Two entries of XPOL are equal:\n" ); fprintf ( stderr, " XPOL(%d) = %g.\n", i, xpol[i] ); fprintf ( stderr, " XPOL(%d) = %g.\n", j, xpol[j] ); exit ( 1 ); } } } /* Evaluate the polynomial. */ *pval = 1.0; for ( i = 0; i < npol; i++ ) { if ( i != ipol ) { *pval = *pval * ( xval - xpol[i] ) / ( xpol[ipol] - xpol[i] ); } } /* Evaluate the derivative, which can be found by summing up the result of differentiating one factor at a time, successively. */ *dpdx = 0.0; for ( i = 0; i < npol; i++ ) { if ( i != ipol ) { p2 = 1.0; for ( j = 0; j < npol; j++ ) { if ( j == i ) { p2 = p2 / ( xpol[ipol] - xpol[j] ); } else if ( j != ipol ) { p2 = p2 * ( xval - xpol[j] ) / ( xpol[ipol] - xpol[j] ); } } *dpdx = *dpdx + p2; } } return 0; } /******************************************************************************/ void r8poly_multiply ( int na, double a[], int nb, double b[], double c[] ) /******************************************************************************/ /* Purpose: r8poly_multiply() computes the product of two polynomials. Discussion: The polynomials are in power sum form. The power sum form is: p(x) = a(0) + a(1)*x + ... + a(n-1)*x^(n-1) + a(n)*x^(n) Licensing: This code is distributed under the MIT license. Modified: 28 March 2009 Author: John Burkardt Input: int NA, the dimension of A. double A[NA+1], the coefficients of the first polynomial factor. int NB, the dimension of B. double B[NB+1], the coefficients of the second polynomial factor. Output: double C[NA+NB+1], the coefficients of A * B. */ { double *d; int i; int j; d = ( double * ) malloc ( ( na + nb + 1 ) * sizeof ( double ) ); for ( i = 0; i <= na + nb; i++ ) { d[i] = 0.0; } for ( i = 0; i <= na; i++ ) { for ( j = 0; j <= nb; j++ ) { d[i+j] = d[i+j] + a[i] * b[j]; } } for ( i = 0; i <= na + nb; i++ ) { c[i] = d[i]; } free ( d ); return; } /******************************************************************************/ void r8poly_power ( int na, double a[], int p, double b[] ) /******************************************************************************/ /* Purpose: r8poly_power() computes a positive integer power of a polynomial. Discussion: The power sum form is: p(x) = a(0) + a(1)*x + ... + a(n-1)*x^(n-1) + a(n)*x^(n) Licensing: This code is distributed under the MIT license. Modified: 13 June 2015 Author: John Burkardt Input: int NA, the dimension of A. double A[NA+1], the polynomial to be raised to the power. int P, the nonnegative power to which A is raised. Output: double B[P*NA+1], the power of the polynomial. */ { int i; int j; int nonzer; /* Zero out B. */ for ( i = 0; i <= p * na; i++ ) { b[i] = 0.0; } /* Search for the first nonzero element in A. */ nonzer = -1; for ( i = 0; i <= na; i++ ) { if ( a[i] != 0.0 ) { nonzer = i; break; } } if ( nonzer == -1 ) { return; } b[0] = pow ( a[nonzer], p ); for ( i = 1; i <= p * ( na - nonzer ); i++ ) { if ( i + nonzer <= na ) { b[i] = ( double ) ( i * p ) * b[0] * a[i+nonzer]; } else { b[i] = 0.0; } for ( j = 1; j <= i-1; j++ ) { if ( j+nonzer <= na ) { b[i] = b[i] - ( double ) ( i - j ) * a[j+nonzer] * b[i-j]; } if ( i-j+nonzer <= na ) { b[i] = b[i] + ( double ) ( i - j ) * ( double ) p * b[j] * a[i-j+nonzer]; } } b[i] = b[i] / ( ( double ) i * a[nonzer] ); } /* Shift B up. */ for ( i = p*nonzer; i <= p*na; i++ ) { b[i] = b[i-p*nonzer]; } for ( i = 0; i <= p * nonzer-1; i++ ) { b[i] = 0.0; } return; } /******************************************************************************/ void r8poly_print ( int n, double a[], char *title ) /******************************************************************************/ /* Purpose: r8poly_print() prints a polynomial. Licensing: This code is distributed under the MIT license. Modified: 15 July 2015 Author: John Burkardt Input: int N, the dimension of A. double A[N+1], the polynomial coefficients. A(0) is the constant term and A(N) is the coefficient of X**N. char *TITLE, a title. */ { int i; double mag; char plus_minus; if ( 0 < s_len_trim ( title ) ) { printf ( "\n" ); printf ( "%s\n", title ); } printf ( "\n" ); if ( n < 0 ) { printf ( " p(x) = 0\n" ); return; } if ( a[n] < 0.0 ) { plus_minus = '-'; } else { plus_minus = ' '; } mag = fabs ( a[n] ); if ( 2 <= n ) { printf ( " p(x) = %c%g * x^%d\n", plus_minus, mag, n ); } else if ( n == 1 ) { printf ( " p(x) = %c%g * x\n", plus_minus, mag ); } else if ( n == 0 ) { printf ( " p(x) = %c%g\n", plus_minus, mag ); } for ( i = n - 1; 0 <= i; i-- ) { if ( a[i] < 0.0 ) { plus_minus = '-'; } else { plus_minus = '+'; } mag = fabs ( a[i] ); if ( mag != 0.0 ) { if ( 2 <= i ) { printf ( " %c%g * x^%d\n", plus_minus, mag, i ); } else if ( i == 1 ) { printf ( " %c%g * x\n", plus_minus, mag ); } else if ( i == 0 ) { printf ( " %c%g\n", plus_minus, mag ); } } } 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) * X^N + ... + 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) * Z^N + ... + 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 MIT license. Modified: 25 August 2011 Reference: Press, Flannery, Teukolsky, Vetterling, Numerical Recipes: The Art of Scientific Computing, Cambridge University Press. Input: double SHIFT, SCALE, the shift and scale applied to X, so that Z = SCALE * X + SHIFT. int N, the number of coefficients. double POLY_COF[N+1]. the coefficient array in terms of the X variable. Output: double POLY_COF[N+1]: the coefficient array in terms of the Z variable. */ { int i; int j; for ( i = 1; i <= n; i++ ) { for ( j = i; j <= n; j++ ) { poly_cof[j] = poly_cof[j] / scale; } } for ( i = 0; i <= n - 1; i++ ) { for ( j = n - 1; i <= j; j-- ) { poly_cof[j] = poly_cof[j] - shift * poly_cof[j+1]; } } return; } /******************************************************************************/ double r8poly_value ( int m, double c[], double x ) /******************************************************************************/ /* Purpose: r8poly_value() evaluates a polynomial using a naive method. Discussion: The polynomial p(x) = c0 + c1 * x + c2 * x^2 + ... + cm * x^m is to be evaluated at the value X. Licensing: This code is distributed under the MIT license. Modified: 22 August 2017 Author: John Burkardt Input: int M, the degree of the polynomial. double C[0:M], the coefficients of the polynomial. C[0] is the constant term. double X, the point at which the polynomial is to be evaluated. Output: double R8POLY_VALUE, the value of the polynomial at X. */ { int i; double value; double xi; value = c[0]; xi = 1.0; for ( i = 1; i <= m; i++ ) { xi = xi * x; value = value + c[i] * xi; } return value; } /******************************************************************************/ double r8poly_value_horner ( int m, double c[], double x ) /******************************************************************************/ /* Purpose: r8poly_value_horner() evaluates a polynomial using Horner's method. Discussion: The polynomial p(x) = c0 + c1 * x + c2 * x^2 + ... + cm * x^m is to be evaluated at the value X. Licensing: This code is distributed under the MIT license. Modified: 02 January 2015 Author: John Burkardt Input: int M, the degree of the polynomial. double C[M+1], the coefficients of the polynomial. C[0] is the constant term. Output: double X, the point at which the polynomial is to be evaluated. double R8POLY_VALUE_HORNER, the value of the polynomial at X. */ { int i; double value; value = c[m]; for ( i = m - 1; 0 <= i; i-- ) { value = value * x + c[i]; } return value; } /******************************************************************************/ double *r8poly_values_horner ( int m, double c[], int n, double x[] ) /******************************************************************************/ /* Purpose: r8poly_values_horner() evaluates a polynomial using Horner's method. Discussion: The polynomial p(x) = c0 + c1 * x + c2 * x^2 + ... + cm * x^m is to be evaluated at the vector of values X. Licensing: This code is distributed under the MIT license. Modified: 03 December 2013 Author: John Burkardt Input: int M, the degree. double C[M+1], the polynomial coefficients. C[I] is the coefficient of X^I. int N, the number of evaluation points. double X[N], the evaluation points. Output: double R8POLY_VALUES_HORNER[N], the polynomial values. */ { int i; int j; double *p; p = ( double * ) malloc ( n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { p[j] = c[m]; } for ( i = m - 1; 0 <= i; i-- ) { for ( j = 0; j < n; j++ ) { p[j] = p[j] * x[j] + c[i]; } } return p; } /******************************************************************************/ double *r8poly_value_2d ( int m, double c[], int n, double x[], double y[] ) /******************************************************************************/ /* Purpose: r8poly_value_2d() evaluates a polynomial in 2 variables, X and Y. Discussion: We assume the polynomial is of total degree M, and has the form: p(x,y) = c00 + c10 * x + c01 * y + c20 * x^2 + c11 * xy + c02 * y^2 + ... + cm0 * x^(m) + ... + c0m * y^m. Licensing: This code is distributed under the MIT license. Modified: 23 September 2012 Author: John Burkardt Input: int M, the degree of the polynomial. double C[T(M+1)], the polynomial coefficients. C[0] is the constant term. T(M+1) is the M+1-th triangular number. The coefficients are stored consistent with the following ordering of monomials: 1, X, Y, X^2, XY, Y^2, X^3, X^2Y, XY^2, Y^3, X^4, ... int N, the number of evaluation points. double X[N], Y[N], the evaluation points. Output: double R8POLY_VALUE_2D[N], the value of the polynomial at the evaluation points. */ { int ex; int ey; int i; int j; double *p; int s; p = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { p[i] = 0.0; } j = 0; for ( s = 0; s <= m; s++ ) { for ( ex = s; 0 <= ex; ex-- ) { ey = s - ex; for ( i = 0; i < n; i++ ) { p[i] = p[i] + c[j] * pow ( x[i], ex ) * pow ( y[i], ey ); } j = j + 1; } } return p; } /******************************************************************************/ double r8poly2_discriminant ( double a, double b, double c ) /******************************************************************************/ /* Purpose: r8poly2_discriminant() returns the discriminant of a quadratic polynomial. Discussion: The polynomial has the form: A * X^2 + B * X + C = 0 If the discriminant is zero, then the polynomial has a multiple root. Licensing: This code is distributed under the MIT license. Modified: 14 July 2025 Author: John Burkardt Input: real A, B, C, the coefficients of the polynomial. Output: real VALUE, the discriminant. */ { double value; value = b * b - 4.0 * a * c; return value; } /******************************************************************************/ int r8poly2_ex ( double x1, double y1, double x2, double y2, double x3, double y3, double *x, double *y ) /******************************************************************************/ /* Purpose: r8poly2_ex() finds the extremal point of a parabola determined by three points. Licensing: This code is distributed under the MIT license. Modified: 27 August 2011 Author: John Burkardt Input: double X1, Y1, X2, Y2, X3, Y3, the coordinates of three points on the parabola. X1, X2 and X3 must be distinct. Output: double *X, *Y, the X coordinate of the extremal point of the parabola, and the value of the parabola at that point. int R8POLY2_EX, error flag. 0, no error. 1, two of the X values are equal. 2, the data lies on a straight line; there is no finite extremal point. 3, the data lies on a horizontal line; every point is "extremal". */ { double bot; *x = 0.0; *y = 0.0; if ( x1 == x2 || x2 == x3 || x3 == x1 ) { return 1; } if ( y1 == y2 && y2 == y3 && y3 == y1 ) { *x = x1; *y = y1; return 3; } bot = ( x2 - x3 ) * y1 + ( x3 - x1 ) * y2 + ( x1 - x2 ) * y3; if ( bot == 0.0 ) { return 2; } *x = 0.5 * ( x1 * x1 * ( y3 - y2 ) + x2 * x2 * ( y1 - y3 ) + x3 * x3 * ( y2 - y1 ) ) / ( ( x2 - x3 ) * y1 + ( x3 - x1 ) * y2 + ( x1 - x2 ) * y3 ); *y = - ( ( *x - x2 ) * ( *x - x3 ) * ( x2 - x3 ) * y1 + ( *x - x1 ) * ( *x - x3 ) * ( x3 - x1 ) * y2 + ( *x - x1 ) * ( *x - x2 ) * ( x1 - x2 ) * y3 ) / ( ( x1 - x2 ) * ( x2 - x3 ) * ( x3 - x1 ) ); return 0; } /******************************************************************************/ int r8poly2_ex2 ( double x1, double y1, double x2, double y2, double x3, double y3, double *x, double *y, double *a, double *b, double *c ) /******************************************************************************/ /* Purpose: r8poly2_ex2() finds the extremal point of a parabola determined by three points. Licensing: This code is distributed under the MIT license. Modified: 27 August 2011 Author: John Burkardt Input: double X1, Y1, X2, Y2, X3, Y3, the coordinates of three points on the parabola. X1, X2 and X3 must be distinct. Output: double *X, *Y, the X coordinate of the extremal point of the parabola, and the value of the parabola at that point. double *A, *B, *C, the coefficients that define the parabola: P(X) = A * X**2 + B * X + C. int R8POLY2_EX2, error flag. 0, no error. 1, two of the X values are equal. 2, the data lies on a straight line; there is no finite extremal point. 3, the data lies on a horizontal line; any point is an "extremal point". */ { double v[3*3]; double *w; *a = 0.0; *b = 0.0; *c = 0.0; *x = 0.0; *y = 0.0; if ( x1 == x2 || x2 == x3 || x3 == x1 ) { return 1; } if ( y1 == y2 && y2 == y3 && y3 == y1 ) { *x = x1; *y = y1; return 3; } /* Set up the Vandermonde matrix. */ v[0+0*3] = 1.0; v[0+1*3] = x1; v[0+2*3] = x1 * x1; v[1+0*3] = 1.0; v[1+1*3] = x2; v[1+2*3] = x2 * x2; v[2+0*3] = 1.0; v[2+1*3] = x3; v[2+2*3] = x3 * x3; /* Get the inverse. */ w = r8mat_inverse_3d ( v ); /* Compute the parabolic coefficients. */ *c = w[0+0*3] * y1 + w[0+1*3] * y2 + w[0+2*3] * y3; *b = w[1+0*3] * y1 + w[1+1*3] * y2 + w[1+2*3] * y3; *a = w[2+0*3] * y1 + w[2+1*3] * y2 + w[2+2*3] * y3; /* Determine the extremal point. */ if ( *a == 0.0 ) { return 2; } *x = - *b / ( 2.0 * *a ); *y = *a * *x * *x + *b * *x + *c; return 0; } /******************************************************************************/ void r8poly2_root ( double a, double b, double c, double complex *r1, double complex *r2 ) /******************************************************************************/ /* Purpose: r8poly2_root() returns the two roots of a quadratic polynomial. Discussion: The polynomial has the form: A * X * X + B * X + C = 0 Licensing: This code is distributed under the MIT license. Modified: 08 August 2018 Author: John Burkardt Input: double A, B, C, the coefficients of the polynomial. A must not be zero. Output: double complex *R1, *R2, the roots of the polynomial, which might be real and distinct, real and equal, or complex conjugates. */ { double complex disc; double complex q; if ( a == 0.0 ) { printf ( "\n" ); printf ( "R8POLY2_ROOT(): Fatal error!\n" ); printf ( " The coefficient A is zero.\n" ); exit ( 1 ); } disc = b * b - 4.0 * a * c; q = - 0.5 * ( b + r8_sign ( b ) * csqrt ( disc ) ); *r1 = q / a; *r2 = c / q; return; } /******************************************************************************/ void r8poly2_rroot ( double a, double b, double c, double *r1, double *r2 ) /******************************************************************************/ /* Purpose: r8poly2_rroot() returns the real parts of the roots of a quadratic polynomial. Example: A B C roots R1 R2 -- -- -- ------------------ -- -- 1 -4 3 1 3 1 3 1 0 4 2*i - 2*i 0 0 1 -6 10 3 + i 3 - i 3 3 Licensing: This code is distributed under the MIT license. Modified: 01 December 2016 Author: John Burkardt Input: double A, B, C, the coefficients of the quadratic polynomial A * X^2 + B * X + C = 0 whose roots are desired. A must not be zero. Output: double *R1, *R2, the real parts of the roots of the polynomial. */ { double disc; double q; if ( a == 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8POLY2_RROOT(): Fatal error!\n" ); fprintf ( stderr, " The coefficient A is zero.\n" ); exit ( 1 ); } disc = b * b - 4.0 * a * c; if ( 0.0 <= disc ) { q = ( b + r8_sign ( b ) * sqrt ( disc ) ); *r1 = -0.5 * q / a; *r2 = -2.0 * c / q; } else { *r1 = b / 2.0 / a; *r2 = b / 2.0 / a; } return; } /******************************************************************************/ void r8poly2_val ( double x1, double y1, double x2, double y2, double x3, double y3, double x, double *y, double *yp, double *ypp ) /******************************************************************************/ /* Purpose: r8poly2_val() evaluates a parabola defined by three data values. Licensing: This code is distributed under the MIT license. Modified: 27 August 2011 Author: John Burkardt Input: double X1, Y1, X2, Y2, X3, Y3, three pairs of data values. If the X values are distinct, then all the Y values represent actual values of the parabola. Three special cases are allowed: X1 = X2 =/= X3: Y2 is the derivative at X1; X1 =/= X2 = X3: Y3 is the derivative at X3; X1 = X2 = X3: Y2 is the derivative at X1, and Y3 is the second derivative at X1. double X, an abscissa at which the parabola is to be evaluated. Output: double *Y, *YP, *YPP, the values of the parabola and its first and second derivatives at X. */ { int distinct; double dif1; double dif2; double temp; /* If any X's are equal, put them and the Y data first. */ if ( x1 == x2 && x2 == x3 ) { distinct = 1; } else if ( x1 == x2 ) { distinct = 2; } else if ( x1 == x3 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8POLY2_VAL(): Fatal error!\n" ); fprintf ( stderr, " X1 = X3 =/= X2.\n" ); return; } else if ( x2 == x3 ) { distinct = 2; temp = x1; x1 = x3; x3 = temp; temp = y1; y1 = y2; y2 = y3; y3 = y1; } else { distinct = 3; } /* Set up the coefficients. */ if ( distinct == 1 ) { dif1 = y2; dif2 = 0.5 * y3; } else if ( distinct == 2 ) { dif1 = y2; dif2 = ( ( y3 - y1 ) / ( x3 - x1 ) - y2 ) / ( x3 - x2 ); } else if ( distinct == 3 ) { dif1 = ( y2 - y1 ) / ( x2 - x1 ); dif2 = ( ( y3 - y1 ) / ( x3 - x1 ) - ( y2 - y1 ) / ( x2 - x1 ) ) / ( x3 - x2 ); } /* Evaluate. */ *y = y1 + ( x - x1 ) * dif1 + ( x - x1 ) * ( x - x2 ) * dif2; *yp = dif1 + ( 2.0 * x - x1 - x2 ) * dif2; *ypp = 2.0 * dif2; return; } /******************************************************************************/ void r8poly2_val2 ( int ndata, double tdata[], double ydata[], int left, double tval, double *yval ) /******************************************************************************/ /* Purpose: r8poly2_val2() evaluates a parabolic function through 3 points in a table. Discussion: This routine is a utility routine used by OVERHAUSER_SPLINE_VAL. It constructs the parabolic interpolant through the data in 3 consecutive entries of a table and evaluates this interpolant at a given abscissa value. Licensing: This code is distributed under the MIT license. Modified: 27 August 2011 Author: John Burkardt Input: int NDATA, the number of data points. NDATA must be at least 3. double TDATA[NDATA], the abscissas of the data points. The values in TDATA must be in strictly ascending order. double YDATA[NDATA], the data points corresponding to the abscissas. int LEFT, the location of the first of the three consecutive data points through which the parabolic interpolant must pass. 0 <= LEFT <= NDATA - 3. double TVAL, the value of T at which the parabolic interpolant is to be evaluated. Normally, TDATA[0] <= TVAL <= T[NDATA-1], and the data will be interpolated. For TVAL outside this range, extrapolation will be used. Output: double *YVAL, the value of the parabolic interpolant at TVAL. */ { double dif1; double dif2; double t1; double t2; double t3; double y1; double y2; double y3; /* Check. */ if ( left < 0 || ndata-3 < left ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "RPOLY2_VAL2(): Fatal error!\n" ); fprintf ( stderr, " LEFT < 0 or NDATA-3 < LEFT.\n" ); exit ( 1 ); } /* Copy out the three abscissas. */ t1 = tdata[left]; t2 = tdata[left+1]; t3 = tdata[left+2]; if ( t2 <= t1 || t3 <= t2 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "RPOLY2_VAL2(): Fatal error!\n" ); fprintf ( stderr, " T2 <= T1 or T3 <= T2.\n" ); fprintf ( stderr, " T1 = %g\n", t1 ); fprintf ( stderr, " T2 = %g\n", t2 ); fprintf ( stderr, " T3 = %g\n", t3 ); exit ( 1 ); } /* Construct and evaluate a parabolic interpolant for the data. */ y1 = ydata[left]; y2 = ydata[left+1]; y3 = ydata[left+2]; dif1 = ( y2 - y1 ) / ( t2 - t1 ); dif2 = ( ( y3 - y1 ) / ( t3 - t1 ) - ( y2 - y1 ) / ( t2 - t1 ) ) / ( t3 - t2 ); *yval = y1 + ( tval - t1 ) * ( dif1 + ( tval - t2 ) * dif2 ); return; } /******************************************************************************/ double r8poly3_discriminant ( double a, double b, double c, double d ) /******************************************************************************/ /* Purpose: r8poly3_discriminant() returns the discriminant of a cubic polynomial. Discussion: The polynomial has the form A * X^3 + B * X^2 + C * X + D = 0 If the discriminant is zero, then the polynomial has a multiple root. Licensing: This code is distributed under the MIT license. Modified: 14 July 2025 Author: John Burkardt Reference: John D Cook, When a cubic or quartic has a double root, https://www.johndcook.com/blog/2022/11/16/cubic-quartic-double-root/ Posted 16 November 2022. Input: real A, B, C, D, the coefficients of the polynomial. A must not be zero. Output: real VALUE, the discriminant. */ { double value; value = 18.0 * a * b * c * d - 4.0 * pow ( b, 3 ) * d + pow ( b, 2 ) * pow ( c, 2 ) - 4.0 * a * pow ( c, 3 ) - 27.0 * pow ( a, 2 ) * pow ( d, 2 ); return value; } /******************************************************************************/ void r8poly3_root ( double a, double b, double c, double d, double complex *r1, double complex *r2, double complex *r3 ) /******************************************************************************/ /* Purpose: r8poly3_root() returns the three roots of a cubic polynomial. Discussion: The polynomial has the form A * X^3 + B * X^2 + C * X + D = 0 Licensing: This code is distributed under the MIT license. Modified: 08 August 2018 Input: double A, B, C, D, the coefficients of the polynomial. A must not be zero. Output: double complex *R1, *R2, *R3, the roots of the polynomial, which will include at least one real root. */ { double complex i; double q; double r; const double r8_pi = 3.141592653589793; double s1; double s2; double temp; double theta; if ( a == 0.0 ) { printf ( "\n" ); printf ( "r8poly3_root(): Fatal error!\n" ); printf ( " A must not be zero!\n" ); exit ( 1 ); } i = I; q = ( ( b / a ) * ( b / a ) - 3.0 * ( c / a ) ) / 9.0; r = ( 2.0 * ( b / a ) * ( b / a ) * ( b / a ) - 9.0 * ( b / a ) * ( c / a ) + 27.0 * ( d / a ) ) / 54.0; if ( r * r < q * q * q ) { theta = acos ( r / sqrt ( q * q * q ) ); *r1 = - 2.0 * csqrt ( q ) * cos ( theta / 3.0 ); *r2 = - 2.0 * csqrt ( q ) * cos ( ( theta + 2.0 * r8_pi ) / 3.0 ); *r3 = - 2.0 * csqrt ( q ) * cos ( ( theta + 4.0 * r8_pi ) / 3.0 ); } else if ( q * q * q <= r * r ) { temp = - r + sqrt ( r * r - q * q * q ); s1 = r8_sign ( temp ) * pow ( fabs ( temp ), 1.0 / 3.0 ); temp = - r - sqrt ( r * r - q * q * q ); s2 = r8_sign ( temp ) * pow ( fabs ( temp ), 1.0 / 3.0 ); *r1 = s1 + s2; *r2 = - 0.5 * ( s1 + s2 ) + i * 0.5 * sqrt ( 3.0 ) * ( s1 - s2 ); *r3 = - 0.5 * ( s1 + s2 ) - i * 0.5 * sqrt ( 3.0 ) * ( s1 - s2 ); } *r1 = *r1 - b / ( 3.0 * a ); *r2 = *r2 - b / ( 3.0 * a ); *r3 = *r3 - b / ( 3.0 * a ); return; } /******************************************************************************/ double r8poly4_discriminant ( double a, double b, double c, double d, double e, double *p, double *r, double *d0, double *d1 ) /******************************************************************************/ /* Purpose: r8poly4_discriminant() returns the discriminant of a quartic polynomial. Discussion: The polynomial has the form a * x^4 + b * x^3 + c * x^2 + d * x + e = 0 If the discriminant is negative, then the polynomial has two distinct real roots and two distinct complex conjugate roots. If the discriminant is positive, then all four roots are real, or all four roots are strictly complex. If the discriminant is positive, let P=8*a*c-3*b^2 and D = 64*a^3*e - 16*a^2*c^2 + 16*a*b^2*c - 16*a^2*b*d - 3*b^4. If P<0 and D<0, then all four roots are real and distinct. If P>0 or D>0, there are two pairs of non-real complex conjugate roots Licensing: This code is distributed under the MIT license. Modified: 14 July 2025 Author: John Burkardt Reference: John D Cook, When a cubic or quartic has a double root, https://www.johndcook.com/blog/2022/11/16/cubic-quartic-double-root/ Posted 16 November 2022. Input: real A, B, C, D, E: the coefficients of the polynomial. A must not be zero. Output: real disc: the discriminant. real p, r, d0, d1: additional discriminant information. */ { double disc; disc = 256.0 * pow ( a, 3 ) * pow ( e, 3 ) - 192.0 * pow ( a, 2 ) * b * d * pow ( e, 2 ) - 128.0 * pow ( a, 2 ) * pow ( c, 2 ) * pow ( e, 2 ) + 144.0 * pow ( a, 2 ) * c * pow ( d, 2 ) * e - 27.0 * pow ( a, 2 ) * pow ( d, 4 ) + 144.0 * a * pow ( b, 2 ) * c * pow ( e, 2 ) - 6.0 * a * pow ( b, 2 ) * pow ( d, 2 ) * e - 80.0 * a * b * pow ( c, 2 ) * d * e + 18.0 * a * b * c * pow ( d, 3 ) + 16.0 * a * pow ( c, 4 ) * e - 4.0 * a * pow ( c, 3 ) * pow ( d, 2 ) - 27.0 * pow ( b, 4 ) * pow ( e, 2 ) + 18.0 * pow ( b, 3 ) * c * d * e - 4.0 * pow ( b, 3 ) * pow ( d, 3 ) - 4.0 * pow ( b, 2 ) * pow ( c, 3 ) * e + pow ( b, 2 ) * pow ( c, 2 ) * pow ( d, 2 ); *p = 8.0 * a * c - 3.0 * pow ( b, 2 ); *r = pow ( b, 3 ) + 8.0 * d * pow ( a, 2 ) - 4.0 * a * b * c; *d0 = pow ( c, 2 ) - 3.0 * b * d + 12.0 * a * e; *d1 = 64.0 * pow ( a, 3 ) * e - 16.0 * pow ( a, 2 ) * pow ( c, 2 ) + 16.0 * a * pow ( b, 2 ) * c - 16.0 * pow ( a, 2 ) * b * d - 3.0 * pow ( b, 4 ); return disc; } /******************************************************************************/ void r8poly4_root ( double a, double b, double c, double d, double e, double complex *r1, double complex *r2, double complex *r3, double complex *r4 ) /******************************************************************************/ /* Purpose: r8poly4_root() returns the four roots of a quartic polynomial. Discussion: The polynomial has the form: A * X^4 + B * X^3 + C * X^2 + D * X + E = 0 Licensing: This code is distributed under the MIT license. Modified: 08 August 2018 Input: double A, B, C, D, E, the coefficients of the polynomial. A must not be zero. Output: double complex *R1, *R2, *R3, *R4, the roots of the polynomial. */ { double a3; double a4; double b3; double b4; double c3; double c4; double d3; double d4; double complex p; double complex q; double complex r; double complex zero; zero = 0.0; if ( a == 0.0 ) { printf ( "\n" ); printf ( "R8POLY4_ROOT(): Fatal error!\n" ); printf ( " A must not be zero!\n" ); exit ( 1 ); } a4 = b / a; b4 = c / a; c4 = d / a; d4 = e / a; /* Set the coefficients of the resolvent cubic equation. */ a3 = 1.0; b3 = - b4; c3 = a4 * c4 - 4.0 * d4; d3 = - a4 * a4 * d4 + 4.0 * b4 * d4 - c4 * c4; /* Find the roots of the resolvent cubic. */ r8poly3_root ( a3, b3, c3, d3, r1, r2, r3 ); /* Choose one root of the cubic, here R1. Set R = sqrt ( 0.25 * A4 ^ 2 - B4 + R1 ) */ r = csqrt ( 0.25 * a4 * a4 - b4 + *r1 ); if ( r != zero ) { p = csqrt ( 0.75 * a4 * a4 - r * r - 2.0 * b4 + 0.25 * ( 4.0 * a4 * b4 - 8.0 * c4 - a4 * a4 * a4 ) / r ); q = csqrt ( 0.75 * a4 * a4 - r * r - 2.0 * b4 - 0.25 * ( 4.0 * a4 * b4 - 8.0 * c4 - a4 * a4 * a4 ) / r ); } else { p = csqrt ( 0.75 * a4 * a4 - 2.0 * b4 + 2.0 * sqrt ( *r1 * *r1 - 4.0 * d4 ) ); q = csqrt ( 0.75 * a4 * a4 - 2.0 * b4 - 2.0 * sqrt ( *r1 * *r1 - 4.0 * d4 ) ); } /* Set the roots. */ *r1 = -0.25 * a4 + 0.5 * r + 0.5 * p; *r2 = -0.25 * a4 + 0.5 * r - 0.5 * p; *r3 = -0.25 * a4 - 0.5 * r + 0.5 * q; *r4 = -0.25 * a4 - 0.5 * r - 0.5 * q; return; } /******************************************************************************/ void r8r8_print ( double a1, double a2, char *title ) /******************************************************************************/ /* Purpose: r8r8_print() prints an R8R8. Discussion: An R8R8 is a pair of R8 values, regarded as a single item. A format is used which suggests a coordinate pair: Example: Center : ( 1.23, 7.45 ) Licensing: This code is distributed under the MIT license. Modified: 15 May 2010 Author: John Burkardt Input: double A1, A2, the coordinates of the vector. char *TITLE, a title. */ { fprintf ( stdout, "%s: ( %g, %g )\n", title, a1, a2 ); return; } /******************************************************************************/ double *r8vec_copy_new ( int n, double a1[] ) /******************************************************************************/ /* Purpose: r8vec_copy_new() copies an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 26 August 2008 Author: John Burkardt Input: int N, the number of entries in the vectors. double A1[N], the vector to be copied. Output: double R8VEC_COPY_NEW[N], the copy of A1. */ { double *a2; int i; a2 = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { a2[i] = a1[i]; } return a2; } /******************************************************************************/ double *r8vec_even_new ( int n, double alo, double ahi ) /******************************************************************************/ /* Purpose: r8vec_even_new() returns an R8VEC of values evenly spaced between ALO and AHI. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 17 February 2004 Author: John Burkardt Input: int N, the number of values. double ALO, AHI, the low and high values. Output: double R8VEC_EVEN_NEW[N], N evenly spaced values. */ { double *a; int i; a = ( double * ) malloc ( n * sizeof ( double ) ); if ( n == 1 ) { a[0] = 0.5 * ( alo + ahi ); } else { for ( i = 1; i <= n; i++ ) { a[i-1] = ( ( double ) ( n - i ) * alo + ( double ) ( i - 1 ) * ahi ) / ( double ) ( n - 1 ); } } return a; } /******************************************************************************/ double r8vec_even_select ( int n, double xlo, double xhi, int ival ) /******************************************************************************/ /* Purpose: r8vec_even_select() returns the I-th of N evenly spaced values in [ XLO, XHI ]. Discussion: An R8VEC is a vector of R8's. XVAL = ( (N-IVAL) * XLO + (IVAL-1) * XHI ) / ( N - 1 ) Licensing: This code is distributed under the MIT license. Modified: 20 August 2010 Author: John Burkardt Input: int N, the number of values. double XLO, XHI, the low and high values. int IVAL, the index of the desired point. IVAL is normally between 1 and N, but may be any integer value. Output: double R8VEC_EVEN_SELECT, the IVAL-th of N evenly spaced values between XLO and XHI. Unless N = 1, X(1) = XLO and X(N) = XHI. If N = 1, then X(1) = 0.5*(XLO+XHI). */ { double xval; if ( n == 1 ) { xval = 0.5 * ( xlo + xhi ); } else { xval = ( ( double ) ( n - ival ) * xlo + ( double ) ( ival - 1 ) * xhi ) / ( double ) ( n - 1 ); } return xval; } /******************************************************************************/ double *r8vec_indicator1_new ( int n ) /******************************************************************************/ /* Purpose: r8vec_indicator1_new() sets an R8VEC to the indicator vector {1,2,3,...}. Licensing: This code is distributed under the MIT license. Modified: 27 September 2014 Author: John Burkardt Input: int N, the number of elements of A. Output: double R8VEC_INDICATOR1_NEW[N], the array. */ { double *a; int i; a = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { a[i] = ( double ) ( i + 1 ); } return a; } /******************************************************************************/ bool r8vec_is_distinct ( int n, double x[] ) /******************************************************************************/ /* Purpose: r8vec_is_distinct() is true if the entries in an R8VEC are distinct. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 31 March 2018 Author: John Burkardt Input: int N, the number of entries in the vector. double X[N], the vector to be checked. Output: bool R8VEC_IS_DISTINCT is true if all N elements of X are distinct. */ { int i; int j; bool value; value = true; for ( i = 1; i < n; i++ ) { for ( j = 0; j < i; j++ ) { if ( x[i] == x[j] ) { value = false; break; } } } return value; } /******************************************************************************/ double *r8vec_linspace_new ( int n, double a, double b ) /******************************************************************************/ /* Purpose: r8vec_linspace_new() creates a vector of linearly spaced values. Discussion: An R8VEC is a vector of R8's. 4 points evenly spaced between 0 and 12 will yield 0, 4, 8, 12. In other words, the interval is divided into N-1 even subintervals, and the endpoints of intervals are used as the points. Licensing: This code is distributed under the MIT license. Modified: 29 March 2011 Author: John Burkardt Input: int N, the number of entries in the vector. double A, B, the first and last entries. Output: double R8VEC_LINSPACE_NEW[N], a vector of linearly spaced data. */ { int i; double *x; x = ( double * ) malloc ( n * sizeof ( double ) ); if ( n == 1 ) { x[0] = ( a + b ) / 2.0; } else { for ( i = 0; i < n; i++ ) { x[i] = ( ( double ) ( n - 1 - i ) * a + ( double ) ( i ) * b ) / ( double ) ( n - 1 ); } } return x; } /******************************************************************************/ void r8vec_print ( int n, double a[], char *title ) /******************************************************************************/ /* Purpose: r8vec_print() prints an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 08 April 2009 Author: John Burkardt Input: int N, the number of components of the vector. double A[N], the vector to be printed. char *TITLE, a title. */ { int i; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); fprintf ( stdout, "\n" ); for ( i = 0; i < n; i++ ) { fprintf ( stdout, " %8d: %14g\n", i, a[i] ); } return; } /******************************************************************************/ void r8vec_transpose_print ( int n, double a[], char *title ) /******************************************************************************/ /* Purpose: r8vec_transpose_print() prints an R8VEC "transposed". Discussion: An R8VEC is a vector of R8's. Example: A = (/ 1.0, 2.1, 3.2, 4.3, 5.4, 6.5, 7.6, 8.7, 9.8, 10.9, 11.0 /) TITLE = 'My vector: ' My vector: 1.0 2.1 3.2 4.3 5.4 6.5 7.6 8.7 9.8 10.9 11.0 Licensing: This code is distributed under the MIT license. Modified: 12 May 2014 Author: John Burkardt Input: int N, the number of components of the vector. double A[N], the vector to be printed. char *TITLE, a title. */ { int i; int ihi; int ilo; int title_length; title_length = s_len_trim ( title ); for ( ilo = 0; ilo < n; ilo = ilo + 5 ) { if ( ilo == 0 ) { printf ( "%s", title ); } else { for ( i = 0; i < title_length; i++ ) { printf ( " " ); } } printf ( " " ); ihi = i4_min ( ilo + 5, n ); for ( i = ilo; i < ihi; i++ ) { printf ( " %12g", a[i] ); } printf ( "\n" ); } return; } /******************************************************************************/ double *r8vec_uniform_01_new ( int n, int *seed ) /******************************************************************************/ /* Purpose: r8vec_uniform_01_new() returns a unit pseudorandom R8VEC. Discussion: This routine implements the recursion seed = 16807 * seed mod ( 2^31 - 1 ) unif = seed / ( 2^31 - 1 ) The integer arithmetic never requires more than 32 bits, including a sign bit. Licensing: This code is distributed under the MIT license. Modified: 19 August 2004 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Second Edition, Springer, 1987, ISBN: 0387964673, LC: QA76.9.C65.B73. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, December 1986, pages 362-376. Pierre L'Ecuyer, Random Number Generation, in Handbook of Simulation, edited by Jerry Banks, Wiley, 1998, ISBN: 0471134031, LC: T57.62.H37. Peter Lewis, Allen Goodman, James Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, Number 2, 1969, pages 136-143. Input: int N, the number of entries in the vector. int *SEED, a seed for the random number generator. Output: int *SEED: the updated seed. double R8VEC_UNIFORM_01_NEW[N], the vector of pseudorandom values. */ { int i; const int i4_huge = 2147483647; int k; double *r; if ( *seed == 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8VEC_UNIFORM_01_NEW(): Fatal error!\n" ); fprintf ( stderr, " Input value of SEED = 0.\n" ); exit ( 1 ); } r = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + i4_huge; } r[i] = ( double ) ( *seed ) * 4.656612875E-10; } return r; } /******************************************************************************/ double *r8vec_zeros_new ( int n ) /******************************************************************************/ /* Purpose: r8vec_zeros_new() creates and zeroes an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 25 March 2009 Author: John Burkardt Input: int N, the number of entries in the vector. Output: double R8VEC_ZEROS_NEW[N], a vector of zeroes. */ { double *a; int i; a = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { a[i] = 0.0; } return a; } /******************************************************************************/ void r8vec2_print ( int n, double a1[], double a2[], char *title ) /******************************************************************************/ /* Purpose: r8vec2_print() prints an R8VEC2. Discussion: An R8VEC2 is a dataset consisting of N pairs of real values, stored as two separate vectors A1 and A2. Licensing: This code is distributed under the MIT license. Modified: 26 March 2009 Author: John Burkardt Input: int N, the number of components of the vector. double A1[N], double A2[N], the vectors to be printed. char *TITLE, a title. */ { int i; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); fprintf ( stdout, "\n" ); for ( i = 0; i < n; i++ ) { fprintf ( stdout, " %4d: %14g %14g\n", i, a1[i], a2[i] ); } return; } /******************************************************************************/ double *roots_to_r8poly ( int n, double x[] ) /******************************************************************************/ /* Purpose: roots_to_r8poly() converts polynomial roots to polynomial coefficients. Licensing: This code is distributed under the MIT license. Modified: 19 July 2010 Author: John Burkardt Input: int N, the number of roots specified. double X[N], the roots. Output: double ROOTS_TO_R8POLY[N+1], the coefficients of the polynomial. */ { double *c; int i; int j; c = r8vec_zeros_new ( n + 1 ); /* Initialize C to (0, 0, ..., 0, 1). Essentially, we are setting up a divided difference table. */ c[n] = 1.0; /* Convert to standard polynomial form by shifting the abscissas of the divided difference table to 0. */ for ( j = 1; j <= n; j++ ) { for ( i = 1; i <= n + 1 - j; i++ ) { c[n-i] = c[n-i] - x[n+1-i-j] * c[n-i+1]; } } return c; } /******************************************************************************/ int s_len_trim ( char *s ) /******************************************************************************/ /* Purpose: s_len_trim() returns the length of a string to the last nonblank. Licensing: This code is distributed under the MIT license. Modified: 26 April 2003 Author: John Burkardt 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 != ' ' ) { return n; } t--; n--; } return n; }