# include # include # include # include "polynomial_conversion.h" /******************************************************************************/ double *bernstein_to_legendre01 ( int n, double bcoef[] ) /******************************************************************************/ /* Purpose: bernstein_to_legendre01() converts from Bernstein to Legendre01 form. Licensing: This code is distributed under the MIT license. Modified: 24 February 2024 Author: John Burkardt. Input: integer n: the order of the polynomial. real bcoef(1:n+1): the Bernstein coefficients of the polynomial. Output: real lcoef(1:n+1): the Legendre01 coefficients of the polynomial. */ { double *A; double *lcoef; if ( n < 0 ) { return NULL; } /* Get the matrix. */ A = bernstein_to_legendre01_matrix ( n ); /* Apply the transformation. */ lcoef = r8mat_mv_new ( n + 1, n + 1, A, bcoef ); free ( A ); return lcoef; } /******************************************************************************/ double *bernstein_to_legendre01_matrix ( int n ) /******************************************************************************/ /* Purpose: bernstein_to_legendre01_matrix() returns the Bernstein-to-Legendre01 matrix. Discussion: The Legendre polynomials are often defined on [-1,+1], while the Bernstein polynomials are defined on [0,1]. For this function, the Legendre polynomials have been shifted to share the [0,1] interval of definition. Example: bernstein_to_legendre01_matrix ( 5 ): 0.1667 0.1667 0.1667 0.1667 0.1667 0.1667 -0.3571 -0.2143 -0.0714 0.0714 0.2143 0.3571 0.2976 -0.0595 -0.2381 -0.2381 -0.0595 0.2976 -0.1389 0.1944 0.1111 -0.1111 -0.1944 0.1389 0.0357 -0.1071 0.0714 0.0714 -0.1071 0.0357 -0.0040 0.0198 -0.0397 0.0397 -0.0198 0.0040 Licensing: This code is distributed under the MIT license. Modified: 24 February 2024 Author: John Burkardt Input: int N, the maximum degree of the polynomials. Output: double A[(N+1)*(N+1)]: the matrix. */ { double *A; int i; int j; int k; A = ( double * ) malloc ( ( n + 1 ) * ( n + 1 ) * sizeof ( double ) ); for ( j = 0; j <= n; j++ ) { for ( i = 0; i <= n; i++ ) { A[i+j*(n+1)] = 0.0; } } for ( i = 0; i <= n; i++ ) { for ( j = 0; j <= n; j++ ) { for ( k = 0; k <= i; k++ ) { A[i+j*(n+1)] = A[i+j*(n+1)] + r8_mop ( i + k ) * pow ( r8_choose ( i, k ), 2 ) / r8_choose ( n + i, j + k ); } A[i+j*(n+1)] = A[i+j*(n+1)] * r8_choose ( n, j ) * ( double ) ( 2 * i + 1 ) / ( double ) ( n + i + 1 ); } } return A; } /******************************************************************************/ double *bernstein_to_monomial ( int n, double bcoef[] ) /******************************************************************************/ /* Purpose: bernstein_to_monomial() converts from Bernstein to monomial form. Licensing: This code is distributed under the MIT license. Modified: 20 February 2024 Author: John Burkardt. Input: integer n: the order of the polynomial. real bcoef(1:n+1): the Bernstein coefficients of the polynomial. Output: real mcoef(1:n+1): the monomial coefficients of the polynomial. */ { double *A; double *mcoef; if ( n < 0 ) { return NULL; } /* Get the matrix. */ A = bernstein_to_monomial_matrix ( n + 1 ); /* Apply the transformation. */ mcoef = r8mat_mv_new ( n + 1, n + 1, A, bcoef ); free ( A ); return mcoef; } /******************************************************************************/ double *bernstein_to_monomial_matrix ( int n ) /******************************************************************************/ /* Purpose: bernstein_to_monomial_matrix() converts Bernstein to monomial format. Discussion: The Bernstein matrix of order N is an NxN matrix A which can be used to transform a vector of power basis coefficients C representing a polynomial P(X) to a corresponding Bernstein basis coefficient vector B: B = A * C The N power basis vectors are ordered as (1,X,X^2,...X^(N-1)) and the N Bernstein basis vectors as ((1-X)^(N-1), X*(1_X)^(N-2),...,X^(N-1)). Example: N = 5 1 -4 6 -4 1 0 4 -12 12 -4 0 0 6 -12 6 0 0 0 4 -4 0 0 0 0 1 Licensing: This code is distributed under the MIT license. Modified: 10 July 2011 Author: John Burkardt Input: int N, the order of the matrix. Output: double A[N*N], the matrix. */ { double *A; int i; int j; 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; } } for ( j = 0; j < n; j++ ) { for ( i = 0; i <= j; i++ ) { A[i+j*n] = r8_mop ( j - i ) * r8_choose ( n - 1 - i, j - i ) * r8_choose ( n - 1, i ); } } return A; } /******************************************************************************/ double *chebyshev_to_monomial ( int n, double ccoef[] ) /******************************************************************************/ /* Purpose: chebyshev_to_monomial() converts from Chebyshev to monomial form. Licensing: This code is distributed under the MIT license. Modified: 18 February 2024 Author: Original Fortran77 version by Fred Krogh. This version by John Burkardt. Input: int n: the order of the polynomial. double ccoef(0:n): the Chebyshev coefficients of the polynomial. Output: double mcoef(0:n): the monomial coefficients of the polynomial. */ { int i; int j; double *mcoef; double tp; if ( n < 0 ) { return NULL; } mcoef = ( double * ) malloc ( ( n + 1 ) * sizeof ( double ) ); for ( i = 0; i <= n; i++ ) { mcoef[i] = ccoef[i]; } tp = 1.0; for ( j = 0; j <= n - 2; j++ ) { for ( i = n - 2; j <= i; i-- ) { mcoef[i] = mcoef[i] - mcoef[i+2]; } mcoef[j+1] = 0.5 * mcoef[j+1]; mcoef[j] = tp * mcoef[j]; tp = 2.0 * tp; } mcoef[n] = tp * mcoef[n]; if ( 0 < n ) { mcoef[n-1] = tp * mcoef[n-1]; } return mcoef; } /******************************************************************************/ double *chebyshev_to_monomial_matrix ( int n ) /******************************************************************************/ /* Purpose: chebyshev_to_monomial_matrix() converts Chebyshev T to monomial. 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 ( 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 = 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 ( int n, double alpha, double gcoef[] ) /******************************************************************************/ /* Purpose: gegenbauer_to_monomial() converts from Gegenbauer to monomial form. Licensing: This code is distributed under the MIT license. Modified: 25 April 2024 Author: John Burkardt. Input: integer n: the order of the polynomial. double alpha: the Gegenbauer parameter. real gcoef(1:n+1): the gegenbauer coefficients of the polynomial. Output: real mcoef(1:n+1): the monomial coefficients of the polynomial. */ { double *A; double *mcoef; if ( n < 0 ) { return NULL; } /* Get the matrix. */ A = gegenbauer_to_monomial_matrix ( n + 1, alpha ); /* Apply the transformation. */ mcoef = r8mat_mv_new ( n + 1, n + 1, A, gcoef ); free ( A ); return mcoef; } /******************************************************************************/ 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: 01 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 ) ); if ( n <= 0 ) { return A; } 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 * 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 ( int n, double hcoef[] ) /******************************************************************************/ /* Purpose: hermite_to_monomial() converts from Hermite to monomial form. Licensing: This code is distributed under the MIT license. Modified: 21 February 2024 Author: John Burkardt. Input: integer n: the order of the polynomial. real hcoef(1:n+1): the Hermite coefficients of the polynomial. Output: real mcoef(1:n+1): the monomial coefficients of the polynomial. */ { double *A; double *mcoef; if ( n < 0 ) { return NULL; } /* Get the matrix. */ A = hermite_to_monomial_matrix ( n + 1 ); /* Apply the transformation. */ mcoef = r8mat_mv_new ( n + 1, n + 1, A, hcoef ); free ( A ); return mcoef; } /******************************************************************************/ double *hermite_to_monomial_matrix ( int n ) /******************************************************************************/ /* Purpose: hermite_to_monomial_matrix() converts Hermite to monomial form. Example: 1 . -2 . 12 . -120 . 1680 . -30240 . 2 . -12 . 120 . -1680 . 30240 . . . 4 . -48 . 720 . -13440 . 302400 . . . 8 . -160 . 3360 . -80640 . . . . . 16 . -480 . 13440 . -403200 . . . . . 32 . -1344 . 48384 . . . . . . . 64 . -3584 . 161280 . . . . . . . 128 . -9216 . . . . . . . . . 256 . -23040 . . . . . . . . . 512 . . . . . . . . . . . 1024 Licensing: This code is distributed under the MIT license. Modified: 21 February 2024 Author: John Burkardt Input: int N, the order of the matrix. Output: double A[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 * A[i-1+(j-1)*n] - 2.0 * ( double ) ( j - 1 ) * A[i+(j-2)*n]; } } } return A; } /******************************************************************************/ 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 *laguerre_to_monomial ( int n, double lcoef[] ) /******************************************************************************/ /* Purpose: laguerre_to_monomial() converts from Laguerre to monomial form. Licensing: This code is distributed under the MIT license. Modified: 22 February 2024 Author: John Burkardt. Input: integer n: the order of the polynomial. real bcoef(1:n+1): the Laguerre coefficients of the polynomial. Output: real mcoef(1:n+1): the monomial coefficients of the polynomial. */ { double *A; double *mcoef; if ( n < 0 ) { return NULL; } /* Get the matrix. */ A = laguerre_to_monomial_matrix ( n + 1 ); /* Apply the transformation. */ mcoef = r8mat_mv_new ( n + 1, n + 1, A, lcoef ); free ( A ); return mcoef; } /******************************************************************************/ double *laguerre_to_monomial_matrix ( int n ) /******************************************************************************/ /* Purpose: laguerre_to_monomial_matrix() converts Laguerre to monomial form. 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 Licensing: This code is distributed under the MIT license. Modified: 22 February 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 A[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 ( int n, double lcoef[] ) /******************************************************************************/ /* Purpose: legendre_to_monomial() converts from Legendre to monomial form. Licensing: This code is distributed under the MIT license. Modified: 19 February 2024 Author: John Burkardt. Input: integer n: the order of the polynomial. real lcoef(1:n+1): the Legendre coefficients of the polynomial. Output: real mcoef(1:n+1): the monomial coefficients of the polynomial. */ { double *A; double *mcoef; if ( n < 0 ) { return NULL; } /* Get the matrix. */ A = legendre_to_monomial_matrix ( n + 1 ); /* Apply the transformation. */ mcoef = r8mat_mv_new ( n + 1, n + 1, A, lcoef ); free ( A ); return mcoef; } /******************************************************************************/ double *legendre_to_monomial_matrix ( int n ) /******************************************************************************/ /* Purpose: legendre_to_monomial_matrix() converts Legendre to monomial form. 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 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 A[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; } /******************************************************************************/ double *legendre01_to_bernstein ( int n, double lcoef[] ) /******************************************************************************/ /* Purpose: legendre01_to_bernsteinl() converts from Legendre01 to Bernstein form. Licensing: This code is distributed under the MIT license. Modified: 24 February 2024 Author: John Burkardt. Input: integer n: the order of the polynomial. real lcoef(1:n+1): the Legendre01 coefficients of the polynomial. Output: real bcoef(1:n+1): the Bernstein coefficients of the polynomial. */ { double *A; double *bcoef; if ( n < 0 ) { return NULL; } /* Get the matrix. */ A = legendre01_to_bernstein_matrix ( n ); /* Apply the transformation. */ bcoef = r8mat_mv_new ( n + 1, n + 1, A, lcoef ); free ( A ); return bcoef; } /******************************************************************************/ double *legendre01_to_bernstein_matrix ( int n ) /******************************************************************************/ /* Purpose: legendre01_to_bernstein_matrix() returns the Legendre-to-Bernstein matrix. Discussion: The Legendre polynomials are often defined on [-1,+1], while the Bernstein polynomials are defined on [0,1]. For this function, the Legendre polynomials have been shifted to share the [0,1] interval of definition. Example: legendre01_to_bernstein_matrix(5): 1.0000 -1.0000 1.0000 -1.0000 1.0000 -1.0000 1.0000 -0.6000 -0.2000 1.4000 -3.0000 5.0000 1.0000 -0.2000 -0.8000 0.8000 2.0000 -10.0000 1.0000 0.2000 -0.8000 -0.8000 2.0000 10.0000 1.0000 0.6000 -0.2000 -1.4000 -3.0000 -5.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 Licensing: This code is distributed under the MIT license. Modified: 24 February 2024 Author: John Burkardt Input: int N, the maximum degree of the polynomials. Output: double A[(N+1)*(N+1)], the matrix. */ { double *A; int i; int j; int k; A = ( double * ) malloc ( ( n + 1 ) * ( n + 1 ) * sizeof ( double ) ); for ( j = 0; j <= n; j++ ) { for ( i = 0; i <= n; i++ ) { A[i+j*(n+1)] = 0.0; } } for ( i = 0; i <= n; i++ ) { for ( j = 0; j <= n; j++ ) { for ( k = i4_max ( 0, i + j - n ); k <= i4_min ( i, j ); k++ ) { A[i+j*(n+1)] = A[i+j*(n+1)] + r8_mop ( j + k ) * pow ( r8_choose ( j, k ), 2 ) * r8_choose ( n - j, i - k ); } A[i+j*(n+1)] = A[i+j*(n+1)] / r8_choose ( n, i ); } } return A; } /******************************************************************************/ double *monomial_to_bernstein ( int n, double mcoef[] ) /******************************************************************************/ /* Purpose: monomial_to_bernstein() converts from monomial to Bernstein form. Licensing: This code is distributed under the MIT license. Modified: 20 February 2024 Author: John Burkardt. Input: integer n: the order of the polynomial. real mcoef(1:n+1): the monomial coefficients of the polynomial. Output: real bcoef(1:n+1): the Bernstein coefficients of the polynomial. */ { double *A; double *bcoef; if ( n < 0 ) { return NULL; } /* Get the matrix. */ A = monomial_to_bernstein_matrix ( n + 1 ); /* Apply the transformation. */ bcoef = r8mat_mv_new ( n + 1, n + 1, A, mcoef ); free ( A ); return bcoef; } /******************************************************************************/ double *monomial_to_bernstein_matrix ( int n ) /******************************************************************************/ /* Purpose: monomial_to_bernstein_matrix() converts monomial to Bernstein form. Discussion: The inverse Bernstein matrix of order N is an NxN matrix A which can be used to transform a vector of Bernstein basis coefficients B representing a polynomial P(X) to a corresponding power basis coefficient vector C: C = A * B The N power basis vectors are ordered as (1,X,X^2,...X^(N-1)) and the N Bernstein basis vectors as ((1-X)^(N-1), X*(1_X)^(N-2),...,X^(N-1)). Example: N = 5 1.0000 1.0000 1.0000 1.0000 1.0000 0 0.2500 0.5000 0.7500 1.0000 0 0 0.1667 0.5000 1.0000 0 0 0 0.2500 1.0000 0 0 0 0 1.0000 Licensing: This code is distributed under the MIT license. Modified: 10 July 2011 Author: John Burkardt Input: int N, the order of the matrix. Output: double A[N*N], the matrix. */ { double *A; int i; int j; 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; } } for ( j = 0; j < n; j++ ) { for ( i = 0; i <= j; i++ ) { A[i+j*n] = r8_choose ( j, i ) / r8_choose ( n - 1, i ); } } return A; } /******************************************************************************/ double *monomial_to_chebyshev ( int n, double mcoef[] ) /******************************************************************************/ /* Purpose: monomial_to_chebyshev() converts from monomial to Chebyshev form. Licensing: This code is distributed under the MIT license. Modified: 18 February 2024 Author: Original Fortran77 version by Fred Krogh. This version by John Burkardt. Input: int n: the order of the polynomial. double mcoef(0:n): the monomial coefficients of the polynomial. Output: double ccoef(0:n): the Chebyshev coefficients of the polynomial. */ { double *ccoef; int i; int j; double tp; if ( n < 0 ) { return NULL; } ccoef = ( double * ) malloc ( ( n + 1 ) * sizeof ( double ) ); for ( i = 0; i <= n; i++ ) { ccoef[i] = mcoef[i]; } tp = pow ( 0.5, n - 1 ); ccoef[n] = tp * ccoef[n]; if ( n == 0 ) { return ccoef; } ccoef[n-1] = tp * ccoef[n-1]; for ( j = n - 2; 0 <= j; j-- ) { tp = 2.0 * tp; ccoef[j] = tp * ccoef[j]; ccoef[j+1] = 2.0 * ccoef[j+1]; for ( i = j; i <= n - 2; i++ ) { ccoef[i] = ccoef[i] + ccoef[i+2]; } } return ccoef; } /******************************************************************************/ double *monomial_to_chebyshev_matrix ( int n ) /******************************************************************************/ /* Purpose: monomial_to_chebyshev_matrix() converts from monomial to Chebyshev T form. Example: N = 11 Each column must be divided by the divisor below it. 1 . 1 . 3 . 10 . 35 . 126 . 1 . 3 . 10 . 35 . 126 . . . 1 . 4 . 15 . 56 . 210 . . . 1 . 5 . 21 . 84 . . . . . 1 . 6 . 28 . 120 . . . . . 1 . 7 . 36 . . . . . . . 1 . 8 . 45 . . . . . . . 1 . 9 . . . . . . . . . 1 . 10 . . . . . . . . . 1 . . . . . . . . . . . 1 /1 /1 /2 /4 /8 /16 /32 /64 /128 /256 /512 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_INVERSE[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 = 2; j < n; j++ ) { for ( i = 0; i < n; i++ ) { if ( i == 0 ) { A[i+j*n] = A[i+1+(j-1)*n] / 2.0; } else if ( i == 1 ) { A[i+j*n] = ( 2.0 * A[i-1+(j-1)*n] + A[i+1+(j-1)*n] ) / 2.0; } else if ( i < n - 1 ) { A[i+j*n] = ( A[i-1+(j-1)*n] + A[i+1+(j-1)*n] ) / 2.0; } else { A[i+j*n] = A[i-1+(j-1)*n] / 2.0; } } } return A; } /******************************************************************************/ double *monomial_to_gegenbauer ( int n, double alpha, double mcoef[] ) /******************************************************************************/ /* Purpose: monomial_to_gegenbauer() converts from monomial to Gegenbauer form. Licensing: This code is distributed under the MIT license. Modified: 25 April 2024 Author: John Burkardt. Input: integer n: the order of the polynomial. real alpha: the parameter. real mcoef(1:n+1): the monomial coefficients of the polynomial. Output: real gcoef(1:n+1): the Gegenbauer coefficients of the polynomial. */ { double *A; double *gcoef; if ( n < 0 ) { return NULL; } /* Get the matrix. */ A = monomial_to_gegenbauer_matrix ( n + 1, alpha ); /* Appy the transformation. */ gcoef = r8mat_mv_new ( n + 1, n + 1, A, mcoef ); free ( A ); return gcoef; } /******************************************************************************/ double *monomial_to_gegenbauer_matrix ( int n, double alpha ) /******************************************************************************/ /* Purpose: monomial_to_gegenbauer_matrix(): monomial to Gegenbauer conversion matrix. Licensing: This code is distributed under the MIT license. Modified: 24 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 bot; int i; int ilo; int j; double top; 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; } for ( j = 0; j < n; j++ ) { ilo = ( j % 2 ); for ( i = ilo; i <= j; i = i + 2 ) { top = ( i + alpha ) * r8_factorial ( j ) * tgamma ( alpha ); bot = pow ( 2.0, j ) * r8_factorial ( ( j - i ) / 2 ) * tgamma ( ( j + i ) / 2 + alpha + 1.0 ); A[i+j*n] = top / bot; } } return A; } /******************************************************************************/ double *monomial_to_hermite ( int n, double mcoef[] ) /******************************************************************************/ /* Purpose: monomial_to_hermite() converts from monomial to Hermite form. Licensing: This code is distributed under the MIT license. Modified: 21 February 2024 Author: John Burkardt. Input: integer n: the order of the polynomial. real mcoef(1:n+1): the monomial coefficients of the polynomial. Output: real hcoef(1:n+1): the Hermite coefficients of the polynomial. */ { double *A; double *hcoef; if ( n < 0 ) { return NULL; } /* Get the matrix. */ A = monomial_to_hermite_matrix ( n + 1 ); /* Appy the transformation. */ hcoef = r8mat_mv_new ( n + 1, n + 1, A, mcoef ); free ( A ); return hcoef; } /******************************************************************************/ double *monomial_to_hermite_matrix ( int n ) /******************************************************************************/ /* Purpose: monomial_to_hermite_matrix() converts monomial to Hermite form. 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 Licensing: This code is distributed under the MIT license. Modified: 21 February 2024 Author: John Burkardt Input: int N, the order of the matrix. Output: double A[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] = 0.5; if ( n == 2 ) { return A; } for ( j = 2; j < n; j++ ) { for ( i = 0; i < n; i++ ) { if ( i == 0 ) { A[i+j*n] = ( ( double ) ( j - 1 ) * A[i+(j-2)*n] ) / 2.0; } else { A[i+j*n] = ( ( double ) ( j - 1 ) * A[i+(j-2)*n] + A[i-1+(j-1)*n] ) / 2.0; } } } return A; } /******************************************************************************/ double *monomial_to_laguerre ( int n, double mcoef[] ) /******************************************************************************/ /* Purpose: monomial_to_laguerre() converts from monomial to Laguerre form. Licensing: This code is distributed under the MIT license. Modified: 22 February 2024 Author: John Burkardt. Input: integer n: the order of the polynomial. real mcoef(1:n+1): the monomial coefficients of the polynomial. Output: real lcoef(1:n+1): the Laguerre coefficients of the polynomial. */ { double *A; double *lcoef; if ( n < 0 ) { return NULL; } /* Get the matrix. */ A = monomial_to_laguerre_matrix ( n + 1 ); /* Apply the transformation. */ lcoef = r8mat_mv_new ( n + 1, n + 1, A, mcoef ); free ( A ); return lcoef; } /******************************************************************************/ double *monomial_to_laguerre_matrix ( int n ) /******************************************************************************/ /* Purpose: monomial_to_laguerre_matrix() converts from monomial to Laguerre form. Example: N = 9 1 1 2 6 24 120 720 5040 40320 . -1 -4 -18 -96 -600 -4320 -35280 -322560 . . 2 18 144 1200 10800 105840 1128960 . . . -6 -96 -1200 -14400 -176400 -2257920 . . . . 24 600 10800 176400 2822400 . . . . . -120 -4320 -105840 -2257920 . . . . . . 720 35280 1128960 . . . . . . . -5040 -322560 . . . . . . . . 40320 Licensing: This code is distributed under the MIT license. Modified: 24 February 2024 Author: John Burkardt Input: int N, the order of the matrix. Output: double A[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++ ) { for ( i = 0; i < n; i++ ) { if ( i == 0 ) { A[i+j*n] = j * A[i+(j-1)*n]; } else { A[i+j*n] = j * ( A[i+(j-1)*n] - A[i-1+(j-1)*n] ); } } } return A; } /******************************************************************************/ double *monomial_to_legendre ( int n, double mcoef[] ) /******************************************************************************/ /* Purpose: monomial_to_legendre() converts from monomial to Legendre form. Licensing: This code is distributed under the MIT license. Modified: 19 February 2024 Author: John Burkardt. Input: integer n: the order of the polynomial. real mcoef(1:n+1): the monomial coefficients of the polynomial. Output: real lcoef(1:n+1): the Legendre coefficients of the polynomial. */ { double *A; double *lcoef; if ( n < 0 ) { return NULL; } /* Get the matrix. */ A = monomial_to_legendre_matrix ( n + 1 ); /* Apply the transformation. */ lcoef = r8mat_mv_new ( n + 1, n + 1, A, mcoef ); free ( A ); return lcoef; } /******************************************************************************/ double *monomial_to_legendre_matrix ( int n ) /******************************************************************************/ /* Purpose: monomial_to_legendre_matrix() converts from monomial to Legendre form. Example: N = 11 (each column must be divided by the underlying factor). 1 . 1 . 7 . 33 . 715 . 4199 . 1 . 3 . 27 . 143 . 3315 . . . 2 . 20 . 110 . 2600 .16150 . . . 2 . 28 . 182 . 4760 . . . . . 8 . 72 . 2160 .15504 . . . . . 8 . 88 . 2992 . . . . . . . 16 . 832 . 7904 . . . . . . . 16 . 960 . . . . . . . . . 128 . 2176 . . . . . . . . . 128 . . . . . . . . . . . 256 /1 /1 /3 /5 /35 /63 /231 /429 /6435/12155/46189 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 A[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 = 2; j <= n - 1; j++ ) { for ( i = 0; i <= n - 1; i++ ) { if ( i == 0 ) { A[i+j*n] = ( i + 1 ) * A[i+1+(j-1)*n] / ( 2 * i + 3 ); } else if ( i < n - 1 ) { A[i+j*n] = ( i ) * A[i-1+(j-1)*n] / ( 2 * i - 1 ) + ( i + 1 ) * A[i+1+(j-1)*n] / ( 2 * i + 3 ); } else { A[i+j*n] = ( i ) * A[i-1+(j-1)*n] / ( 2 * i - 1 ); } } } return A; } /******************************************************************************/ double r8_choose ( int n, int k ) /******************************************************************************/ /* Purpose: r8_choose() computes the binomial coefficient C(N,K) as an R8. Discussion: The value is calculated in such a way as to avoid overflow and roundoff. The calculation is done in R8 arithmetic. The formula used is: C(N,K) = N / ( K * (N-K) ) Licensing: This code is distributed under the MIT license. Modified: 01 July 2008 Author: John Burkardt Reference: ML Wolfson, HV Wright, Algorithm 160: Combinatorial of M Things Taken N at a Time, Communications of the ACM, Volume 6, Number 4, April 1963, page 161. Input: int N, K, the values of N and K. Output: double R8_CHOOSE, the number of combinations of N things taken K at a time. */ { int i; int mn; int mx; int value; mn = i4_min ( k, n - k ); if ( mn < 0 ) { value = 0.0; } else if ( mn == 0 ) { value = 1.0; } else { mx = i4_max ( k, n - k ); value = ( double ) ( mx + 1 ); for ( i = 2; i <= mn; i++ ) { value = ( value * ( double ) ( mx + i ) ) / ( double ) i; } } return value; } /******************************************************************************/ double r8_factorial ( int n ) /******************************************************************************/ /* Purpose: r8_factorial() computes the factorial of N. Discussion: factorial ( N ) = product ( 1 <= I <= N ) I Licensing: This code is distributed under the MIT license. Modified: 26 June 2008 Author: John Burkardt Input: int N: the argument of the factorial function. If N is less than 1, the function value is returned as 1. Output: double R8_FACTORIAL: the factorial of N. */ { int i; double value; value = 1.0; for ( i = 1; i <= n; i++ ) { value = value * ( double ) ( i ); } return value; } /******************************************************************************/ double r8_mop ( int i ) /******************************************************************************/ /* Purpose: r8_mop() returns the I-th power of -1 as an R8 value. Discussion: An R8 is an double value. Licensing: This code is distributed under the MIT license. Modified: 01 July 2008 Author: John Burkardt Input: int I, the power of -1. Output: double R8_MOP, the I-th power of -1. */ { double value; if ( ( i % 2 ) == 0 ) { value = 1.0; } else { value = -1.0; } return value; } /******************************************************************************/ 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; }