# include # include # include # include # include # include "sparse_grid_hw.h" /******************************************************************************/ int cce_order ( int l ) /******************************************************************************/ /* Purpose: CCE_ORDER: order of a Clenshaw-Curtis Exponential rule from the level. Licensing: This code is distributed under the MIT license. Modified: 22 February 2014 Author: John Burkardt. Parameters: Input, int L, the level of the rule. 1 <= L. Output, int CCE_ORDER, the order of the rule. */ { int n; if ( l < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "CCE_ORDER - Fatal error!\n" ); fprintf ( stderr, " 1 <= L required.\n" ); fprintf ( stderr, " Input L = %d\n", l ); exit ( 1 ); } else if ( l == 1 ) { n = 1; } else { n = i4_power ( 2, l - 1 ) + 1; } return n; } /******************************************************************************/ int ccl_order ( int l ) /******************************************************************************/ /* Purpose: CCL_ORDER computes the order of a CCL rule from the level. Licensing: This code is distributed under the MIT license. Modified: 19 February 2014 Author: John Burkardt. Parameters: Input, int L, the level of the rule. 1 <= L. Output, int CCL_ORDER, the order of the rule. */ { int n; if ( l < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "CCL_ORDER - Fatal error!\n" ); fprintf ( stderr, " 1 <= L required.\n" ); fprintf ( stderr, " Input L = %d\n", l ); exit ( 1 ); } n = 2 * l - 1; return n; } /******************************************************************************/ int ccs_order ( int l ) /******************************************************************************/ /* Purpose: CCS_ORDER: order of a "slow growth" Clenshaw Curtis quadrature rule. Discussion: Our convention is that the abscissas are numbered from left to right. The rule is defined on [0,1]. The integral to approximate: Integral ( 0 <= X <= 1 ) F(X) dX The quadrature rule: Sum ( 1 <= I <= N ) W(I) * F ( X(I) ) The input value L requests a rule of precision at least 2*L-1. In order to preserve nestedness, this function returns the order of a rule which is the smallest value of the form 1+2^E which is greater than or equal to 2*L-1. L 2*L-1 N -- ----- -- 1 1 1 2 3 3 3 5 5 4 7 9 5 9 9 6 11 17 7 13 17 8 15 17 9 17 17 10 19 33 Licensing: This code is distributed under the MIT license. Modified: 08 December 2012 Author: John Burkardt Parameters: Input, int L, the level of the rule. 1 <= L. Output, int CCS_ORDER, the appropriate order. */ { int n; if ( l < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "CCS_ORDER - Fatal error!\n" ); fprintf ( stderr, " Illegal value of L = %d\n", l ); exit ( 1 ); } /* Find the order N that satisfies the precision requirement. */ if ( l == 1 ) { n = 1; } else { n = 3; while ( n < 2 * l - 1 ) { n = 2 * n - 1; } } return n; } /******************************************************************************/ void cc ( int n, double x[], double w[] ) /******************************************************************************/ /* Purpose: CC computes a Clenshaw Curtis quadrature rule based on order. Discussion: Our convention is that the abscissas are numbered from left to right. The rule is defined on [0,1]. The integral to approximate: Integral ( 0 <= X <= 1 ) F(X) dX The quadrature rule: Sum ( 1 <= I <= N ) W(I) * F ( X(I) ) Licensing: This code is distributed under the MIT license. Modified: 08 December 2012 Author: John Burkardt Parameters: Input, int N, the order of the rule. 1 <= N. Output, double X[N], the abscissas. Output, double W[N], the weights. */ { double b; int i; int j; double pi = 3.141592653589793; double theta; if ( n < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "CC - Fatal error!\n" ); fprintf ( stderr, " Illegal value of N = %d\n", n ); exit ( 1 ); } if ( n == 1 ) { x[0] = 0.0; w[0] = 2.0; } else { for ( i = 0; i < n; i++ ) { x[i] = cos ( ( double ) ( n - 1 - i ) * pi / ( double ) ( n - 1 ) ); } x[0] = -1.0; if ( ( n % 2 ) == 1 ) { x[(n+1)/2-1] = 0.0; } x[n-1] = +1.0; for ( i = 0; i < n; i++ ) { theta = ( double ) ( i ) * pi / ( double ) ( n - 1 ); w[i] = 1.0; for ( j = 1; j <= ( n - 1 ) / 2; j++ ) { if ( 2 * j == ( n - 1 ) ) { b = 1.0; } else { b = 2.0; } w[i] = w[i] - b * cos ( 2.0 * ( double ) ( j ) * theta ) / ( double ) ( 4 * j * j - 1 ); } } w[0] = w[0] / ( double ) ( n - 1 ); for ( j = 1; j < n - 1; j++ ) { w[j] = 2.0 * w[j] / ( double ) ( n - 1 ); } w[n-1] = w[n-1] / ( double ) ( n - 1 ); } /* Transform from [-1,+1] to [0,1]. */ rule_adjust ( -1.0, +1.0, 0.0, +1.0, n, x, w ); return; } /******************************************************************************/ double cpu_time ( ) /******************************************************************************/ /* Purpose: CPU_TIME returns the current reading on the CPU clock. Discussion: The CPU time measurements available through this routine are often not very accurate. In some cases, the accuracy is no better than a hundredth of a second. Licensing: This code is distributed under the MIT license. Modified: 06 June 2005 Author: John Burkardt Parameters: Output, double CPU_TIME, the current reading of the CPU clock, in seconds. */ { double value; value = ( double ) clock ( ) / ( double ) CLOCKS_PER_SEC; return value; } /******************************************************************************/ double fn_integral ( int d ) /******************************************************************************/ /* Purpose: FN_INTEGRAL is the integral of the Hermite test function. Licensing: This code is distributed under the MIT license. Modified: 08 December 2012 Author: John Burkardt. Parameters: Input, int D, the spatial dimension. Output, double FN_INTEGRAL, the integral value. */ { int exponent = 6; double value; value = ( double ) ( i4_factorial2 ( exponent - 1 ) ); return value; } /******************************************************************************/ double *fn_value ( int d, int n, double x[] ) /******************************************************************************/ /* Purpose: FN_VALUE is a Hermite test function. Licensing: This code is distributed under the MIT license. Modified: 08 May 2012 Author: John Burkardt. Parameters: Input, int D, the spatial dimension. Input, int N, the number of points. Input, double X[D*N], the points. Output, double FN_VALUE[N], the function values. */ { int exponent = 6; double *fx; int i; fx = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { fx[i] = pow ( x[0+i*d], exponent ); } return fx; } /******************************************************************************/ double fu_integral ( int d ) /******************************************************************************/ /* Purpose: FU_INTEGRAL is the integral of the test function for the [0,1]^D interval. Discussion: The same function, integrated over [-1,+1]^D, has an integral that is 2^D times larger. Licensing: This code is distributed under the MIT license. Modified: 08 December 2012 Author: John Burkardt. Parameters: Input, int D, the spatial dimension. Output, double FU_INTEGRAL, the integral value. */ { double value; value = pow ( 0.5 * erf ( 0.5 / sqrt ( 2.0 ) ), d ); return value; } /******************************************************************************/ double *fu_value ( int d, int n, double x[] ) /******************************************************************************/ /* Purpose: FU_VALUE is a sample function for the [0,1]^D interval. Licensing: This code is distributed under the MIT license. Modified: 08 December 2012 Author: John Burkardt. Parameters: Input, int D, the spatial dimension. Input, int N, the number of points. Input, double X[D*N], the points. Output, double FU_VALUE[N], the function values. */ { double *fx; int i; int j; double pi = 3.141592653589793; fx = ( double * ) malloc ( n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { fx[j] = 1.0; for ( i = 0; i < d; i++ ) { fx[j] = fx[j] * exp ( - pow ( x[i+j*d] / 2.0, 2 ) / 2.0 ) / 2.0 / sqrt ( 2.0 * pi ); } } return fx; } /******************************************************************************/ int *get_seq ( int d, int norm, int seq_num ) /******************************************************************************/ /* Purpose: GET_SEQ generates all positive integer D-vectors that sum to NORM. Discussion: This function computes a list, in reverse dictionary order, of all D-vectors of positive values that sum to NORM. For example, call get_seq ( 3, 5, 6, fs ) returns 3 1 1 2 2 1 2 1 2 1 3 1 1 2 2 1 1 3 Licensing: This code is distributed under the MIT license. Modified: 08 December 2012 Author: Original MATLAB version by Florian Heiss, Viktor Winschel. C version by John Burkardt. Reference: Florian Heiss, Viktor Winschel, Likelihood approximation by numerical integration on sparse grids, Journal of Econometrics, Volume 144, 2008, pages 62-80. Parameters: Input, int D, the dimension. 1 <= D. Input, int NORM, the value that each row must sum to. D <= NORM. Input, int SEQ_NUM, the number of rows of FS. Output, int GET_SEQ[SEQ_NUM*D]. Each row of FS represents one vector with all elements positive and summing to NORM. */ { int a; int c; int *fs; int i; int row; int *seq; seq = ( int * ) malloc ( d * sizeof ( int ) ); fs = ( int * ) malloc ( seq_num * d * sizeof ( int ) ); if ( norm < d ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GET_SEQ - Fatal error!\n" ); fprintf ( stderr, " NORM = %d < D = %d\n", norm, d ); exit ( 1 ) ; } for ( i = 0; i < d; i++ ) { seq[i] = 0; } /* The algorithm is written to work with vectors whose minimum value is allowed to be zero. So we subtract D from NORM at the beginning and then increment the result vectors by 1 at the end! */ a = norm - d; seq[0] = a; row = 0; for ( i = 0; i < d; i++ ) { fs[row+i*seq_num] = seq[i] + 1; } c = 0; while ( seq[d-1] < a ) { if ( c == d - 1 ) { for ( i = c - 1; 0 <= i; i-- ) { c = i; if ( seq[i] != 0 ) { break; } } } seq[c] = seq[c] - 1; c = c + 1; seq[c] = a; for ( i = 0; i < c; i++ ) { seq[c] = seq[c] - seq[i]; } if ( c < d - 1 ) { for ( i = c + 1; i < d; i++ ) { seq[i] = 0; } } row = row + 1; for ( i = 0; i < d; i++ ) { fs[row+i*seq_num] = seq[i] + 1; } } free ( seq ); return fs; } /******************************************************************************/ void gqn ( int n, double x[], double w[] ) /******************************************************************************/ /* Purpose: GQN provides data for Gauss quadrature with a normal weight. Discussion: This data assumes integration over the interval (-oo,+oo) with weight function w(x) = exp(-x*x/2)/sqrt(2*pi). Licensing: This code is distributed under the MIT license. Modified: 10 December 2012 Author: Original MATLAB version by Florian Heiss, Viktor Winschel. C version by John Burkardt. Reference: Florian Heiss, Viktor Winschel, Likelihood approximation by numerical integration on sparse grids, Journal of Econometrics, Volume 144, 2008, pages 62-80. Parameters: Input, int N, the number of points and weights. 1 <= N <= 25. Output, double X[N], the nodes. Output, double W[N], the weights. */ { double x01[1] = { 0.00000000000000000 }; double w01[1] = { 1.00000000000000000 }; double x02[2] = { -1.00000000000000000, 1.00000000000000000 }; double w02[2] = { 0.50000000000000000, 0.50000000000000000 }; double x03[3] = { -1.73205080756887719, 0.00000000000000000, 1.73205080756887719 }; double w03[3] = { 0.166666666666666741, 0.66666666666666663, 0.166666666666666741 }; double x04[4] = { -2.33441421833897733, -0.741963784302725915, 0.741963784302725915, 2.33441421833897733 }; double w04[4] = { 0.0458758547680684983, 0.454124145231931453, 0.454124145231931453, 0.0458758547680684983 }; double x05[5] = { -2.85697001387280558, -1.35562617997426593, 0.00000000000000000, 1.35562617997426593, 2.85697001387280558 }; double w05[5] = { 0.011257411327720691, 0.22207592200561263, 0.533333333333333437, 0.22207592200561263, 0.011257411327720691 }; double x06[6] = { -3.32425743355211933, -1.88917587775371087, -0.616706590192594217, 0.616706590192594217, 1.88917587775371087, 3.32425743355211933 }; double w06[6] = { 0.00255578440205624308, 0.0886157460419145226, 0.408828469556029195, 0.408828469556029195, 0.0886157460419145226, 0.00255578440205624308 }; double x07[7] = { -3.75043971772574247, -2.36675941073454155, -1.15440539473996817, 0.00000000000000000, 1.15440539473996817, 2.36675941073454155, 3.75043971772574247 }; double w07[7] = { 0.000548268855972218754, 0.0307571239675864909, 0.240123178605012505, 0.457142857142857573, 0.240123178605012505, 0.0307571239675864909, 0.000548268855972218754 }; double x08[8] = { -4.14454718612589446, -2.80248586128754162, -1.63651904243510815, -0.539079811351375171, 0.539079811351375171, 1.63651904243510815, 2.80248586128754162, 4.14454718612589446 }; double w08[8] = { 0.000112614538375367836, 0.00963522012078826297, 0.117239907661758971, 0.373012257679077364, 0.373012257679077364, 0.117239907661758971, 0.00963522012078826297, 0.000112614538375367836 }; double x09[9] = { -4.51274586339978256, -3.20542900285647026, -2.07684797867783022, -1.02325566378913257, 0.00000000000000000, 1.02325566378913257, 2.07684797867783022, 3.20542900285647026, 4.51274586339978256 }; double w09[9] = { 2.23458440077465626E-05, 0.0027891413212317675, 0.0499164067652179688, 0.244097502894939089, 0.406349206349206848, 0.244097502894939089, 0.0499164067652179688, 0.0027891413212317675, 2.23458440077465626E-05 }; double x10[10] = { -4.85946282833231269, -3.58182348355192692, -2.48432584163895465, -1.46598909439115821, -0.484935707515497638, 0.484935707515497638, 1.46598909439115821, 2.48432584163895465, 3.58182348355192692, 4.85946282833231269 }; double w10[10] = { 4.3106526307183106E-06, 0.000758070934312219725, 0.0191115805007703171, 0.135483702980267295, 0.344642334932019401, 0.344642334932019401, 0.135483702980267295, 0.0191115805007703171, 0.000758070934312219725, 4.3106526307183106E-06 }; double x11[11] = { -5.18800122437487143, -3.93616660712997746, -2.8651231606436447, -1.87603502015484591, -0.928868997381063877, 0.00000000000000000, 0.928868997381063877, 1.87603502015484591, 2.8651231606436447, 3.93616660712997746, 5.18800122437487143 }; double w11[11] = { 8.12184979021490357E-07, 0.000195671930271223244, 0.0067202852355372697, 0.0661387460710576441, 0.242240299873970027, 0.369408369408369575, 0.242240299873970027, 0.0661387460710576441, 0.0067202852355372697, 0.000195671930271223244, 8.12184979021490357E-07 }; double x12[12] = { -5.50090170446774795, -4.27182584793228148, -3.22370982877009737, -2.25946445100079929, -1.34037519715161668, -0.444403001944139009, 0.444403001944139009, 1.34037519715161668, 2.25946445100079929, 3.22370982877009737, 4.27182584793228148, 5.50090170446774795 }; double w12[12] = { 1.49992716763715968E-07, 4.83718492259060763E-05, 0.00220338068753318491, 0.0291166879123641378, 0.146967048045329951, 0.321664361512830066, 0.321664361512830066, 0.146967048045329951, 0.0291166879123641378, 0.00220338068753318491, 4.83718492259060763E-05, 1.49992716763715968E-07 }; double x13[13] = { -5.8001672523865011, -4.59139844893652072, -3.56344438028163468, -2.62068997343221488, -1.7254183795882394, -0.856679493519450053, 0.00000000000000000, 0.856679493519450053, 1.7254183795882394, 2.62068997343221488, 3.56344438028163468, 4.59139844893652072, 5.8001672523865011 }; double w13[13] = { 2.72262764280590389E-08, 1.15265965273338848E-05, 0.000681236350442926191, 0.0117705605059965426, 0.0791689558604501409, 0.237871522964135884, 0.340992340992341492, 0.237871522964135884, 0.0791689558604501409, 0.0117705605059965426, 0.000681236350442926191, 1.15265965273338848E-05, 2.72262764280590389E-08 }; double x14[14] = { -6.08740954690129144, -4.89693639734556463, -3.88692457505976963, -2.96303657983866753, -2.08834474570194439, -1.24268895548546432, -0.412590457954601808, 0.412590457954601808, 1.24268895548546432, 2.08834474570194439, 2.96303657983866753, 3.88692457505976963, 4.89693639734556463, 6.08740954690129144 }; double w14[14] = { 4.86816125774838718E-09, 2.66099134406763342E-06, 0.00020033955376074381, 0.00442891910694740657, 0.0386501088242534319, 0.154083339842513656, 0.302634626813019447, 0.302634626813019447, 0.154083339842513656, 0.0386501088242534319, 0.00442891910694740657, 0.00020033955376074381, 2.66099134406763342E-06, 4.86816125774838718E-09 }; double x15[15] = { -6.36394788882983775, -5.19009359130478209, -4.19620771126901548, -3.28908242439876641, -2.43243682700975805, -1.60671006902873015, -0.799129068324548109, 0.00000000000000000, 0.799129068324548109, 1.60671006902873015, 2.43243682700975805, 3.28908242439876641, 4.19620771126901548, 5.19009359130478209, 6.36394788882983775 }; double w15[15] = { 8.58964989963318053E-10, 5.97541959792059611E-07, 5.64214640518901565E-05, 0.00156735750354995707, 0.0173657744921376159, 0.0894177953998444436, 0.232462293609732223, 0.318259518259518204, 0.232462293609732223, 0.0894177953998444436, 0.0173657744921376159, 0.00156735750354995707, 5.64214640518901565E-05, 5.97541959792059611E-07, 8.58964989963318053E-10 }; double x16[16] = { -6.63087819839312953, -5.47222570594934332, -4.49295530252001196, -3.60087362417154866, -2.76024504763070189, -1.95198034571633361, -1.1638291005549648, -0.386760604500557381, 0.386760604500557381, 1.1638291005549648, 1.95198034571633361, 2.76024504763070189, 3.60087362417154866, 4.49295530252001196, 5.47222570594934332, 6.63087819839312953 }; double w16[16] = { 1.49781472316183141E-10, 1.30947321628682029E-07, 1.53000321624872858E-05, 0.000525984926573909786, 0.0072669376011847411, 0.0472847523540140674, 0.158338372750949252, 0.286568521238012408, 0.286568521238012408, 0.158338372750949252, 0.0472847523540140674, 0.0072669376011847411, 0.000525984926573909786, 1.53000321624872858E-05, 1.30947321628682029E-07, 1.49781472316183141E-10 }; double x17[17] = { -6.88912243989533302, -5.74446007865940711, -4.77853158962998403, -3.90006571719801043, -3.07379717532819408, -2.28101944025298886, -1.50988330779674085, -0.751842600703896302, 0.00000000000000000, 0.751842600703896302, 1.50988330779674085, 2.28101944025298886, 3.07379717532819408, 3.90006571719801043, 4.77853158962998403, 5.74446007865940711, 6.88912243989533302 }; double w17[17] = { 2.58431491937491514E-11, 2.80801611793057831E-08, 4.0126794479798725E-06, 0.000168491431551339447, 0.00285894606228464989, 0.023086657025711152, 0.0974063711627180806, 0.226706308468978768, 0.299538370126607556, 0.226706308468978768, 0.0974063711627180806, 0.023086657025711152, 0.00285894606228464989, 0.000168491431551339447, 4.0126794479798725E-06, 2.80801611793057831E-08, 2.58431491937491514E-11 }; double x18[18] = { -7.13946484914647961, -6.00774591135959746, -5.05407268544274046, -4.1880202316294044, -3.37473653577809074, -2.59583368891124033, -1.83977992150864567, -1.09839551809150127, -0.365245755507697667, 0.365245755507697667, 1.09839551809150127, 1.83977992150864567, 2.59583368891124033, 3.37473653577809074, 4.1880202316294044, 5.05407268544274046, 6.00774591135959746, 7.13946484914647961 }; double w18[18] = { 4.41658876935870775E-12, 5.90548847883654844E-09, 1.02155239763698159E-06, 5.17989614411619621E-05, 0.00106548479629164959, 0.0105165177519413525, 0.0548966324802226541, 0.160685303893512627, 0.272783234654287887, 0.272783234654287887, 0.160685303893512627, 0.0548966324802226541, 0.0105165177519413525, 0.00106548479629164959, 5.17989614411619621E-05, 1.02155239763698159E-06, 5.90548847883654844E-09, 4.41658876935870775E-12 }; double x19[19] = { -7.38257902403043165, -6.2628911565132519, -5.32053637733603857, -4.46587262683103159, -3.66441654745063827, -2.89805127651575356, -2.15550276131693508, -1.4288766760783731, -0.712085044042379933, 0.00000000000000000, 0.712085044042379933, 1.4288766760783731, 2.15550276131693508, 2.89805127651575356, 3.66441654745063827, 4.46587262683103159, 5.32053637733603857, 6.2628911565132519, 7.38257902403043165 }; double w19[19] = { 7.4828300540572308E-13, 1.22037084844747862E-09, 2.53222003209286807E-07, 1.53511459546667444E-05, 0.000378502109414267593, 0.00450723542034203555, 0.0286666910301184956, 0.103603657276143998, 0.220941712199143658, 0.283773192751521075, 0.220941712199143658, 0.103603657276143998, 0.0286666910301184956, 0.00450723542034203555, 0.000378502109414267593, 1.53511459546667444E-05, 2.53222003209286807E-07, 1.22037084844747862E-09, 7.4828300540572308E-13 }; double x20[20] = { -7.61904854167975909, -6.51059015701365507, -5.57873880589320148, -4.73458133404605519, -3.9439673506573163, -3.18901481655339003, -2.45866361117236787, -1.74524732081412703, -1.04294534880275092, -0.346964157081355917, 0.346964157081355917, 1.04294534880275092, 1.74524732081412703, 2.45866361117236787, 3.18901481655339003, 3.9439673506573163, 4.73458133404605519, 5.57873880589320148, 6.51059015701365507, 7.61904854167975909 }; double w20[20] = { 1.25780067243793047E-13, 2.4820623623151838E-10, 6.12749025998295973E-08, 4.40212109023086457E-06, 0.000128826279961928981, 0.00183010313108049175, 0.0139978374471010428, 0.0615063720639760295, 0.161739333984000255, 0.26079306344955544, 0.26079306344955544, 0.161739333984000255, 0.0615063720639760295, 0.0139978374471010428, 0.00183010313108049175, 0.000128826279961928981, 4.40212109023086457E-06, 6.12749025998295973E-08, 2.4820623623151838E-10, 1.25780067243793047E-13 }; double x21[21] = { -7.84938289511382248, -6.75144471871746088, -5.82938200730447065, -4.99496394478202532, -4.21434398168842161, -3.46984669047537642, -2.75059298105237326, -2.0491024682571628, -1.35976582321123041, -0.678045692440644054, 0.00000000000000000, 0.678045692440644054, 1.35976582321123041, 2.0491024682571628, 2.75059298105237326, 3.46984669047537642, 4.21434398168842161, 4.99496394478202532, 5.82938200730447065, 6.75144471871746088, 7.84938289511382248 }; double w21[21] = { 2.09899121956566525E-14, 4.97536860412174643E-11, 1.45066128449307397E-08, 1.22535483614825217E-06, 4.21923474255158655E-05, 0.000708047795481537364, 0.00643969705140877684, 0.0339527297865428387, 0.10839228562641938, 0.215333715695059824, 0.270260183572877066, 0.215333715695059824, 0.10839228562641938, 0.0339527297865428387, 0.00643969705140877684, 0.000708047795481537364, 4.21923474255158655E-05, 1.22535483614825217E-06, 1.45066128449307397E-08, 4.97536860412174643E-11, 2.09899121956566525E-14 }; double x22[22] = { -8.07402998402171157, -6.98598042401881525, -6.07307495112289786, -5.2477244337144251, -4.47636197731086849, -3.74149635026651772, -3.03240422783167629, -2.34175999628770803, -1.6641248391179071, -0.995162422271215541, -0.331179315715273814, 0.331179315715273814, 0.995162422271215541, 1.6641248391179071, 2.34175999628770803, 3.03240422783167629, 3.74149635026651772, 4.47636197731086849, 5.2477244337144251, 6.07307495112289786, 6.98598042401881525, 8.07402998402171157 }; double w22[22] = { 3.47946064787714279E-15, 9.84137898234601051E-12, 3.36651415945821088E-09, 3.31985374981400429E-07, 1.33459771268087124E-05, 0.000262283303255964159, 0.00280876104757721073, 0.0175690728808057736, 0.0671963114288898905, 0.161906293413675378, 0.250243596586935013, 0.250243596586935013, 0.161906293413675378, 0.0671963114288898905, 0.0175690728808057736, 0.00280876104757721073, 0.000262283303255964159, 1.33459771268087124E-05, 3.31985374981400429E-07, 3.36651415945821088E-09, 9.84137898234601051E-12, 3.47946064787714279E-15 }; double x23[23] = { -8.29338602741735365, -7.21465943505186225, -6.31034985444839958, -5.49347398647179475, -4.73072419745147332, -4.00477532173330442, -3.30504002175296518, -2.62432363405918201, -1.9573275529334242, -1.29987646830397896, -0.648471153534495803, 0.00000000000000000, 0.648471153534495803, 1.29987646830397896, 1.9573275529334242, 2.62432363405918201, 3.30504002175296518, 4.00477532173330442, 4.73072419745147332, 5.49347398647179475, 6.31034985444839958, 7.21465943505186225, 8.29338602741735365 }; double w23[23] = { 5.73238316780208728E-16, 1.92293531156779128E-12, 7.67088886239990765E-10, 8.77506248386171607E-08, 4.08997724499215494E-06, 9.34081860903129835E-05, 0.00116762863749786134, 0.00857967839146566401, 0.0388671837034809467, 0.112073382602620911, 0.209959669577542613, 0.258509740808839039, 0.209959669577542613, 0.112073382602620911, 0.0388671837034809467, 0.00857967839146566401, 0.00116762863749786134, 9.34081860903129835E-05, 4.08997724499215494E-06, 8.77506248386171607E-08, 7.67088886239990765E-10, 1.92293531156779128E-12, 5.73238316780208728E-16 }; double x24[24] = { -8.50780351919525835, -7.43789066602166304, -6.54167500509863409, -5.73274717525120092, -4.97804137463912078, -4.26038360501990532, -3.56930676407356096, -2.89772864322331403, -2.24046785169175244, -1.59348042981642024, -0.953421922932109256, -0.317370096629452314, 0.317370096629452314, 0.953421922932109256, 1.59348042981642024, 2.24046785169175244, 2.89772864322331403, 3.56930676407356096, 4.26038360501990532, 4.97804137463912078, 5.73274717525120092, 6.54167500509863409, 7.43789066602166304, 8.50780351919525835 }; double w24[24] = { 9.39019368904192022E-17, 3.71497415276241595E-13, 1.71866492796486901E-10, 2.26746167348046514E-08, 1.21765974544258296E-06, 3.20950056527459886E-05, 0.000464718718779397633, 0.00397660892918131129, 0.0211263444089670287, 0.0720693640171784361, 0.161459512867000249, 0.240870115546640562, 0.240870115546640562, 0.161459512867000249, 0.0720693640171784361, 0.0211263444089670287, 0.00397660892918131129, 0.000464718718779397633, 3.20950056527459886E-05, 1.21765974544258296E-06, 2.26746167348046514E-08, 1.71866492796486901E-10, 3.71497415276241595E-13, 9.39019368904192022E-17 }; double x25[25] = { -8.71759767839958855, -7.65603795539307619, -6.76746496380971685, -5.96601469060670198, -5.21884809364427937, -4.50892992296728501, -3.82590056997249173, -3.16277567938819271, -2.51447330395220581, -1.8770583699478387, -1.24731197561678919, -0.622462279186076106, 0.00000000000000000, 0.622462279186076106, 1.24731197561678919, 1.8770583699478387, 2.51447330395220581, 3.16277567938819271, 3.82590056997249173, 4.50892992296728501, 5.21884809364427937, 5.96601469060670198, 6.76746496380971685, 7.65603795539307619, 8.71759767839958855 }; double w25[25] = { 1.53003899799868247E-17, 7.10210303700392527E-14, 3.79115000047718706E-11, 5.7380238688993763E-09, 3.53015256024549785E-07, 1.06721949052025363E-05, 0.0001777669069265266, 0.00175785040526379608, 0.0108567559914623159, 0.0433799701676449712, 0.114880924303951637, 0.204851025650340413, 0.248169351176485475, 0.204851025650340413, 0.114880924303951637, 0.0433799701676449712, 0.0108567559914623159, 0.00175785040526379608, 0.0001777669069265266, 1.06721949052025363E-05, 3.53015256024549785E-07, 5.7380238688993763E-09, 3.79115000047718706E-11, 7.10210303700392527E-14, 1.53003899799868247E-17 }; if ( n == 1 ) { r8vec_copy ( n, x01, x ); r8vec_copy ( n, w01, w ); } else if ( n == 2 ) { r8vec_copy ( n, x02, x ); r8vec_copy ( n, w02, w ); } else if ( n == 3 ) { r8vec_copy ( n, x03, x ); r8vec_copy ( n, w03, w ); } else if ( n == 4 ) { r8vec_copy ( n, x04, x ); r8vec_copy ( n, w04, w ); } else if ( n == 5 ) { r8vec_copy ( n, x05, x ); r8vec_copy ( n, w05, w ); } else if ( n == 6 ) { r8vec_copy ( n, x06, x ); r8vec_copy ( n, w06, w ); } else if ( n == 7 ) { r8vec_copy ( n, x07, x ); r8vec_copy ( n, w07, w ); } else if ( n == 8 ) { r8vec_copy ( n, x08, x ); r8vec_copy ( n, w08, w ); } else if ( n == 9 ) { r8vec_copy ( n, x09, x ); r8vec_copy ( n, w09, w ); } else if ( n == 10 ) { r8vec_copy ( n, x10, x ); r8vec_copy ( n, w10, w ); } else if ( n == 11 ) { r8vec_copy ( n, x11, x ); r8vec_copy ( n, w11, w ); } else if ( n == 12 ) { r8vec_copy ( n, x12, x ); r8vec_copy ( n, w12, w ); } else if ( n == 13 ) { r8vec_copy ( n, x13, x ); r8vec_copy ( n, w13, w ); } else if ( n == 14 ) { r8vec_copy ( n, x14, x ); r8vec_copy ( n, w14, w ); } else if ( n == 15 ) { r8vec_copy ( n, x15, x ); r8vec_copy ( n, w15, w ); } else if ( n == 16 ) { r8vec_copy ( n, x16, x ); r8vec_copy ( n, w16, w ); } else if ( n == 17 ) { r8vec_copy ( n, x17, x ); r8vec_copy ( n, w17, w ); } else if ( n == 18 ) { r8vec_copy ( n, x18, x ); r8vec_copy ( n, w18, w ); } else if ( n == 19 ) { r8vec_copy ( n, x19, x ); r8vec_copy ( n, w19, w ); } else if ( n == 20 ) { r8vec_copy ( n, x20, x ); r8vec_copy ( n, w20, w ); } else if ( n == 21 ) { r8vec_copy ( n, x21, x ); r8vec_copy ( n, w21, w ); } else if ( n == 22 ) { r8vec_copy ( n, x22, x ); r8vec_copy ( n, w22, w ); } else if ( n == 23 ) { r8vec_copy ( n, x23, x ); r8vec_copy ( n, w23, w ); } else if ( n == 24 ) { r8vec_copy ( n, x24, x ); r8vec_copy ( n, w24, w ); } else if ( n == 25 ) { r8vec_copy ( n, x25, x ); r8vec_copy ( n, w25, w ); } else { fprintf ( stderr, "\n" ); fprintf ( stderr, "GQN - Fatal error!\n" ); fprintf ( stderr, " Value of N must be between 1 and 25.\n" ); exit ( 1 ); } return; } /******************************************************************************/ int gqn_order ( int l ) /******************************************************************************/ /* Purpose: GQN_ORDER computes the order of a GQN rule from the level. Licensing: This code is distributed under the MIT license. Modified: 08 December 2012 Author: John Burkardt. Parameters: Input, int L, the level of the rule. 1 <= L. Output, int GQN_ORDER, the order of the rule. */ { int n; if ( l < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GQN_ORDER - Fatal error!\n" ); fprintf ( stderr, " 1 <= L required.\n" ); fprintf ( stderr, " Input L = %d\n", l ); exit ( 1 ); } else if ( l <= 25 ) { n = l; } else { fprintf ( stderr, "\n" ); fprintf ( stderr, "GQN_ORDER - Fatal error!\n" ); fprintf ( stderr, " L <= 25 required.\n" ); fprintf ( stderr, " Input L = %d\n", l ); exit ( 1 ); } return n; } /******************************************************************************/ int gqn2_order ( int l ) /******************************************************************************/ /* Purpose: GQN2_ORDER computes the order of a GQN rule from the level. Discussion: For this version of the order routine, we have n = 2 * l - 1 Licensing: This code is distributed under the MIT license. Modified: 07 February 2014 Author: John Burkardt. Parameters: Input, int L, the level of the rule. 1 <= L. Output, int GQN2_ORDER, the order of the rule. */ { int n; if ( l < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GQN2_ORDER - Fatal error!\n" ); fprintf ( stderr, " 1 <= L required.\n" ); fprintf ( stderr, " Input L = %d\n", l ); exit ( 1 ); } else if ( l <= 13 ) { n = 2 * l - 1; } else { fprintf ( stderr, "\n" ); fprintf ( stderr, "GQN2_ORDER - Fatal error!\n" ); fprintf ( stderr, " L <= 13 required.\n" ); fprintf ( stderr, " Input L = %d\n", l ); exit ( 1 ); } return n; } /******************************************************************************/ void gqu ( int n, double x[], double w[] ) /******************************************************************************/ /* Purpose: GQU provides data for Gauss quadrature with a uniform weight. Discussion: This data assumes integration over the interval [0,1] with weight function w(x) = 1. Licensing: This code is distributed under the MIT license. Modified: 10 December 2012 Author: John Burkardt. Reference: Florian Heiss, Viktor Winschel, Likelihood approximation by numerical integration on sparse grids, Journal of Econometrics, Volume 144, 2008, pages 62-80. Parameters: Input, int N, the number of points and weights. 1 <= N <= 25. Output, double X[N], the nodes. Output, double W[N], the weights. */ { double x01[1] = { 0.500000000000000000 }; double w01[1] = { 1.000000000000000000 }; double x02[2] = { 0.211324865405187134, 0.788675134594812866 }; double w02[2] = { 0.500000000000000000, 0.500000000000000000 }; double x03[3] = { 0.112701665379258298, 0.500000000000000000, 0.887298334620741702 }; double w03[3] = { 0.277777777777777124, 0.444444444444445697, 0.277777777777777124 }; double x04[4] = { 0.0694318442029737692, 0.330009478207571871, 0.669990521792428129, 0.930568155797026231 }; double w04[4] = { 0.173927422568724843, 0.326072577431275157, 0.326072577431275157, 0.173927422568724843 }; double x05[5] = { 0.0469100770306680737, 0.230765344947158502, 0.500000000000000000, 0.769234655052841498, 0.953089922969331926 }; double w05[5] = { 0.118463442528091739, 0.239314335249685012, 0.284444444444446554, 0.239314335249685012, 0.118463442528091739 }; double x06[6] = { 0.033765242898423975, 0.169395306766867648, 0.380690406958401506, 0.619309593041598494, 0.830604693233132352, 0.966234757101576025 }; double w06[6] = { 0.0856622461895818338, 0.180380786524070719, 0.233956967286347461, 0.233956967286347461, 0.180380786524070719, 0.0856622461895818338 }; double x07[7] = { 0.0254460438286208124, 0.12923440720030277, 0.297077424311301463, 0.500000000000000000, 0.702922575688698537, 0.87076559279969723, 0.974553956171379188 }; double w07[7] = { 0.0647424830844317012, 0.139852695744639349, 0.190915025252560905, 0.2089795918367362, 0.190915025252560905, 0.139852695744639349, 0.0647424830844317012 }; double x08[8] = { 0.0198550717512319119, 0.101666761293186525, 0.237233795041835505, 0.408282678752175054, 0.591717321247824946, 0.762766204958164495, 0.898333238706813475, 0.980144928248768088 }; double w08[8] = { 0.0506142681451851803, 0.111190517226687935, 0.156853322938944689, 0.181341891689182133, 0.181341891689182133, 0.156853322938944689, 0.111190517226687935, 0.0506142681451851803 }; double x09[9] = { 0.0159198802461869571, 0.0819844463366821152, 0.193314283649704821, 0.337873288298095487, 0.500000000000000000, 0.662126711701904513, 0.806685716350295179, 0.918015553663317885, 0.984080119753813043 }; double w09[9] = { 0.0406371941807845832, 0.0903240803474292531, 0.130305348201468441, 0.156173538520002264, 0.165119677500630752, 0.156173538520002264, 0.130305348201468441, 0.0903240803474292531, 0.0406371941807845832 }; double x10[10] = { 0.0130467357414141283, 0.0674683166555076763, 0.160295215850487782, 0.283302302935376393, 0.42556283050918442, 0.57443716949081558, 0.716697697064623607, 0.839704784149512218, 0.932531683344492324, 0.986953264258585872 }; double w10[10] = { 0.0333356721543420012, 0.0747256745752905988, 0.109543181257991576, 0.13463335965499873, 0.147762112357377129, 0.147762112357377129, 0.13463335965499873, 0.109543181257991576, 0.0747256745752905988, 0.0333356721543420012 }; double x11[11] = { 0.010885670926971569, 0.0564687001159522861, 0.134923997212975322, 0.240451935396594152, 0.365228422023827548, 0.500000000000000000, 0.634771577976172452, 0.759548064603405848, 0.865076002787024678, 0.943531299884047714, 0.989114329073028431 }; double w11[11] = { 0.0278342835580849164, 0.0627901847324526252, 0.0931451054638675197, 0.11659688229599563, 0.131402272255123881, 0.136462543388950863, 0.131402272255123881, 0.11659688229599563, 0.0931451054638675197, 0.0627901847324526252, 0.0278342835580849164 }; double x12[12] = { 0.00921968287664043373, 0.0479413718147625456, 0.115048662902847654, 0.206341022856691314, 0.316084250500909936, 0.437383295744265488, 0.562616704255734512, 0.683915749499090064, 0.793658977143308686, 0.884951337097152346, 0.952058628185237454, 0.990780317123359566 }; double w12[12] = { 0.0235876681932543145, 0.0534696629976592758, 0.0800391642716734436, 0.101583713361533282, 0.116746268269177805, 0.124573522906701886, 0.124573522906701886, 0.116746268269177805, 0.101583713361533282, 0.0800391642716734436, 0.0534696629976592758, 0.0235876681932543145 }; double x13[13] = { 0.00790847264070593248, 0.0412008003885109275, 0.0992109546333450609, 0.178825330279829942, 0.275753624481776649, 0.384770842022432613, 0.500000000000000000, 0.615229157977567387, 0.724246375518223351, 0.821174669720170058, 0.900789045366654939, 0.958799199611489072, 0.992091527359294068 }; double w13[13] = { 0.0202420023826562281, 0.0460607499188643785, 0.0694367551098938746, 0.0890729903809732021, 0.103908023768444616, 0.113141590131449032, 0.116275776615437407, 0.113141590131449032, 0.103908023768444616, 0.0890729903809732021, 0.0694367551098938746, 0.0460607499188643785, 0.0202420023826562281 }; double x14[14] = { 0.00685809565159378742, 0.0357825581682131855, 0.0863993424651174902, 0.156353547594157316, 0.24237568182092295, 0.340443815536055183, 0.445972525646328166, 0.554027474353671834, 0.659556184463944817, 0.75762431817907705, 0.843646452405842684, 0.91360065753488251, 0.964217441831786815, 0.993141904348406213 }; double w14[14] = { 0.0175597301658745736, 0.0400790435798802913, 0.0607592853439517105, 0.0786015835790969952, 0.0927691987389691608, 0.102599231860648107, 0.107631926731579161, 0.107631926731579161, 0.102599231860648107, 0.0927691987389691608, 0.0786015835790969952, 0.0607592853439517105, 0.0400790435798802913, 0.0175597301658745736 }; double x15[15] = { 0.0060037409897573113, 0.0313633037996470243, 0.0758967082947863414, 0.137791134319914965, 0.214513913695730585, 0.302924326461218252, 0.399402953001282701, 0.500000000000000000, 0.600597046998717299, 0.697075673538781748, 0.785486086304269415, 0.862208865680085035, 0.924103291705213659, 0.968636696200352976, 0.993996259010242689 }; double w15[15] = { 0.0153766209980574341, 0.0351830237440541593, 0.0535796102335861571, 0.0697853389630773147, 0.083134602908497196, 0.0930805000077812861, 0.0992157426635560391, 0.101289120962780907, 0.0992157426635560391, 0.0930805000077812861, 0.083134602908497196, 0.0697853389630773147, 0.0535796102335861571, 0.0351830237440541593, 0.0153766209980574341 }; double x16[16] = { 0.00529953250417503074, 0.0277124884633836999, 0.0671843988060840669, 0.122297795822498445, 0.191061877798678115, 0.270991611171386371, 0.359198224610370542, 0.452493745081181231, 0.547506254918818769, 0.640801775389629458, 0.729008388828613629, 0.808938122201321885, 0.877702204177501555, 0.932815601193915933, 0.9722875115366163, 0.994700467495824969 }; double w16[16] = { 0.0135762297058759553, 0.0311267619693239538, 0.0475792558412465455, 0.062314485627767105, 0.0747979944082885623, 0.0845782596975014622, 0.0913017075224620001, 0.0947253052275344315, 0.0947253052275344315, 0.0913017075224620001, 0.0845782596975014622, 0.0747979944082885623, 0.062314485627767105, 0.0475792558412465455, 0.0311267619693239538, 0.0135762297058759553 }; double x17[17] = { 0.00471226234279131795, 0.024662239115616158, 0.0598804231365070994, 0.109242998051599316, 0.171164420391654692, 0.243654731456761531, 0.324384118273061794, 0.410757909252076114, 0.500000000000000000, 0.589242090747923886, 0.675615881726938206, 0.756345268543238469, 0.828835579608345308, 0.890757001948400684, 0.940119576863492901, 0.975337760884383842, 0.995287737657208682 }; double w17[17] = { 0.01207415143427314, 0.0277297646869936118, 0.0425180741585896443, 0.0559419235967020534, 0.0675681842342628902, 0.0770228805384053083, 0.0840020510782251428, 0.0882813526834964474, 0.0897232351781034193, 0.0882813526834964474, 0.0840020510782251428, 0.0770228805384053083, 0.0675681842342628902, 0.0559419235967020534, 0.0425180741585896443, 0.0277297646869936118, 0.01207415143427314 }; double x18[18] = { 0.00421741578953449547, 0.022088025214301199, 0.0536987667512220934, 0.0981475205137384288, 0.154156478469823388, 0.22011458446302623, 0.294124419268578685, 0.374056887154247231, 0.457612493479132354, 0.542387506520867646, 0.625943112845752769, 0.705875580731421315, 0.77988541553697377, 0.845843521530176612, 0.901852479486261571, 0.946301233248777907, 0.977911974785698801, 0.995782584210465505 }; double w18[18] = { 0.0108080067632407191, 0.0248572744474849679, 0.0382128651274446646, 0.0504710220531437159, 0.061277603355739306, 0.0703214573353254518, 0.0773423375631328014, 0.0821382418729165037, 0.084571191481571939, 0.084571191481571939, 0.0821382418729165037, 0.0773423375631328014, 0.0703214573353254518, 0.061277603355739306, 0.0504710220531437159, 0.0382128651274446646, 0.0248572744474849679, 0.0108080067632407191 }; double x19[19] = { 0.00379657807820787951, 0.0198959239325849913, 0.048422048192590994, 0.0886426717314285906, 0.139516911332385307, 0.19972734766915945, 0.267714629312019503, 0.341717950018185057, 0.419820677179887358, 0.500000000000000000, 0.580179322820112642, 0.658282049981814943, 0.732285370687980497, 0.80027265233084055, 0.860483088667614693, 0.911357328268571409, 0.951577951807409006, 0.980104076067415009, 0.99620342192179212 }; double w19[19] = { 0.00973089411486243415, 0.0224071133828498206, 0.0345222713688206687, 0.0457450108112251244, 0.0557833227736671128, 0.064376981269668232, 0.0713033510868034126, 0.0763830210329299597, 0.0794844216969773365, 0.0805272249243919463, 0.0794844216969773365, 0.0763830210329299597, 0.0713033510868034126, 0.064376981269668232, 0.0557833227736671128, 0.0457450108112251244, 0.0345222713688206687, 0.0224071133828498206, 0.00973089411486243415 }; double x20[20] = { 0.00343570040745255767, 0.0180140363610430398, 0.0438827858743370269, 0.0804415140888905533, 0.126834046769924602, 0.181973159636742432, 0.244566499024586381, 0.313146955642290226, 0.386107074429177466, 0.461736739433251331, 0.538263260566748669, 0.613892925570822534, 0.686853044357709774, 0.755433500975413619, 0.818026840363257568, 0.873165953230075398, 0.919558485911109447, 0.956117214125662973, 0.98198596363895696, 0.996564299592547442 }; double w20[20] = { 0.0088070035695753026, 0.0203007149001935561, 0.0313360241670545686, 0.0416383707883524329, 0.0509650599086203179, 0.0590972659807592476, 0.0658443192245883463, 0.0710480546591911871, 0.0745864932363019956, 0.0763766935653631129, 0.0763766935653631129, 0.0745864932363019956, 0.0710480546591911871, 0.0658443192245883463, 0.0590972659807592476, 0.0509650599086203179, 0.0416383707883524329, 0.0313360241670545686, 0.0203007149001935561, 0.0088070035695753026 }; double x21[21] = { 0.00312391468980521836, 0.0163865807168468436, 0.0399503329247996586, 0.0733183177083414073, 0.115780018262161111, 0.166430597901293886, 0.224190582056390086, 0.287828939896280556, 0.355989341598799469, 0.427219072919552412, 0.500000000000000000, 0.572780927080447588, 0.644010658401200531, 0.712171060103719444, 0.775809417943609914, 0.833569402098706114, 0.884219981737838889, 0.926681682291658593, 0.960049667075200341, 0.983613419283153156, 0.996876085310194782 }; double w21[21] = { 0.00800861412888644909, 0.018476894885426285, 0.0285672127134286406, 0.0380500568141897075, 0.0467222117280169935, 0.0543986495835743558, 0.06091570802686435, 0.0661344693166688452, 0.0699436973955366581, 0.0722622019949851341, 0.073040566824845346, 0.0722622019949851341, 0.0699436973955366581, 0.0661344693166688452, 0.06091570802686435, 0.0543986495835743558, 0.0467222117280169935, 0.0380500568141897075, 0.0285672127134286406, 0.018476894885426285,0.00800861412888644909 }; double x22[22] = { 0.0028527072588003799, 0.0149697510822857094, 0.036521613906413064, 0.0670937111398499653, 0.106091597010395944, 0.152756368406658627, 0.206179798246544199, 0.265322081006621469, 0.329032089553957907, 0.396069786655889322, 0.465130363340138908, 0.534869636659861092, 0.603930213344110678, 0.670967910446042093, 0.734677918993378531, 0.793820201753455801, 0.847243631593341373, 0.893908402989604056, 0.932906288860150035, 0.963478386093586936, 0.985030248917714291, 0.99714729274119962 }; double w22[22] = { 0.00731399764913532799, 0.0168874507924071104, 0.0261466675763416916, 0.0348982342122602998, 0.0429708031085339753, 0.0502070722214406004, 0.0564661480402697119, 0.0616261884052562506, 0.0655867523935313168, 0.0682707491730076971, 0.0696259364278161291, 0.0696259364278161291, 0.0682707491730076971, 0.0655867523935313168, 0.0616261884052562506, 0.0564661480402697119, 0.0502070722214406004, 0.0429708031085339753, 0.0348982342122602998, 0.0261466675763416916, 0.0168874507924071104, 0.00731399764913532799 }; double x23[23] = { 0.00261533250122392147, 0.013728764390942394, 0.0335144565869919253, 0.061623820864779244, 0.0975557991905799948, 0.140669318434024859, 0.190195062118176939, 0.245249261076996294, 0.304849480984854537, 0.367932159514827495, 0.433371587850766904, 0.500000000000000000, 0.566628412149233096, 0.632067840485172505, 0.695150519015145463, 0.754750738923003706, 0.809804937881823061, 0.859330681565975141, 0.902444200809420005, 0.938376179135220756, 0.966485543413008075, 0.986271235609057606, 0.997384667498776079 }; double w23[23] = { 0.0067059297435702412, 0.015494002928489686, 0.0240188358655423692, 0.0321162107042629943, 0.0396407058883595509, 0.0464578830300175633, 0.052446045732270824, 0.057498320111205814, 0.0615245421533648154, 0.06445286109404115, 0.0662310197023484037, 0.0668272860930531759, 0.0662310197023484037, 0.06445286109404115, 0.0615245421533648154, 0.057498320111205814, 0.052446045732270824, 0.0464578830300175633, 0.0396407058883595509, 0.0321162107042629943, 0.0240188358655423692, 0.015494002928489686, 0.0067059297435702412 }; double x24[24] = { 0.00240639000148934468, 0.0126357220143452631, 0.0308627239986336566, 0.0567922364977995198, 0.0899990070130485265, 0.129937904210722821, 0.175953174031512227, 0.227289264305580163, 0.283103246186977464, 0.342478660151918302, 0.404440566263191803, 0.467971553568697241, 0.532028446431302759, 0.595559433736808197, 0.657521339848081698, 0.716896753813022536, 0.772710735694419837, 0.824046825968487773, 0.870062095789277179, 0.910000992986951474, 0.94320776350220048, 0.969137276001366343, 0.987364277985654737, 0.997593609998510655 }; double w24[24] = { 0.00617061489999283508, 0.014265694314466934, 0.0221387194087098796, 0.0296492924577183847, 0.0366732407055402054, 0.0430950807659766927, 0.0488093260520570393, 0.0537221350579829143, 0.0577528340268628829, 0.0608352364639017928, 0.0629187281734143178, 0.0639690976733762462, 0.0639690976733762462, 0.0629187281734143178, 0.0608352364639017928, 0.0577528340268628829, 0.0537221350579829143, 0.0488093260520570393, 0.0430950807659766927, 0.0366732407055402054, 0.0296492924577183847, 0.0221387194087098796, 0.014265694314466934,0.00617061489999283508 }; double x25[25] = { 0.00222151510475088187, 0.0116680392702412927, 0.0285127143855128384, 0.052504001060862393, 0.0832786856195830705, 0.120370368481321099, 0.163216815763265854, 0.211168534879388581, 0.263498634277142485, 0.319413847095306069, 0.378066558139505737, 0.438567653694644788, 0.500000000000000000, 0.561432346305355212, 0.621933441860494263, 0.680586152904693931, 0.736501365722857515, 0.788831465120611419, 0.836783184236734146, 0.879629631518678901, 0.91672131438041693, 0.947495998939137607, 0.971487285614487162, 0.988331960729758707, 0.997778484895249118 }; double w25[25] = { 0.00569689925051255347, 0.0131774933075161083, 0.0204695783506531476, 0.0274523479879176906, 0.0340191669061785454, 0.0400703501675005319, 0.0455141309914819034, 0.0502679745335253628, 0.0542598122371318672, 0.0574291295728558623, 0.0597278817678924615, 0.0611212214951551217, 0.061588026863357799, 0.0611212214951551217, 0.0597278817678924615, 0.0574291295728558623, 0.0542598122371318672, 0.0502679745335253628, 0.0455141309914819034, 0.0400703501675005319, 0.0340191669061785454, 0.0274523479879176906, 0.0204695783506531476, 0.0131774933075161083, 0.00569689925051255347 }; if ( n == 1 ) { r8vec_copy ( n, x01, x ); r8vec_copy ( n, w01, w ); } else if ( n == 2 ) { r8vec_copy ( n, x02, x ); r8vec_copy ( n, w02, w ); } else if ( n == 3 ) { r8vec_copy ( n, x03, x ); r8vec_copy ( n, w03, w ); } else if ( n == 4 ) { r8vec_copy ( n, x04, x ); r8vec_copy ( n, w04, w ); } else if ( n == 5 ) { r8vec_copy ( n, x05, x ); r8vec_copy ( n, w05, w ); } else if ( n == 6 ) { r8vec_copy ( n, x06, x ); r8vec_copy ( n, w06, w ); } else if ( n == 7 ) { r8vec_copy ( n, x07, x ); r8vec_copy ( n, w07, w ); } else if ( n == 8 ) { r8vec_copy ( n, x08, x ); r8vec_copy ( n, w08, w ); } else if ( n == 9 ) { r8vec_copy ( n, x09, x ); r8vec_copy ( n, w09, w ); } else if ( n == 10 ) { r8vec_copy ( n, x10, x ); r8vec_copy ( n, w10, w ); } else if ( n == 11 ) { r8vec_copy ( n, x11, x ); r8vec_copy ( n, w11, w ); } else if ( n == 12 ) { r8vec_copy ( n, x12, x ); r8vec_copy ( n, w12, w ); } else if ( n == 13 ) { r8vec_copy ( n, x13, x ); r8vec_copy ( n, w13, w ); } else if ( n == 14 ) { r8vec_copy ( n, x14, x ); r8vec_copy ( n, w14, w ); } else if ( n == 15 ) { r8vec_copy ( n, x15, x ); r8vec_copy ( n, w15, w ); } else if ( n == 16 ) { r8vec_copy ( n, x16, x ); r8vec_copy ( n, w16, w ); } else if ( n == 17 ) { r8vec_copy ( n, x17, x ); r8vec_copy ( n, w17, w ); } else if ( n == 18 ) { r8vec_copy ( n, x18, x ); r8vec_copy ( n, w18, w ); } else if ( n == 19 ) { r8vec_copy ( n, x19, x ); r8vec_copy ( n, w19, w ); } else if ( n == 20 ) { r8vec_copy ( n, x20, x ); r8vec_copy ( n, w20, w ); } else if ( n == 21 ) { r8vec_copy ( n, x21, x ); r8vec_copy ( n, w21, w ); } else if ( n == 22 ) { r8vec_copy ( n, x22, x ); r8vec_copy ( n, w22, w ); } else if ( n == 23 ) { r8vec_copy ( n, x23, x ); r8vec_copy ( n, w23, w ); } else if ( n == 24 ) { r8vec_copy ( n, x24, x ); r8vec_copy ( n, w24, w ); } else if ( n == 25 ) { r8vec_copy ( n, x25, x ); r8vec_copy ( n, w25, w ); } else { fprintf ( stderr, "\n" ); fprintf ( stderr, "GQU - Fatal error!\n" ); fprintf ( stderr, " Value of N must be between 1 and 25.\n" ); exit ( 1 ); } return; } /******************************************************************************/ int gqu_order ( int l ) /******************************************************************************/ /* Purpose: GQU_ORDER computes the order of a GQU rule from the level. Licensing: This code is distributed under the MIT license. Modified: 08 December 2012 Author: John Burkardt. Parameters: Input, int L, the level of the rule. 1 <= L <= 25. Output, int GQU_ORDER, the order of the rule. */ { int n; if ( l < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GQU_ORDER - Fatal error!\n" ); fprintf ( stderr, " L <= 25 required.\n" ); fprintf ( stderr, " Input L = %d\n", l ); exit ( 1 ); } else if ( l <= 25 ) { n = l; } else { fprintf ( stderr, "\n" ); fprintf ( stderr, "GQU_ORDER - Fatal error!\n" ); fprintf ( stderr, " L <= 25 required.\n" ); fprintf ( stderr, " Input L = %d\n", l ); exit ( 1 ); } return n; } /******************************************************************************/ int i4_choose ( int n, int k ) /******************************************************************************/ /* Purpose: I4_CHOOSE computes the binomial coefficient C(N,K). Discussion: The value is calculated in such a way as to avoid overflow and roundoff. The calculation is done in integer arithmetic. The formula used is: C(N,K) = N! / ( K! * (N-K)! ) Licensing: This code is distributed under the MIT license. Modified: 03 May 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. Parameters: Input, int N, K, are the values of N and K. Output, int I4_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; } else if ( mn == 0 ) { value = 1; } else { mx = i4_max ( k, n - k ); value = mx + 1; for ( i = 2; i <= mn; i++ ) { value = ( value * ( mx + i ) ) / i; } } return value; } /******************************************************************************/ int i4_factorial2 ( int n ) /******************************************************************************/ /* Purpose: I4_FACTORIAL2 computes the double factorial function. Discussion: FACTORIAL2( N ) = Product ( N * (N-2) * (N-4) * ... * 2 ) (N even) = Product ( N * (N-2) * (N-4) * ... * 1 ) (N odd) Licensing: This code is distributed under the MIT license. Modified: 04 March 2008 Author: John Burkardt Parameters: Input, int N, the argument of the double factorial function. If N is less than 1, I4_FACTORIAL2 is returned as 1. Output, int I4_FACTORIAL2, the value of the double factorial function. */ { int n_copy; int value; if ( n < 1 ) { value = 1; return value; } n_copy = n; value = 1; while ( 1 < n_copy ) { value = value * n_copy; n_copy = n_copy - 2; } return value; } /******************************************************************************/ 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 Parameters: 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 Parameters: 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; } /******************************************************************************/ int i4_mop ( int i ) /******************************************************************************/ /* Purpose: I4_MOP returns the I-th power of -1 as an I4 value. Discussion: An I4 is an int value. Licensing: This code is distributed under the MIT license. Modified: 16 November 2007 Author: John Burkardt Parameters: Input, int I, the power of -1. Output, int I4_MOP, the I-th power of -1. */ { int value; if ( ( i % 2 ) == 0 ) { value = 1; } else { value = -1; } return value; } /******************************************************************************/ int i4_power ( int i, int j ) /******************************************************************************/ /* Purpose: I4_POWER returns the value of I^J. Licensing: This code is distributed under the MIT license. Modified: 23 October 2007 Author: John Burkardt Parameters: Input, int I, J, the base and the power. J should be nonnegative. Output, int I4_POWER, the value of I^J. */ { int k; int value; if ( j < 0 ) { if ( i == 1 ) { value = 1; } else if ( i == 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4_POWER - Fatal error!\n" ); fprintf ( stderr, " I^J requested, with I = 0 and J negative.\n" ); exit ( 1 ); } else { value = 0; } } else if ( j == 0 ) { if ( i == 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4_POWER - Fatal error!\n" ); fprintf ( stderr, " I^J requested, with I = 0 and J = 0.\n" ); exit ( 1 ); } else { value = 1; } } else if ( j == 1 ) { value = i; } else { value = 1; for ( k = 1; k <= j; k++ ) { value = value * i; } } return value; } /******************************************************************************/ void i4mat_print ( int m, int n, int a[], char *title ) /******************************************************************************/ /* Purpose: I4MAT_PRINT prints an I4MAT. Discussion: An I4MAT is an MxN array of I4's, stored by (I,J) -> [I+J*M]. Licensing: This code is distributed under the MIT license. Modified: 28 May 2008 Author: John Burkardt Parameters: Input, int M, the number of rows in A. Input, int N, the number of columns in A. Input, int A[M*N], the M by N matrix. Input, char *TITLE, a title. */ { i4mat_print_some ( m, n, a, 1, 1, m, n, title ); return; } /******************************************************************************/ void i4mat_print_some ( int m, int n, int a[], int ilo, int jlo, int ihi, int jhi, char *title ) /******************************************************************************/ /* Purpose: I4MAT_PRINT_SOME prints some of an I4MAT. Discussion: An I4MAT is an MxN array of I4's, stored by (I,J) -> [I+J*M]. Licensing: This code is distributed under the MIT license. Modified: 20 August 2010 Author: John Burkardt Parameters: Input, int M, the number of rows of the matrix. M must be positive. Input, int N, the number of columns of the matrix. N must be positive. Input, int A[M*N], the matrix. Input, int ILO, JLO, IHI, JHI, designate the first row and column, and the last row and column to be printed. Input, char *TITLE, a title. */ { # define INCX 10 int i; int i2hi; int i2lo; int j; int j2hi; int j2lo; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); if ( m <= 0 || n <= 0 ) { fprintf ( stdout, "\n" ); fprintf ( stdout, " (None)\n" ); return; } /* Print the columns of the matrix, in strips of INCX. */ for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX ) { j2hi = j2lo + INCX - 1; j2hi = i4_min ( j2hi, n ); j2hi = i4_min ( j2hi, jhi ); fprintf ( stdout, "\n" ); /* For each column J in the current range... Write the header. */ fprintf ( stdout, " Col:" ); for ( j = j2lo; j <= j2hi; j++ ) { fprintf ( stdout, " %6d", j - 1 ); } fprintf ( stdout, "\n" ); fprintf ( stdout, " Row\n" ); fprintf ( stdout, "\n" ); /* Determine the range of the rows in this strip. */ i2lo = i4_max ( ilo, 1 ); i2hi = i4_min ( ihi, m ); for ( i = i2lo; i <= i2hi; i++ ) { /* Print out (up to INCX) entries in row I, that lie in the current strip. */ fprintf ( stdout, "%5d:", i - 1 ); for ( j = j2lo; j <= j2hi; j++ ) { fprintf ( stdout, " %6d", a[i-1+(j-1)*m] ); } fprintf ( stdout, "\n" ); } } return; # undef INCX } /******************************************************************************/ int *i4vec_cum0_new ( int n, int a[] ) /******************************************************************************/ /* Purpose: I4VEC_CUM0_NEW computes the cumulutive sum of the entries of an I4VEC. Discussion: An I4VEC is a vector of I4's. This routine returns a vector of length N+1, with the first value being 0. Example: Input: A = (/ 1, 2, 3, 4 /) Output: A_CUM = (/ 0, 1, 3, 6, 10 /) Licensing: This code is distributed under the MIT license. Modified: 22 December 2010 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_CUM0_NEW[N+1], the cumulative sum of the entries of A. */ { int *a_cum; int i; a_cum = ( int * ) malloc ( ( n + 1 ) * sizeof ( int ) ); a_cum[0] = 0; for ( i = 1; i <= n; i++ ) { a_cum[i] = a_cum[i-1] + a[i-1]; } return a_cum; } /******************************************************************************/ void i4vec_print ( int n, int a[], char *title ) /******************************************************************************/ /* Purpose: I4VEC_PRINT prints an I4VEC. Discussion: An I4VEC is a vector of I4's. Licensing: This code is distributed under the MIT license. Modified: 14 November 2003 Author: John Burkardt Parameters: Input, int N, the number of components of the vector. Input, int A[N], the vector to be printed. Input, char *TITLE, a title. */ { int i; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); fprintf ( stdout, "\n" ); for ( i = 0; i < n; i++ ) { fprintf ( stdout, " %6d: %8d\n", i, a[i] ); } return; } /******************************************************************************/ int i4vec_product ( int n, int a[] ) /******************************************************************************/ /* Purpose: I4VEC_PRODUCT multiplies the entries of an I4VEC. Discussion: An I4VEC is a vector of I4's. Example: Input: A = ( 1, 2, 3, 4 ) Output: I4VEC_PRODUCT = 24 Licensing: This code is distributed under the MIT license. Modified: 26 August 2008 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Input, int A[N], the vector Output, int I4VEC_PRODUCT, the product of the entries of A. */ { int i; int product; product = 1; for ( i = 0; i < n; i++ ) { product = product * a[i]; } return product; } /******************************************************************************/ 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; } /******************************************************************************/ void i4vec_transpose_print ( int n, int a[], char *title ) /******************************************************************************/ /* Purpose: I4VEC_TRANSPOSE_PRINT prints an I4VEC "transposed". Discussion: An I4VEC is a vector of I4's. Example: A = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 } TITLE = "My vector: " My vector: 1 2 3 4 5 6 7 8 9 10 11 Licensing: This code is distributed under the MIT license. Modified: 25 July 2010 Author: John Burkardt Parameters: Input, int N, the number of components of the vector. Input, int A[N], the vector to be printed. Input, char *TITLE, a title. */ { int i; int ihi; int ilo; int title_len; title_len = strlen ( title ); for ( ilo = 1; ilo <= n; ilo = ilo + 5 ) { ihi = i4_min ( ilo + 5 - 1, n ); if ( ilo == 1 ) { printf ( "%s", title ); } else { for ( i = 1; i <= title_len; i++ ) { printf ( " " ); } } for ( i = ilo; i <= ihi; i++ ) { printf ( "%12d", a[i-1] ); } printf ( "\n" ); } return; } /******************************************************************************/ void kpn ( int n, double x[], double w[] ) /******************************************************************************/ /* Purpose: KPN provides data for Kronrod-Patterson quadrature with a normal weight. Discussion: This data assumes integration over the interval (-oo,+oo) with weight function w(x) = exp(-x*x/2)/sqrt(2*pi). Licensing: This code is distributed under the MIT license. Modified: 10 December 2012 Author: Original MATLAB version by Florian Heiss, Viktor Winschel. C version by John Burkardt. Reference: Florian Heiss, Viktor Winschel, Likelihood approximation by numerical integration on sparse grids, Journal of Econometrics, Volume 144, 2008, pages 62-80. Alan Genz, Bradley Keister, Fully symmetric interpolatory rules for multiple integrals over infinite regions with Gaussian weight, Journal of Computational and Applied Mathematics, Volume 71, 1996, pages 299-309. Thomas Patterson, The optimal addition of points to quadrature formulae, Mathematics of Computation, Volume 22, Number 104, October 1968, pages 847-856. Parameters: Input, int N, the order of the rule. Output, double X[N], the nodes. Output, double W[N], the weights. */ { double x01[1] = { 0.0000000000000000 }; double w01[1] = { 1.0000000000000000}; double x03[3] = { -1.73205080756887719, 0.000000000000000000, 1.73205080756887719 }; double w03[3] = { 0.166666666666666657, 0.66666666666666663, 0.166666666666666657}; double x07[7] = { -4.18495601767273229, -1.73205080756887719, -0.741095349994540853, 0.00000000000000000, 0.741095349994540853, 1.73205080756887719, 4.18495601767273229 }; double w07[7] = { 0.000695684158369139867, 0.138553274729749237, 0.13137860698313561, 0.458744868257491889, 0.13137860698313561, 0.138553274729749237, 0.000695684158369139867 }; double x09[9] = { -4.18495601767273229, -2.86127957605705818, -1.73205080756887719, -0.741095349994540853, 0.00000000000000000, 0.741095349994540853, 1.73205080756887719, 2.86127957605705818, 4.18495601767273229 }; double w09[9] = { 9.42694575565174701E-05, 0.00799632547089352934, 0.0948509485094851251, 0.270074329577937755, 0.253968253968254065, 0.270074329577937755, 0.0948509485094851251, 0.00799632547089352934, 9.42694575565174701E-05 }; double x17[17] = { -6.36339449433636961, -5.18701603991365623, -4.18495601767273229, -2.86127957605705818, -2.59608311504920231, -1.73205080756887719, -1.23042363402730603, -0.741095349994540853, 0.00000000000000000, 0.741095349994540853, 1.23042363402730603, 1.73205080756887719, 2.59608311504920231, 2.86127957605705818, 4.18495601767273229, 5.18701603991365623, 6.36339449433636961 }; double w17[17] = { 2.11364995054242569E-08, -8.20492075415092169E-07, 0.000105637836154169414, 0.00703348023782790748, 0.0019656770938777492, 0.0886810021520280101, 0.0141926548264493645, 0.254561232041712215, 0.266922230335053023, 0.254561232041712215, 0.0141926548264493645, 0.0886810021520280101, 0.0019656770938777492, 0.00703348023782790748, 0.000105637836154169414, -8.20492075415092169E-07, 2.11364995054242569E-08 }; double x19[19] = { -6.36339449433636961, -5.18701603991365623, -4.18495601767273229, -3.20533379449919442, -2.86127957605705818, -2.59608311504920231, -1.73205080756887719, -1.23042363402730603, -0.741095349994540853, 0.0000000000000000, 0.741095349994540853, 1.23042363402730603, 1.73205080756887719, 2.59608311504920231, 2.86127957605705818, 3.20533379449919442, 4.18495601767273229, 5.18701603991365623, 6.36339449433636961 }; double w19[19] = { 8.62968460222986318E-10, 6.09480873146898402E-07, 6.01233694598479965E-05, 0.00288488043650675591, -0.00633722479337375712, 0.0180852342547984622, 0.0640960546868076103, 0.0611517301252477163, 0.208324991649608771, 0.303467199854206227, 0.208324991649608771, 0.0611517301252477163, 0.0640960546868076103, 0.0180852342547984622, -0.00633722479337375712, 0.00288488043650675591, 6.01233694598479965E-05, 6.09480873146898402E-07, 8.62968460222986318E-10 }; double x31[31] = { -9.0169397898903032, -7.98077179859056063, -7.12210670080461661, -6.36339449433636961, -5.18701603991365623, -4.18495601767273229, -3.63531851903727832, -3.20533379449919442, -2.86127957605705818, -2.59608311504920231, -2.23362606167694189, -1.73205080756887719, -1.23042363402730603, -0.741095349994540853, -0.248992297579960609, 0.00000000000000000, 0.248992297579960609, 0.741095349994540853, 1.23042363402730603, 1.73205080756887719, 2.23362606167694189, 2.59608311504920231, 2.86127957605705818, 3.20533379449919442, 3.63531851903727832, 4.18495601767273229, 5.18701603991365623, 6.36339449433636961, 7.12210670080461661, 7.98077179859056063, 9.0169397898903032 }; double w31[31] = { 1.26184642808151181E-15, -1.4840835740298868E-13, 5.11580531055042083E-12, 7.92982678648693382E-10, 6.14358432326179133E-07, 5.94749611639316215E-05, 1.50442053909142189E-05, 0.00272984304673340016, -0.00556100630683581572, 0.0165924926989360101, 0.00176084755813180017, 0.0617185325658671791, 0.0654173928360925611, 0.199688635117345498, 0.0281281015400331666, 0.25890005324151566, 0.0281281015400331666, 0.199688635117345498, 0.0654173928360925611, 0.0617185325658671791, 0.00176084755813180017, 0.0165924926989360101, -0.00556100630683581572, 0.00272984304673340016, 1.50442053909142189E-05, 5.94749611639316215E-05, 6.14358432326179133E-07, 7.92982678648693382E-10, 5.11580531055042083E-12, -1.4840835740298868E-13, 1.26184642808151181E-15 }; double x33[33] = { -9.0169397898903032, -7.98077179859056063, -7.12210670080461661, -6.36339449433636961, -5.69817776848810986, -5.18701603991365623, -4.18495601767273229, -3.63531851903727832, -3.20533379449919442, -2.86127957605705818, -2.59608311504920231, -2.23362606167694189, -1.73205080756887719, -1.23042363402730603, -0.741095349994540853, -0.248992297579960609, 0.00000000000000000, 0.248992297579960609, 0.741095349994540853, 1.23042363402730603, 1.73205080756887719, 2.23362606167694189, 2.59608311504920231, 2.86127957605705818, 3.20533379449919442, 3.63531851903727832, 4.18495601767273229, 5.18701603991365623, 5.69817776848810986, 6.36339449433636961, 7.12210670080461661, 7.98077179859056063, 9.0169397898903032 }; double w33[33] = { -9.93139132868224651E-16, 2.66406251662316506E-13, -1.93413050008809555E-11, 1.5542195992782658E-09, -1.34860173485429301E-08, 6.90862611791137378E-07, 5.56911589810814793E-05, 8.32360452957667447E-05, 0.00212022595595963252, -0.00277121890077892431, 0.01152924706539879, 0.00735301102049550764, 0.0546775561434630422, 0.0774436027462994808, 0.176075987415714591, 0.103876871255742839, 0.139110222363380387, 0.103876871255742839, 0.176075987415714591, 0.0774436027462994808, 0.0546775561434630422, 0.00735301102049550764, 0.01152924706539879, -0.00277121890077892431, 0.00212022595595963252, 8.32360452957667447E-05, 5.56911589810814793E-05, 6.90862611791137378E-07, -1.34860173485429301E-08, 1.5542195992782658E-09, -1.93413050008809555E-11, 2.66406251662316506E-13, -9.93139132868224651E-16 }; double x35[35] = { -9.0169397898903032, -7.98077179859056063, -7.12210670080461661, -6.36339449433636961, -5.69817776848810986, -5.18701603991365623, -4.73643308595229673, -4.18495601767273229, -3.63531851903727832, -3.20533379449919442, -2.86127957605705818, -2.59608311504920231, -2.23362606167694189, -1.73205080756887719, -1.23042363402730603, -0.741095349994540853, -0.248992297579960609, 0.00000000000000000, 0.248992297579960609, 0.741095349994540853, 1.23042363402730603, 1.73205080756887719, 2.23362606167694189, 2.59608311504920231, 2.86127957605705818, 3.20533379449919442, 3.63531851903727832, 4.18495601767273229, 4.73643308595229673, 5.18701603991365623, 5.69817776848810986, 6.36339449433636961, 7.12210670080461661, 7.98077179859056063, 9.0169397898903032 }; double w35[35] = { 1.05413265823340136E-18, 5.45004126506381281E-15, 3.09722235760629949E-12, 4.60117603486559168E-10, 2.13941944795610622E-08, 2.46764213457981401E-07, 2.73422068011878881E-06, 3.57293481989753322E-05, 0.000275242141167851312, 0.000818953927502267349, 0.00231134524035220713, 0.00315544626918755127, 0.015673473751851151, 0.0452736854651503914, 0.0923647267169863534, 0.148070831155215854, 0.191760115888044341, 0.000514894508069213769, 0.191760115888044341, 0.148070831155215854, 0.0923647267169863534, 0.0452736854651503914, 0.015673473751851151, 0.00315544626918755127, 0.00231134524035220713, 0.000818953927502267349, 0.000275242141167851312, 3.57293481989753322E-05, 2.73422068011878881E-06, 2.46764213457981401E-07, 2.13941944795610622E-08, 4.60117603486559168E-10, 3.09722235760629949E-12, 5.45004126506381281E-15, 1.05413265823340136E-18 }; if ( n == 1 ) { r8vec_copy ( n, x01, x ); r8vec_copy ( n, w01, w ); } else if ( n == 3 ) { r8vec_copy ( n, x03, x ); r8vec_copy ( n, w03, w ); } else if ( n == 7 ) { r8vec_copy ( n, x07, x ); r8vec_copy ( n, w07, w ); } else if ( n == 9 ) { r8vec_copy ( n, x09, x ); r8vec_copy ( n, w09, w ); } else if ( n == 17 ) { r8vec_copy ( n, x17, x ); r8vec_copy ( n, w17, w ); } else if ( n == 19 ) { r8vec_copy ( n, x19, x ); r8vec_copy ( n, w19, w ); } else if ( n == 31 ) { r8vec_copy ( n, x31, x ); r8vec_copy ( n, w31, w ); } else if ( n == 33 ) { r8vec_copy ( n, x33, x ); r8vec_copy ( n, w33, w ); } else if ( n == 35 ) { r8vec_copy ( n, x35, x ); r8vec_copy ( n, w35, w ); } else { fprintf ( stderr, "\n" ); fprintf ( stderr, "KPN - Fatal error!\n" ); fprintf ( stderr, " Illegal value of N.\n" ); exit ( 1 ); } return; } /******************************************************************************/ int kpn_order ( int l ) /******************************************************************************/ /* Purpose: KPN_ORDER computes the order of a KPN rule from the level. Licensing: This code is distributed under the MIT license. Modified: 08 December 2012 Author: John Burkardt. Parameters: Input, int L, the level of the rule. 1 <= L <= 25 Output, int KPN_ORDER, the order of the rule. */ { int n; if ( l < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "KPN_ORDER - Fatal error!\n" ); fprintf ( stderr, " 1 <= L <= 25 required.\n" ); fprintf ( stderr, " Input L = %d\n", l ); exit ( 1 ); } else if ( l == 1 ) { n = 1; } else if ( l <= 3 ) { n = 3; } else if ( l == 4 ) { n = 7; } else if ( l <= 8 ) { n = 9; } else if ( l == 9 ) { n = 17; } else if ( l <= 15 ) { n = 19; } else if ( l == 16 ) { n = 31; } else if ( l == 17 ) { n = 33; } else if ( l <= 25 ) { n = 35; } else { fprintf ( stderr, "\n" ); fprintf ( stderr, "KPN_ORDER - Fatal error!\n" ); fprintf ( stderr, " 1 <= L <= 25 required.\n" ); fprintf ( stderr, " Input L = %d\n", l ); exit ( 1 ); } return n; } /******************************************************************************/ void kpu ( int n, double x[], double w[] ) /******************************************************************************/ /* Purpose: KPU provides data for Kronrod-Patterson quadrature with a uniform weight. Discussion: This data assumes integration over the interval [0,1] with weight function w(x) = 1. This data was originally supplied with only 7 digit accuracy. It has been replaced by higher accuracy data, which is defined over [-1,+1], but adjusted to the interval [0,1] before return. Licensing: This code is distributed under the MIT license. Modified: 11 December 2012 Author: John Burkardt. Reference: Florian Heiss, Viktor Winschel, Likelihood approximation by numerical integration on sparse grids, Journal of Econometrics, Volume 144, 2008, pages 62-80. Alan Genz, Bradley Keister, Fully symmetric interpolatory rules for multiple integrals over infinite regions with Gaussian weight, Journal of Computational and Applied Mathematics, Volume 71, 1996, pages 299-309. Thomas Patterson, The optimal addition of points to quadrature formulae, Mathematics of Computation, Volume 22, Number 104, October 1968, pages 847-856. Parameters: Input, int N, the order of the rule. Only 1, 3, 7, 15, 31 and 63 are legal input values for N. Output, double X[N], the nodes. Output, double W[N], the weights. */ { double x01[1] = { 0.0000000 }; double w01[1] = { 2.0000000 }; double x03[3] = { -0.77459666924148337704, 0.0, 0.77459666924148337704 }; double w03[3] = { 0.555555555555555555556, 0.888888888888888888889, 0.555555555555555555556 }; double x07[7] = { -0.96049126870802028342, -0.77459666924148337704, -0.43424374934680255800, 0.0, 0.43424374934680255800, 0.77459666924148337704, 0.96049126870802028342 }; double w07[7] = { 0.104656226026467265194, 0.268488089868333440729, 0.401397414775962222905, 0.450916538658474142345, 0.401397414775962222905, 0.268488089868333440729, 0.104656226026467265194 }; double x15[15] = { -0.99383196321275502221, -0.96049126870802028342, -0.88845923287225699889, -0.77459666924148337704, -0.62110294673722640294, -0.43424374934680255800, -0.22338668642896688163, 0.0, 0.22338668642896688163, 0.43424374934680255800, 0.62110294673722640294, 0.77459666924148337704, 0.88845923287225699889, 0.96049126870802028342, 0.99383196321275502221 }; double w15[15] = { 0.0170017196299402603390, 0.0516032829970797396969, 0.0929271953151245376859, 0.134415255243784220360, 0.171511909136391380787, 0.200628529376989021034, 0.219156858401587496404, 0.225510499798206687386, 0.219156858401587496404, 0.200628529376989021034, 0.171511909136391380787, 0.134415255243784220360, 0.0929271953151245376859, 0.0516032829970797396969, 0.0170017196299402603390 }; double x31[31] = { -0.99909812496766759766, -0.99383196321275502221, -0.98153114955374010687, -0.96049126870802028342, -0.92965485742974005667, -0.88845923287225699889, -0.83672593816886873550, -0.77459666924148337704, -0.70249620649152707861, -0.62110294673722640294, -0.53131974364437562397, -0.43424374934680255800, -0.33113539325797683309, -0.22338668642896688163, -0.11248894313318662575, 0.0, 0.11248894313318662575, 0.22338668642896688163, 0.33113539325797683309, 0.43424374934680255800, 0.53131974364437562397, 0.62110294673722640294, 0.70249620649152707861, 0.77459666924148337704, 0.83672593816886873550, 0.88845923287225699889, 0.92965485742974005667, 0.96049126870802028342, 0.98153114955374010687, 0.99383196321275502221, 0.99909812496766759766 }; double w31[31] = { 0.00254478079156187441540, 0.00843456573932110624631, 0.0164460498543878109338, 0.0258075980961766535646, 0.0359571033071293220968, 0.0464628932617579865414, 0.0569795094941233574122, 0.0672077542959907035404, 0.0768796204990035310427, 0.0857559200499903511542, 0.0936271099812644736167, 0.100314278611795578771, 0.105669893580234809744, 0.109578421055924638237, 0.111956873020953456880, 0.112755256720768691607, 0.111956873020953456880, 0.109578421055924638237, 0.105669893580234809744, 0.100314278611795578771, 0.0936271099812644736167, 0.0857559200499903511542, 0.0768796204990035310427, 0.0672077542959907035404, 0.0569795094941233574122, 0.0464628932617579865414, 0.0359571033071293220968, 0.0258075980961766535646, 0.0164460498543878109338, 0.00843456573932110624631, 0.00254478079156187441540 }; double x63[63] = { -0.99987288812035761194, -0.99909812496766759766, -0.99720625937222195908, -0.99383196321275502221, -0.98868475754742947994, -0.98153114955374010687, -0.97218287474858179658, -0.96049126870802028342, -0.94634285837340290515, -0.92965485742974005667, -0.91037115695700429250, -0.88845923287225699889, -0.86390793819369047715, -0.83672593816886873550, -0.80694053195021761186, -0.77459666924148337704, -0.73975604435269475868, -0.70249620649152707861, -0.66290966002478059546, -0.62110294673722640294, -0.57719571005204581484, -0.53131974364437562397, -0.48361802694584102756, -0.43424374934680255800, -0.38335932419873034692, -0.33113539325797683309, -0.27774982202182431507, -0.22338668642896688163, -0.16823525155220746498, -0.11248894313318662575, -0.056344313046592789972, 0.0, 0.056344313046592789972, 0.11248894313318662575, 0.16823525155220746498, 0.22338668642896688163, 0.27774982202182431507, 0.33113539325797683309, 0.38335932419873034692, 0.43424374934680255800, 0.48361802694584102756, 0.53131974364437562397, 0.57719571005204581484, 0.62110294673722640294, 0.66290966002478059546, 0.70249620649152707861, 0.73975604435269475868, 0.77459666924148337704, 0.80694053195021761186, 0.83672593816886873550, 0.86390793819369047715, 0.88845923287225699889, 0.91037115695700429250, 0.92965485742974005667, 0.94634285837340290515, 0.96049126870802028342, 0.97218287474858179658, 0.98153114955374010687, 0.98868475754742947994, 0.99383196321275502221, 0.99720625937222195908, 0.99909812496766759766, 0.99987288812035761194 }; double w63[63] = { 0.000363221481845530659694, 0.00126515655623006801137, 0.00257904979468568827243, 0.00421763044155885483908, 0.00611550682211724633968, 0.00822300795723592966926, 0.0104982469096213218983, 0.0129038001003512656260, 0.0154067504665594978021, 0.0179785515681282703329, 0.0205942339159127111492, 0.0232314466399102694433, 0.0258696793272147469108, 0.0284897547458335486125, 0.0310735511116879648799, 0.0336038771482077305417, 0.0360644327807825726401, 0.0384398102494555320386, 0.0407155101169443189339, 0.0428779600250077344929, 0.0449145316536321974143, 0.0468135549906280124026, 0.0485643304066731987159, 0.0501571393058995374137, 0.0515832539520484587768, 0.0528349467901165198621, 0.0539054993352660639269, 0.0547892105279628650322, 0.0554814043565593639878, 0.0559784365104763194076, 0.0562776998312543012726, 0.0563776283603847173877, 0.0562776998312543012726, 0.0559784365104763194076, 0.0554814043565593639878, 0.0547892105279628650322, 0.0539054993352660639269, 0.0528349467901165198621, 0.0515832539520484587768, 0.0501571393058995374137, 0.0485643304066731987159, 0.0468135549906280124026, 0.0449145316536321974143, 0.0428779600250077344929, 0.0407155101169443189339, 0.0384398102494555320386, 0.0360644327807825726401, 0.0336038771482077305417, 0.0310735511116879648799, 0.0284897547458335486125, 0.0258696793272147469108, 0.0232314466399102694433, 0.0205942339159127111492, 0.0179785515681282703329, 0.0154067504665594978021, 0.0129038001003512656260, 0.0104982469096213218983, 0.00822300795723592966926, 0.00611550682211724633968, 0.00421763044155885483908, 0.00257904979468568827243, 0.00126515655623006801137, 0.000363221481845530659694 }; if ( n == 1 ) { r8vec_copy ( n, x01, x ); r8vec_copy ( n, w01, w ); } else if ( n == 3 ) { r8vec_copy ( n, x03, x ); r8vec_copy ( n, w03, w ); } else if ( n == 7 ) { r8vec_copy ( n, x07, x ); r8vec_copy ( n, w07, w ); } else if ( n == 15 ) { r8vec_copy ( n, x15, x ); r8vec_copy ( n, w15, w ); } else if ( n == 31 ) { r8vec_copy ( n, x31, x ); r8vec_copy ( n, w31, w ); } else if ( n == 63 ) { r8vec_copy ( n, x63, x ); r8vec_copy ( n, w63, w ); } else { fprintf ( stderr, "\n" ); fprintf ( stderr, "KPU - Fatal error!\n" ); fprintf ( stderr, " Illegal value of N.\n" ); exit ( 1 ); } /* The rule as stored is for the interval [-1,+1]. Adjust it to the interval [0,1]. */ rule_adjust ( -1.0, +1.0, 0.0, +1.0, n, x, w ); return; } /******************************************************************************/ int kpu_order ( int l ) /******************************************************************************/ /* Purpose: KPU_ORDER computes the order of a KPU rule from the level. Licensing: This code is distributed under the MIT license. Modified: 08 December 2012 Author: John Burkardt. Parameters: Input, int L, the level of the rule. 1 <= L <= 25 Output, int KPU_ORDER, the order of the rule. */ { int n; if ( l < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "KPU_ORDER - Fatal error!\n" ); fprintf ( stderr, " 1 <= L <= 25 required.\n" ); fprintf ( stderr, " Input L = %d\n", l ); exit ( 1 ); } else if ( l == 1 ) { n = 1; } else if ( l <= 3 ) { n = 3; } else if ( l <= 6 ) { n = 7; } else if ( l <= 12 ) { n = 15; } else if ( l <= 24 ) { n = 31; } else if ( l <= 25 ) { n = 63; } else { fprintf ( stderr, "\n" ); fprintf ( stderr, "KPU_ORDER - Fatal error!\n" ); fprintf ( stderr, " 1 <= L <= 25 required.\n" ); fprintf ( stderr, " Input L = %d\n", l ); exit ( 1 ); } return n; } /******************************************************************************/ int num_seq ( int n, int k ) /******************************************************************************/ /* Purpose: NUM_SEQ returns the number of compositions of the integer N into K parts. Discussion: A composition of the integer N into K parts is an ordered sequence of K nonnegative integers which sum to N. The compositions (1,2,1) and (1,1,2) are considered to be distinct. The 28 compositions of 6 into three parts are: 6 0 0, 5 1 0, 5 0 1, 4 2 0, 4 1 1, 4 0 2, 3 3 0, 3 2 1, 3 1 2, 3 0 3, 2 4 0, 2 3 1, 2 2 2, 2 1 3, 2 0 4, 1 5 0, 1 4 1, 1 3 2, 1 2 3, 1 1 4, 1 0 5, 0 6 0, 0 5 1, 0 4 2, 0 3 3, 0 2 4, 0 1 5, 0 0 6. The formula for the number of compositions of N into K parts is Number = ( N + K - 1 )! / ( N! * ( K - 1 )! ) Describe the composition using N '1's and K-1 dividing lines '|'. The number of distinct permutations of these symbols is the number of compositions. This is equal to the number of permutations of N+K-1 things, with N identical of one kind and K-1 identical of another. Thus, for the above example, we have: Number = ( 6 + 3 - 1 )! / ( 6! * (3-1)! ) = 28 Licensing: This code is distributed under the MIT license. Modified: 08 December 2012 Author: John Burkardt Reference: Albert Nijenhuis, Herbert Wilf, Combinatorial Algorithms for Computers and Calculators, Second Edition, Academic Press, 1978, ISBN: 0-12-519260-6, LC: QA164.N54. Parameters: Input, int N, the integer whose compositions are desired. Input, int K, the number of parts in the composition. Output, int NUM_SEQ, the number of compositions of N into K parts. */ { int value; value = i4_choose ( n + k - 1, n ); return value; } /******************************************************************************/ void nwspgr ( void rule ( int n, double x[], double w[] ), int rule_order ( int l ), int dim, int k, int r_size, int *s_size, double nodes[], double weights[] ) /******************************************************************************/ /* Purpose: NWSPGR generates nodes and weights for sparse grid integration. Licensing: This code is distributed under the MIT license. Modified: 19 April 2013 Author: Original MATLAB version by Florian Heiss, Viktor Winschel. C version by John Burkardt. Reference: Florian Heiss, Viktor Winschel, Likelihood approximation by numerical integration on sparse grids, Journal of Econometrics, Volume 144, 2008, pages 62-80. Parameters: Input, void RULE ( int n, double x[], double w[] ), the name of a function which is given the order N and returns the points X and weights W of the corresponding 1D quadrature rule. Input, int RULE_ORDER ( int l ), the name of a function which is given the level L and returns the order N of the corresponding 1D rule. Input, int DIM, the spatial dimension. Input, int K, the level of the sparse rule. Input, int R_SIZE, the "size" of the sparse rule. Output, int *S_SIZE, the size of the sparse rule, after duplicate points have been merged. Output, double NODES[DIM*R_SIZE], the nodes of the sparse rule. Output, double WEIGHTS[R_SIZE], the weights of the sparse rule. */ { int bq; int equal; int i; int *is; int j; int j2; int level; int *lr; int maxq; int minq; int n; int *n1d; int n1d_total; int nc; int np; int *nr; int q; int r; int *roff; int *rq; int seq_num; double t; double *w; double *w1d; int *w1d_off; double *wc; double *wp; double *wr; double *x; double *x1d; int *x1d_off; double *xc; double *xp; double *xr; for ( j = 0; j < r_size; j++ ) { for ( i = 0; i < dim; i++ ) { nodes[i+j*dim] = 0.0; } } for ( j = 0; j < r_size; j++ ) { weights[j] = 0.0; } /* Create cell arrays that will contain the points and weights for levels 1 through K. */ n1d = ( int * ) malloc ( k * sizeof ( int ) ); x1d_off = ( int * ) malloc ( ( k + 1 ) * sizeof ( int ) ); w1d_off = ( int * ) malloc ( ( k + 1 ) * sizeof ( int ) ); x1d_off[0] = 0; w1d_off[0] = 0; for ( level = 1; level <= k; level++ ) { n = rule_order ( level ); n1d[level-1] = n; x1d_off[level] = x1d_off[level-1] + n; w1d_off[level] = w1d_off[level-1] + n; } n1d_total = x1d_off[k]; /* Calculate all the 1D rules needed. */ x1d = ( double * ) malloc ( n1d_total * sizeof ( double ) ); w1d = ( double * ) malloc ( n1d_total * sizeof ( double ) ); for ( level = 1; level <= k; level++ ) { n = n1d[level-1]; x = ( double * ) malloc ( n * sizeof ( double ) ); w = ( double * ) malloc ( n * sizeof ( double ) ); rule ( n, x, w ); r8cvv_rset ( n1d_total, x1d, k, x1d_off, level - 1, x ); r8cvv_rset ( n1d_total, w1d, k, w1d_off, level - 1, w ); free ( x ); free ( w ); } /* Construct the sparse grid. */ minq = i4_max ( 0, k - dim ); maxq = k - 1; /* Q is the level total. */ lr = ( int * ) malloc ( dim * sizeof ( int ) ); nr = ( int * ) malloc ( dim * sizeof ( int ) ); r = 0; for ( q = minq; q <= maxq; q++ ) { /* BQ is the combinatorial coefficient applied to the component product rules which have level Q. */ bq = i4_mop ( maxq - q ) * i4_choose ( dim - 1, dim + q - k ); /* Compute the D-dimensional row vectors that sum to DIM+Q. */ seq_num = num_seq ( q, dim ); is = get_seq ( dim, q + dim, seq_num ); /* Allocate new rows for nodes and weights. */ rq = ( int * ) malloc ( seq_num * sizeof ( int ) ); for ( j = 0; j < seq_num; j++ ) { rq[j] = 1; for ( i = 0; i < dim; i++ ) { level = is[j+i*seq_num] - 1; rq[j] = rq[j] * n1d[level]; } } /* Generate every product rule whose level total is Q. */ for ( j2 = 0; j2 < seq_num; j2++ ) { for ( i = 0; i < dim; i++ ) { lr[i] = is[j2+i*seq_num]; } for ( i = 0; i < dim; i++ ) { nr[i] = rule_order ( lr[i] ); } roff = r8cvv_offset ( dim, nr ); nc = i4vec_sum ( dim, nr ); wc = ( double * ) malloc ( nc * sizeof ( double ) ); xc = ( double * ) malloc ( nc * sizeof ( double ) ); /* Corrected first dimension to N1D_TOTAL in following two calls to R8CVV_GET_NEW, 19 April 2013. */ for ( i = 0; i < dim; i++ ) { xr = r8cvv_rget_new ( n1d_total, x1d, k, x1d_off, lr[i] - 1 ); wr = r8cvv_rget_new ( n1d_total, w1d, k, w1d_off, lr[i] - 1 ); r8cvv_rset ( nc, xc, dim, roff, i, xr ); r8cvv_rset ( nc, wc, dim, roff, i, wr ); free ( wr ); free ( xr ); } np = rq[j2]; wp = ( double * ) malloc ( np * sizeof ( double ) ); xp = ( double * ) malloc ( dim * np * sizeof ( double ) ); tensor_product_cell ( nc, xc, wc, dim, nr, roff, np, xp, wp ); /* Append the new nodes and weights to the arrays. */ for ( j = 0; j < np; j++ ) { for ( i = 0; i < dim; i++ ) { nodes[i+(r+j)*dim] = xp[i+j*dim]; } } for ( j = 0; j < np; j++ ) { weights[r+j] = bq * wp[j]; } /* Update the count. */ r = r + rq[j2]; free ( roff ); free ( wc ); free ( wp ); free ( xc ); free ( xp ); } free ( is ); free ( rq ); } free ( lr ); free ( n1d ); free ( nr ); free ( w1d ); free ( w1d_off ); free ( x1d ); free ( x1d_off ); /* Reorder the rule so the points are in ascending lexicographic order. */ rule_sort ( dim, r_size, nodes, weights ); /* Suppress duplicate points and merge weights. */ r = 0; for ( j = 1; j < r_size; j++ ) { equal = 1; for ( i = 0; i < dim; i++ ) { if ( nodes[i+r*dim] != nodes[i+j*dim] ) { equal = 0; break; } } if ( equal ) { weights[r] = weights[r] + weights[j]; } else { r = r + 1; weights[r] = weights[j]; for ( i = 0; i < dim; i++ ) { nodes[i+r*dim] = nodes[i+j*dim]; } } } r = r + 1; *s_size = r; /* Zero out unneeded entries. */ for ( j = r; j < r_size; j++ ) { for ( i = 0; i < dim; i++ ) { nodes[i+j*dim] = 0.0; } } for ( j = r; j < r_size; j++ ) { weights[j] = 0.0; } /* Normalize the weights to sum to 1. */ t = r8vec_sum ( r, weights ); for ( j = 0; j < r; j++ ) { weights[j] = weights[j] / t; } return; } /******************************************************************************/ int nwspgr_size ( int rule_order ( int l ), int dim, int k ) /******************************************************************************/ /* Purpose: NWSPGR_SIZE determines the size of a sparse grid rule. Discussion: This routine does a "raw" count, that is, it does not notice that many points may be repeated, in which case, the size of the rule could be reduced by merging repeated points and combining the corresponding weights. Licensing: This code is distributed under the MIT license. Modified: 10 December 2012 Author: John Burkardt. Reference: Florian Heiss, Viktor Winschel, Likelihood approximation by numerical integration on sparse grids, Journal of Econometrics, Volume 144, 2008, pages 62-80. Parameters: Input, int RULE_ORDER ( int l ), the name of a function which is given the level L and returns the order N of the corresponding rule. Input, int DIM, the dimension of the integration problem. Input, int K, the level. When using the built in 1D rules, the resulting sparse grid will be exact for polynomials up to total order 2*K-1. When using the built-in quadrature rules, the maximum value of K that is available is 25. Output, int NWSPGR_SIZE, the "size" of the rule, that is, the number of weights and multidimensional quadrature points that will be needed. The size of the rule will be reduced when duplicate points are merged. */ { int i; int *is; int j; int level; int maxq; int minq; int n; int *n1d; int q; int r_size; int *rq; int seq_num; /* Determine the size of each 1D rule. */ n1d = ( int * ) malloc ( k * sizeof ( int ) ); for ( level = 1; level <= k; level++ ) { n = rule_order ( level ); n1d[level-1] = n; } /* Go through the motions of generating the rules. */ minq = i4_max ( 0, k - dim ); maxq = k - 1; r_size = 0; for ( q = minq; q <= maxq; q++ ) { /* Compute the D-dimensional vectors that sum to Q+DIM. */ seq_num = num_seq ( q, dim ); is = get_seq ( dim, q + dim, seq_num ); /* Determine the size of each rule. */ rq = ( int * ) malloc ( seq_num * sizeof ( int ) ); for ( j = 0; j < seq_num; j++ ) { rq[j] = 1; for ( i = 0; i < dim; i++ ) { rq[j] = rq[j] * n1d[is[j+i*seq_num]-1]; } } /* Add the sizes to the total. */ r_size = r_size + i4vec_sum ( seq_num, rq ); free ( is ); free ( rq ); } free ( n1d ); return r_size; } /******************************************************************************/ void quad_rule_print ( int m, int n, double x[], double w[], char *title ) /******************************************************************************/ /* Purpose: QUAD_RULE_PRINT prints a multidimensional quadrature rule. Licensing: This code is distributed under the MIT license. Modified: 08 December 2012 Author: John Burkardt Parameters: Input, int M, the spatial dimension. Input, int N, the number of points. Input, double X[M*N], the abscissas. Input, double W[N], the weights. Input, char *TITLE, a title for the rule. */ { int i; int j; printf ( "\n" ); printf ( "%s\n", title ); printf ( "\n" ); for ( j = 0; j < n; j++ ) { printf ( " %2d %10.6f * f (", j, w[j] ); for ( i = 0; i < m; i++ ) { printf ( "%10.6f", x[i+j*m] ); if ( i < m - 1 ) { printf ( "," ); } else { printf ( " )\n" ); } } } return; } /******************************************************************************/ double r8_uniform_01 ( int *seed ) /******************************************************************************/ /* Purpose: R8_UNIFORM_01 returns a pseudorandom R8 scaled to [0,1]. Discussion: This routine implements the recursion seed = 16807 * seed mod ( 2^31 - 1 ) r8_uniform_01 = seed / ( 2^31 - 1 ) The integer arithmetic never requires more than 32 bits, including a sign bit. If the initial seed is 12345, then the first three computations are Input Output R8_UNIFORM_01 SEED SEED 12345 207482415 0.096616 207482415 1790989824 0.833995 1790989824 2035175616 0.947702 Licensing: This code is distributed under the MIT license. Modified: 11 August 2004 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Springer Verlag, pages 201-202, 1983. Pierre L'Ecuyer, Random Number Generation, in Handbook of Simulation edited by Jerry Banks, Wiley Interscience, page 95, 1998. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, pages 362-376, 1986. P A Lewis, A S Goodman, J M Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, pages 136-143, 1969. Parameters: Input/output, int *SEED, the "seed" value. Normally, this value should not be 0. On output, SEED has been updated. Output, double R8_UNIFORM_01, a new pseudorandom variate, strictly between 0 and 1. */ { int k; double r; k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + 2147483647; } r = ( ( double ) ( *seed ) ) * 4.656612875E-10; return r; } /******************************************************************************/ int *r8cvv_offset ( int m, int nr[] ) /******************************************************************************/ /* Purpose: R8CVV_OFFSET determines the row offsets of an R8CVV. Discussion: An R8CVV is a "vector of vectors" of R8's. Licensing: This code is distributed under the MIT license. Modified: 16 November 2012 Author: John Burkardt Parameters: Input, int M, the number of rows in the array. Input, int NR(M), the row sizes. Output, int R8CVV_OFFSET[M+1], the row offsets. */ { int i; int *roff; roff = ( int * ) malloc ( ( m + 1 ) * sizeof ( int ) ); roff[0] = 0; for ( i = 0; i < m; i++ ) { roff[i+1] = roff[i] + nr[i]; } return roff; } /******************************************************************************/ void r8cvv_print ( int mn, double a[], int m, int roff[], char *title ) /******************************************************************************/ /* Purpose: R8CVV_PRINT prints an R8CVV. Discussion: An R8CVV is a "vector of vectors" of R8's. Licensing: This code is distributed under the MIT license. Modified: 16 November 2012 Author: John Burkardt Parameters: Input, int MN, the size of the cell array. Input, double A[MN], the cell array. Input, int M, the number of rows in the array. Input, int ROFF[M+1], the row offsets. Input, char *TITLE, a title. */ { int i; int k; int k1; int k2; int khi; int klo; printf ( "\n" ); printf ( "%s\n", title ); printf ( "\n" ); for ( i = 0; i < m; i++ ) { k1 = roff[i]; k2 = roff[i+1]; for ( klo = k1; klo < k2; klo = klo + 5 ) { khi = i4_min ( klo + 5, k2 ); if ( klo == k1 ) { printf ( "%5d", i ); } else { printf ( " " ); } printf ( " " ); for ( k = klo; k < khi; k++ ) { printf ( "%14g", a[k] ); } printf ( "\n" ); } } return; } /******************************************************************************/ double *r8cvv_rget_new ( int mn, double a[], int m, int roff[], int i ) /******************************************************************************/ /* Purpose: R8CVV_RGET_NEW gets row I from an R8CVV. Discussion: An R8CVV is a "vector of vectors" of R8's. Licensing: This code is distributed under the MIT license. Modified: 16 November 2012 Author: John Burkardt Parameters: Input, int MN, the size of the cell array. Input, double A[MN], the cell array. Input, int M, the number of rows in the array. Input, int ROFF[M+1], the row offsets. Input, int I, the row. 0 <= I < M. Output, double R8CVV_RGET_NEW[NR[I]], the value of A(I,*). */ { double *ai; int j; int k1; int k2; int nv; k1 = roff[i]; k2 = roff[i+1]; nv = k2 - k1; ai = ( double * ) malloc ( nv * sizeof ( double ) ); for ( j = 0; j < nv; j++ ) { ai[j] = a[k1+j]; } return ai; } /******************************************************************************/ void r8cvv_rset ( int mn, double a[], int m, int roff[], int i, double ai[] ) /******************************************************************************/ /* Purpose: R8CVV_RSET sets row I from an R8CVV. Discussion: An R8CVV is a "vector of vectors" of R8's. Licensing: This code is distributed under the MIT license. Modified: 16 November 2012 Author: John Burkardt Parameters: Input, int MN, the size of the cell array. Input/output, double A[MN], the cell array. Input, int M, the number of rows in the array. Input, int ROFF[M+1], the row offsets. Input, int I, the row. 0 <= I < M. Input, double AI[NR[I]], the new value of A(I,*). */ { int j; int k1; int k2; int nv; k1 = roff[i]; k2 = roff[i+1]; nv = k2 - k1; for ( j = 0; j < nv; j++ ) { a[k1+j] = ai[j]; } return; } /******************************************************************************/ double *r8mat_normal_01_new ( int m, int n, int *seed ) /******************************************************************************/ /* Purpose: R8MAT_NORMAL_01_NEW returns a unit pseudonormal R8MAT. Licensing: This code is distributed under the MIT license. Modified: 03 October 2005 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Springer Verlag, pages 201-202, 1983. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, pages 362-376, 1986. Peter Lewis, Allen Goodman, James Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, pages 136-143, 1969. Parameters: Input, int M, N, the number of rows and columns in the array. Input/output, int *SEED, the "seed" value, which should NOT be 0. On output, SEED has been updated. Output, double R8MAT_NORMAL_01_NEW[M*N], the array of pseudonormal values. */ { double *r; r = r8vec_normal_01_new ( m * n, seed ); return r; } /******************************************************************************/ void r8mat_transpose_print ( int m, int n, double a[], char *title ) /******************************************************************************/ /* Purpose: R8MAT_TRANSPOSE_PRINT prints an R8MAT, transposed. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. Licensing: This code is distributed under the MIT license. Modified: 28 May 2008 Author: John Burkardt Parameters: Input, int M, N, the number of rows and columns. Input, double A[M*N], an M by N matrix to be printed. Input, char *TITLE, a title. */ { r8mat_transpose_print_some ( m, n, a, 1, 1, m, n, title ); return; } /******************************************************************************/ void r8mat_transpose_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi, int jhi, char *title ) /******************************************************************************/ /* Purpose: R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT, transposed. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. Licensing: This code is distributed under the MIT license. Modified: 20 August 2010 Author: John Burkardt Parameters: Input, int M, N, the number of rows and columns. Input, double A[M*N], an M by N matrix to be printed. Input, int ILO, JLO, the first row and column to print. Input, int IHI, JHI, the last row and column to print. Input, char *TITLE, a title. */ { # define INCX 5 int i; int i2; int i2hi; int i2lo; int inc; int j; int j2hi; int j2lo; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); if ( m <= 0 || n <= 0 ) { fprintf ( stdout, "\n" ); fprintf ( stdout, " (None)\n" ); return; } for ( i2lo = i4_max ( ilo, 1 ); i2lo <= i4_min ( ihi, m ); i2lo = i2lo + INCX ) { i2hi = i2lo + INCX - 1; i2hi = i4_min ( i2hi, m ); i2hi = i4_min ( i2hi, ihi ); inc = i2hi + 1 - i2lo; fprintf ( stdout, "\n" ); fprintf ( stdout, " Row:" ); for ( i = i2lo; i <= i2hi; i++ ) { fprintf ( stdout, " %7d ", i - 1 ); } fprintf ( stdout, "\n" ); fprintf ( stdout, " Col\n" ); fprintf ( stdout, "\n" ); j2lo = i4_max ( jlo, 1 ); j2hi = i4_min ( jhi, n ); for ( j = j2lo; j <= j2hi; j++ ) { fprintf ( stdout, "%5d:", j - 1 ); for ( i2 = 1; i2 <= inc; i2++ ) { i = i2lo - 1 + i2; fprintf ( stdout, " %14f", a[(i-1)+(j-1)*m] ); } fprintf ( stdout, "\n" ); } } return; # undef INCX } /******************************************************************************/ double *r8mat_uniform_01_new ( int m, int n, int *seed ) /******************************************************************************/ /* Purpose: R8MAT_UNIFORM_01_NEW fills an R8MAT with pseudorandom values scaled to [0,1]. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. This routine implements the recursion seed = 16807 * seed mod ( 2^31 - 1 ) unif = seed / ( 2^31 - 1 ) The integer arithmetic never requires more than 32 bits, including a sign bit. Licensing: This code is distributed under the MIT license. Modified: 30 June 2009 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Springer Verlag, pages 201-202, 1983. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, pages 362-376, 1986. Philip Lewis, Allen Goodman, James Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, pages 136-143, 1969. Parameters: Input, int M, N, the number of rows and columns. Input/output, int *SEED, the "seed" value. Normally, this value should not be 0, otherwise the output value of SEED will still be 0, and R8_UNIFORM will be 0. On output, SEED has been updated. Output, double R8MAT_UNIFORM_01_NEW[M*N], a matrix of pseudorandom values. */ { int i; int j; int k; double *r; r = ( double * ) malloc ( m * n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + 2147483647; } r[i+j*m] = ( double ) ( *seed ) * 4.656612875E-10; } } return r; } /******************************************************************************/ void r8vec_copy ( int n, double a1[], double a2[] ) /******************************************************************************/ /* Purpose: R8VEC_COPY copies an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 03 July 2005 Author: John Burkardt Parameters: Input, int N, the number of entries in the vectors. Input, double A1[N], the vector to be copied. Input, double A2[N], the copy of A1. */ { int i; for ( i = 0; i < n; i++ ) { a2[i] = a1[i]; } return; } /******************************************************************************/ void r8vec_direct_product ( int factor_index, int factor_order, double factor_value[], int factor_num, int point_num, double x[] ) /******************************************************************************/ /* Purpose: R8VEC_DIRECT_PRODUCT creates a direct product of R8VEC's. Discussion: An R8VEC is a vector of R8's. To explain what is going on here, suppose we had to construct a multidimensional quadrature rule as the product of K rules for 1D quadrature. The product rule will be represented as a list of points and weights. The J-th item in the product rule will be associated with item J1 of 1D rule 1, item J2 of 1D rule 2, ..., item JK of 1D rule K. In particular, X(J) = ( X(1,J1), X(2,J2), ..., X(K,JK)) and W(J) = W(1,J1) * W(2,J2) * ... * W(K,JK) So we can construct the quadrature rule if we can properly distribute the information in the 1D quadrature rules. This routine carries out that task. Another way to do this would be to compute, one by one, the set of all possible indices (J1,J2,...,JK), and then index the appropriate information. An advantage of the method shown here is that you can process the K-th set of information and then discard it. Example: Rule 1: Order = 4 X(1:4) = ( 1, 2, 3, 4 ) Rule 2: Order = 3 X(1:3) = ( 10, 20, 30 ) Rule 3: Order = 2 X(1:2) = ( 100, 200 ) Product Rule: Order = 24 X(1:24) = ( 1, 10, 100 ) ( 2, 10, 100 ) ( 3, 10, 100 ) ( 4, 10, 100 ) ( 1, 20, 100 ) ( 2, 20, 100 ) ( 3, 20, 100 ) ( 4, 20, 100 ) ( 1, 30, 100 ) ( 2, 30, 100 ) ( 3, 30, 100 ) ( 4, 30, 100 ) ( 1, 10, 200 ) ( 2, 10, 200 ) ( 3, 10, 200 ) ( 4, 10, 200 ) ( 1, 20, 200 ) ( 2, 20, 200 ) ( 3, 20, 200 ) ( 4, 20, 200 ) ( 1, 30, 200 ) ( 2, 30, 200 ) ( 3, 30, 200 ) ( 4, 30, 200 ) Licensing: This code is distributed under the MIT license. Modified: 31 May 2009 Author: John Burkardt Parameters: Input, int FACTOR_INDEX, the index of the factor being processed. The first factor processed must be factor 0. Input, int FACTOR_ORDER, the order of the factor. Input, double FACTOR_VALUE[FACTOR_ORDER], the factor values for factor FACTOR_INDEX. Input, int FACTOR_NUM, the number of factors. Input, int POINT_NUM, the number of elements in the direct product. Input/output, double X[FACTOR_NUM*POINT_NUM], the elements of the direct product, which are built up gradually. Local Parameters: Local, int START, the first location of a block of values to set. Local, int CONTIG, the number of consecutive values to set. Local, int SKIP, the distance from the current value of START to the next location of a block of values to set. Local, int REP, the number of blocks of values to set. */ { static int contig = 0; int i; int j; int k; static int rep = 0; static int skip = 0; int start; if ( factor_index == 0 ) { contig = 1; skip = 1; rep = point_num; for ( j = 0; j < point_num; j++ ) { for ( i = 0; i < factor_num; i++ ) { x[i+j*factor_num] = 0.0; } } } rep = rep / factor_order; skip = skip * factor_order; for ( i = 0; i < factor_order; i++ ) { start = 0 + i * contig; for ( k = 1; k <= rep; k++ ) { for ( j = start; j < start + contig; j++ ) { x[factor_index+j*factor_num] = factor_value[i]; } start = start + skip; } } contig = contig * factor_order; return; } /******************************************************************************/ void r8vec_direct_product2 ( int factor_index, int factor_order, double factor_value[], int factor_num, int point_num, double w[] ) /******************************************************************************/ /* Purpose: R8VEC_DIRECT_PRODUCT2 creates a direct product of R8VEC's. Discussion: An R8VEC is a vector of R8's. To explain what is going on here, suppose we had to construct a multidimensional quadrature rule as the product of K rules for 1D quadrature. The product rule will be represented as a list of points and weights. The J-th item in the product rule will be associated with item J1 of 1D rule 1, item J2 of 1D rule 2, ..., item JK of 1D rule K. In particular, X(J) = ( X(1,J1), X(2,J2), ..., X(K,JK)) and W(J) = W(1,J1) * W(2,J2) * ... * W(K,JK) So we can construct the quadrature rule if we can properly distribute the information in the 1D quadrature rules. This routine carries out that task for the weights W. Another way to do this would be to compute, one by one, the set of all possible indices (J1,J2,...,JK), and then index the appropriate information. An advantage of the method shown here is that you can process the K-th set of information and then discard it. Example: Rule 1: Order = 4 W(1:4) = ( 2, 3, 5, 7 ) Rule 2: Order = 3 W(1:3) = ( 11, 13, 17 ) Rule 3: Order = 2 W(1:2) = ( 19, 23 ) Product Rule: Order = 24 W(1:24) = ( 2 * 11 * 19 ) ( 3 * 11 * 19 ) ( 4 * 11 * 19 ) ( 7 * 11 * 19 ) ( 2 * 13 * 19 ) ( 3 * 13 * 19 ) ( 5 * 13 * 19 ) ( 7 * 13 * 19 ) ( 2 * 17 * 19 ) ( 3 * 17 * 19 ) ( 5 * 17 * 19 ) ( 7 * 17 * 19 ) ( 2 * 11 * 23 ) ( 3 * 11 * 23 ) ( 5 * 11 * 23 ) ( 7 * 11 * 23 ) ( 2 * 13 * 23 ) ( 3 * 13 * 23 ) ( 5 * 13 * 23 ) ( 7 * 13 * 23 ) ( 2 * 17 * 23 ) ( 3 * 17 * 23 ) ( 5 * 17 * 23 ) ( 7 * 17 * 23 ) Licensing: This code is distributed under the MIT license. Modified: 31 May 2009 Author: John Burkardt Parameters: Input, int FACTOR_INDEX, the index of the factor being processed. The first factor processed must be factor 0. Input, int FACTOR_ORDER, the order of the factor. Input, double FACTOR_VALUE[FACTOR_ORDER], the factor values for factor FACTOR_INDEX. Input, int FACTOR_NUM, the number of factors. Input, int POINT_NUM, the number of elements in the direct product. Input/output, double W[POINT_NUM], the elements of the direct product, which are built up gradually. Local Parameters: Local, integer START, the first location of a block of values to set. Local, integer CONTIG, the number of consecutive values to set. Local, integer SKIP, the distance from the current value of START to the next location of a block of values to set. Local, integer REP, the number of blocks of values to set. */ { static int contig = 0; int i; int j; int k; static int rep = 0; static int skip = 0; int start; if ( factor_index == 0 ) { contig = 1; skip = 1; rep = point_num; for ( i = 0; i < point_num; i++ ) { w[i] = 1.0; } } rep = rep / factor_order; skip = skip * factor_order; for ( j = 0; j < factor_order; j++ ) { start = 0 + j * contig; for ( k = 1; k <= rep; k++ ) { for ( i = start; i < start + contig; i++ ) { w[i] = w[i] * factor_value[j]; } start = start + skip; } } contig = contig * factor_order; return; } /******************************************************************************/ double r8vec_dot_product ( int n, double a1[], double a2[] ) /******************************************************************************/ /* Purpose: R8VEC_DOT_PRODUCT computes the dot product of a pair of R8VEC's. Licensing: This code is distributed under the MIT license. Modified: 26 July 2007 Author: John Burkardt Parameters: Input, int N, the number of entries in the vectors. Input, double A1[N], A2[N], the two vectors to be considered. Output, double R8VEC_DOT_PRODUCT, the dot product of the vectors. */ { int i; double value; value = 0.0; for ( i = 0; i < n; i++ ) { value = value + a1[i] * a2[i]; } return value; } /******************************************************************************/ double *r8vec_normal_01_new ( int n, int *seed ) /******************************************************************************/ /* Purpose: R8VEC_NORMAL_01_NEW returns a unit pseudonormal R8VEC. Discussion: The standard normal probability distribution function (PDF) has mean 0 and standard deviation 1. This routine can generate a vector of values on one call. It has the feature that it should provide the same results in the same order no matter how we break up the task. Before calling this routine, the user may call RANDOM_SEED in order to set the seed of the random number generator. The Box-Muller method is used, which is efficient, but generates an even number of values each time. On any call to this routine, an even number of new values are generated. Depending on the situation, one value may be left over. In that case, it is saved for the next call. Licensing: This code is distributed under the MIT license. Modified: 18 February 2012 Author: John Burkardt Parameters: Input, int N, the number of values desired. If N is negative, then the code will flush its internal memory; in particular, if there is a saved value to be used on the next call, it is instead discarded. This is useful if the user has reset the random number seed, for instance. Input/output, int *SEED, a seed for the random number generator. Output, double R8VEC_NORMAL_01_NEW[N], a sample of the standard normal PDF. Local parameters: Local, int MADE, records the number of values that have been computed. On input with negative N, this value overwrites the return value of N, so the user can get an accounting of how much work has been done. Local, double R[N+1], is used to store some uniform random values. Its dimension is N+1, but really it is only needed to be the smallest even number greater than or equal to N. Local, int SAVED, is 0 or 1 depending on whether there is a single saved value left over from the previous call. Local, int X_LO, X_HI, records the range of entries of X that we need to compute. This starts off as 1:N, but is adjusted if we have a saved value that can be immediately stored in X(1), and so on. Local, double Y, the value saved from the previous call, if SAVED is 1. */ { # define R8_PI 3.141592653589793 int i; int m; static int made = 0; double *r; static int saved = 0; double *x; int x_hi; int x_lo; static double y = 0.0; x = ( double * ) malloc ( n * sizeof ( double ) ); /* I'd like to allow the user to reset the internal data. But this won't work properly if we have a saved value Y. I'm making a crock option that allows the user to signal explicitly that any internal memory should be flushed, by passing in a negative value for N. */ if ( n < 0 ) { made = 0; saved = 0; y = 0.0; return NULL; } else if ( n == 0 ) { return NULL; } /* Record the range of X we need to fill in. */ x_lo = 1; x_hi = n; /* Use up the old value, if we have it. */ if ( saved == 1 ) { x[0] = y; saved = 0; x_lo = 2; } /* Maybe we don't need any more values. */ if ( x_hi - x_lo + 1 == 0 ) { } /* If we need just one new value, do that here to avoid null arrays. */ else if ( x_hi - x_lo + 1 == 1 ) { r = r8vec_uniform_01_new ( 2, seed ); x[x_hi-1] = sqrt ( - 2.0 * log ( r[0] ) ) * cos ( 2.0 * R8_PI * r[1] ); y = sqrt ( - 2.0 * log ( r[0] ) ) * sin ( 2.0 * R8_PI * r[1] ); saved = 1; made = made + 2; free ( r ); } /* If we require an even number of values, that's easy. */ else if ( ( x_hi - x_lo + 1 ) % 2 == 0 ) { m = ( x_hi - x_lo + 1 ) / 2; r = r8vec_uniform_01_new ( 2*m, seed ); for ( i = 0; i <= 2*m-2; i = i + 2 ) { x[x_lo+i-1] = sqrt ( - 2.0 * log ( r[i] ) ) * cos ( 2.0 * R8_PI * r[i+1] ); x[x_lo+i ] = sqrt ( - 2.0 * log ( r[i] ) ) * sin ( 2.0 * R8_PI * r[i+1] ); } made = made + x_hi - x_lo + 1; free ( r ); } /* If we require an odd number of values, we generate an even number, and handle the last pair specially, storing one in X(N), and saving the other for later. */ else { x_hi = x_hi - 1; m = ( x_hi - x_lo + 1 ) / 2 + 1; r = r8vec_uniform_01_new ( 2*m, seed ); for ( i = 0; i <= 2*m-4; i = i + 2 ) { x[x_lo+i-1] = sqrt ( - 2.0 * log ( r[i] ) ) * cos ( 2.0 * R8_PI * r[i+1] ); x[x_lo+i ] = sqrt ( - 2.0 * log ( r[i] ) ) * sin ( 2.0 * R8_PI * r[i+1] ); } i = 2*m - 2; x[x_lo+i-1] = sqrt ( - 2.0 * log ( r[i] ) ) * cos ( 2.0 * R8_PI * r[i+1] ); y = sqrt ( - 2.0 * log ( r[i] ) ) * sin ( 2.0 * R8_PI * r[i+1] ); saved = 1; made = made + x_hi - x_lo + 2; free ( r ); } return x; # undef R8_PI } /******************************************************************************/ void r8vec_print ( int n, double a[], char *title ) /******************************************************************************/ /* Purpose: R8VEC_PRINT prints an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 08 April 2009 Author: John Burkardt Parameters: Input, int N, the number of components of the vector. Input, double A[N], the vector to be printed. Input, char *TITLE, a title. */ { int i; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); fprintf ( stdout, "\n" ); for ( i = 0; i < n; i++ ) { fprintf ( stdout, " %8d: %14f\n", i, a[i] ); } return; } /******************************************************************************/ double r8vec_sum ( int n, double a[] ) /******************************************************************************/ /* Purpose: R8VEC_SUM returns the sum of an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 26 August 2008 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Input, double A[N], the vector. Output, double R8VEC_SUM, the sum of the vector. */ { int i; double value; value = 0.0; for ( i = 0; i < n; i++ ) { value = value + a[i]; } return value; } /******************************************************************************/ void r8vec_transpose_print ( int n, double a[], char *title ) /******************************************************************************/ /* Purpose: R8VEC_TRANSPOSE_PRINT prints an R8VEC "transposed". Discussion: An R8VEC is a vector of R8's. Example: A = (/ 1.0, 2.1, 3.2, 4.3, 5.4, 6.5, 7.6, 8.7, 9.8, 10.9, 11.0 /) TITLE = 'My vector: ' My vector: 1.0 2.1 3.2 4.3 5.4 6.5 7.6 8.7 9.8 10.9 11.0 Licensing: This code is distributed under the MIT license. Modified: 12 November 2010 Author: John Burkardt Parameters: Input, int N, the number of components of the vector. Input, double A[N], the vector to be printed. Input, char *TITLE, a title. */ { int i; int ihi; int ilo; printf ( "\n" ); printf ( "%s\n", title ); printf ( "\n" ); if ( n <= 0 ) { printf ( " (Empty)\n" ); return; } for ( ilo = 0; ilo < n; ilo = ilo + 5 ) { ihi = i4_min ( ilo + 5, n ); for ( i = ilo; i < ihi; i++ ) { printf ( " %12f", a[i] ); } printf ( "\n" ); } return; } /******************************************************************************/ double *r8vec_uniform_01_new ( int n, int *seed ) /******************************************************************************/ /* Purpose: R8VEC_UNIFORM_01_NEW returns a unit pseudorandom R8VEC. Discussion: This routine implements the recursion seed = 16807 * seed mod ( 2^31 - 1 ) unif = seed / ( 2^31 - 1 ) The integer arithmetic never requires more than 32 bits, including a sign bit. Licensing: This code is distributed under the MIT license. Modified: 19 August 2004 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Second Edition, Springer, 1987, ISBN: 0387964673, LC: QA76.9.C65.B73. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, December 1986, pages 362-376. Pierre L'Ecuyer, Random Number Generation, in Handbook of Simulation, edited by Jerry Banks, Wiley, 1998, ISBN: 0471134031, LC: T57.62.H37. Peter Lewis, Allen Goodman, James Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, Number 2, 1969, pages 136-143. Parameters: Input, int N, the number of entries in the vector. Input/output, int *SEED, a seed for the random number generator. Output, double R8VEC_UNIFORM_01_NEW[N], the vector of pseudorandom values. */ { int i; int i4_huge = 2147483647; int k; double *r; if ( *seed == 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8VEC_UNIFORM_01_NEW - Fatal error!\n" ); fprintf ( stderr, " Input value of SEED = 0.\n" ); exit ( 1 ); } r = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + i4_huge; } r[i] = ( double ) ( *seed ) * 4.656612875E-10; } return r; } /******************************************************************************/ void rule_adjust ( double a, double b, double c, double d, int n, double x[], double w[] ) /******************************************************************************/ /* Purposse: RULE_ADJUST adjusts a 1D quadrature rule from [A,B] to [C,D]. Discussion: This function is only appropriate for cases involving finite intervals and a uniform weighting function. Licensing: This code is distributed under the MIT license. Modified: 21 February 2014 Author: John Burkardt Parameters: Input, double A, B, the left and right endpoints of the original interval. Input, double C, D, the left and right endpoints of the new interval. Input, int N, the order of the rule. 1 <= N. Input/output, double X[N], the abscissas. Input/output, double W[N], the weights. */ { int i; double s; for ( i = 0; i < n; i++ ) { x[i] = ( ( b - x[i] ) * c + ( x[i] - a ) * d ) / ( b - a ); } s = ( d - c ) / ( b - a ); for ( i = 0; i < n; i++ ) { w[i] = s * w[i]; } return; } /******************************************************************************/ void rule_sort ( int m, int n, double x[], double w[] ) /******************************************************************************/ /* Purpose: RULE_SORT sorts a multidimensional quadrature rule. Discussion: This routine reorders the items in the rule so that the points occur in increasing lexicographic order. Licensing: This code is distributed under the MIT license. Modified: 08 December 2012 Author: John Burkardt Parameters: Input, int M, N, the number of rows and columns. Input/output, double X[M*N], W[N]. The points and weights. */ { int i; int indx; int isgn; int j1; int j2; double t; double ww; if ( m <= 0 ) { return; } if ( n <= 1 ) { return; } /* Initialize. */ indx = 0; isgn = 0; j1 = 0; j2 = 0; /* Call the external heap sorter. */ while ( 1 ) { sort_heap_external ( n, &indx, &j1, &j2, isgn ); /* Interchange columns J1 and J2. */ if ( 0 < indx ) { for ( i = 0; i < m; i++ ) { t = x[i+(j1-1)*m]; x[i+(j1-1)*m] = x[i+(j2-1)*m]; x[i+(j2-1)*m] = t; } ww = w[j1-1]; w[j1-1] = w[j2-1]; w[j2-1] = ww; } /* Compare columns J1 and J2. */ else if ( indx < 0 ) { isgn = 0; for ( i = 0; i < m; i++ ) { if ( x[i+(j1-1)*m] < x[i+(j2-1)*m] ) { isgn = -1; break; } else if ( x[i+(j2-1)*m] < x[i+(j1-1)*m] ) { isgn = +1; break; } } } /* The columns are sorted. */ else if ( indx == 0 ) { break; } } return; } /******************************************************************************/ void sort_heap_external ( int n, int *indx, int *i, int *j, int isgn ) /******************************************************************************/ /* Purpose: SORT_HEAP_EXTERNAL externally sorts a list of items into ascending order. Discussion: The actual list is not passed to the routine. Hence it may consist of integers, reals, numbers, names, etc. The user, after each return from the routine, will be asked to compare or interchange two items. The current version of this code mimics the FORTRAN version, so the values of I and J, in particular, are FORTRAN indices. Licensing: This code is distributed under the MIT license. Modified: 05 February 2004 Author: Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf. C version by John Burkardt. Reference: Albert Nijenhuis, Herbert Wilf, Combinatorial Algorithms, Academic Press, 1978, second edition, ISBN 0-12-519260-6. Parameters: Input, int N, the length of the input list. Input/output, int *INDX. The user must set INDX to 0 before the first call. On return, if INDX is greater than 0, the user must interchange items I and J and recall the routine. If INDX is less than 0, the user is to compare items I and J and return in ISGN a negative value if I is to precede J, and a positive value otherwise. If INDX is 0, the sorting is done. Output, int *I, *J. On return with INDX positive, elements I and J of the user's list should be interchanged. On return with INDX negative, elements I and J are to be compared by the user. Input, int ISGN. On return with INDX negative, the user should compare elements I and J of the list. If item I is to precede item J, set ISGN negative, otherwise set ISGN positive. */ { static int i_save = 0; static int j_save = 0; static int k = 0; static int k1 = 0; static int n1 = 0; /* INDX = 0: This is the first call. */ if ( *indx == 0 ) { i_save = 0; j_save = 0; k = n / 2; k1 = k; n1 = n; } /* INDX < 0: The user is returning the results of a comparison. */ else if ( *indx < 0 ) { if ( *indx == -2 ) { if ( isgn < 0 ) { i_save = i_save + 1; } j_save = k1; k1 = i_save; *indx = -1; *i = i_save; *j = j_save; return; } if ( 0 < isgn ) { *indx = 2; *i = i_save; *j = j_save; return; } if ( k <= 1 ) { if ( n1 == 1 ) { i_save = 0; j_save = 0; *indx = 0; } else { i_save = n1; j_save = 1; n1 = n1 - 1; *indx = 1; } *i = i_save; *j = j_save; return; } k = k - 1; k1 = k; } /* 0 < INDX: the user was asked to make an interchange. */ else if ( *indx == 1 ) { k1 = k; } for ( ; ; ) { i_save = 2 * k1; if ( i_save == n1 ) { j_save = k1; k1 = i_save; *indx = -1; *i = i_save; *j = j_save; return; } else if ( i_save <= n1 ) { j_save = i_save + 1; *indx = -2; *i = i_save; *j = j_save; return; } if ( k <= 1 ) { break; } k = k - 1; k1 = k; } if ( n1 == 1 ) { i_save = 0; j_save = 0; *indx = 0; *i = i_save; *j = j_save; } else { i_save = n1; j_save = 1; n1 = n1 - 1; *indx = 1; *i = i_save; *j = j_save; } return; } /******************************************************************************/ int symmetric_sparse_size ( int nr, int dim, double nodes[], double x0 ) /******************************************************************************/ /* Purpose: SYMMETRIC_SPARSE_SIZE sizes a symmetric sparse rule. Licensing: This code is distributed under the MIT license. Modified: 08 December 2012 Author: John Burkardt. Reference: Florian Heiss, Viktor Winschel, Likelihood approximation by numerical integration on sparse grids, Journal of Econometrics, Volume 144, 2008, pages 62-80. Parameters: Input, int DIM, the dimension. Input, int NR, the dimension of the rule in the positive orthant. Input, double NODES[NR*DIM], the nodes for the positive orthant. Input, double X0, the point of symmetry for the 1D rule, typically 0. Output, int SYMMETRIC_SPARSE_SIZE, the dimension of the rule when "unfolded" to the full space. */ { int count; int j; int nr2; int r; /* Count the size of the full rule. */ nr2 = 0; for ( r = 0; r < nr; r++ ) { count = 1; for ( j = 0; j < dim; j++ ) { if ( nodes[r+j*nr] != x0 ) { count = 2 * count; } } nr2 = nr2 + count; } return nr2; } /******************************************************************************/ void tensor_product ( int d, int order1d[], int n1d, double x1d[], double w1d[], int n, double xnd[], double wnd[] ) /******************************************************************************/ /* Purpose: TENSOR_PRODUCT generates a tensor product quadrature rule. Discussion: The Kronecker product of an K by L matrix A and an M by N matrix B is the K*M by L*N matrix formed by a(1,1) * B, a(1,2) * B, ..., a(1,l) * B a(2,1) * B, a(2,2) * B, ..., a(2,l) * B .......... .......... .... .......... a(k,1) * B, a(k,2) * B, ..., a(k,l) * B Licensing: This code is distributed under the MIT license. Modified: 09 December 2012 Author: Original MATLAB version by Florian Heiss, Viktor Winschel. C version by John Burkardt. Reference: Florian Heiss, Viktor Winschel, Likelihood approximation by numerical integration on sparse grids, Journal of Econometrics, Volume 144, 2008, pages 62-80. Parameters: Input, int D, the spatial dimension. Input, int ORDER1D[D], the order of each 1D rule. Input, int N1D, the number of 1D items. Input, double X1D[N1D], the 1D nodes. Input, double W1D[N1D], the 1D weights. Input, int N, the number of N-dimensional items. Output, double XND[D*N], the nodes. Output, double WND[N], the weights. */ { int i; int i1; int i2; /* Compute the weights. */ i2 = -1; for ( i = 0; i < d; i++ ) { i1 = i2 + 1; i2 = i2 + order1d[i]; r8vec_direct_product2 ( i, order1d[i], w1d+i1, d, n, wnd ); } /* Compute the points. */ i2 = -1; for ( i = 0; i < d; i++ ) { i1 = i2 + 1; i2 = i2 + order1d[i]; r8vec_direct_product ( i, order1d[i], x1d+i1, d, n, xnd ); } return; } /******************************************************************************/ void tensor_product_cell ( int nc, double xc[], double wc[], int dim, int nr[], int roff[], int np, double xp[], double wp[] ) /******************************************************************************/ /* Purpose: TENSOR_PRODUCT_CELL generates a tensor product quadrature rule. Discussion: The Kronecker product of an K by L matrix A and an M by N matrix B is the K*M by L*N matrix formed by a(1,1) * B, a(1,2) * B, ..., a(1,l) * B a(2,1) * B, a(2,2) * B, ..., a(2,l) * B .......... .......... .... .......... a(k,1) * B, a(k,2) * B, ..., a(k,l) * B The 1D factors are stored in a kind of cell array structure. Licensing: This code is distributed under the MIT license. Modified: 10 December 2012 Author: John Burkardt. Parameters: Input, int NC, the number of items in the cell arrays. Input, double XC[NC], a cell array containing points for 1D rules. Input, double WC[NC], a cell array containing weights for 1D rules. Input, int DIM, the spatial dimension. Input, int NR[DIM], the length of each row of the cell array. Input, int ROFF[DIM+1], offsets for the cell arrays. Input, int NP, the number of points in the product rule. Output, double XP[DIM*NP], the nodes. Output, double WP[NP], the weights. */ { int i; int n1d; double *w1d; double *x1d; /* Compute the weights. */ for ( i = 0; i < dim; i++ ) { n1d = nr[i]; w1d = r8cvv_rget_new ( nc, wc, dim, roff, i ); r8vec_direct_product2 ( i, n1d, w1d, dim, np, wp ); free ( w1d ); } /* Compute the points. */ for ( i = 0; i < dim; i++ ) { n1d = nr[i]; x1d = r8cvv_rget_new ( nc, xc, dim, roff, i ); r8vec_direct_product ( i, n1d, x1d, dim, np, xp ); free ( x1d ); } return; } /******************************************************************************/ void timestamp ( void ) /******************************************************************************/ /* Purpose: TIMESTAMP prints the current YMDHMS date as a time stamp. Example: 31 May 2001 09:45:54 AM Licensing: This code is distributed under the MIT license. Modified: 24 September 2003 Author: John Burkardt Parameters: None */ { # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct tm *tm; time_t now; now = time ( NULL ); tm = localtime ( &now ); strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm ); fprintf ( stdout, "%s\n", time_buffer ); return; # undef TIME_SIZE }