# include # include # include # include # include # include "companion_matrix.h" # include "eigs.h" # include int main ( ); void companion_chebyshev_test ( ); void companion_gegenbauer_test ( ); void companion_hermite_test ( ); void companion_laguerre_test ( ); void companion_legendre_test ( ); void companion_monomial_test ( ); double *chebyshev_to_monomial_matrix ( int n ); double *gegenbauer_to_monomial_matrix ( int n, double alpha ); double *hermite_to_monomial_matrix ( int n ); double *laguerre_to_monomial_matrix ( int n ); double *legendre_to_monomial_matrix ( int n ); void polynomial_chebyshev_print ( int d, double p[], char *label ); void polynomial_gegenbauer_print ( int d, double p[], char *label ); void polynomial_hermite_print ( int d, double p[], char *label ); void polynomial_laguerre_print ( int d, double p[], char *label ); void polynomial_legendre_print ( int d, double p[], char *label ); void polynomial_monomial_print ( int d, double p[], char *label ); void c8vec_print ( int n, double complex a[], char *title ); double *r8mat_mv_new ( int m, int n, double a[], double x[] ); void r8mat_print ( int m, int n, double a[], char *title ); void r8mat_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi, int jhi, char *title ); double complex *r8poly_roots ( int n, double p[] ); void timestamp ( ); /******************************************************************************/ int main ( ) /******************************************************************************/ /* Purpose: companion_matrix_test() tests companion_matrix(). Licensing: This code is distributed under the MIT license. Modified: 11 June 2024 Author: John Burkardt */ { timestamp ( ); printf ( "\n" ); printf ( "companion_matrix_test():\n" ); printf ( " C version \n" ); printf ( " companion_matrix() computes the companion matrix\n" ); printf ( " of a polynomial, in various bases.\n" ); companion_chebyshev_test ( ); companion_gegenbauer_test ( ); companion_hermite_test ( ); companion_laguerre_test ( ); companion_legendre_test ( ); companion_monomial_test ( ); /* Terminate. */ printf ( "\n" ); printf ( "companion_matrix_test():\n" ); printf ( " Normal end of execution.\n" ); printf ( "\n" ); timestamp ( ); return 0; } /******************************************************************************/ void companion_chebyshev_test ( ) /******************************************************************************/ /* Purpose: companion_chebyshev_test() tests companion_chebyshev(). Licensing: This code is distributed under the MIT license. Modified: 10 June 2024 Author: John Burkardt */ { double *A; double *C; int n = 5; double p[6] = { 6.0, 5.0, 4.0, 3.0, 2.0, 1.0 }; double *q; double complex *r; double complex *s; printf ( "\n" ); printf ( "companion_chebyshev_test():\n" ); printf ( " companion_chebyshev() computes the companion matrix\n" ); printf ( " of a polynomial p(x) in the Chebyshev basis.\n" ); polynomial_chebyshev_print ( n, p, " p(x)" ); /* Compute monomial form so we can get roots directly. */ A = chebyshev_to_monomial_matrix ( n + 1 ); q = r8mat_mv_new ( n + 1, n + 1, A, p ); polynomial_monomial_print ( n, q, " Monomial q(x) =" ); r = r8poly_roots ( n, q ); c8vec_print ( n, r, " Roots of q(x):" ); C = companion_chebyshev ( n, p ); r8mat_print ( n, n, C, " Chebyshev companion matrix C(p):" ); s = eigs ( n, C ); c8vec_print ( n, s, " Eigenvalues of C(p):" ); free ( A ); free ( C ); free ( q ); free ( r ); free ( s ); return; } /******************************************************************************/ void companion_gegenbauer_test ( ) /******************************************************************************/ /* Purpose: companion_gegenbauer_test() tests companion_gegenbauer(). Licensing: This code is distributed under the MIT license. Modified: 11 June 2024 Author: John Burkardt */ { double *A; double alpha; double *C; int i; int n = 5; double p[6] = { 6.0, 5.0, 4.0, 3.0, 2.0, 1.0 }; double *q; double complex *r; double complex *s; printf ( "\n" ); printf ( "companion_gegenbauer_test():\n" ); printf ( " companion_gegenbauer() computes the companion matrix\n" ); printf ( " of a polynomial p(x) in the Gegenbauer basis.\n" ); polynomial_chebyshev_print ( n, p, " p(x)" ); /* Try several parameter values. */ for ( i = 1; i <= 2; i++ ) { alpha = i / 2.0; printf ( "\n" ); printf ( " Gegenbauer parameter alpha will be %g\n", alpha ); /* Compute monomial form so we can get roots directly. */ A = gegenbauer_to_monomial_matrix ( n + 1, alpha ); q = r8mat_mv_new ( n + 1, n + 1, A, p ); polynomial_monomial_print ( n, q, " Monomial q(x) =" ); r = r8poly_roots ( n, q ); c8vec_print ( n, r, " Roots of q(x):" ); C = companion_gegenbauer ( n, p, alpha ); r8mat_print ( n, n, C, " Gegenbauer companion matrix C(p):" ); s = eigs ( n, C ); c8vec_print ( n, s, " Eigenvalues of C(p):" ); free ( A ); free ( C ); free ( q ); free ( r ); free ( s ); } return; } /******************************************************************************/ void companion_hermite_test ( ) /******************************************************************************/ /* Purpose: companion_hermite_test() tests companion_hermite(). Licensing: This code is distributed under the MIT license. Modified: 11 June 2024 Author: John Burkardt */ { double *A; double *C; int n = 5; double p[6] = { 6.0, 5.0, 4.0, 3.0, 2.0, 1.0 }; double *q; double complex *r; double complex *s; printf ( "\n" ); printf ( "companion_hermite_test():\n" ); printf ( " companion_hermite() computes the companion matrix\n" ); printf ( " of a polynomial p(x) in the Hermite basis.\n" ); polynomial_hermite_print ( n, p, " p(x)" ); /* Compute monomial form so we can get roots directly. */ A = hermite_to_monomial_matrix ( n + 1 ); q = r8mat_mv_new ( n + 1, n + 1, A, p ); polynomial_monomial_print ( n, q, " Monomial q(x) =" ); r = r8poly_roots ( n, q ); c8vec_print ( n, r, " Roots of q(x):" ); C = companion_hermite ( n, p ); r8mat_print ( n, n, C, " Hermite companion matrix C(p):" ); s = eigs ( n, C ); c8vec_print ( n, s, " Eigenvalues of C(p):" ); free ( A ); free ( C ); free ( q ); free ( r ); free ( s ); return; } /******************************************************************************/ void companion_laguerre_test ( ) /******************************************************************************/ /* Purpose: companion_laguerre_test() tests companion_laguerre(). Licensing: This code is distributed under the MIT license. Modified: 11 June 2024 Author: John Burkardt */ { double *A; double *C; int n = 5; double p[6] = { 6.0, 5.0, 4.0, 3.0, 2.0, 1.0 }; double *q; double complex *r; double complex *s; printf ( "\n" ); printf ( "companion_laguerre_test():\n" ); printf ( " companion_laguerre() computes the companion matrix\n" ); printf ( " of a polynomial p(x) in the Laguerre basis.\n" ); polynomial_laguerre_print ( n, p, " p(x)" ); /* Compute monomial form so we can get roots directly. */ A = laguerre_to_monomial_matrix ( n + 1 ); q = r8mat_mv_new ( n + 1, n + 1, A, p ); polynomial_monomial_print ( n, q, " Monomial q(x) =" ); r = r8poly_roots ( n, q ); c8vec_print ( n, r, " Roots of q(x):" ); C = companion_laguerre ( n, p ); r8mat_print ( n, n, C, " Laguerre companion matrix C(p):" ); s = eigs ( n, C ); c8vec_print ( n, s, " Eigenvalues of C(p):" ); free ( A ); free ( C ); free ( q ); free ( r ); free ( s ); return; } /******************************************************************************/ void companion_legendre_test ( ) /******************************************************************************/ /* Purpose: companion_legendre_test() tests companion_legendre(). Licensing: This code is distributed under the MIT license. Modified: 11 June 2024 Author: John Burkardt */ { double *A; double *C; int n = 5; double p[6] = { 6.0, 5.0, 4.0, 3.0, 2.0, 1.0 }; double *q; double complex *r; double complex *s; printf ( "\n" ); printf ( "companion_legendre_test():\n" ); printf ( " companion_legendre() computes the companion matrix\n" ); printf ( " of a polynomial p(x) in the Legendre basis.\n" ); polynomial_legendre_print ( n, p, " p(x)" ); /* Compute monomial form so we can get roots directly. */ A = legendre_to_monomial_matrix ( n + 1 ); q = r8mat_mv_new ( n + 1, n + 1, A, p ); polynomial_monomial_print ( n, q, " Monomial q(x) =" ); r = r8poly_roots ( n, q ); c8vec_print ( n, r, " Roots of q(x):" ); C = companion_legendre ( n, p ); r8mat_print ( n, n, C, " Legendre companion matrix C(p):" ); s = eigs ( n, C ); c8vec_print ( n, s, " Eigenvalues of C(p):" ); free ( A ); free ( C ); free ( q ); free ( r ); free ( s ); return; } /******************************************************************************/ void companion_monomial_test ( ) /******************************************************************************/ /* Purpose: companion_monomial_test() tests companion_monomial(). Licensing: This code is distributed under the MIT license. Modified: 10 June 2024 Author: John Burkardt */ { double *C; int n = 5; double p[6] = { 6.0, 5.0, 4.0, 3.0, 2.0, 1.0 }; double complex *r; double complex *s; printf ( "\n" ); printf ( "companion_monomial_test():\n" ); printf ( " companion_monomial() computes the companion matrix\n" ); printf ( " of a polynomial p(x) in the monomial basis.\n" ); polynomial_monomial_print ( n, p, " p(x) =" ); /* Get roots directly. */ r = r8poly_roots ( n, p ); c8vec_print ( n, r, " Roots of q(x):" ); C = companion_monomial ( n, p ); r8mat_print ( n, n, C, " Monomial companion matrix C(p):" ); s = eigs ( n, C ); c8vec_print ( n, s, " Eigenvalues of C(p):" ); free ( C ); free ( r ); free ( s ); return; } /******************************************************************************/ double *chebyshev_to_monomial_matrix ( int n ) /******************************************************************************/ /* Purpose: chebyshev_to_monomial_matrix() returns the chebyshev_to_monomial matrix. Discussion This is the Chebyshev T matrix, associated with the Chebyshev "T" polynomials, or Chebyshev polynomials of the first kind. Example: N = 11 1 . -1 . 1 . -1 . 1 . -1 . 1 . -3 . 5 . -7 . 9 . . . 2 . -8 . 18 . -32 . 50 . . . 4 . -20 . 56 . -120 . . . . . 8 . -48 . 160 . -400 . . . . . 16 . -112 . 432 . . . . . . . 32 . -256 . 1120 . . . . . . . 64 . -576 . . . . . . . . . 128 . -1280 . . . . . . . . . 256 . . . . . . . . . . . 512 Properties: A is generally not symmetric: A' /= A. A is integral, therefore det ( A ) is integral, and det ( A ) * inverse ( A ) is integral. A is reducible. A is lower triangular. Each row of A sums to 1. det ( A ) = 2^( (N-1) * (N-2) / 2 ) A is not normal: A' * A /= A * A'. For I = 1: LAMBDA(1) = 1 For 1 < I LAMBDA(I) = 2^(I-2) The family of matrices is nested as a function of N. Licensing: This code is distributed under the MIT license. Modified: 04 April 2024 Author: John Burkardt Input: int N, the order of the matrix. Output: double CHEBY_T[N*N]: the matrix. */ { double *A; int i; int j; if ( n <= 0 ) { return NULL; } A = ( double * ) malloc ( n * n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { for ( j = 0; j < n; j++ ) { A[i+j*n] = 0.0; } } A[0+0*n] = 1.0; if ( n == 1 ) { return A; } A[1+1*n] = 1.0; if ( n == 2 ) { return A; } for ( j = 2; j < n; j++ ) { i = 0; A[i+j*n] = - A[i+(j-2)*n]; for ( i = 1; i < n; i++ ) { A[i+j*n] = 2.0 * A[i-1+(j-1)*n] - A[i+(j-2)*n]; } } return A; } /******************************************************************************/ double *gegenbauer_to_monomial_matrix ( int n, double alpha ) /******************************************************************************/ /* Purpose: gegenbauer_to_monomial_matrix(): Gegenbauer to monomial conversion matrix. Licensing: This code is distributed under the MIT license. Modified: 05 April 2024 Author: John Burkardt Input: int N: the order of A. double alpha: the parameter. Output: double A(N,N): the matrix. */ { double *A; double c1; double c2; int i; int j; int nn; A = ( double * ) malloc ( n * n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { for ( i = 0; i < n; i++ ) { A[i+j*n] = 0.0; } } if ( n <= 0 ) { return A; } A[0+0*n] = 1.0; if ( n == 1 ) { return A; } A[1+1*n] = 2.0 * alpha; /* Perform convex sum. Translating "(n+1) C(n+1) = 2 (n+alpha) x C(n) - ( n + 2 alpha - 1 ) C(n-1)" drove me nuts, between indexing at 1 rather than 0, and dealing with the interpretation of "n+1", because I now face the rare "off by 2" error! */ for ( j = 2; j < n; j++ ) { nn = j - 1; c1 = ( 2.0 * nn + 2.0 * alpha ) / ( nn + 1 ); c2 = ( - nn - 2.0 * alpha + 1.0 ) / ( nn + 1 ); for ( i = 1; i <= j; i++ ) { A[i+j*n] = c1 * A[i-1+(j-1)*n]; } for ( i = 0; i <= j - 2; i++ ) { A[i+j*n] = A[i+j*n] + c2 * A[i+(j-2)*n]; } } return A; } /******************************************************************************/ double *hermite_to_monomial_matrix ( int n ) /******************************************************************************/ /* Purpose: hermite_to_monomial_matrix() returns the HERMITE matrix. Example: N = 11 (Each column must be divided by the factor below it). 1 . 2 . 12 . 120 . 1680 . 30240 . 1 . 6 . 60 . 840 . 15120 . . . 1 . 12 . 180 . 3360 . 75600 . . . 1 . 20 . 420 . 10080 . . . . . 1 . 30 . 840 . 25200 . . . . . 1 . 42 . 1512 . . . . . . . 1 . 56 . 2520 . . . . . . . 1 . 72 . . . . . . . . . 1 . 90 . . . . . . . . . 1 . . . . . . . . . . . 1 /1 /2 /4 /8 /16 /32 /64 /128 /256 /512 /1024 Properties: A is generally not symmetric: A' /= A. A is lower triangular. det ( A ) = 2^((N*(N-1))/2) LAMBDA(I) = 2^(I-1). A is integral, therefore det ( A ) is integral, and det ( A ) * inverse ( A ) is integral. A is reducible. Successive diagonals are zero, positive, zero, negative. A is generally not normal: A' * A /= A * A'. The family of matrices is nested as a function of N. Licensing: This code is distributed under the MIT license. Modified: 04 April 2024 Author: John Burkardt Input: int N, the order of the matrix. Output: double HERMITE[N*N], the matrix. */ { double *A; int i; int j; if ( n <= 0 ) { return NULL; } A = ( double * ) malloc ( n * n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { for ( i = 0; i < n; i++ ) { A[i+j*n] = 0.0; } } A[0+0*n] = 1.0; if ( n == 1 ) { return A; } A[1+1*n] = 2.0; if ( n == 2 ) { return A; } for ( j = 2; j < n; j++ ) { for ( i = 0; i < n; i++ ) { if ( i == 0 ) { A[i+j*n] = - 2.0 * ( double ) ( j - 1 ) * A[i+(j-2)*n]; } else { A[i+j*n] = - 2.0 * ( double ) ( j - 1 ) * A[i+(j-2)*n] + 2.0 * A[i-1+(j-1)*n]; } } } return A; } /******************************************************************************/ double *laguerre_to_monomial_matrix ( int n ) /******************************************************************************/ /* Purpose: laguerre_to_monomial_matrix() returns the LAGUERRE matrix. Example: N = 8 (each column must be divided by the factor below it.) 1 1 2 6 24 120 720 5040 . -1 -4 -18 -96 -600 -4320 -35280 . . 1 9 72 600 5400 52920 . . . 1 -16 -200 -2400 -29400 . . . . 1 25 450 7350 . . . . . -1 -36 -882 . . . . . . 1 49 . . . . . . . -1 /1 /1 /2 /6 /24 /120 /720 /5040 Properties: A is generally not symmetric: A' /= A. A is lower triangular. The columns of A alternate in sign. A(I,1) = 1 A(I,I) = (-1)^(I-1) / (I-1). LAMBDA(I) = (-1)^(I-1) / (I-1). det ( A ) = product ( 1 <= I <= N ) (-1)^(I-1) / (I-1) The I-th row contains the coefficients of the Laguerre polynomial of degree I-1. The family of matrices is nested as a function of N. Licensing: This code is distributed under the MIT license. Modified: 04 April 2024 Author: John Burkardt Reference: Milton Abramowitz, Irene Stegun, Handbook of Mathematical Functions, US Department of Commerce, 1964. Input: int N, the order of the matrix. Output: double LAGUERRE[N*N], the matrix. */ { double *A; int i; int j; if ( n <= 0 ) { return NULL; } A = ( double * ) malloc ( n * n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { for ( i = 0; i < n; i++ ) { A[i+j*n] = 0.0; } } A[0+0*n] = 1.0; if ( n == 1 ) { return A; } A[0+1*n] = 1.0; A[1+1*n] = - 1.0; if ( n == 2 ) { return A; } for ( j = 2; j < n; j++ ) { A[0+j*n] = ( ( double ) ( 2 * j - 1 ) * A[0+(j-1)*n] + ( double ) ( - j + 1 ) * A[0+(j-2)*n] ) / ( double ) ( j ); for ( i = 1; i < n; i++ ) { A[i+j*n] = ( ( double ) ( 2 * j - 1 ) * A[i+(j-1)*n] - ( double ) ( 1 ) * A[i-1+(j-1)*n] + ( double ) ( - j + 1 ) * A[i+(j-2)*n] ) / ( double ) ( j ); } } return A; } /******************************************************************************/ double *legendre_to_monomial_matrix ( int n ) /******************************************************************************/ /* Purpose: legendre_to_monomial_matrix() returns the LEGENDRE matrix. Example: N = 11 (each column must be divided by factor at bottom) 1 . -1 . 3 . -5 . 35 . -63 . 1 . -3 . 15 . -25 . 315 . . . 3 . -30 . 105 . -1260 . 3465 . . . 5 . -70 . 315 . -4620 . . . . . 35 . -315 . 6930 .-30030 . . . . . 63 . -693 . 18018 . . . . . . . 231 . -12012 . 90090 . . . . . . . 429 .-25740 . . . . . . . . . 6435 -109395 . . . . . . . . . 12155 . . . . . . . . . . . 46189 /1 /1 /2 /2 /8 /8 /16 /16 /128 /128 /256 Properties: A is generally not symmetric: A' /= A. A is lower triangular. The elements of each row sum to 1. Because it has a constant row sum of 1, A has an eigenvalue of 1, and a (right) eigenvector of ( 1, 1, 1, ..., 1 ). A is reducible. The diagonals form a pattern of zero, positive, zero, negative. The family of matrices is nested as a function of N. Licensing: This code is distributed under the MIT license. Modified: 19 February 2024 Author: John Burkardt Input: int N, the order of the matrix. Output: double LEGENDRE_MATRIX[N*N], the matrix. */ { double *A; int i; int j; if ( n <= 0 ) { return NULL; } A = ( double * ) malloc ( n * n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { for ( i = 0; i < n; i++ ) { A[i+j*n] = 0.0; } } A[0+0*n] = 1.0; if ( n == 1 ) { return A; } A[1+1*n] = 1.0; if ( n == 2 ) { return A; } for ( j = 3; j <= n; j++ ) { for ( i = 1; i <= n; i++ ) { if ( i == 1 ) { A[i-1+(j-1)*n] = - ( j - 2 ) * A[i-1+(j-3)*n] / ( j - 1 ); } else { A[i-1+(j-1)*n] = ( ( 2 * j - 3 ) * A[i-2+(j-2)*n] + ( - j + 2 ) * A[i-1+(j-3)*n] ) / ( j - 1 ); } } } return A; } /******************************************************************************/ void polynomial_chebyshev_print ( int d, double p[], char *label ) /******************************************************************************/ /* Purpose: polynomial_chebyshev_print() prints a Chebyshev polynomial. Discussion: The power sum form is: p(x) = a(0)*T0(x) + a(1) * T1(x) + ... + a(n-1) * T(n-1)(x) + a(n) * Tn(x) Licensing: This code is distributed under the MIT license. Modified: 10 June 2024 Author: John Burkardt Input: integer d: the degree of the polynomial. real p(0:d), the polynomial coefficients. p(0) is the constant term and p(N) is the coefficient of X^N. char *label: a title. */ { int i; double mag; char plus_minus; printf ( "\n" ); printf ( "%s\n", label ); if ( d < 0 ) { printf ( " 0\n" ); return; } for ( i = d; 0 <= i; i-- ) { if ( p[i] < 0.0 ) { plus_minus = '-'; } else { plus_minus = '+'; } mag = fabs ( p[i] ); if ( mag != 0.0 ) { printf ( " %c%g * T%d(x)\n", plus_minus, mag, i ); } } return; } /******************************************************************************/ void polynomial_gegenbauer_print ( int d, double p[], char *label ) /******************************************************************************/ /* Purpose: polynomial_gegenbauer_print() prints a Gegenbauer polynomial. Discussion: The power sum form is: p(x) = a(0)*C0(x) + a(1) * C1(x) + ... + a(n-1) * C(n-1)(x) + a(n) * Cn(x) Licensing: This code is distributed under the MIT license. Modified: 10 June 2024 Author: John Burkardt Input: integer d: the degree of the polynomial. real p(0:d), the polynomial coefficients. p(0) is the constant term and p(N) is the coefficient of X^N. char *label: a title. */ { int i; double mag; char plus_minus; printf ( "\n" ); printf ( "%s\n", label ); if ( d < 0 ) { printf ( " 0\n" ); return; } for ( i = d; 0 <= i; i-- ) { if ( p[i] < 0.0 ) { plus_minus = '-'; } else { plus_minus = '+'; } mag = fabs ( p[i] ); if ( mag != 0.0 ) { printf ( " %c%g * C%d(x)\n", plus_minus, mag, i ); } } return; } /******************************************************************************/ void polynomial_hermite_print ( int d, double p[], char *label ) /******************************************************************************/ /* Purpose: polynomial_hermite_print() prints a Hermite polynomial. Discussion: The power sum form is: p(x) = a(0)*H0(x) + a(1) * H1(x) + ... + a(n-1) * H(n-1)(x) + a(n) * Hn(x) Licensing: This code is distributed under the MIT license. Modified: 10 June 2024 Author: John Burkardt Input: integer d: the degree of the polynomial. real p(0:d), the polynomial coefficients. p(0) is the constant term and p(N) is the coefficient of X^N. char *label: a title. */ { int i; double mag; char plus_minus; printf ( "\n" ); printf ( "%s\n", label ); if ( d < 0 ) { printf ( " 0\n" ); return; } for ( i = d; 0 <= i; i-- ) { if ( p[i] < 0.0 ) { plus_minus = '-'; } else { plus_minus = '+'; } mag = fabs ( p[i] ); if ( mag != 0.0 ) { printf ( " %c%g * H%d(x)\n", plus_minus, mag, i ); } } return; } /******************************************************************************/ void polynomial_laguerre_print ( int d, double p[], char *label ) /******************************************************************************/ /* Purpose: polynomial_laguerre_print() prints a Laguerre polynomial. Discussion: The power sum form is: p(x) = a(0)*L0(x) + a(1) * L1(x) + ... + a(n-1) * L(n-1)(x) + a(n) * Ln(x) Licensing: This code is distributed under the MIT license. Modified: 10 June 2024 Author: John Burkardt Input: integer d: the degree of the polynomial. real p(0:d), the polynomial coefficients. p(0) is the constant term and p(N) is the coefficient of X^N. char *label: a title. */ { int i; double mag; char plus_minus; printf ( "\n" ); printf ( "%s\n", label ); if ( d < 0 ) { printf ( " 0\n" ); return; } for ( i = d; 0 <= i; i-- ) { if ( p[i] < 0.0 ) { plus_minus = '-'; } else { plus_minus = '+'; } mag = fabs ( p[i] ); if ( mag != 0.0 ) { printf ( " %c%g * L%d(x)\n", plus_minus, mag, i ); } } return; } /******************************************************************************/ void polynomial_legendre_print ( int d, double p[], char *label ) /******************************************************************************/ /* Purpose: polynomial_legendre_print() prints a Legendre polynomial. Discussion: The power sum form is: p(x) = a(0)*P0(x) + a(1) * P1(x) + ... + a(n-1) * P(n-1)(x) + a(n) * Pn(x) Licensing: This code is distributed under the MIT license. Modified: 10 June 2024 Author: John Burkardt Input: integer d: the degree of the polynomial. real p(0:d), the polynomial coefficients. p(0) is the constant term and p(N) is the coefficient of X^N. char *label: a title. */ { int i; double mag; char plus_minus; printf ( "\n" ); printf ( "%s\n", label ); if ( d < 0 ) { printf ( " 0\n" ); return; } for ( i = d; 0 <= i; i-- ) { if ( p[i] < 0.0 ) { plus_minus = '-'; } else { plus_minus = '+'; } mag = fabs ( p[i] ); if ( mag != 0.0 ) { printf ( " %c%g * P%d(x)\n", plus_minus, mag, i ); } } return; } /******************************************************************************/ void polynomial_monomial_print ( int d, double p[], char *label ) /******************************************************************************/ /* Purpose: polynomial_monomial_print() prints 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: 10 June 2024 Author: John Burkardt Input: integer d: the degree of the polynomial. real p(0:d), the polynomial coefficients. p(0) is the constant term and p(N) is the coefficient of X^N. char *label: a title. */ { int i; double mag; char plus_minus; printf ( "\n" ); printf ( "%s\n", label ); if ( d < 0 ) { printf ( " 0\n" ); return; } for ( i = d; 0 <= i; i-- ) { if ( p[i] < 0.0 ) { plus_minus = '-'; } else { plus_minus = '+'; } mag = fabs ( p[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 c8vec_print ( int n, double complex a[], char *title ) /******************************************************************************/ /* Purpose: c8vec_print() prints a C8VEC. Discussion: A C8VEC is a vector of double complex values. Licensing: This code is distributed under the MIT license. Modified: 27 January 2013 Author: John Burkardt Input: int N, the number of components of the vector. double complex 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: %14f %14f\n", i, creal ( a[i] ), cimag ( a[i] ) ); } return; } /******************************************************************************/ double *r8mat_mv_new ( int m, int n, double a[], double x[] ) /******************************************************************************/ /* Purpose: r8mat_mv_new() multiplies a matrix times a vector. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. For this routine, the result is returned as the function value. Licensing: This code is distributed under the MIT license. Modified: 15 April 2014 Author: John Burkardt Input: int M, N, the number of rows and columns of the matrix. double A[M,N], the M by N matrix. double X[N], the vector to be multiplied by A. Output: double R8MAT_MV_NEW[M], the product A*X. */ { int i; int j; double *y; y = ( double * ) malloc ( m * sizeof ( double ) ); for ( i = 0; i < m; i++ ) { y[i] = 0.0; for ( j = 0; j < n; j++ ) { y[i] = y[i] + a[i+j*m] * x[j]; } } return y; } /******************************************************************************/ 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 complex *r8poly_roots ( int n, double p[] ) /******************************************************************************/ /* Purpose: r8poly_roots() returns the roots of a polynomial. Discussion: The monomial form of a polynomial is: p(x) = c(0)*x^0 + c(1)*x + ... + c(n)*x^(n) Licensing: This code is distributed under the MIT license. Modified: 10 June 2024 Author: John Burkardt Input: integer n: the polynomial degree. real p(n+1): the polynomial coefficients, constant first. Output: double complex r(n): the polynomial roots. */ { int i; double complex *r; double *z; /* Use GSL to compute roots. */ z = ( double * ) malloc ( 2 * n * sizeof ( double ) ); gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc ( n + 1 ); gsl_poly_complex_solve ( p, n + 1, w, z ); gsl_poly_complex_workspace_free ( w ); /* Transfer pairs of real, imaginary parts to double complex vector. */ r = ( double complex * ) malloc ( n * sizeof ( double complex ) ); for ( i = 0; i < n; i++ ) { r[i] = z[2*i] + z[2*i+1] * I; } free ( z ); return r; } /******************************************************************************/ void timestamp ( ) /******************************************************************************/ /* Purpose: timestamp() prints the current YMDHMS date as a time stamp. Example: 17 June 2014 09:45:54 AM Licensing: This code is distributed under the MIT license. Modified: 01 May 2021 Author: John Burkardt */ { # 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 ); printf ( "%s\n", time_buffer ); return; # undef TIME_SIZE }