rule.c

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

  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];
    if( g[dim - 1] > gd ) {
      count i, ix, j = dim, dx = dim - 1;
      for( i = 0; i < --j; ++i ) {
        creal tmp = g[i];
        g[i] = g[j];
        g[j] = tmp;
        if( tmp <= gd ) --dx;
        if( g[i] > gd ) ix = i;
      }
      if( g[dx] <= gd ) dx = ix;
      g[dim] = g[dx];
      g[dx] = 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 Sample(cRule *rule, void *voidregion, cint flags)
{
  TYPEDEFREGION;
  TYPEDEFSET;

  Region *const region = (Region *)voidregion;
  creal vol = ldexp(1., -region->div);

  real *x = rule->x, *f = rule->f;
  Set *first = (Set *)rule->first, *last = (Set *)rule->last, *s;
  creal *errcoeff = rule->errcoeff;
  creal ratio = Sq(first[2].gen[0]/first[1].gen[0]);

  ccount offset = 2*ndim_*ncomp_;
  count dim, comp, rul, n, maxdim = 0;
  real maxrange = 0;

  for( dim = 0; dim < ndim_; ++dim ) {
    cBounds *b = &region->bounds[dim];
    creal range = b->upper - b->lower;
    if( range > maxrange ) {
      maxrange = range;
      maxdim = dim;
    }
  }

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

  DoSample(rule->n, rule->x, f);

  for( comp = 0; comp < ncomp_; ++comp ) {
    Result *r = &region->result[comp];
    real sum[nrules];
    creal *f1 = f;
    creal base = *f1*2*(1 - ratio);
    real maxdiff = 0;
    count bisectdim = maxdim;

    for( dim = 0; dim < ndim_; ++dim ) {
      creal *fp = f1 + ncomp_;
      creal *fm = fp + ncomp_;
      creal fourthdiff = fabs(base +
        ratio*(fp[0] + fm[0]) - (fp[offset] + fm[offset]));
      f1 = fm;
      if( fourthdiff > maxdiff ) {
        maxdiff = fourthdiff;
        bisectdim = dim;
      }
    }
    r->bisectdim = bisectdim;

    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;
    }

    r->avg = vol*sum[0];
    r->err = 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]) );
  }

  if( VERBOSE > 2 ) {
    char s[64*NDIM + 128*NCOMP], *p = s;

    for( dim = 0; dim < ndim_; ++dim ) {
      cBounds *b = &region->bounds[dim];
      p += sprintf(p,
        (dim == 0) ? "\nRegion (" REALF ") - (" REALF ")" :
                     "\n       (" REALF ") - (" REALF ")",
        b->lower, b->upper);
    }

    for( comp = 0; comp < ncomp_; ++comp ) {
      cResult *r = &region->result[comp];
      p += sprintf(p, "\n[" COUNT "] "
        REAL " +- " REAL, comp + 1, r->avg, r->err);
    }

    Print(s);
  }
}