rule.c

The CUBA library provides new implementation of four general-purpose multidimensional integration al

     .89567365764160676508,    -.36479356986049146661,
     .0035417564516782676826,  -.072609367395893679605,
     .10557491625218991012,     .0021486025550098687713,
    -.032268563892953949998,    .010636783990231217481,
     .014689102496143490175,    .51134708346467591431,
     .45976448120806344646,     .18239678493024573331,
    -.04508628929435784076,     .21415883524352793401,
    -.027351546526545644722,    .054941067048711234101,
     .11937596202570775297,     .65089519391920250593,
     .14744939829434460168,     .057693384490973483573,
     .034999626602143583822,  -1.3868627719278281436,
    -.2386668732575008879,      .015532417276607053264,
     .0035328099607090870236,   .09231719987444221619,
     .02254314464717892038,     .013675773263272822361,
    -.32544759695960125297,     .0017708782258391338413,
     .0010743012775049343856,   .25150011495314791996 };

  static creal g[] = {
     .47795365790226950619,     .20302858736911986780,
     .44762735462617812882,     .125,
     .34303789878087814570 };

  enum { nsets = 9 };

  TYPEDEFSET;

  ccount ndim = ndim_;
  ccount twondim = 1 << ndim;
  count dim, n, r;
  Set *first, *last, *s, *t;

  Alloc(first, nsets);
  Clear(first, nsets);

  last = first;
  n = last->n = 1;
  last->weight[0] = ndim*(ndim*(ndim*w[0] + w[1]) + w[2]) + w[3];
  last->weight[1] = ndim*(ndim*(ndim*w[4] + w[5]) + w[6]) - w[7];
  last->weight[2] = ndim*w[8] - last->weight[1];
  last->weight[3] = ndim*(ndim*w[9] + w[10]) - 1 + last->weight[0];
  last->weight[4] = ndim*w[11] + 1 - last->weight[0];

  ++last;
  n += last->n = 2*ndim;
  last->weight[0] = ndim*(ndim*w[12] + w[13]) + w[14];
  last->weight[1] = ndim*(ndim*w[15] + w[16]) + w[17];
  last->weight[2] = w[18] - last->weight[1];
  last->weight[3] = ndim*w[19] + w[20] + last->weight[0];
  last->weight[4] = w[21] - last->weight[0];
  last->gen[0] = g[0];

  ++last;
  n += last->n = 2*ndim;
  last->weight[0] = ndim*w[22] + w[23];
  last->weight[1] = ndim*w[24] + w[25];
  last->weight[2] = w[26] - last->weight[1];
  last->weight[3] = ndim*w[27] + w[28];
  last->weight[4] = -last->weight[0];
  last->gen[0] = g[1];

  ++last;
  n += last->n = 2*ndim;
  last->weight[0] = w[29];
  last->weight[1] = w[30];
  last->weight[2] = -w[29];
  last->weight[3] = w[31];
  last->weight[4] = -w[29];
  last->gen[0] = g[2];

  ++last;
  n += last->n = 2*ndim;
  last->weight[2] = w[32];
  last->gen[0] = g[3];

  ++last;
  n += last->n = 2*ndim*(ndim - 1);
  last->weight[0] = w[33] - ndim*w[12];
  last->weight[1] = w[34] - ndim*w[15];
  last->weight[2] = -last->weight[1];
  last->weight[3] = w[35] + last->weight[0];
  last->weight[4] = -last->weight[0];
  last->gen[0] = g[0];
  last->gen[1] = g[0];

  ++last;
  n += last->n = 4*ndim*(ndim - 1);
  last->weight[0] = w[36];
  last->weight[1] = w[37];
  last->weight[2] = -w[37];
  last->weight[3] = w[38];
  last->weight[4] = -w[36];
  last->gen[0] = g[0];
  last->gen[1] = g[1];

  ++last;
  n += last->n = 4*ndim*(ndim - 1)*(ndim - 2)/3;
  last->weight[0] = w[39];
  last->weight[1] = w[40];
  last->weight[2] = -w[40];
  last->weight[3] = w[39];
  last->weight[4] = -w[39];
  last->gen[0] = g[0];
  last->gen[1] = g[0];
  last->gen[2] = g[0];

  ++last;
  n += last->n = twondim;
  last->weight[0] = w[41]/twondim;
  last->weight[1] = w[7]/twondim;
  last->weight[2] = -last->weight[1];
  last->weight[3] = last->weight[0];
  last->weight[4] = -last->weight[0];
  for( dim = 0; dim < ndim; ++dim )
    last->gen[dim] = g[4];

  rule->first = first;
  rule->last = last;
  rule->errcoeff[0] = 5;
  rule->errcoeff[1] = 1;
  rule->errcoeff[2] = 5;
  rule->n = n;

  for( s = first; s <= last; ++s )
    for( r = 1; r < nrules - 1; ++r ) {
      creal scale = (s->weight[r] == 0) ? 100 :
        -s->weight[r + 1]/s->weight[r];
      real sum = 0;
      for( t = first; t <= last; ++t )
        sum += t->n*fabs(t->weight[r + 1] + scale*t->weight[r]);
      s->scale[r] = scale;
      s->norm[r] = 1/sum;
    }
}

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

