# include # include # include # include # include # include # include "test_interp_nd.h" # include "r8lib.h" /******************************************************************************/ double csevl ( double x, double a[], int n ) /******************************************************************************/ /* Purpose: CSEVL evaluates a Chebyshev series. Licensing: This code is distributed under the MIT license. Modified: 17 January 2012 Author: C version by John Burkardt. Reference: Roger Broucke, Algorithm 446: Ten Subroutines for the Manipulation of Chebyshev Series, Communications of the ACM, Volume 16, Number 4, April 1973, pages 254-256. Parameters: Input, double X, the evaluation point. Input, double CS[N], the Chebyshev coefficients. Input, int N, the number of Chebyshev coefficients. Output, double CSEVL, the Chebyshev series evaluated at X. */ { double b0; double b1; double b2; int i; double twox; double value; if ( n < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "CSEVL - Fatal error!\n" ); fprintf ( stderr, " Number of terms <= 0.\n" ); exit ( 1 ); } if ( 1000 < n ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "CSEVL - Fatal error!\n" ); fprintf ( stderr, " Number of terms greater than 1000.\n" ); exit ( 1 ); } if ( x < -1.1 || 1.1 < x ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "CSEVL - Fatal error!\n" ); fprintf ( stderr, " X outside (-1,+1).\n" ); exit ( 1 ); } twox = 2.0 * x; b1 = 0.0; b0 = 0.0; for ( i = n - 1; 0 <= i; i-- ) { b2 = b1; b1 = b0; b0 = twox * b1 - b2 + a[i]; } value = 0.5 * ( b0 - b2 ); return value; } /******************************************************************************/ int i4vec_sum ( int n, int a[] ) /******************************************************************************/ /* Purpose: I4VEC_SUM sums the entries of an I4VEC. Discussion: An I4VEC is a vector of I4's. Example: Input: A = ( 1, 2, 3, 4 ) Output: I4VEC_SUM = 10 Licensing: This code is distributed under the MIT license. Modified: 29 May 2003 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Input, int A[N], the vector to be summed. Output, int I4VEC_SUM, the sum of the entries of A. */ { int i; int sum; sum = 0; for ( i = 0; i < n; i++ ) { sum = sum + a[i]; } return sum; } /******************************************************************************/ int inits ( double dos[], int nos, double eta ) /******************************************************************************/ /* Purpose: INITS initializes a Chebyshev series. Licensing: This code is distributed under the MIT license. Modified: 17 January 2012 Author: C version by John Burkardt. Reference: Roger Broucke, Algorithm 446: Ten Subroutines for the Manipulation of Chebyshev Series, Communications of the ACM, Volume 16, Number 4, April 1973, pages 254-256. Parameters: Input, double DOS[NOS], the Chebyshev coefficients. Input, int NOS, the number of coefficients. Input, double ETA, the desired accuracy. Output, int INITS, the number of terms of the series needed to ensure the requested accuracy. */ { double err; int i; int value; if ( nos < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "INITS - Fatal error!\n" ); fprintf ( stderr, " Number of coefficients < 1.\n" ); exit ( 1 ); } err = 0.0; for ( i = nos - 1; 0 <= i; i-- ) { err = err + fabs ( dos[i] ); if ( eta < err ) { value = i + 1; return value; } } value = i; fprintf ( stderr, "\n" ); fprintf ( stderr, "INITS - Warning!\n" ); fprintf ( stderr, " ETA may be too small.\n" ); return value; } /******************************************************************************/ double *p00_c ( int prob, int m, int *seed ) /******************************************************************************/ /* Purpose: P00_CW computes a random C parameter vector for any problem. Licensing: This code is distributed under the MIT license. Modified: 27 August 2012 Author: John Burkardt Parameters: Input, int PROB, the problem number. Input, int M, the spatial dimension. Input/output, int *SEED, a seed for the random number generator. Output, double P00_C[M], the parameter vector. */ { double b[6] = { 1.5, 0.0, 1.85, 7.03, 20.4, 4.3 }; double *c; double c_sum; int i; b[2-1] = ( double ) ( m ); c = r8vec_uniform_01_new ( m, seed ); c_sum = r8vec_sum ( m, c ); for ( i = 0; i < m; i++ ) { c[i] = b[prob-1] * c[i] / c_sum; } return c; } /******************************************************************************/ double *p00_d ( int prob, int m, int id, double c[], double w[], int n, double x[] ) /******************************************************************************/ /* Purpose: P00_D returns a derivative component of any function. Licensing: This code is distributed under the MIT license. Modified: 27 August 2012 Author: John Burkardt Parameters: Input, int PROB, the index of the function. Input, int M, the spatial dimension. Input, int ID, the spatial coordinate to differentiate. Input, double C[M], W[M], the problem parameters. Input, int N, the number of evaluation points. Input, double X[M*N], the evalution points. Output, double P00_D[N], the ID-th derivative component. */ { double *d; if ( id < 0 || m < id ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "P00_D - Fatal error!\n" ); fprintf ( stderr, " Illegal spatial coordinate ID = %d\n", id ); exit ( 1 ); } if ( prob == 1 ) { d = p01_d ( m, id, c, w, n, x ); } else if ( prob == 2 ) { d = p02_d ( m, id, c, w, n, x ); } else if ( prob == 3 ) { d = p03_d ( m, id, c, w, n, x ); } else if ( prob == 4 ) { d = p04_d ( m, id, c, w, n, x ); } else if ( prob == 5 ) { d = p05_d ( m, id, c, w, n, x ); } else if ( prob == 6 ) { d = p06_d ( m, id, c, w, n, x ); } else { fprintf ( stderr, "\n" ); fprintf ( stderr, "P00_D - Fatal error!\n" ); fprintf ( stderr, " Illegal function index PROB = %d\n", prob ); exit ( 1 ); } return d; } /******************************************************************************/ double *p00_f ( int prob, int m, double c[], double w[], int n, double x[] ) /******************************************************************************/ /* Purpose: P00_F returns the value of any function. Licensing: This code is distributed under the MIT license. Modified: 27 August 2012 Author: John Burkardt Parameters: Input, int PROB, the index of the function. Input, int M, the spatial dimension. Input, double C[M], W[M], the problem parameters. Input, int N, the number of evaluation points. Input, double X[M*N], the evalution points. Output, double P00_F[N], the function values. */ { double *f; if ( prob == 1 ) { f = p01_f ( m, c, w, n, x ); } else if ( prob == 2 ) { f = p02_f ( m, c, w, n, x ); } else if ( prob == 3 ) { f = p03_f ( m, c, w, n, x ); } else if ( prob == 4 ) { f = p04_f ( m, c, w, n, x ); } else if ( prob == 5 ) { f = p05_f ( m, c, w, n, x ); } else if ( prob == 6 ) { f = p06_f ( m, c, w, n, x ); } else { fprintf ( stderr, "\n" ); fprintf ( stderr, "P00_F - Fatal error!\n" ); fprintf ( stderr, " Illegal function index PROB = %d\n", prob ); exit ( 1 ); } return f; } /******************************************************************************/ int p00_prob_num ( ) /******************************************************************************/ /* Purpose: P00_PROB_NUM returns the number of test functions available. Licensing: This code is distributed under the MIT license. Modified: 27 August 2012 Author: John Burkardt Parameters: Output, int P00_PROB_NUM, the number of test functions. */ { int prob_num; prob_num = 6; return prob_num; } /******************************************************************************/ double p00_q ( int prob, int m, double c[], double w[] ) /******************************************************************************/ /* Purpose: P00_Q returns the integral of any function. Licensing: This code is distributed under the MIT license. Modified: 29 August 2012 Author: John Burkardt Parameters: Input, int PROB, the index of the function. Input, int M, the spatial dimension. Input, double C[M], W[M], the problem parameters. Output, double P00_Q, the integral. */ { double q; if ( prob == 1 ) { q = p01_q ( m, c, w ); } else if ( prob == 2 ) { q = p02_q ( m, c, w ); } else if ( prob == 3 ) { q = p03_q ( m, c, w ); } else if ( prob == 4 ) { q = p04_q ( m, c, w ); } else if ( prob == 5 ) { q = p05_q ( m, c, w ); } else if ( prob == 6 ) { q = p06_q ( m, c, w ); } else { fprintf ( stderr, "\n" ); fprintf ( stderr, "P00_Q - Fatal error!\n" ); fprintf ( stderr, " Illegal function index PROB = %d\n", prob ); exit ( 1 ); } return q; } /******************************************************************************/ char *p00_title ( int prob ) /******************************************************************************/ /* Purpose: P00_TITLE returns the title for any function. Licensing: This code is distributed under the MIT license. Modified: 27 August 2012 Author: John Burkardt Parameters: Input, int PROB, the index of the function. Output, char *P00_TITLE, the function title. */ { char *title; if ( prob == 1 ) { title = p01_title ( ); } else if ( prob == 2 ) { title = p02_title ( ); } else if ( prob == 3 ) { title = p03_title ( ); } else if ( prob == 4 ) { title = p04_title ( ); } else if ( prob == 5 ) { title = p05_title ( ); } else if ( prob == 6 ) { title = p06_title ( ); } else { fprintf ( stderr, "\n" ); fprintf ( stderr, "P00_TITLE - Fatal error!\n" ); fprintf ( stderr, " Illegal function index PROB = %d\n", prob ); exit ( 1 ); } return title; } /******************************************************************************/ double *p00_w ( int prob, int m, int *seed ) /******************************************************************************/ /* Purpose: P00_W computes a random W parameter vector for any problem. Licensing: This code is distributed under the MIT license. Modified: 27 August 2012 Author: John Burkardt Parameters: Input, int PROB, the problem number. Input, int M, the spatial dimension. Input/output, int *SEED, a seed for the random number generator. Output, double P00_W[M], the parameter vector. */ { double *w; w = r8vec_uniform_01_new ( m, seed ); return w; } /******************************************************************************/ double *p01_d ( int m, int id, double c[], double w[], int n, double x[] ) /******************************************************************************/ /* Purpose: P01_D evaluates any derivative component for problem p01. Discussion: f(x) = cos ( 2 * pi * w(1) + sum ( c(1:m) * x(1:m) ) ) Default values are: c(1:m) = 1/m w(1) = 0.3 Licensing: This code is distributed under the MIT license. Modified: 27 August 2012 Author: John Burkardt Reference: Alan Genz, A Package for Testing Multiple Integration Subroutines, in Numerical Integration: Recent Developments, Software and Applications, edited by Patrick Keast and Graeme Fairweather, Reidel, 1987, pages 337-340, ISBN: 9027725144, LC: QA299.3.N38. Parameters: Input, int M, the dimension of the argument. Input, int ID, the spatial coordinate to differentiate. Input, double C[M], W[M], the problem parameters. Input, int N, the number of points. Input, double X[M*N], the evaluation points. Output, double P01_D[N], the ID-th derivative component. */ { double *d; int i; int j; double pi = 3.141592653589793; d = ( double * ) malloc ( n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { d[j] = 2.0 * pi * w[0]; } for ( i = 0; i < m; i++ ) { for ( j = 0; j < n; j++ ) { d[j] = d[j] + c[i] * x[i+j*m]; } } for ( j = 0; j < n; j++ ) { d[j] = - c[id] * sin ( d[j] ); } return d; } /******************************************************************************/ double *p01_f ( int m, double c[], double w[], int n, double x[] ) /******************************************************************************/ /* Purpose: P01_F evaluates the function for problem p01. Discussion: f(x) = cos ( 2 * pi * w(1) + sum ( c(1:m) * x(1:m) ) ) Default values are: c(1:m) = 1/m w(1) = 0.3 Licensing: This code is distributed under the MIT license. Modified: 27 August 2012 Author: John Burkardt Reference: Alan Genz, A Package for Testing Multiple Integration Subroutines, in Numerical Integration: Recent Developments, Software and Applications, edited by Patrick Keast and Graeme Fairweather, Reidel, 1987, pages 337-340, ISBN: 9027725144, LC: QA299.3.N38. Parameters: Input, int M, the dimension of the argument. Input, double C[M], W[M], the problem parameters. Input, int N, the number of points. Input, double X[M*N], the evaluation points. Output, double P01_F[N], the function values. */ { double *f; int i; int j; double pi = 3.141592653589793; f = ( double * ) malloc ( n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { f[j] = 2.0 * pi * w[0]; } for ( i = 0; i < m; i++ ) { for ( j = 0; j < n; j++ ) { f[j] = f[j] + c[i] * x[i+j*m]; } } for ( j = 0; j < n; j++ ) { f[j] = cos ( f[j] ); } return f; } /******************************************************************************/ double p01_q ( int m, double c[], double w[] ) /******************************************************************************/ /* Purpose: P01_Q evaluates the integral for problem p01. Licensing: This code is distributed under the MIT license. Modified: 29 August 2012 Author: John Burkardt Reference: Alan Genz, A Package for Testing Multiple Integration Subroutines, in Numerical Integration: Recent Developments, Software and Applications, edited by Patrick Keast and Graeme Fairweather, Reidel, 1987, pages 337-340, ISBN: 9027725144, LC: QA299.3.N38. Parameters: Input, int M, the dimension of the argument. Input, double C[M], W[M], the problem parameters. Output, double P01_Q, the integral. */ { double c_prod; double c_sum; int i; double pi = 3.141592653589793; double q; c_sum = r8vec_sum ( m, c ); c_prod = 1.0; for ( i = 0; i < m; i++ ) { c_prod = c_prod * sin ( 0.5 * c[i] ) / c[i]; } q = pow ( 2.0, m ) * cos ( 2.0 * pi * w[0] + 0.5 * c_sum ) * c_prod; return q; } /******************************************************************************/ char *p01_title ( ) /******************************************************************************/ /* Purpose: P01_TITLE returns the name of problem p01. Licensing: This code is distributed under the MIT license. Modified: 27 August 2012 Author: John Burkardt Parameters: Output, char *P01_TITLE, the title of the problem. */ { char *title; title = ( char * ) malloc ( 80 * sizeof ( char ) ); strcpy ( title, "Oscillatory" ); return title; } /******************************************************************************/ double *p02_d ( int m, int id, double c[], double w[], int n, double x[] ) /******************************************************************************/ /* Purpose: P02_D evaluates an derivative component for problem p02. Discussion: f(x) = 1 / product ( c(1:m)^(-2) + ( x(1:m) - w(1:m) )^2 ) Default values are: c(1:m) = 1 w(1:m) = 0.5 Licensing: This code is distributed under the MIT license. Modified: 27 August 2012 Author: John Burkardt Reference: Alan Genz, A Package for Testing Multiple Integration Subroutines, in Numerical Integration: Recent Developments, Software and Applications, edited by Patrick Keast and Graeme Fairweather, Reidel, 1987, pages 337-340, ISBN: 9027725144, LC: QA299.3.N38. Parameters: Input, int M, the dimension of the argument. Input, int ID, the spatial coordinate to differentiate. Input, double C[M], W[M], the problem parameters. Input, int N, the number of points. Input, double X[M*N], the evaluation points. Output, double P02_D[N], the ID-th derivative component. */ { double *d; int i; int j; d = ( double * ) malloc ( n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { d[j] = 1.0; } for ( i = 0; i < m; i++ ) { for ( j = 0; j < n; j++ ) { d[j] = d[j] * ( pow ( c[i], - 2 ) + pow ( x[i+j*m] - w[i], 2 ) ); } } for ( j = 0; j < n; j++ ) { d[j] = - 2.0 / d[j] * ( x[id+j*m] - w[id] ) / ( pow ( c[id], - 2 ) + pow ( x[id+j*m] - w[id], 2 ) ); } return d; } /******************************************************************************/ double *p02_f ( int m, double c[], double w[], int n, double x[] ) /******************************************************************************/ /* Purpose: P02_F evaluates the function for problem p02. Discussion: f(x) = 1 / product ( c(1:m)^(-2) + ( x(1:m) - w(1:m) )^2 ) Default values are: c(1:m) = 1 w(1:m) = 0.5 Licensing: This code is distributed under the MIT license. Modified: 27 August 2012 Author: John Burkardt Reference: Alan Genz, A Package for Testing Multiple Integration Subroutines, in Numerical Integration: Recent Developments, Software and Applications, edited by Patrick Keast and Graeme Fairweather, Reidel, 1987, pages 337-340, ISBN: 9027725144, LC: QA299.3.N38. Parameters: Input, int M, the dimension of the argument. Input, double C[M], W[M], the problem parameters. Input, int N, the number of points. Input, double X[M*N], the evaluation points. Output, double P02_F[N], the function values. */ { double *f; int i; int j; f = ( double * ) malloc ( n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { f[j] = 1.0; } for ( i = 0; i < m; i++ ) { for ( j = 0; j < n; j++ ) { f[j] = f[j] * ( pow ( c[i], - 2 ) + pow ( x[i+j*m] - w[i], 2 ) ); } } for ( j = 0; j < n; j++ ) { f[j] = 1.0 / f[j]; } return f; } /******************************************************************************/ double p02_q ( int m, double c[], double w[] ) /******************************************************************************/ /* Purpose: P02_Q evaluates the integral for problem p02. Licensing: This code is distributed under the MIT license. Modified: 29 August 2012 Author: John Burkardt Reference: Alan Genz, A Package for Testing Multiple Integration Subroutines, in Numerical Integration: Recent Developments, Software and Applications, edited by Patrick Keast and Graeme Fairweather, Reidel, 1987, pages 337-340, ISBN: 9027725144, LC: QA299.3.N38. Parameters: Input, int M, the dimension of the argument. Input, double C[M], W[M], the problem parameters. Output, double P02_Q, the integral. */ { int i; double q; q = 1.0; for ( i = 0; i < m; i++ ) { q = q * ( atan ( ( 1.0 - w[i] ) * c[i] ) + atan ( w[i] * c[i] ) ) * c[i]; } return q; } /******************************************************************************/ char *p02_title ( ) /******************************************************************************/ /* Purpose: P02_TITLE returns the title of problem p02. Licensing: This code is distributed under the MIT license. Modified: 27 August 2012 Author: John Burkardt Parameters: Output, char *P02_TITLE, the title of the problem. */ { char *title; title = ( char * ) malloc ( 80 * sizeof ( char ) ); strcpy ( title, "Product Peak" ); return title; } /******************************************************************************/ double *p03_d ( int m, int id, double c[], double w[], int n, double x[] ) /******************************************************************************/ /* Purpose: P03_D evaluates any derivative component for problem p03. Discussion: f(x) = 1 / ( 1 + sum ( c(1:m) * x(1:m) ) ) ^ ( m + 1 ) Default values are: c(1:m) = 1/m Licensing: This code is distributed under the MIT license. Modified: 27 August 2012 Author: John Burkardt Reference: Alan Genz, A Package for Testing Multiple Integration Subroutines, in Numerical Integration: Recent Developments, Software and Applications, edited by Patrick Keast and Graeme Fairweather, Reidel, 1987, pages 337-340, ISBN: 9027725144, LC: QA299.3.N38. Parameters: Input, int M, the dimension of the argument. Input, int ID, the spatial coordinate to differentiate. Input, double C[M], W[M], the problem parameters. Input, int N, the number of points. Input, double X[M*N], the evaluation points. Output, double P03_D[N], the ID-th derivative component. */ { double *d; int i; int j; d = ( double * ) malloc ( n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { d[j] = 1.0; } for ( i = 0; i < m; i++ ) { for ( j = 0; j < n; j++ ) { d[j] = d[j] + c[i] * x[i+j*m]; } } for ( j = 0; j < n; j++ ) { d[j] = - c[id] * ( double ) ( m + 1 ) / pow ( d[j], m + 2 ); } return d; } /******************************************************************************/ double *p03_f ( int m, double c[], double w[], int n, double x[] ) /******************************************************************************/ /* Purpose: P03_F evaluates the function for problem p03. Discussion: f(x) = 1 / ( 1 + sum ( c(1:m) * x(1:m) ) ) ^ ( m + 1 ) Default values are: c(1:m) = 1/m Licensing: This code is distributed under the MIT license. Modified: 27 August 2012 Author: John Burkardt Reference: Alan Genz, A Package for Testing Multiple Integration Subroutines, in Numerical Integration: Recent Developments, Software and Applications, edited by Patrick Keast and Graeme Fairweather, Reidel, 1987, pages 337-340, ISBN: 9027725144, LC: QA299.3.N38. Parameters: Input, int M, the dimension of the argument. Input, double C[M], W[M], the problem parameters. Input, int N, the number of points. Input, double X[M*N], the evaluation points. Output, double P03_F[N], the function values. */ { double *f; int i; int j; f = ( double * ) malloc ( n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { f[j] = 1.0; } for ( i = 0; i < m; i++ ) { for ( j = 0; j < n; j++ ) { f[j] = f[j] + c[i] * x[i+j*m]; } } for ( j = 0; j < n; j++ ) { f[j] = 1.0 / pow ( f[j], m + 1 ); } return f; } /******************************************************************************/ double p03_q ( int m, double c[], double w[] ) /******************************************************************************/ /* Purpose: P03_Q evaluates the integral for problem p03. Licensing: This code is distributed under the MIT license. Modified: 29 August 2012 Author: John Burkardt Reference: Alan Genz, A Package for Testing Multiple Integration Subroutines, in Numerical Integration: Recent Developments, Software and Applications, edited by Patrick Keast and Graeme Fairweather, Reidel, 1987, pages 337-340, ISBN: 9027725144, LC: QA299.3.N38. Parameters: Input, int M, the dimension of the argument. Input, double C[M], W[M], the problem parameters. Output, double P03_Q, the integral. */ { int *ivec; double q; int rank; int s; /* Here, we need to generate all possible DIM_NUM tuples with values of 0 or 1. */ ivec = ( int * ) malloc ( m * sizeof ( int ) ); q = 0.0; rank = 0; while ( 1 ) { tuple_next ( 0, 1, m, &rank, ivec ); if ( rank == 0 ) { break; } s = i4vec_sum ( m, ivec ); q = q + r8_mop ( s ) / ( 1.0 + r8vec_i4vec_dot_product ( m, c, ivec ) ); } q = q / ( r8_factorial ( m ) * r8vec_product ( m, c ) ); free ( ivec ); return q; } /******************************************************************************/ char *p03_title ( ) /******************************************************************************/ /* Purpose: P03_TITLE returns the title of problem p03. Licensing: This code is distributed under the MIT license. Modified: 27 August 2012 Author: John Burkardt Parameters: Output, char *P03_TITLE, the title of the problem. */ { char *title; title = ( char * ) malloc ( 80 * sizeof ( char ) ); strcpy ( title, "Corner Peak" ); return title; } /******************************************************************************/ double *p04_d ( int m, int id, double c[], double w[], int n, double x[] ) /******************************************************************************/ /* Purpose: P04_D evaluates any derivative component for problem p04. Discussion: f(x) = exp ( - sum ( c(1:m)^2 * ( x(1:m) - w(1:m) )^2 ) Default values are: c(1:m) = 1 / m w(1:m) = 0.5 Licensing: This code is distributed under the MIT license. Modified: 27 August 2012 Author: John Burkardt Reference: Alan Genz, A Package for Testing Multiple Integration Subroutines, in Numerical Integration: Recent Developments, Software and Applications, edited by Patrick Keast and Graeme Fairweather, Reidel, 1987, pages 337-340, ISBN: 9027725144, LC: QA299.3.N38. Parameters: Input, int M, the dimension of the argument. Input, int ID, the spatial coordinate to differentiate. Input, double C[M], W[M], the problem parameters. Input, int N, the number of points. Input, double X[M*N], the evaluation points. Output, double P04_D[N], the ID-th derivative component. */ { double *d; int i; int j; d = ( double * ) malloc ( n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { d[j] = 0.0; } for ( i = 0; i < m; i++ ) { for ( j = 0; j < n; j++ ) { d[j] = d[j] + pow ( c[i] * ( x[i+j*m] - w[i] ), 2 ); } } for ( j = 0; j < n; j++ ) { d[j] = exp ( - d[j] ) * pow ( c[id], 2 ) * ( - 2.0 ) * ( x[id+j*m] - w[id] ); } return d; } /******************************************************************************/ double *p04_f ( int m, double c[], double w[], int n, double x[] ) /******************************************************************************/ /* Purpose: P04_F evaluates the function for problem p04. Discussion: f(x) = exp ( - sum ( c(1:m)^2 * ( x(1:m) - w(1:m) )^2 ) Default values are: c(1:m) = 1 / m w(1:m) = 0.5 Licensing: This code is distributed under the MIT license. Modified: 27 August 2012 Author: John Burkardt Reference: Alan Genz, A Package for Testing Multiple Integration Subroutines, in Numerical Integration: Recent Developments, Software and Applications, edited by Patrick Keast and Graeme Fairweather, Reidel, 1987, pages 337-340, ISBN: 9027725144, LC: QA299.3.N38. Parameters: Input, int M, the dimension of the argument. Input, double C[M], W[M], the problem parameters. Input, int N, the number of points. Input, double X[M*N], the evaluation points. Output, double P04_F[N], the function values. */ { double *f; int i; int j; f = ( double * ) malloc ( n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { f[j] = 0.0; } for ( i = 0; i < m; i++ ) { for ( j = 0; j < n; j++ ) { f[j] = f[j] + pow ( c[i] * ( x[i+j*m] - w[i] ), 2 ); } } for ( j = 0; j < n; j++ ) { f[j] = exp ( - f[j] ); } return f; } /******************************************************************************/ double p04_q ( int m, double c[], double w[] ) /******************************************************************************/ /* Purpose: P04_Q evaluates the integral for problem p04. Licensing: This code is distributed under the MIT license. Modified: 29 August 2012 Author: John Burkardt Reference: Alan Genz, A Package for Testing Multiple Integration Subroutines, in Numerical Integration: Recent Developments, Software and Applications, edited by Patrick Keast and Graeme Fairweather, Reidel, 1987, pages 337-340, ISBN: 9027725144, LC: QA299.3.N38. Parameters: Input, int M, the dimension of the argument. Input, double C[M], W[M], the problem parameters. Output, double P04_Q, the integral. */ { int i; double pi = 3.141592653589793; double q; q = 1.0; for ( i = 0; i < m; i++ ) { q = q * sqrt ( pi ) * ( r8_error ( c[i] * ( 1.0 - w[i] ) ) + r8_error ( c[i] * w[i] ) ) / ( 2.0 * c[i] ); } return q; } /******************************************************************************/ char *p04_title ( ) /******************************************************************************/ /* Purpose: P04_TITLE returns the title of problem p04. Licensing: This code is distributed under the MIT license. Modified: 27 August 2012 Author: John Burkardt Parameters: Output, char *P04_TITLE, the title of the problem. */ { char *title; title = ( char * ) malloc ( 80 * sizeof ( char ) ); strcpy ( title, "Gaussian" ); return title; } /******************************************************************************/ double *p05_d ( int m, int id, double c[], double w[], int n, double x[] ) /******************************************************************************/ /* Purpose: P05_D evaluates any derivative component for problem p05. Discussion: f(x) = exp ( - sum ( c(1:m) * abs ( x(1:m) - w(1:m) ) ) ) Default values are: c(1:m) = 2.0 w(1:m) = 0.5 Licensing: This code is distributed under the MIT license. Modified: 27 August 2012 Author: John Burkardt Reference: Alan Genz, A Package for Testing Multiple Integration Subroutines, in Numerical Integration: Recent Developments, Software and Applications, edited by Patrick Keast and Graeme Fairweather, Reidel, 1987, pages 337-340, ISBN: 9027725144, LC: QA299.3.N38. Parameters: Input, int M, the dimension of the argument. Input, int ID, the spatial coordinate to differentiate. Input, double C[M], W[M], the problem parameters. Input, int N, the number of points. Input, double X[M*N], the evaluation points. Output, double P05_D[N], the ID-th derivative component. */ { double *d; int i; int j; d = ( double * ) malloc ( n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { d[j] = 0.0; } for ( i = 0; i < m; i++ ) { for ( j = 0; j < n; j++ ) { d[j] = d[j] + c[i] * fabs ( x[i+j*m] - w[i] ); } } for ( j = 0; j < n; j++ ) { d[j] = exp ( - d[j] ); } for ( j = 0; j < n; j++ ) { if ( x[id+j*m] - w[id] <= 0.0 ) { d[j] = d[j] * c[id]; } else { d[j] = - d[j] * c[id]; } } return d; } /******************************************************************************/ double *p05_f ( int m, double c[], double w[], int n, double x[] ) /******************************************************************************/ /* Purpose: P05_F evaluates the function for problem p05. Discussion: f(x) = exp ( - sum ( c(1:m) * abs ( x(1:m) - w(1:m) ) ) ) Default values are: c(1:m) = 2.0 w(1:m) = 0.5 Licensing: This code is distributed under the MIT license. Modified: 27 August 2012 Author: John Burkardt Reference: Alan Genz, A Package for Testing Multiple Integration Subroutines, in Numerical Integration: Recent Developments, Software and Applications, edited by Patrick Keast and Graeme Fairweather, Reidel, 1987, pages 337-340, ISBN: 9027725144, LC: QA299.3.N38. Parameters: Input, int M, the dimension of the argument. Input, double C[M], W[M], the problem parameters. Input, int N, the number of points. Input, double X[M*N], the evaluation points. Output, double P05_F[N], the function values. */ { double *f; int i; int j; f = ( double * ) malloc ( n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { f[j] = 0.0; } for ( i = 0; i < m; i++ ) { for ( j = 0; j < n; j++ ) { f[j] = f[j] + c[i] * fabs ( x[i+j*m] - w[i] ); } } for ( j = 0; j < n; j++ ) { f[j] = exp ( - f[j] ); } return f; } /******************************************************************************/ double p05_q ( int m, double c[], double w[] ) /******************************************************************************/ /* Purpose: P05_Q evaluates the integral for problem p05. Licensing: This code is distributed under the MIT license. Modified: 29 August 2012 Author: John Burkardt Reference: Alan Genz, A Package for Testing Multiple Integration Subroutines, in Numerical Integration: Recent Developments, Software and Applications, edited by Patrick Keast and Graeme Fairweather, Reidel, 1987, pages 337-340, ISBN: 9027725144, LC: QA299.3.N38. Parameters: Input, int M, the dimension of the argument. Input, double C[M], W[M], the problem parameters. Output, double P05_Q, the integral. */ { int i; double q; q = 1.0; for ( i = 0; i < m; i++ ) { /* W < 0 < 1 | X - W | = X - W from 0 to 1. */ if ( w[i] < 0.0 ) { q = q * ( exp ( - c[i] * ( - w[i] ) ) - exp ( - c[i] * ( 1.0 - w[i] ) ) ) / c[i]; } /* 0 < W < 1 | X - W | = W - X from 0 to Z, = X - W from Z to 1. */ else if ( w[i] < 1.0 ) { q = q * ( 2.0 - exp ( - c[i] * ( w[i] ) ) - exp ( - c[i] * ( 1.0 - w[i] ) ) ) / c[i]; } /* 0 < 1 < W | X - W | = W - X from 0 to 1. */ else { q = q * ( exp ( - c[i] * ( w[i] - 1.0 ) ) - exp ( - c[i] * ( w[i] ) ) ) / c[i]; } } return q; } /******************************************************************************/ char *p05_title ( ) /******************************************************************************/ /* Purpose: P05_TITLE returns the title of problem p05. Licensing: This code is distributed under the MIT license. Modified: 27 August 2012 Author: John Burkardt Parameters: Output, char *P05_TITLE, the title of the problem. */ { char *title; title = ( char * ) malloc ( 80 * sizeof ( char ) ); strcpy ( title, "Continuous" ); return title; } /******************************************************************************/ double *p06_d ( int m, int id, double c[], double w[], int n, double x[] ) /******************************************************************************/ /* Purpose: P06_D evaluates any derivative component for problem p06. Discussion: f(x) = exp ( c(1:m) * x(1:m) ) if x(1) <= w(1) and x(2) <= w(2). 0 otherwise Default values are: c(1:m) = 0.5^(1/m) w(1:2) = 0.5^(1/m) Licensing: This code is distributed under the MIT license. Modified: 27 August 2012 Author: John Burkardt Reference: Alan Genz, A Package for Testing Multiple Integration Subroutines, in Numerical Integration: Recent Developments, Software and Applications, edited by Patrick Keast and Graeme Fairweather, Reidel, 1987, pages 337-340, ISBN: 9027725144, LC: QA299.3.N38. Parameters: Input, int M, the dimension of the argument. Input, int ID, the spatial coordinate to differentiate. Input, double C[M], W[M], the problem parameters. Input, int N, the number of points. Input, double X[M*N], the evaluation points. Output, double P06_D[N], the ID-th derivative component. */ { double *d; int i; int j; d = ( double * ) malloc ( n * sizeof ( double ) ); if ( m == 1 ) { for ( j = 0; j < n; j++ ) { d[j] = c[0] * exp ( c[0] * x[0+j*m] ); } for ( j = 0; j < n; j++ ) { if ( w[0] < x[0+j*m] ) { d[j] = 0.0; } } } else { for ( j = 0; j < n; j++ ) { d[j] = 0.0; } for ( i = 0; i < m; i++ ) { for ( j = 0; j < n; j++ ) { d[j] = d[j] + c[i] * x[i+j*m]; } } for ( j = 0; j < n; j++ ) { d[j] = c[id] * exp ( d[j] ); } for ( j = 0; j < n; j++ ) { if ( w[0] < x[0+j*m] || w[1] < x[1+j*m] ) { d[j] = 0.0; } } } return d; } /******************************************************************************/ double *p06_f ( int m, double c[], double w[], int n, double x[] ) /******************************************************************************/ /* Purpose: P06_F evaluates the function for problem p06. Discussion: f(x) = exp ( c(1:m) * x(1:m) ) if x(1) <= w(1) and x(2) <= w(2). 0 otherwise Default values are: c(1:m) = 0.5^(1/m) w(1:2) = 0.5^(1/m) Licensing: This code is distributed under the MIT license. Modified: 27 August 2012 Author: John Burkardt Reference: Alan Genz, A Package for Testing Multiple Integration Subroutines, in Numerical Integration: Recent Developments, Software and Applications, edited by Patrick Keast and Graeme Fairweather, Reidel, 1987, pages 337-340, ISBN: 9027725144, LC: QA299.3.N38. Parameters: Input, int M, the dimension of the argument. Input, double C[M], W[M], the problem parameters. Input, int N, the number of points. Input, double X[M*N], the evaluation points. Output, double P06_F[N], the function values. */ { double *f; int i; int j; f = ( double * ) malloc ( n * sizeof ( double ) ); if ( m == 1 ) { for ( j = 0; j < n; j++ ) { f[j] = exp ( c[0] * x[0+j*m] ); } for ( j = 0; j < n; j++ ) { if ( w[0] < x[0+j*m] ) { f[j] = 0.0; } } } else { for ( j = 0; j < n; j++ ) { f[j] = 0.0; } for ( i = 0; i < m; i++ ) { for ( j = 0; j < n; j++ ) { f[j] = f[j] + c[i] * x[i+j*m]; } } for ( j = 0; j < n; j++ ) { f[j] = exp ( f[j] ); } for ( j = 0; j < n; j++ ) { if ( w[0] < x[0+j*m] || w[1] < x[1+j*m] ) { f[j] = 0.0; } } } return f; } /******************************************************************************/ double p06_q ( int m, double c[], double w[] ) /******************************************************************************/ /* Purpose: P06_Q evaluates the integral for problem p06. Licensing: This code is distributed under the MIT license. Modified: 29 August 2012 Author: John Burkardt Reference: Alan Genz, A Package for Testing Multiple Integration Subroutines, in Numerical Integration: Recent Developments, Software and Applications, edited by Patrick Keast and Graeme Fairweather, Reidel, 1987, pages 337-340, ISBN: 9027725144, LC: QA299.3.N38. Parameters: Input, int M, the dimension of the argument. Input, double C[M], W[M], the problem parameters. Output, double P06_Q, the integral. */ { int i; double q; /* To simplify the calculation, force W(3:M) to be at least 1.0. */ for ( i = 2; i < m; i++ ) { w[i] = 1.0; } q = 1.0; for ( i = 0; i < m; i++ ) { if ( w[i] <= 0.0 ) { q = q * 0.0; } else if ( w[i] <= 1.0 ) { if ( c[i] == 0.0 ) { q = q * w[i]; } else { q = q * ( exp ( c[i] * w[i] ) - 1.0 ) / c[i]; } } else { if ( c[i] != 0.0 ) { q = q * ( exp ( c[i] * w[i] ) - 1.0 ) / c[i]; } } } return q; } /******************************************************************************/ char *p06_title ( ) /******************************************************************************/ /* Purpose: P06_TITLE returns the title of problem p06. Licensing: This code is distributed under the MIT license. Modified: 27 August 2012 Author: John Burkardt Parameters: Output, char *P06_TITLE, the title of the problem. */ { char *title; title = ( char * ) malloc ( 80 * sizeof ( char ) ); strcpy ( title, "Discontinuous" ); return title; } /******************************************************************************/ double r8_error ( double x ) /******************************************************************************/ /* Purpose: R8_ERROR evaluates the error function of an R8 argument. Licensing: This code is distributed under the MIT license. Modified: 17 January 2012 Author: Original FORTRAN77 version by Wayne Fullerton. C version by John Burkardt. Reference: Wayne Fullerton, Portable Special Function Routines, in Portability of Numerical Software, edited by Wayne Cowell, Lecture Notes in Computer Science, Volume 57, Springer 1977, ISBN: 978-3-540-08446-4, LC: QA297.W65. Parameters: Input, double X, the argument. Output, double R8_ERF, the error function of X. */ { static double erfcs[21] = { -0.49046121234691808039984544033376E-01, -0.14226120510371364237824741899631, +0.10035582187599795575754676712933E-01, -0.57687646997674847650827025509167E-03, +0.27419931252196061034422160791471E-04, -0.11043175507344507604135381295905E-05, +0.38488755420345036949961311498174E-07, -0.11808582533875466969631751801581E-08, +0.32334215826050909646402930953354E-10, -0.79910159470045487581607374708595E-12, +0.17990725113961455611967245486634E-13, -0.37186354878186926382316828209493E-15, +0.71035990037142529711689908394666E-17, -0.12612455119155225832495424853333E-18, +0.20916406941769294369170500266666E-20, -0.32539731029314072982364160000000E-22, +0.47668672097976748332373333333333E-24, -0.65980120782851343155199999999999E-26, +0.86550114699637626197333333333333E-28, -0.10788925177498064213333333333333E-29, +0.12811883993017002666666666666666E-31 }; static int nterf = 0; static double sqeps = 0.0; static double sqrtpi = 1.77245385090551602729816748334115; double value; static double xbig = 0.0; double y; if ( nterf == 0 ) { nterf = inits ( erfcs, 21, 0.1 * r8_mach ( 3 ) ); xbig = sqrt ( - log ( sqrtpi * r8_mach ( 3 ) ) ); sqeps = sqrt ( 2.0 * r8_mach ( 3 ) ); } y = fabs ( x ); if ( y <= sqeps ) { value = 2.0 * x / sqrtpi; } else if ( y <= 1.0 ) { value = x * ( 1.0 + csevl ( 2.0 * x * x - 1.0, erfcs, nterf ) ); } else if ( y <= xbig ) { value = 1.0 - r8_errorc ( y ); if ( x < 0.0 ) { value = - value; } } else { value = 1.0; if ( x < 0.0 ) { value = - value; } } return value; } /******************************************************************************/ double r8_errorc ( double x ) /******************************************************************************/ /* Purpose: R8_ERRORC evaluates the co-error function of an R8 argument. Licensing: This code is distributed under the MIT license. Modified: 17 January 2012 Author: Original FORTRAN77 version by Wayne Fullerton. C version by John Burkardt. Reference: Wayne Fullerton, Portable Special Function Routines, in Portability of Numerical Software, edited by Wayne Cowell, Lecture Notes in Computer Science, Volume 57, Springer 1977, ISBN: 978-3-540-08446-4, LC: QA297.W65. Parameters: Input, double X, the argument. Output, double R8_ERRORC, the co-error function of X. */ { static double erc2cs[49] = { -0.6960134660230950112739150826197E-01, -0.4110133936262089348982212084666E-01, +0.3914495866689626881561143705244E-02, -0.4906395650548979161280935450774E-03, +0.7157479001377036380760894141825E-04, -0.1153071634131232833808232847912E-04, +0.1994670590201997635052314867709E-05, -0.3642666471599222873936118430711E-06, +0.6944372610005012589931277214633E-07, -0.1371220902104366019534605141210E-07, +0.2788389661007137131963860348087E-08, -0.5814164724331161551864791050316E-09, +0.1238920491752753181180168817950E-09, -0.2690639145306743432390424937889E-10, +0.5942614350847910982444709683840E-11, -0.1332386735758119579287754420570E-11, +0.3028046806177132017173697243304E-12, -0.6966648814941032588795867588954E-13, +0.1620854541053922969812893227628E-13, -0.3809934465250491999876913057729E-14, +0.9040487815978831149368971012975E-15, -0.2164006195089607347809812047003E-15, +0.5222102233995854984607980244172E-16, -0.1269729602364555336372415527780E-16, +0.3109145504276197583836227412951E-17, -0.7663762920320385524009566714811E-18, +0.1900819251362745202536929733290E-18, -0.4742207279069039545225655999965E-19, +0.1189649200076528382880683078451E-19, -0.3000035590325780256845271313066E-20, +0.7602993453043246173019385277098E-21, -0.1935909447606872881569811049130E-21, +0.4951399124773337881000042386773E-22, -0.1271807481336371879608621989888E-22, +0.3280049600469513043315841652053E-23, -0.8492320176822896568924792422399E-24, +0.2206917892807560223519879987199E-24, -0.5755617245696528498312819507199E-25, +0.1506191533639234250354144051199E-25, -0.3954502959018796953104285695999E-26, +0.1041529704151500979984645051733E-26, -0.2751487795278765079450178901333E-27, +0.7290058205497557408997703680000E-28, -0.1936939645915947804077501098666E-28, +0.5160357112051487298370054826666E-29, -0.1378419322193094099389644800000E-29, +0.3691326793107069042251093333333E-30, -0.9909389590624365420653226666666E-31, +0.2666491705195388413323946666666E-31 }; static double erfccs[59] = { +0.715179310202924774503697709496E-01, -0.265324343376067157558893386681E-01, +0.171115397792085588332699194606E-02, -0.163751663458517884163746404749E-03, +0.198712935005520364995974806758E-04, -0.284371241276655508750175183152E-05, +0.460616130896313036969379968464E-06, -0.822775302587920842057766536366E-07, +0.159214187277090112989358340826E-07, -0.329507136225284321486631665072E-08, +0.722343976040055546581261153890E-09, -0.166485581339872959344695966886E-09, +0.401039258823766482077671768814E-10, -0.100481621442573113272170176283E-10, +0.260827591330033380859341009439E-11, -0.699111056040402486557697812476E-12, +0.192949233326170708624205749803E-12, -0.547013118875433106490125085271E-13, +0.158966330976269744839084032762E-13, -0.472689398019755483920369584290E-14, +0.143587337678498478672873997840E-14, -0.444951056181735839417250062829E-15, +0.140481088476823343737305537466E-15, -0.451381838776421089625963281623E-16, +0.147452154104513307787018713262E-16, -0.489262140694577615436841552532E-17, +0.164761214141064673895301522827E-17, -0.562681717632940809299928521323E-18, +0.194744338223207851429197867821E-18, -0.682630564294842072956664144723E-19, +0.242198888729864924018301125438E-19, -0.869341413350307042563800861857E-20, +0.315518034622808557122363401262E-20, -0.115737232404960874261239486742E-20, +0.428894716160565394623737097442E-21, -0.160503074205761685005737770964E-21, +0.606329875745380264495069923027E-22, -0.231140425169795849098840801367E-22, +0.888877854066188552554702955697E-23, -0.344726057665137652230718495566E-23, +0.134786546020696506827582774181E-23, -0.531179407112502173645873201807E-24, +0.210934105861978316828954734537E-24, -0.843836558792378911598133256738E-25, +0.339998252494520890627359576337E-25, -0.137945238807324209002238377110E-25, +0.563449031183325261513392634811E-26, -0.231649043447706544823427752700E-26, +0.958446284460181015263158381226E-27, -0.399072288033010972624224850193E-27, +0.167212922594447736017228709669E-27, -0.704599152276601385638803782587E-28, +0.297976840286420635412357989444E-28, -0.126252246646061929722422632994E-28, +0.539543870454248793985299653154E-29, -0.238099288253145918675346190062E-29, +0.109905283010276157359726683750E-29, -0.486771374164496572732518677435E-30, +0.152587726411035756763200828211E-30 }; static double erfcs[21] = { -0.49046121234691808039984544033376E-01, -0.14226120510371364237824741899631, +0.10035582187599795575754676712933E-01, -0.57687646997674847650827025509167E-03, +0.27419931252196061034422160791471E-04, -0.11043175507344507604135381295905E-05, +0.38488755420345036949961311498174E-07, -0.11808582533875466969631751801581E-08, +0.32334215826050909646402930953354E-10, -0.79910159470045487581607374708595E-12, +0.17990725113961455611967245486634E-13, -0.37186354878186926382316828209493E-15, +0.71035990037142529711689908394666E-17, -0.12612455119155225832495424853333E-18, +0.20916406941769294369170500266666E-20, -0.32539731029314072982364160000000E-22, +0.47668672097976748332373333333333E-24, -0.65980120782851343155199999999999E-26, +0.86550114699637626197333333333333E-28, -0.10788925177498064213333333333333E-29, +0.12811883993017002666666666666666E-31 }; double eta; static int nterc2 = 0; static int nterf = 0; static int nterfc = 0; static double sqeps = 0.0; static double sqrtpi = 1.77245385090551602729816748334115; double value; static double xmax = 0.0; static double xsml = 0.0; double y; if ( nterf == 0 ) { eta = 0.1 * r8_mach ( 3 ); nterf = inits ( erfcs, 21, eta ); nterfc = inits ( erfccs, 59, eta ); nterc2 = inits ( erc2cs, 49, eta ); xsml = - sqrt ( - log ( sqrtpi * r8_mach ( 3 ) ) ); xmax = sqrt (- log ( sqrtpi * r8_mach ( 1 ) ) ); xmax = xmax - 0.5 * log ( xmax ) / xmax - 0.01; sqeps = sqrt ( 2.0 * r8_mach ( 3 ) ); } if ( x <= xsml ) { value = 2.0; return value; } if ( xmax < x ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8_ERRORC - Warning!\n" ); fprintf ( stderr, " X so big that ERFC underflows.\n" ); value = 0.0; return value; } y = fabs ( x ); if ( y < sqeps ) { value = 1.0 - 2.0 * x / sqrtpi; return value; } else if ( y <= 1.0 ) { value = 1.0 - x * ( 1.0 + csevl ( 2.0 * x * x - 1.0, erfcs, nterf ) ); return value; } y = y * y; if ( y <= 4.0 ) { value = exp ( - y ) / fabs ( x ) * ( 0.5 + csevl ( ( 8.0 / y - 5.0 ) / 3.0, erc2cs, nterc2 ) ); } else { value = exp ( - y ) / fabs ( x ) * ( 0.5 + csevl ( 8.0 / y - 1.0, erfccs, nterfc ) ); } if ( x < 0.0 ) { value = 2.0 - value; } return value; } /******************************************************************************/ double r8_mach ( int i ) /******************************************************************************/ /* Purpose: R8_MACH returns double precision real machine constants. Discussion: Assuming that the internal representation of a double precision real number is in base B, with T the number of base-B digits in the mantissa, and EMIN the smallest possible exponent and EMAX the largest possible exponent, then R8_MACH(1) = B^(EMIN-1), the smallest positive magnitude. R8_MACH(2) = B^EMAX*(1-B^(-T)), the largest magnitude. R8_MACH(3) = B^(-T), the smallest relative spacing. R8_MACH(4) = B^(1-T), the largest relative spacing. R8_MACH(5) = log10(B). Licensing: This code is distributed under the MIT license. Modified: 24 April 2007 Author: Original FORTRAN77 version by Phyllis Fox, Andrew Hall, Norman Schryer. C version by John Burkardt. Reference: Phyllis Fox, Andrew Hall, Norman Schryer, Algorithm 528: Framework for a Portable Library, ACM Transactions on Mathematical Software, Volume 4, Number 2, June 1978, page 176-188. Parameters: Input, int I, chooses the parameter to be returned. 1 <= I <= 5. Output, double R8_MACH, the value of the chosen parameter. */ { double value; if ( i < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8_MACH - Fatal error!\n" ); fprintf ( stderr, " The input argument I is out of bounds.\n" ); fprintf ( stderr, " Legal values satisfy 1 <= I <= 5.\n" ); fprintf ( stderr, " I = %d\n", i ); value = 0.0; exit ( 1 ); } else if ( i == 1 ) { value = 4.450147717014403E-308; } else if ( i == 2 ) { value = 8.988465674311579E+307; } else if ( i == 3 ) { value = 1.110223024625157E-016; } else if ( i == 4 ) { value = 2.220446049250313E-016; } else if ( i == 5 ) { value = 0.301029995663981E+000; } else if ( 5 < i ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8_MACH - Fatal error!\n" ); fprintf ( stderr, " The input argument I is out of bounds.\n" ); fprintf ( stderr, " Legal values satisfy 1 <= I <= 5.\n" ); fprintf ( stderr, " I = %d\n", i ); value = 0.0; exit ( 1 ); } return value; } /******************************************************************************/ void tuple_next ( int m1, int m2, int n, int *rank, int x[] ) /******************************************************************************/ /* Purpose: TUPLE_NEXT computes the next element of a tuple space. Discussion: The elements are N vectors. Each entry is constrained to lie between M1 and M2. The elements are produced one at a time. The first element is (M1,M1,...,M1), the second element is (M1,M1,...,M1+1), and the last element is (M2,M2,...,M2) Intermediate elements are produced in lexicographic order. Example: N = 2, M1 = 1, M2 = 3 INPUT OUTPUT ------- ------- Rank X Rank X ---- --- ----- --- 0 * * 1 1 1 1 1 1 2 1 2 2 1 2 3 1 3 3 1 3 4 2 1 4 2 1 5 2 2 5 2 2 6 2 3 6 2 3 7 3 1 7 3 1 8 3 2 8 3 2 9 3 3 9 3 3 0 0 0 Licensing: This code is distributed under the MIT license. Modified: 29 April 2003 Author: John Burkardt Parameters: Input, int M1, M2, the minimum and maximum entries. Input, int N, the number of components. Input/output, int *RANK, counts the elements. On first call, set RANK to 0. Thereafter, the output value of RANK will indicate the order of the element returned. When there are no more elements, RANK will be returned as 0. Input/output, int X[N], on input the previous tuple. On output, the next tuple. */ { int i; int j; if ( m2 < m1 ) { *rank = 0; return; } if ( *rank <= 0 ) { for ( i = 0; i < n; i++ ) { x[i] = m1; } *rank = 1; } else { *rank = *rank + 1; i = n - 1; for ( ; ; ) { if ( x[i] < m2 ) { x[i] = x[i] + 1; break; } x[i] = m1; if ( i == 0 ) { *rank = 0; for ( j = 0; j < n; j++ ) { x[j] = m1; } break; } i = i - 1; } } return; }