C 代码 计算伯恩斯坦多项式.rar

# include <math.h> # include <stdio.h> # include <stdlib.h> # include <string.h> # include <time.h> # include "bernstein_polynomial.h" /******************************************************************************/ double *bernstein_matrix ( int n ) /******************************************************************************/ /* Purpose: BERNSTEIN_MATRIX returns the Bernstein matrix. 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)). For N = 5, the matrix has the form: 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 and the numbers in each column represent the coefficients in the power series expansion of a Bernstein polynomial, so that B(5,4) = - 4 x^4 + 12 x^3 - 12 x^2 + 4 x Licensing: This code is distributed under the GNU LGPL license. Modified: 11 February 2012 Author: John Burkardt Parameters: Input, int N, the order of the matrix. Output, double BERNSTEIN_MATRIX[N*N], the Bernstein matrix. */ { double *a; int i; int j; a = ( double * ) malloc ( n * n * sizeof ( double ) ); 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 ); } for ( i = j + 1; i < n; i++ ) { a[i+j*n] = 0.0; } } return a; } /******************************************************************************/ double bernstein_matrix_determinant ( int n ) /******************************************************************************/ /* Purpose: BERNSTEIN_MATRIX_DETERMINANT returns the determinant of the BERNSTEIN matrix. Licensing: This code is distributed under the GNU LGPL license. Modified: 14 March 2015 Author: John Burkardt Parameters: Input, int N, the order of the matrix. Output, double BERNSTEIN_MATRIX_DETERMINANT, the determinant. */ { int i; double value; value = 1.0; for ( i = 0; i < n; i++ ) { value = value * r8_choose ( n - 1, i ); } return value; } /******************************************************************************/ double *bernstein_matrix_inverse ( int n ) /******************************************************************************/ /* Purpose: BERNSTEIN_MATRIX_INVERSE returns the inverse Bernstein matrix. 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)). For N = 5, the matrix has the form: 1 1 1 1 1 0 1/4 1/2 3/4 1 0 0 1/6 1/2 1 0 0 0 1/4 1 0 0 0 0 1 Licensing: This code is distributed under the GNU LGPL license. Modified: 11 February 2012 Author: John Burkardt Parameters: Input, int N, the order of the matrix. Output, double BERNSTEIN_MATRIX_INVERSE[N*N], the inverse Bernstein matrix. */ { double *a; int i; int j; a = ( double * ) malloc ( n * n * sizeof ( double ) ); 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 ); } for ( i = j + 1; i < n; i++ ) { a[i+j*n] = 0.0; } } return a; } /******************************************************************************/ double *bernstein_poly_01 ( int n, double x ) /******************************************************************************/ /* Purpose: BERNSTEIN_POLY_01 evaluates the Bernstein polynomials based in [0,1]. Discussion: The Bernstein polynomials are assumed to be based on [0,1]. The formula is: B(N,I)(X) = [N!/(I!*(N-I)!)] * (1-X)^(N-I) * X^I First values: B(0,0)(X) = 1 B(1,0)(X) = 1-X B(1,1)(X) = X B(2,0)(X) = (1-X)^2 B(2,1)(X) = 2 * (1-X) * X B(2,2)(X) = X^2 B(3,0)(X) = (1-X)^3 B(3,1)(X) = 3 * (1-X)^2 * X B(3,2)(X) = 3 * (1-X) * X^2 B(3,3)(X) = X^3 B(4,0)(X) = (1-X)^4 B(4,1)(X) = 4 * (1-X)^3 * X B(4,2)(X) = 6 * (1-X)^2 * X^2 B(4,3)(X) = 4 * (1-X) * X^3 B(4,4)(X) = X^4 Special values: B(N,I)(X) has a unique maximum value at X = I/N. B(N,I)(X) has an I-fold zero at 0 and and N-I fold zero at 1. B(N,I)(1/2) = C(N,K) / 2^N For a fixed X and N, the polynomials add up to 1: Sum ( 0 <= I <= N ) B(N,I)(X) = 1 Licensing: This code is distributed under the GNU LGPL license. Modified: 11 February 2012 Author: John Burkardt Parameters: Input, int N, the degree of the Bernstein polynomials to be used. For any N, there is a set of N+1 Bernstein polynomials, each of degree N, which form a basis for polynomials on [0,1]. Input, double X, the evaluation point. Output, double BERNSTEIN_POLY_01[N+1], the values of the N+1 Bernstein polynomials at X. */ { double *bern; int i; int j; bern = ( double * ) malloc ( ( n + 1 ) * sizeof ( double ) ); if ( n == 0 ) { bern[0] = 1.0; } else if ( 0 < n ) { bern[0] = 1.0 - x; bern[1] = x; for ( i = 2; i <= n; i++ ) { bern[i] = x * bern[i-1]; for ( j = i - 1; 1 <= j; j-- ) { bern[j] = x * bern[j-1] + ( 1.0 - x ) * bern[j]; } bern[0] = ( 1.0 - x ) * bern[0]; } } return bern; } /******************************************************************************/ double *bernstein_poly_01_matrix ( int m, int n, double x[] ) /******************************************************************************/ /* Purpose: BERNSTEIN_POLY_01_MATRIX evaluates the Bernstein polynomials based in [0,1]. Discussion: The Bernstein polynomials are assumed to be based on [0,1]. The formula is: B(N,I)(X) = [N!/(I!*(N-I)!)] * (1-X)^(N-I) * X^I First values: B(0,0)(X) = 1 B(1,0)(X) = 1-X B(1,1)(X) = X B(2,0)(X) = (1-X)^2 B(2,1)(X) = 2 * (1-X) * X B(2,2)(X) = X^2 B(3,0)(X) = (1-X)^3 B(3,1)(X) = 3 * (1-X)^2 * X B(3,2)(X) = 3 * (1-X) * X^2 B(3,3)(X) = X^3 B(4,0)(X) = (1-X)^4 B(4,1)(X) = 4 * (1-X)^3 * X B(4,2)(X) = 6 * (1-X)^2 * X^2 B(4,3)(X) = 4 * (1-X) * X^3 B(4,4)(X) = X^4 Special values: B(N,I)(X) has a unique maximum value at X = I/N. B(N,I)(X) has an I-fold zero at 0 and and N-I fold zero at 1. B(N,I)(1/2) = C(N,K) / 2^N For a fixed X and N, the polynomials add up to 1: Sum ( 0 <= I <= N ) B(N,I)(X) = 1 Licensing: This code is distributed under the GNU LGPL license. Modified: 27 January 2016 Author: John Burkardt Parameters: Input, int M, the number of evaluation points. Input, int N, the degree of the Bernstein polynomials to be used. For any N, there is a set of N+1 Bernstein polynomials, each of degree N, which form a basis for polynomials on [0,1]. Input, double X[M], the evaluation points. Output, double BERNSTEIN_POLY_01_MATRIX[M*(N+1)], the values of the N+1 Bernstein polynomials at the evaluation points. */ { double *b; int i; int j; int k; b =