static void Rule7Alloc(Rule *rule)
{
  static creal w[] = {
     .019417866674748388428,   -.40385257701150182546,
     .64485668767465982223,     .01177982690775806141,
    -.18041318740733609012,    -.088785828081335044443,
     .056328645808285941374,   -.0097089333373741942142,
    -.99129176779582358138,    -.17757165616267008889,
     .12359398032043233572,     .074978148702033690681,
     .55489147051423559776,     .088041241522692771226,
     .021118358455513385083,   -.0099302203239653333087,
    -.064100053285010904179,    .030381729038221007659,
     .0058899134538790307051,  -.0048544666686870971071,
     .35514331232534017777 };

  static creal g[] = {
     .47795365790226950619,     .20302858736911986780,
     .375,                      .34303789878087814570 };

  enum { nsets = 6 };

  TYPEDEFSET;

  ccount ndim = ndim_;
  ccount twondim = 1 << ndim;
  count dim, n, r;
  Set *first, *last, *s, *t;

  Alloc(first, nsets);
  Clear(first, nsets);

  last = first;
  n = last->n = 1;
  last->weight[0] = ndim*(ndim*w[0] + w[1]) + w[2];
  last->weight[1] = ndim*(ndim*w[3] + w[4]) - w[5];
  last->weight[2] = ndim*w[6] - last->weight[1];
  last->weight[3] = ndim*(ndim*w[7] + w[8]) - w[9];
  last->weight[4] = 1 - last->weight[0];

  ++last;
  n += last->n = 2*ndim;
  last->weight[0] = w[10];
  last->weight[1] = w[11];
  last->weight[2] = -w[10];
  last->weight[3] = w[12];
  last->weight[4] = -w[10];
  last->gen[0] = g[1];

  ++last;
  n += last->n = 2*ndim;
  last->weight[0] = w[13] - ndim*w[0];
  last->weight[1] = w[14] - ndim*w[3];
  last->weight[2] = w[15] - last->weight[1];
  last->weight[3] = w[16] - ndim*w[7];
  last->weight[4] = -last->weight[0];
  last->gen[0] = g[0];

  ++last;
  n += last->n = 2*ndim;
  last->weight[2] = w[17];
  last->gen[0] = g[2];

  ++last;
  n += last->n = 2*ndim*(ndim - 1);
  last->weight[0] = -w[7];
  last->weight[1] = w[18];
  last->weight[2] = -w[18];
  last->weight[3] = w[19];
  last->weight[4] = w[7];
  last->gen[0] = g[0];
  last->gen[1] = g[0];

  ++last;
  n += last->n = twondim;
  last->weight[0] = w[20]/twondim;
  last->weight[1] = w[5]/twondim;
  last->weight[2] = -last->weight[1];
  last->weight[3] = w[9]/twondim;
  last->weight[4] = -last->weight[0];
  for( dim = 0; dim < ndim; ++dim )
    last->gen[dim] = g[3];

  rule->first = first;
  rule->last = last;
  rule->errcoeff[0] = 5;
  rule->errcoeff[1] = 1;
  rule->errcoeff[2] = 5;
  rule->n = n;

  for( s = first; s <= last; ++s )
    for( r = 1; r < nrules - 1; ++r ) {
      creal scale = (s->weight[r] == 0) ? 100 :
        -s->weight[r + 1]/s->weight[r];
      real sum = 0;
      for( t = first; t <= last; ++t )
        sum += t->n*fabs(t->weight[r + 1] + scale*t->weight[r]);
      s->scale[r] = scale;
      s->norm[r] = 1/sum;
    }
}

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

