sparse_grid_hw_jl.c

# include <stdlib.h>
# include <stdio.h>
# include <math.h>
# include <time.h>
# include <string.h>

# include "sparse_grid_hw_jl.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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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.
*/
{
  int i;
  int j;
  int jhi;
  int jlo;

  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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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.
*/
{
  int i;
  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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 *keep;
  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 GNU LGPL 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 n1d_total;
  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;
  }

  n1d_total = i4vec_sum ( k, n1d );
/*
  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 GNU LGPL 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_abs ( double x )

/******************************************************************************/
/*
  Purpose:

    R8_ABS returns the absolute value of an R8.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    07 May 2006

  Author:

    John Burkardt

  Parameters:

    Input, double X, the quantity whose absolute value is desired.

    Output, double R8_ABS, the absolute value of X.
*/
{
  double value;

  if ( 0.0 <= x )
  {
    value = + x;
  }
  else
  {
    value = - x;
  }
  return value;
}
/******************************************************************************/

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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 j;
  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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL license.

  Modified:

    24 September 2003

  Author:

    John Burkardt

  Parameters:

    None
*/
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct tm *tm;
  size_t len;
  time_t now;

  now = time ( NULL );
  tm = localtime ( &now );

  len = strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  fprintf ( stdout, "%s\n", time_buffer );

  return;
# undef TIME_SIZE
}

/******************************************************************************/

int nwspgr_rule_size ( void rule ( int n, double x[], double w[] ), 
  int rule_order ( int l ), int dim, int k, int r_size)

/******************************************************************************/
/*
  Purpose:

    NWSPGR generates size of rule for sparse grid integration.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    21 April 2015

  Author:

    Original MATLAB version by Florian Heiss, Viktor Winschel.
    C version by John Burkardt. Adapted by Alan Crawford.

  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.

*/
{
  int bq;
  int equal;
  int i;
  int *is;
  int j;
  int j2;
  int *keep;
  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;
  double *nodes;
  double *weights;

  nodes = ( double * ) malloc ( dim * r_size * sizeof ( double ) );
  weights = ( double * ) malloc ( r_size * sizeof ( double ) );

  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;

  free(nodes);
  free(weights);

  return r;
}
/********************************************/

int nwspgr_rule_size_wrapper(int rule , int d, int k)

/********************************************/
/* 
Purpose : Calculate size of the Sparse Grid 
*/
{
int n;
int n2;  
  if ( rule==1 ) 
  {
    n = nwspgr_size ( gqu_order, d, k );
    n2 = nwspgr_rule_size ( gqu, gqu_order, d, k, n);
  }
  else if ( rule==2 )
  {
    n = nwspgr_size ( kpu_order, d, k );
    n2 = nwspgr_rule_size ( kpu, kpu_order, d, k, n);
  }
  else if ( rule==3 ) 
  {
    n = nwspgr_size ( gqn_order, d, k );
    n2 = nwspgr_rule_size ( gqn, gqn_order, d, k, n);
  }
  else 
  {
    n = nwspgr_size ( kpn_order, d, k );
    n2 = nwspgr_rule_size ( kpn, kpn_order, d, k, n);
   }

  return n2;
}

/************************************************************/

void nwspgr_wrapper(int rule, int d, int k, double x[], double w[])

/************************************************************/
/* Purpose: Wrapper for computation of Sparse Grid for integration 

Source: www.sparse-grids.de -> 

To approximate a D-variate integral of the form ∫Ωg(x)f(x)dx where Ω⊆R^D and f(x) is a weight function, 
a quadrature rule prescribes a RxD matrix of nodes X and a Rx1 vector of weights w depending on Ω and f.
 The approximation is calculated as g(X)'w.

In all three versions provided here, four quadrature rules for D-dimensional numerical integration are implemented:

Rules for the unweighted integration: Ω=[0,1]^D, f(x)=1
-------------------------------------------------------
GQU: Univariate Gaussian quadrature rules as basis
KPU: Univariate nested quadrature rules as basis - delayed Kronrod-Patterson rules, 
      see Knut Petras (2003): "Smolyak cubature of given polynomial degree with few nodes for increasing dimension." Numerische Mathematik 93, 729-753.
      
Rules for the integration over Ω=R^D with Gaussian weight f(x)=exp(-x'x/2)/sqrt(2π)
----------------------------------------------------------------------------------
GQN: Univariate Gaussian quadrature rules as basis
KPN: Univariate nested quadrature rules as basis, 
      see A. Genz and B. D. Keister (1996): "Fully symmetric interpolatory rules for multiple integrals over infinite regions with Gaussian weight." Journal of Computational and Applied Mathematics 71, 299-309.

Inputs
------
rule = GQU==1, KPU==2, GQN==3, KPN==4
d = dimensions
k = level

Intermediate
------------
n2 is number of sparse points after duplicates removed. 
This is output of nwspgr, as well as x and w

Outputs
-------
n2 is the number of grid points
x = grid of nodes
w = grid of weights

*/
{
double *nodes;
double *weights; 
int n;
int n2;
int j;
  if ( rule==1 ) 
  {
    n = nwspgr_size ( gqu_order, d, k );
    n2 = nwspgr_rule_size_wrapper( rule , d, k );
    nodes = ( double * ) malloc (  d * n * sizeof ( double ) );
    weights = ( double * ) malloc ( n * sizeof ( double ) );
    nwspgr ( gqu, gqu_order, d, k, n, &n2, nodes, weights );
  }
  else if ( rule==2 )
  {
    n = nwspgr_size ( kpu_order, d, k );
    n2 = nwspgr_rule_size_wrapper( rule , d, k );    
    nodes = ( double * ) malloc ( d * n * sizeof ( double ) );
    weights = ( double * ) malloc ( n * sizeof ( double ) );
    nwspgr ( kpu, kpu_order, d, k, n, &n2, nodes, weights );
  }
  else if ( rule==3 ) 
  {
    n = nwspgr_size ( gqn_order, d, k );
    n2 = nwspgr_rule_size_wrapper( rule , d, k );    
    nodes = ( double * ) malloc ( d * n * sizeof ( double ) );
    weights = ( double * ) malloc ( n * sizeof ( double ) );
    nwspgr ( gqn, gqn_order, d, k, n, &n2, nodes, weights);
  }
  else 
  {
    n = nwspgr_size ( kpn_order, d, k );
    n2 = nwspgr_rule_size_wrapper( rule, d, k );        
    nodes = ( double * ) malloc ( d * n * sizeof ( double ) );
    weights = ( double * ) malloc ( n * sizeof ( double ) );
    nwspgr ( kpn, kpn_order, d, k, n, &n2, nodes, weights );
  }
  for ( j = 0; j < d*n2; j++ )
    {
      x[j] = nodes[j];
    }
  for ( j = 0; j < n2; j++ )
    {
      w[j] = weights[j];
  }
  free(nodes);
  free(weights);
}