static real *ExpandFS(cBounds *b, real *g, real *x)
{
  count dim, ndim = ndim_;

next:
  /* Compute centrally symmetric sum for permutation of G */

  for( dim = 0; dim < ndim; ++dim )
    *x++ = (.5 + g[dim])*b[dim].lower + (.5 - g[dim])*b[dim].upper;

  for( dim = 0; dim < ndim; ) {
    g[dim] = -g[dim];
    if( g[dim++] < 0 ) goto next;
  }

  /* Find next distinct permutation of G and loop back for next sum.
     Permutations are generated in reverse lexicographic order. */

  for( dim = 1; dim < ndim; ++dim ) {
    creal gd = g[dim];
    count big = dim - 1;
    if( g[big] > gd ) {
      count i, j = dim, lastbig = big;
      for( i = 0; i < --j; ++i ) {
        creal tmp = g[i];
        g[i] = g[j];
        g[j] = tmp;
        if( tmp <= gd ) --big;
        if( g[i] > gd ) lastbig = i;
      }
      if( g[big] <= gd ) big = lastbig;
      g[dim] = g[big];
      g[big] = gd;
      goto next;
    }
  }

  /* Restore original order to generators */

  for( dim = 0; dim < --ndim; ++dim ) {
    creal tmp = g[dim];
    g[dim] = g[ndim];
    g[ndim] = tmp;
  }

  return x;
}

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

static void SampleRule(cSamples *samples, cBounds *b, creal vol)
{
  TYPEDEFSET;

  real *x = samples->x, *f = samples->f;
  Set *first = (Set *)samples->rule->first;
  Set *last = (Set *)samples->rule->last;
  Set *s;
  creal *errcoeff = samples->rule->errcoeff;
  count comp, rul, n;

  for( s = first; s <= last; ++s )
    if( s->n ) x = ExpandFS(b, s->gen, x);

  DoSample(samples->n, ndim_, samples->x, f);

  for( comp = 0; comp < ncomp_; ++comp ) {
    real sum[nrules];
    creal *f1 = f++;

    Zap(sum);
    for( s = first; s <= last; ++s )
      for( n = s->n; n; --n ) {
        creal fun = *f1;
        f1 += ncomp_;
        for( rul = 0; rul < nrules; ++rul )
          sum[rul] += fun*s->weight[rul];
      }

    /* Search for the null rule, in the linear space spanned by two
       successive null rules in our sequence, which gives the greatest
       error estimate among all normalized (1-norm) null rules in this
       space. */

    for( rul = 1; rul < nrules - 1; ++rul ) {
      real maxerr = 0;
      for( s = first; s <= last; ++s )
        maxerr = Max(maxerr,
          fabs(sum[rul + 1] + s->scale[rul]*sum[rul])*s->norm[rul]);
      sum[rul] = maxerr;
    }

    samples->avg[comp] = vol*sum[0];
    samples->err[comp] = vol*(
      (errcoeff[0]*sum[1] <= sum[2] && errcoeff[0]*sum[2] <= sum[3]) ?
        errcoeff[1]*sum[1] :
        errcoeff[2]*Max(Max(sum[1], sum[2]), sum[3]));
  }
}