diff options
Diffstat (limited to 'src/divonne')
| -rw-r--r-- | src/divonne/Divonne.c | 157 | ||||
| -rw-r--r-- | src/divonne/Divonne.tm | 429 | ||||
| -rw-r--r-- | src/divonne/Explore.c | 160 | ||||
| -rw-r--r-- | src/divonne/FindMinimum.c | 689 | ||||
| -rw-r--r-- | src/divonne/Integrate.c | 452 | ||||
| -rw-r--r-- | src/divonne/Iterate.c | 138 | ||||
| -rw-r--r-- | src/divonne/KorobovCoeff.c | 881 | ||||
| -rw-r--r-- | src/divonne/Rule.c | 685 | ||||
| -rw-r--r-- | src/divonne/Sample.c | 261 | ||||
| -rw-r--r-- | src/divonne/Split.c | 311 | ||||
| -rw-r--r-- | src/divonne/common.c | 34 | ||||
| -rw-r--r-- | src/divonne/decl.h | 127 |
12 files changed, 4324 insertions, 0 deletions
diff --git a/src/divonne/Divonne.c b/src/divonne/Divonne.c new file mode 100644 index 0000000..5fcaada --- /dev/null +++ b/src/divonne/Divonne.c @@ -0,0 +1,157 @@ +/* + Divonne.c + Multidimensional integration by partitioning + originally by J.H. Friedman and M.H. Wright + (CERNLIB subroutine D151) + this version by Thomas Hahn + last modified 15 Nov 11 th +*/ + +#include "decl.h" + +#define Print(s) puts(s); fflush(stdout) + +/*********************************************************************/ + +#define DIVONNE +static void DoSample(This *t, number n, creal *x, real *f, ccount ldx); +static int ExploreParent(This *t, cint iregion); + +/*********************************************************************/ + +static inline count SampleExtra(This *t, cBounds *b) +{ + number n = t->nextra; + t->peakfinder(&t->ndim, b, &n, t->xextra); + DoSample(t, n, t->xextra, t->fextra, t->ldxgiven); + return n; +} + +/*********************************************************************/ + +static inline void AllocGiven(This *t, creal *xgiven) +{ + if( t->ngiven | t->nextra ) { + cnumber nxgiven = t->ngiven*(t->ldxgiven = IMax(t->ldxgiven, t->ndim)); + cnumber nxextra = t->nextra*t->ldxgiven; + cnumber nfgiven = t->ngiven*t->ncomp; + cnumber nfextra = t->nextra*t->ncomp; + + Alloc(t->xgiven, nxgiven + nxextra + nfgiven + nfextra); + t->xextra = t->xgiven + nxgiven; + t->fgiven = t->xextra + nxextra; + t->fextra = t->fgiven + nfgiven; + + if( nxgiven ) { + t->phase = 0; + Copy(t->xgiven, xgiven, nxgiven); + DoSample(t, t->ngiven, t->xgiven, t->fgiven, t->ldxgiven); + } + } +} + +/*********************************************************************/ + +#include "common.c" +#include "DoSample.c" + +Extern void EXPORT(Divonne)(ccount ndim, ccount ncomp, + Integrand integrand, void *userdata, + creal epsrel, creal epsabs, + cint flags, cint seed, + cnumber mineval, cnumber maxeval, + cint key1, cint key2, cint key3, ccount maxpass, + creal border, creal maxchisq, creal mindeviation, + cnumber ngiven, ccount ldxgiven, creal *xgiven, + cnumber nextra, PeakFinder peakfinder, + int *pnregions, number *pneval, int *pfail, + real *integral, real *error, real *prob) +{ + This t; + t.ndim = ndim; + t.ncomp = ncomp; + t.integrand = integrand; + t.userdata = userdata; + t.epsrel = epsrel; + t.epsabs = epsabs; + t.flags = flags; + t.seed = seed; + t.mineval = mineval; + t.maxeval = maxeval; + t.key1 = key1; + t.key2 = key2; + t.key3 = key3; + t.maxpass = maxpass; + t.border.upper = 1 - (t.border.lower = border); + t.maxchisq = maxchisq; + t.mindeviation = mindeviation; + t.ngiven = ngiven; + t.xgiven = NULL; + t.ldxgiven = ldxgiven; + t.nextra = nextra; + t.peakfinder = peakfinder; + t.nregions = 0; + t.neval = 0; + + ForkCores(&t); + AllocGiven(&t, xgiven); + + *pfail = Integrate(&t, integral, error, prob); + *pnregions = t.nregions; + *pneval = t.neval; + + free(t.xgiven); + WaitCores(&t); +} + +/*********************************************************************/ + +Extern void EXPORT(divonne)(ccount *pndim, ccount *pncomp, + Integrand integrand, void *userdata, + creal *pepsrel, creal *pepsabs, + cint *pflags, cint *pseed, + cnumber *pmineval, cnumber *pmaxeval, + cint *pkey1, cint *pkey2, cint *pkey3, ccount *pmaxpass, + creal *pborder, creal *pmaxchisq, creal *pmindeviation, + cnumber *pngiven, ccount *pldxgiven, creal *xgiven, + cnumber *pnextra, PeakFinder peakfinder, + int *pnregions, number *pneval, int *pfail, + real *integral, real *error, real *prob) +{ + This t; + t.ndim = *pndim; + t.ncomp = *pncomp; + t.integrand = integrand; + t.userdata = userdata; + t.epsrel = *pepsrel; + t.epsabs = *pepsabs; + t.flags = *pflags; + t.seed = *pseed; + t.mineval = *pmineval; + t.maxeval = *pmaxeval; + t.key1 = *pkey1; + t.key2 = *pkey2; + t.key3 = *pkey3; + t.maxpass = *pmaxpass; + t.border.upper = 1 - (t.border.lower = *pborder); + t.maxchisq = *pmaxchisq; + t.mindeviation = *pmindeviation; + t.ngiven = *pngiven; + t.xgiven = NULL; + t.ldxgiven = *pldxgiven; + t.nextra = *pnextra; + t.peakfinder = peakfinder; + t.nregions = 0; + t.neval = 0; + + ForkCores(&t); + AllocGiven(&t, xgiven); + + *pfail = Integrate(&t, integral, error, prob); + *pnregions = t.nregions; + *pneval = t.neval; + + free(t.xgiven); + WaitCores(&t); +} + diff --git a/src/divonne/Divonne.tm b/src/divonne/Divonne.tm new file mode 100644 index 0000000..bcb8b88 --- /dev/null +++ b/src/divonne/Divonne.tm @@ -0,0 +1,429 @@ +:Evaluate: BeginPackage["Cuba`"] + +:Evaluate: Divonne::usage = + "Divonne[f, {x, xmin, xmax}..] computes a numerical approximation to the integral of the real scalar or vector function f. + The output is a list with entries of the form {integral, error, chi-square probability} for each component of the integrand." + +:Evaluate: Key1::usage = "Key1 is an option of Divonne. + It determines sampling in the partitioning phase.\n + Special cases:\n + Key1 = 7: use a degree-7 cubature rule,\n + Key1 = 9: use a degree-9 cubature rule,\n + Key1 = 11: use a degree-11 cubature rule (available only in 3 dimensions),\n + Key1 = 13: use a degree-13 cubature rule (available only in 2 dimensions),\n + otherwise a random sample of n1 = Abs[Key1] points is used, where the sign of Key1 determines the type of sample:\n + Key1 > 0: use a Korobov quasi-random sample,\n + Key1 < 0: use a \"standard\" sample." + +:Evaluate: Key2::usage = "Key2 is an option of Divonne. + It determines sampling in the main integration phase.\n + Special cases:\n + Key2 = 7: use a degree-7 cubature rule,\n + Key2 = 9: use a degree-9 cubature rule,\n + Key2 = 11: use a degree-11 cubature rule (available only in 3 dimensions),\n + Key2 = 13: use a degree-13 cubature rule (available only in 2 dimensions),\n + otherwise a random sample is used, where the sign of Key2 determines the type of sample:\n + Key2 > 0: use a Korobov quasi-random sample,\n + Key2 < 0: use a \"standard\" sample,\n + and n2 = Abs[Key2] determines the number of points:\n + n2 >= 40: sample n2 points,\n + n2 < 40: sample n2*nneed points, where nneed is the number of points needed to reach the prescribed accuracy, as estimated by Divonne from the results of the partitioning phase." + +:Evaluate: Key3::usage = "Key3 is an option of Divonne. + It sets the strategy for the refinement phase:\n + Key3 = 0: do not further treat the subregion,\n + Key3 = 1: split the subregion up once more,\n + for other values the region is sampled a third time:\n + Key3 = 7: use a degree-7 cubature rule,\n + Key3 = 9: use a degree-9 cubature rule,\n + Key3 = 11: use a degree-11 cubature rule (available only in 3 dimensions),\n + Key3 = 13: use a degree-13 cubature rule (available only in 2 dimensions),\n + otherwise a random sample is used, where the sign of Key3 determines the type of sample:\n + Key3 > 0: use a Korobov quasi-random sample,\n + Key3 < 0: use a \"standard\" sample,\n + and n3 = Abs[Key3] determines the number of points:\n + n3 >= 40: sample n3 points,\n + n3 < 40: sample n3*nneed points, where nneed is the number of points needed to reach the prescribed accuracy, as estimated by Divonne from the results of the partitioning phase." + +:Evaluate: MaxPass::usage = "MaxPass is an option of Divonne. + It controls the partitioning termination. + The partitioning phase is terminated when the estimated total number of integrand evaluations (partitioning plus main integration) does not decrease for MaxPass successive iterations." + +:Evaluate: Border::usage = "Border is an option of Divonne. + It specifies the width of the border of the integration region. + Points falling into this border region are not sampled directly, but are extrapolated from two samples from the interior. + The border width always refers to the unit hypercube, i.e. it is not rescaled if the integration region is not the unit hypercube." + +:Evaluate: MaxChisq::usage = "MaxChisq is an option of Divonne. + It specifies the maximum chi-square value a single subregion is allowed to have in the main integration phase. + Regions which fail this chi-square test and whose sample averages differ by more than MinDeviation move on to the refinement phase." + +:Evaluate: MinDeviation::usage = "MinDeviation is an option of Divonne. + Regions which fail the chi-square test are not treated further if their sample averages differ by less than MinDeviation. + MinDeviation is specified as the fraction of the requested error of the entire integral." + +:Evaluate: Given::usage = "Given is an option of Divonne. + It provides a list of points where the integrand might have peaks. + Divonne will consider these points when partitioning the integration region." + +:Evaluate: NExtra::usage = "NExtra is an option of Divonne. + It specifies the maximum number of points that will be considered in the output of the PeakFinder function." + +:Evaluate: PeakFinder::usage = "PeakFinder is an option of Divonne. + It specifies the peak-finder function. + This function is called whenever a region is up for subdivision and is supposed to point out possible peaks lying in the region, thus acting as the dynamic counterpart of the static list of points supplied with Given. + It is invoked with two arguments, the multidimensional equivalents of the lower left and upper right corners of the region being investigated, and must return a (possibly empty) list of points." + +:Evaluate: MinPoints::usage = "MinPoints is an option of Divonne. + It specifies the minimum number of points to sample." + +:Evaluate: Final::usage = "Final is an option of Divonne. + It can take the values Last or All which determine whether only the last (largest) or all sets of samples collected on a subregion over the integration phases contribute to the final result." + +:Evaluate: PseudoRandom::usage = "PseudoRandom is an option of Divonne. + It can take the following values: + False for Sobol quasi-random numbers (default), + True or 0 for Mersenne Twister pseudo-random numbers, + any other integer value n for Ranlux pseudo-random numbers of luxury level n." + +:Evaluate: PseudoRandomSeed::usage = "PseudoRandomSeed is an option of Divonne. + It specifies the seed for the pseudo-random number generator." + +:Evaluate: Regions::usage = "Regions is an option of Divonne. + It specifies whether the regions into which the integration region has been cut are returned together with the integration results." + +:Evaluate: Region::usage = "Region[ll, ur, res, df] describes a subregion: + ll and ur are multidimensional equivalents of the region's lower left and upper right corner. + res gives the integration results for the region in a list with entries of the form {integral, error, chi-square} for each component of the integrand. + df is the number of degrees of freedom corresponding to the chi-square values in res." + +:Evaluate: $Phase::usage = "$Phase is a global variable set by Divonne during the evaluation of the integrand to the integration phase:\n + 0 = sampling of the points in xgiven,\n + 1 = partitioning phase,\n + 2 = main integration phase,\n + 3 = refinement phase." + +:Evaluate: MapSample::usage = "MapSample is a function used to map the integrand over the points to be sampled." + + +:Evaluate: Begin["`Divonne`"] + +:Begin: +:Function: Divonne +:Pattern: MLDivonne[ndim_, ncomp_, + epsrel_, epsabs_, flags_, seed_, + mineval_, maxeval_, + key1_, key2_, key3_, maxpass_, + border_, maxchisq_, mindeviation_, + xgiven_, fgiven_, nextra_] +:Arguments: {ndim, ncomp, + epsrel, epsabs, flags, seed, + mineval, maxeval, + key1, key2, key3, maxpass, + border, maxchisq, mindeviation, + xgiven, fgiven, nextra} +:ArgumentTypes: {Integer, Integer, + Real, Real, Integer, Integer, + Integer, Integer, + Integer, Integer, Integer, Integer, + Real, Real, Real, + RealList, RealList, Integer} +:ReturnType: Manual +:End: + +:Evaluate: Attributes[Divonne] = {HoldFirst} + +:Evaluate: Options[Divonne] = {PrecisionGoal -> 3, AccuracyGoal -> 12, + MinPoints -> 0, MaxPoints -> 50000, + Key1 -> 47, Key2 -> 1, Key3 -> 1, MaxPass -> 5, + Border -> 0, MaxChisq -> 10, MinDeviation -> .25, + Given -> {}, NExtra -> 0, PeakFinder -> ({}&), + Verbose -> 1, Final -> All, + PseudoRandom -> False, PseudoRandomSeed -> 5489, + Regions -> False, Compiled -> True} + +:Evaluate: Divonne[f_, v:{_, _, _}.., opt___Rule] := + Block[ {ff = HoldForm[f], ndim = Length[{v}], ncomp, + tags, vars, lower, range, jac, tmp, defs, intT, intX, + rel, abs, mineval, maxeval, key1, key2, key3, maxpass, border, + maxchisq, mindeviation, given, nextra, peakfinder, + final, verbose, level, seed, regions, compiled, + $Phase}, + Message[Divonne::optx, #, Divonne]&/@ + Complement[First/@ {opt}, tags = First/@ Options[Divonne]]; + {rel, abs, mineval, maxeval, key1, key2, key3, maxpass, border, + maxchisq, mindeviation, given, nextra, peakfinder, + verbose, final, level, seed, regions, compiled} = + tags /. {opt} /. Options[Divonne]; + {vars, lower, range} = Transpose[{v}]; + jac = Simplify[Times@@ (range -= lower)]; + tmp = Array[tmpvar, ndim]; + defs = Simplify[lower + range tmp]; + Block[{Set}, define[compiled, tmp, vars, Thread[vars = defs], jac]]; + intT = integrandT[f]; + intX = integrandX[f]; + Block[#, + ncomp = Length[intT@@ RandomReal[1, ndim]]; + MLDivonne[ndim, ncomp, 10.^-rel, 10.^-abs, + Min[Max[verbose, 0], 3] + + If[final === Last, 4, 0] + + If[TrueQ[regions], 128, 0] + + If[IntegerQ[level], 256 level, 0], + If[level =!= False && IntegerQ[seed], seed, 0], + mineval, maxeval, + key1, key2, key3, maxpass, + N[border], N[maxchisq], N[mindeviation], + given, sample[given, 0, intX], nextra] + ]& @ vars + ] + +:Evaluate: tmpvar[n_] := ToExpression["Cuba`Divonne`t" <> ToString[n]] + +:Evaluate: Attributes[foo] = {HoldAll} + +:Evaluate: define[True, tmp_, vars_, defs_, jac_] := ( + TtoX := TtoX = Compile[tmp, defs]; + integrandT[f_] := Compile[tmp, eval[defs, Chop[f jac]//N], + {{_eval, _Real, 1}}]; + integrandX[f_] := Compile[vars, eval[vars, Chop[f jac]//N], + {{_eval, _Real, 1}}] ) + +:Evaluate: define[_, tmp_, vars_, defs_, jac_] := ( + TtoX := TtoX = Function[tmp, defs]; + integrandT[f_] := Function[tmp, eval[defs, Chop[f jac]//N]]; + integrandX[f_] := Function[vars, eval[vars, Chop[f jac]//N]] ) + +:Evaluate: eval[_, f_Real] := {f} + +:Evaluate: eval[_, f:{__Real}] := f + +:Evaluate: eval[x_, _] := (Message[Divonne::badsample, ff, x]; {}) + +:Evaluate: sample[x_, p_, i_:intT] := ( + $Phase = p; + Check[Flatten @ MapSample[i@@ # &, Partition[x, ndim]], {}] ) + +:Evaluate: MapSample = Map + +:Evaluate: findpeak[b_, p_] := Check[Join[#, sample[#, p, intX]]& @ + N[Flatten[peakfinder@@ MapThread[TtoX, Partition[b, 2]]]], {}] + +:Evaluate: region[ll_, ur_, r___] := Region[TtoX@@ ll, TtoX@@ ur, r] + +:Evaluate: Divonne::badsample = "`` is not a real-valued function at ``." + +:Evaluate: Divonne::baddim = "Cannot integrate in `` dimensions." + +:Evaluate: Divonne::badcomp = "Cannot integrate `` components." + +:Evaluate: Divonne::accuracy = + "Desired accuracy was not reached within `` integrand evaluations on `` subregions. + Estimate that MaxPoints needs to be increased by `` for this accuracy." + +:Evaluate: Divonne::success = "Needed `` integrand evaluations on `` subregions." + +:Evaluate: End[] + +:Evaluate: EndPackage[] + + +/* + Divonne.tm + Multidimensional integration by partitioning + originally by J.H. Friedman and M.H. Wright + (CERNLIB subroutine D151) + this version by Thomas Hahn + last modified 12 Oct 11 th +*/ + + +#include "mathlink.h" +#include "decl.h" + +#define ExploreParent Explore + +/*********************************************************************/ + +static void Status(MLCONST char *msg, cint n1, cint n2, cint n3) +{ + MLPutFunction(stdlink, "CompoundExpression", 2); + MLPutFunction(stdlink, "Message", 4); + MLPutFunction(stdlink, "MessageName", 2); + MLPutSymbol(stdlink, "Divonne"); + MLPutString(stdlink, msg); + MLPutInteger(stdlink, n1); + MLPutInteger(stdlink, n2); + MLPutInteger(stdlink, n3); +} + +/*********************************************************************/ + +static void Print(MLCONST char *s) +{ + MLPutFunction(stdlink, "EvaluatePacket", 1); + MLPutFunction(stdlink, "Print", 1); + MLPutString(stdlink, s); + MLEndPacket(stdlink); + + MLNextPacket(stdlink); + MLNewPacket(stdlink); +} + +/*********************************************************************/ + +static void DoSample(This *t, cnumber n, real *x, real *f, ccount ldx) +{ + real *mma_f; + long mma_n; + + if( MLAbort ) longjmp(t->abort, -99); + + MLPutFunction(stdlink, "EvaluatePacket", 1); + MLPutFunction(stdlink, "Cuba`Divonne`sample", 2); + MLPutRealList(stdlink, x, n*t->ndim); + MLPutInteger(stdlink, t->phase); + MLEndPacket(stdlink); + + MLNextPacket(stdlink); + if( !MLGetRealList(stdlink, &mma_f, &mma_n) ) { + MLClearError(stdlink); + MLNewPacket(stdlink); + longjmp(t->abort, -99); + } + + if( mma_n != n*t->ncomp ) { + MLDisownRealList(stdlink, mma_f, mma_n); + longjmp(t->abort, -3); + } + + t->neval += n; + + Copy(f, mma_f, n*t->ncomp); + MLDisownRealList(stdlink, mma_f, mma_n); +} + +/*********************************************************************/ + +static count SampleExtra(This *t, cBounds *b) +{ + count n, nget; + real *mma_f; + long mma_n; + + MLPutFunction(stdlink, "EvaluatePacket", 1); + MLPutFunction(stdlink, "Cuba`Divonne`findpeak", 2); + MLPutRealList(stdlink, (real *)b, 2*t->ndim); + MLPutInteger(stdlink, t->phase); + MLEndPacket(stdlink); + + MLNextPacket(stdlink); + if( !MLGetRealList(stdlink, &mma_f, &mma_n) ) { + MLClearError(stdlink); + MLNewPacket(stdlink); + longjmp(t->abort, -99); + } + + t->neval += nget = mma_n/(t->ndim + t->ncomp); + + n = IMin(nget, t->nextra); + if( n ) { + Copy(t->xextra, mma_f, n*t->ndim); + Copy(t->fextra, mma_f + nget*t->ndim, n*t->ncomp); + } + + MLDisownRealList(stdlink, mma_f, mma_n); + + return n; +} + +/*********************************************************************/ + +#include "common.c" + +static inline void DoIntegrate(This *t) +{ + real integral[NCOMP], error[NCOMP], prob[NCOMP]; + cint fail = Integrate(t, integral, error, prob); + + if( fail < 0 ) { + switch( fail ) { + case -99: + MLPutFunction(stdlink, "Abort", 0); + return; + case -1: + Status("baddim", t->ndim, 0, 0); + break; + case -2: + Status("badcomp", t->ncomp, 0, 0); + break; + } + MLPutSymbol(stdlink, "$Failed"); + } + else { + Status(fail ? "accuracy" : "success", t->neval, t->nregions, fail); + MLPutFunction(stdlink, "Thread", 1); + MLPutFunction(stdlink, "List", 3); + MLPutRealList(stdlink, integral, t->ncomp); + MLPutRealList(stdlink, error, t->ncomp); + MLPutRealList(stdlink, prob, t->ncomp); + } +} + +/*********************************************************************/ + +void Divonne(cint ndim, cint ncomp, + creal epsrel, creal epsabs, + cint flags, cint seed, + cnumber mineval, cnumber maxeval, + cint key1, cint key2, cint key3, cint maxpass, + creal border, creal maxchisq, creal mindeviation, + real *xgiven, clong nxgiven, real *fgiven, clong nfgiven, + cnumber nextra) +{ + This t; + t.ldxgiven = t.ndim = ndim; + t.ncomp = ncomp; + t.epsrel = epsrel; + t.epsabs = epsabs; + t.flags = flags; + t.seed = seed; + t.mineval = mineval; + t.maxeval = maxeval; + t.key1 = key1; + t.key2 = key2; + t.key3 = key3; + t.maxpass = maxpass; + t.border.upper = 1 - (t.border.lower = border); + t.maxchisq = maxchisq; + t.mindeviation = mindeviation; + t.xgiven = NULL; + t.nextra = nextra; + t.nregions = 0; + t.neval = t.ngiven = nxgiven/ndim; + + if( t.ngiven | t.nextra ) { + cnumber nx = nxgiven + nextra*t.ndim; + cnumber nf = nfgiven + nextra*t.ncomp; + + Alloc(t.xgiven, nx + nf); + t.xextra = t.xgiven + nxgiven; + t.fgiven = t.xgiven + nx; + t.fextra = t.fgiven + nfgiven; + + Copy(t.xgiven, xgiven, nxgiven); + Copy(t.fgiven, fgiven, nfgiven); + } + + DoIntegrate(&t); + + free(t.xgiven); + MLEndPacket(stdlink); +} + +/*********************************************************************/ + +int main(int argc, char **argv) +{ + return MLMain(argc, argv); +} + diff --git a/src/divonne/Explore.c b/src/divonne/Explore.c new file mode 100644 index 0000000..d6fba19 --- /dev/null +++ b/src/divonne/Explore.c @@ -0,0 +1,160 @@ +/* + Explore.c + sample region, determine min and max, split if necessary + this file is part of Divonne + last modified 15 Nov 11 th +*/ + + +typedef struct { + real fmin, fmax; + creal *xmin, *xmax; +} Extrema; + +/*********************************************************************/ + +static int Explore(This *t, ccount iregion) +{ + TYPEDEFREGION; + Region *region = RegionPtr(iregion); + cBounds *bounds = region->bounds; + Result *result = region->result; + + count n, dim, comp, maxcomp; + Extrema extrema[NCOMP]; + Result *r; + creal *x; + real *f; + real halfvol, maxerr; + cSamples *samples = &t->samples[region->isamples]; + + /* needed as of gcc 3.3 to make gcc correctly address region #@$&! */ + sizeof(*region); + + for( comp = 0; comp < t->ncomp; ++comp ) { + Extrema *e = &extrema[comp]; + e->fmin = INFTY; + e->fmax = -INFTY; + e->xmin = e->xmax = NULL; + } + + if( region->isamples == 0 ) { /* others already sampled */ + real vol = 1; + for( dim = 0; dim < t->ndim; ++dim ) { + cBounds *b = &bounds[dim]; + vol *= b->upper - b->lower; + } + region->vol = vol; + + for( comp = 0; comp < t->ncomp; ++comp ) { + Result *r = &result[comp]; + r->fmin = INFTY; + r->fmax = -INFTY; + } + + x = t->xgiven; + f = t->fgiven; + n = t->ngiven; + if( t->nextra ) n += SampleExtra(t, bounds); + + for( ; n; --n ) { + for( dim = 0; dim < t->ndim; ++dim ) { + cBounds *b = &bounds[dim]; + if( x[dim] < b->lower || x[dim] > b->upper ) goto skip; + } + for( comp = 0; comp < t->ncomp; ++comp ) { + Extrema *e = &extrema[comp]; + creal y = f[comp]; + if( y < e->fmin ) e->fmin = y, e->xmin = x; + if( y > e->fmax ) e->fmax = y, e->xmax = x; + } +skip: + x += t->ldxgiven; + f += t->ncomp; + } + + samples->sampler(t, iregion); + } + + x = samples->x; + f = samples->f; + for( n = samples->n; n; --n ) { + for( comp = 0; comp < t->ncomp; ++comp ) { + Extrema *e = &extrema[comp]; + creal y = *f++; + if( y < e->fmin ) e->fmin = y, e->xmin = x; + if( y > e->fmax ) e->fmax = y, e->xmax = x; + } + x += t->ndim; + } + t->neval_opt -= t->neval; + + halfvol = .5*region->vol; + maxerr = -INFTY; + maxcomp = -1; + + for( comp = 0; comp < t->ncomp; ++comp ) { + Extrema *e = &extrema[comp]; + Result *r = &result[comp]; + real xtmp[NDIM], ftmp, err; + + if( e->xmin ) { /* not all NaNs */ + t->selectedcomp = comp; + VecCopy(xtmp, e->xmin); + ftmp = FindMinimum(t, bounds, xtmp, e->fmin); + if( ftmp < r->fmin ) { + r->fmin = ftmp; + VecCopy(r->xmin, xtmp); + } + + t->selectedcomp = Tag(comp); + VecCopy(xtmp, e->xmax); + ftmp = -FindMinimum(t, bounds, xtmp, -e->fmax); + if( ftmp > r->fmax ) { + r->fmax = ftmp; + VecCopy(r->xmax, xtmp); + } + } + + r->spread = halfvol*(r->fmax - r->fmin); + err = r->spread/Max(fabs(r->avg), NOTZERO); + if( err > maxerr ) { + maxerr = err; + maxcomp = comp; + } + } + + t->neval_opt += t->neval; + + if( maxcomp == -1 ) { /* all NaNs */ + region->depth = 0; + return -1; + } + + region->cutcomp = maxcomp; + r = ®ion->result[maxcomp]; + if( halfvol*(r->fmin + r->fmax) > r->avg ) { + region->fminor = r->fmin; + region->fmajor = r->fmax; + region->xmajor = r->xmax - (real *)region->result; + } + else { + region->fminor = r->fmax; + region->fmajor = r->fmin; + region->xmajor = r->xmin - (real *)region->result; + } + + if( region->isamples == 0 ) { + if( r->spread < samples->neff*r->err || + r->spread < t->totals[maxcomp].secondspread ) region->depth = 0; + if( region->depth == 0 ) + for( comp = 0; comp < t->ncomp; ++comp ) + t->totals[comp].secondspread = + Max(t->totals[comp].secondspread, result[comp].spread); + } + + if( region->depth ) Split(t, iregion); + + return iregion; +} + diff --git a/src/divonne/FindMinimum.c b/src/divonne/FindMinimum.c new file mode 100644 index 0000000..b476520 --- /dev/null +++ b/src/divonne/FindMinimum.c @@ -0,0 +1,689 @@ +/* + FindMinimum.c + find minimum (maximum) of hyperrectangular region + this file is part of Divonne + last modified 8 Jun 10 th +*/ + + +#define EPS 0x1p-52 +#define RTEPS 0x1p-26 +#define QEPS 0x1p-13 + +#define DELTA 0x1p-16 +#define RTDELTA 0x1p-8 +#define QDELTA 0x1p-4 + +/* +#define DELTA 1e-5 +#define RTDELTA 3.1622776601683791e-3 +#define QDELTA 5.6234132519034912e-2 +*/ + +#define SUFTOL 8*QEPS*QDELTA +#define FTOL 5e-2 +#define GTOL 1e-2 + +#define Hessian(i, j) hessian[(i)*t->ndim + j] + +typedef struct { real dx, f; } Point; + +/*********************************************************************/ + +static inline real Dot(ccount n, creal *a, creal *b) +{ + real sum = 0; + count i; + for( i = 0; i < n; ++i ) sum += a[i]*b[i]; + return sum; +} + +/*********************************************************************/ + +static inline real Length(ccount n, creal *vec) +{ + return sqrt(Dot(n, vec, vec)); +} + +/*********************************************************************/ + +static inline void LinearSolve(cThis *t, ccount n, creal *hessian, + creal *grad, real *p) +{ + int i, j; + real dir; + + for( i = 0; i < n; ++i ) { + dir = -grad[i]; + for( j = 0; j < i; ++j ) + dir -= Hessian(i, j)*p[j]; + p[i] = dir; + } + + while( --i >= 0 ) { + if( Hessian(i, i) <= 0 ) return; + dir = p[i]/Hessian(i, i); + for( j = i + 1; j < n; ++j ) + dir -= Hessian(j, i)*p[j]; + p[i] = dir; + } +} + +/*********************************************************************/ + +static void RenormalizeCholesky(cThis *t, ccount n, real *hessian, + real *z, real alpha) +{ + count i, j; + + for( i = 0; i < n; ++i ) { + creal dir = z[i]; + real beta = alpha*dir; + real gamma = Hessian(i, i); + real gammanew = Hessian(i, i) += beta*dir; + + if( i + 1 >= n || gammanew < 0 || + (gammanew < 1 && gamma > DBL_MAX*gammanew) ) return; + + gamma /= gammanew; + beta /= gammanew; + alpha *= gamma; + + if( gamma < .25 ) { + for( j = i + 1; j < n; ++j ) { + real delta = beta*z[j]; + z[j] -= dir*Hessian(j, i); + Hessian(j, i) = Hessian(j, i)*gamma + delta; + } + } + else { + for( j = i + 1; j < n; ++j ) { + z[j] -= dir*Hessian(j, i); + Hessian(j, i) += beta*z[j]; + } + } + } +} + +/*********************************************************************/ + +static void UpdateCholesky(cThis *t, ccount n, real *hessian, + real *z, real *p) +{ + int i, j; + real gamma = 0; + + for( i = 0; i < n; ++i ) { + real dir = z[i]; + for( j = 0; j < i; ++j ) + dir -= Hessian(i, j)*p[j]; + p[i] = dir; + gamma += Sq(dir)/Hessian(i, i); + } + gamma = Max(fabs(1 - gamma), EPS); + + while( --i >= 0 ) { + creal dir = z[i] = p[i]; + real beta = dir/Hessian(i, i); + creal gammanew = gamma + dir*beta; + Hessian(i, i) *= gamma/gammanew; + beta /= gamma; + gamma = gammanew; + for( j = i + 1; j < n; ++j ) { + creal delta = beta*z[j]; + z[j] += dir*Hessian(j, i); + Hessian(j, i) -= delta; + } + } +} + +/*********************************************************************/ + +static inline void BFGS(cThis *t, ccount n, real *hessian, + creal *gnew, creal *g, real *p, creal dx) +{ + real y[NDIM], c; + count i, j; + + for( i = 0; i < n; ++i ) + y[i] = gnew[i] - g[i]; + c = dx*Dot(n, y, p); + if( c < 1e-10 ) return; + RenormalizeCholesky(t, n, hessian, y, 1/c); + + c = Dot(n, g, p); + if( c >= 0 ) return; + c = 1/sqrt(-c); + for( i = 0; i < n; ++i ) + y[i] = c*g[i]; + UpdateCholesky(t, n, hessian, y, p); + + for( i = 0; i < n - 1; ++i ) + for( j = i + 1; j < n; ++j ) + Hessian(i, j) = Hessian(j, i); +} + +/*********************************************************************/ + +static void Gradient(This *t, ccount nfree, ccount *ifree, + cBounds *b, real *x, creal y, real *grad) +{ + count i; + + for( i = 0; i < nfree; ++i ) { + ccount dim = Untag(ifree[i]); + creal xd = x[dim]; + creal delta = (b[dim].upper - xd < DELTA) ? -DELTA : DELTA; + x[dim] += delta; + grad[i] = (Sample(t, x) - y)/delta; + x[dim] = xd; + } +} + +/*********************************************************************/ + +static Point LineSearch(This *t, ccount nfree, ccount *ifree, + creal *p, creal *xini, real fini, real *x, + real step, creal range, creal grad, + creal ftol, creal xtol, creal gtol) +{ + real tol = ftol, tol2 = tol + tol; + Point cur = {0, fini}; + + VecCopy(x, xini); + + /* don't even try if + a) we'd walk backwards, + b) the range to explore is too small, + c) the gradient is positive, i.e. we'd move uphill */ + + if( step > 0 && range > tol2 && grad <= 0 ) { + creal eps = RTEPS*fabs(range) + ftol; + creal mingrad = -1e-4*grad, maxgrad = -gtol*grad; + + real end = range + eps; + real maxstep = range - eps/(1 + RTEPS); + + Point min = cur, v = cur, w = cur; + Point a = cur, b = {end, 0}; + real a1, b1 = end; + + /* distmin: distance along p from xini to the minimum, + u: second-lowest point, + v: third-lowest point, + a, b: interval in which the minimum is sought. */ + + real distmin = 0, dist, mid, q, r, s; + count i; + int shift; + bool first; + + for( first = true; ; first = false ) { + if( step >= maxstep ) { + step = maxstep; + maxstep = maxstep*(1 + .75*RTEPS) + .75*tol; + } + + cur.dx = (fabs(step) >= tol) ? step : (step > 0) ? tol : -tol; + dist = distmin + cur.dx; + for( i = 0; i < nfree; ++i ) { + ccount dim = ifree[i]; + x[dim] = xini[dim] + dist*p[i]; + } + cur.f = Sample(t, x); + + if( cur.f <= min.f ) { + v = w; + w = min; + min.f = cur.f; + distmin = dist; + + /* shift everything to the new minimum position */ + maxstep -= cur.dx; + v.dx -= cur.dx; + w.dx -= cur.dx; + a.dx -= cur.dx; + b.dx -= cur.dx; + if( cur.dx < 0 ) b = w; + else a = w; + + tol = RTEPS*fabs(distmin) + ftol; + tol2 = tol + tol; + } + else { + if( cur.dx < 0 ) a = cur; + else b = cur; + if( cur.f <= w.f || w.dx == 0 ) v = w, w = cur; + else if( cur.f <= v.f || v.dx == 0 || v.dx == w.dx ) v = cur; + } + + if( distmin + b.dx <= xtol ) break; + if( min.f < fini && + a.f - min.f <= fabs(a.dx)*maxgrad && + (fabs(distmin - range) > tol || maxstep < b.dx) ) break; + + mid = .5*(a.dx + b.dx); + if( fabs(mid) <= tol2 - .5*(b.dx - a.dx) ) break; + + r = q = s = 0; + if( fabs(end) > tol ) { + if( first ) { + creal s1 = w.dx*grad; + creal s2 = w.f - min.f; + s = (s1 - ((distmin == 0) ? 0 : 2*s2))*w.dx; + q = 2*(s2 - s1); + } + else { + creal s1 = w.dx*(v.f - min.f); + creal s2 = v.dx*(w.f - min.f); + s = s1*w.dx - s2*v.dx; + q = 2*(s2 - s1); + } + if( q > 0 ) s = -s; + q = fabs(q); + r = end; + if( step != b1 || b.dx <= maxstep ) end = step; + } + + if( distmin == a.dx ) step = mid; + else if( b.dx > maxstep ) step = (step < b.dx) ? -4*a.dx : maxstep; + else { + real num = a.dx, den = b.dx; + if( fabs(b.dx) <= tol || (w.dx > 0 && fabs(a.dx) > tol) ) + num = b.dx, den = a.dx; + num /= -den; + step = (num < 1) ? .5*den*sqrt(num) : 5/11.*den*(.1 + 1/num); + } + + if( step > 0 ) a1 = a.dx, b1 = step; + else a1 = step, b1 = b.dx; + if( fabs(s) < fabs(.5*q*r) && s > q*a1 && s < q*b1 ) { + step = s/q; + if( step - a.dx < tol2 || b.dx - step < tol2 ) + step = (mid > 0) ? tol : -tol; + } + else end = (mid > 0) ? b.dx : a.dx; + } + + first = true; + if( fabs(distmin - range) < tol ) { + distmin = range; + if( maxstep > b.dx ) first = false; + } + + for( cur.dx = distmin, cur.f = min.f, shift = -1; ; + cur.dx = Max(ldexp(distmin, shift), ftol), shift <<= 1 ) { + for( i = 0; i < nfree; ++i ) { + ccount dim = ifree[i]; + x[dim] = xini[dim] + cur.dx*p[i]; + } + if( !first ) cur.f = Sample(t, x); + + if( cur.dx + b.dx <= xtol ) { + cur.dx = 0; + break; + } + if( fini - cur.f > cur.dx*mingrad ) break; + if( cur.dx <= ftol ) { + cur.dx = 0; + break; + } + first = false; + } + } + + return cur; +} + +/*********************************************************************/ + +static real LocalSearch(This *t, ccount nfree, ccount *ifree, + cBounds *b, creal *x, creal fx, real *z) +{ + real delta, smax, sopp, spmax, snmax; + real y[NDIM], fy, fz, ftest; + real p[NDIM]; + int sign; + count i; + + /* Choose a direction p along which to move away from the + present x. We choose the direction which leads farthest + away from all borders. */ + + smax = INFTY; + for( i = 0; i < nfree; ++i ) { + ccount dim = ifree[i]; + creal sp = b[dim].upper - x[dim]; + creal sn = x[dim] - b[dim].lower; + if( sp < sn ) { + smax = Min(smax, sn); + p[i] = -1; + } + else { + smax = Min(smax, sp); + p[i] = 1; + } + } + smax *= .9; + + /* Move along p until the integrand changes appreciably + or we come close to a border. */ + + VecCopy(y, x); + ftest = SUFTOL*(1 + fabs(fx)); + delta = RTDELTA/5; + do { + delta = Min(5*delta, smax); + for( i = 0; i < nfree; ++i ) { + ccount dim = ifree[i]; + y[dim] = x[dim] + delta*p[i]; + } + fy = Sample(t, y); + if( fabs(fy - fx) > ftest ) break; + } while( delta != smax ); + + /* Construct a second direction p' orthogonal to p, i.e. p.p' = 0. + We let pairs of coordinates cancel in the dot product, + i.e. we choose p'[0] = p[0], p'[1] = -p[1], etc. + (It should really be 1/p and -1/p, but p consists of 1's and -1's.) + For odd nfree, we let the last three components cancel by + choosing p'[nfree - 3] = p[nfree - 3], + p'[nfree - 2] = -1/2 p[nfree - 2], and + p'[nfree - 1] = -1/2 p[nfree - 1]. */ + + sign = (nfree <= 1 && fy > fx) ? 1 : -1; + spmax = snmax = INFTY; + for( i = 0; i < nfree; ++i ) { + ccount dim = ifree[i]; + real sp, sn; + p[i] *= (nfree & 1 && nfree - i <= 2) ? -.5*sign : (sign = -sign); + sp = (b[dim].upper - y[dim])/p[i]; + sn = (y[dim] - b[dim].lower)/p[i]; + if( p[i] > 0 ) { + spmax = Min(spmax, sp); + snmax = Min(snmax, sn); + } + else { + spmax = Min(spmax, -sn); + snmax = Min(snmax, -sp); + } + } + smax = .9*spmax; + sopp = .9*snmax; + + if( nfree > 1 && smax < snmax ) { + real tmp = smax; + smax = sopp; + sopp = tmp; + for( i = 0; i < nfree; ++i ) + p[i] = -p[i]; + } + + /* Move along p' until the integrand changes appreciably + or we come close to a border. */ + + VecCopy(z, y); + ftest = SUFTOL*(1 + fabs(fy)); + delta = RTDELTA/5; + do { + delta = Min(5*delta, smax); + for( i = 0; i < nfree; ++i ) { + ccount dim = ifree[i]; + z[dim] = y[dim] + delta*p[i]; + } + fz = Sample(t, z); + if( fabs(fz - fy) > ftest ) break; + } while( delta != smax ); + + if( fy != fz ) { + real pleneps, grad, range, step; + Point low; + + if( fy > fz ) { + grad = (fz - fy)/delta; + range = smax/.9; + step = Min(delta + delta, smax); + } + else { + grad = (fy - fz)/delta; + range = sopp/.9 + delta; + step = Min(delta + delta, sopp); + VecCopy(y, z); + fy = fz; + for( i = 0; i < nfree; ++i ) + p[i] = -p[i]; + } + + pleneps = Length(nfree, p) + RTEPS; + low = LineSearch(t, nfree, ifree, p, y, fy, z, step, range, grad, + RTEPS/pleneps, 0., RTEPS); + fz = low.f; + } + + if( fz != fx ) { + real pleneps, grad, range, step; + Point low; + + spmax = snmax = INFTY; + for( i = 0; i < nfree; ++i ) { + ccount dim = ifree[i]; + p[i] = z[dim] - x[dim]; + if( p[i] != 0 ) { + creal sp = (b[dim].upper - x[dim])/p[i]; + creal sn = (x[dim] - b[dim].lower)/p[i]; + if( p[i] > 0 ) { + spmax = Min(spmax, sp); + snmax = Min(snmax, sn); + } + else { + spmax = Min(spmax, -sn); + snmax = Min(snmax, -sp); + } + } + } + + grad = fz - fx; + range = spmax; + step = Min(.9*spmax, 2.); + pleneps = Length(nfree, p) + RTEPS; + if( fz > fx ) { + delta = Min(.9*snmax, RTDELTA/pleneps); + for( i = 0; i < nfree; ++i ) { + ccount dim = ifree[i]; + z[dim] = x[dim] - delta*p[i]; + } + fz = Sample(t, z); + if( fz < fx ) { + grad = (fz - fx)/delta; + range = snmax; + step = Min(.9*snmax, delta + delta); + for( i = 0; i < nfree; ++i ) + p[i] = -p[i]; + } + else if( delta < 1 ) grad = (fx - fz)/delta; + } + + low = LineSearch(t, nfree, ifree, p, x, fx, z, step, range, grad, + RTEPS/pleneps, 0., RTEPS); + fz = low.f; + } + + return fz; +} + +/*********************************************************************/ + +static real FindMinimum(This *t, cBounds *b, real *xmin, real fmin) +{ + real hessian[NDIM*NDIM]; + real gfree[NDIM], p[NDIM]; + real tmp[NDIM], ftmp, fini = fmin; + ccount maxeval = t->neval + 50*t->ndim; + count nfree, nfix; + count ifree[NDIM], ifix[NDIM]; + count dim, local; + + Zap(hessian); + for( dim = 0; dim < t->ndim; ++dim ) + Hessian(dim, dim) = 1; + + /* Step 1: - classify the variables as "fixed" (sufficiently close + to a border) and "free", + - if the integrand is flat in the direction of the gradient + w.r.t. the free dimensions, perform a local search. */ + + for( local = 0; local < 2; ++local ) { + bool resample = false; + nfree = nfix = 0; + for( dim = 0; dim < t->ndim; ++dim ) { + if( xmin[dim] < b[dim].lower + (1 + fabs(b[dim].lower))*QEPS ) { + xmin[dim] = b[dim].lower; + ifix[nfix++] = dim; + resample = true; + } + else if( xmin[dim] > b[dim].upper - (1 + fabs(b[dim].upper))*QEPS ) { + xmin[dim] = b[dim].upper; + ifix[nfix++] = Tag(dim); + resample = true; + } + else ifree[nfree++] = dim; + } + + if( resample ) fini = fmin = Sample(t, xmin); + + if( nfree == 0 ) goto releasebounds; + + Gradient(t, nfree, ifree, b, xmin, fmin, gfree); + if( local || Length(nfree, gfree) > GTOL ) break; + + ftmp = LocalSearch(t, nfree, ifree, b, xmin, fmin, tmp); + if( ftmp > fmin - (1 + fabs(fmin))*RTEPS ) + goto releasebounds; + fmin = ftmp; + VecCopy(xmin, tmp); + } + + while( t->neval <= maxeval ) { + + /* Step 2a: perform a quasi-Newton iteration on the free + variables only. */ + + if( nfree > 0 ) { + real plen, pleneps; + real minstep; + count i, mini = 0, minfix = 0; + Point low; + + LinearSolve(t, nfree, hessian, gfree, p); + plen = Length(nfree, p); + pleneps = plen + RTEPS; + + minstep = INFTY; + for( i = 0; i < nfree; ++i ) { + count dim = Untag(ifree[i]); + if( fabs(p[i]) > EPS ) { + real step; + count fix; + if( p[i] < 0 ) { + step = (b[dim].lower - xmin[dim])/p[i]; + fix = dim; + } + else { + step = (b[dim].upper - xmin[dim])/p[i]; + fix = Tag(dim); + } + if( step < minstep ) { + minstep = step; + mini = i; + minfix = fix; + } + } + } + + if( minstep*pleneps <= DELTA ) { +fixbound: + ifix[nfix++] = minfix; + + if( mini < --nfree ) { + creal diag = Hessian(mini, mini); + + Clear(tmp, mini); + for( i = mini; i < nfree; ++i ) + tmp[i] = Hessian(i + 1, mini); + + for( i = mini; i < nfree; ++i ) { + Copy(&Hessian(i, 0), &Hessian(i + 1, 0), i); + Hessian(i, i) = Hessian(i + 1, i + 1); + } + RenormalizeCholesky(t, nfree, hessian, tmp, diag); + + Copy(&ifree[mini], &ifree[mini + 1], nfree - mini); + Copy(&gfree[mini], &gfree[mini + 1], nfree - mini); + } + continue; + } + + low = LineSearch(t, nfree, ifree, p, xmin, fmin, tmp, + Min(minstep, 1.), Min(minstep, 100.), Dot(nfree, gfree, p), + RTEPS/pleneps, DELTA/pleneps, .2); + + if( low.dx > 0 ) { + real fdiff; + + fmin = low.f; + VecCopy(xmin, tmp); + + Gradient(t, nfree, ifree, b, xmin, fmin, tmp); + BFGS(t, nfree, hessian, tmp, gfree, p, low.dx); + VecCopy(gfree, tmp); + + if( fabs(low.dx - minstep) < QEPS*minstep ) goto fixbound; + + fdiff = fini - fmin; + fini = fmin; + if( fdiff > (1 + fabs(fmin))*FTOL || + low.dx*plen > (1 + Length(t->ndim, xmin))*FTOL ) continue; + } + } + + /* Step 2b: check whether freeing any fixed variable will lead + to a reduction in f. */ + +releasebounds: + if( nfix > 0 ) { + real mingrad = INFTY; + count i, mini = 0; + bool repeat = false; + + Gradient(t, nfix, ifix, b, xmin, fmin, tmp); + + for( i = 0; i < nfix; ++i ) { + creal grad = Sign(ifix[i])*tmp[i]; + if( grad < -RTEPS ) { + repeat = true; + if( grad < mingrad ) { + mingrad = grad; + mini = i; + } + } + } + + if( repeat ) { + gfree[nfree] = tmp[mini]; + ifree[nfree] = Untag(ifix[mini]); + Clear(&Hessian(nfree, 0), nfree); + Hessian(nfree, nfree) = 1; + ++nfree; + + --nfix; + Copy(&ifix[mini], &ifix[mini + 1], nfix - mini); + continue; + } + } + + break; + } + + return fmin; +} + diff --git a/src/divonne/Integrate.c b/src/divonne/Integrate.c new file mode 100644 index 0000000..4059a55 --- /dev/null +++ b/src/divonne/Integrate.c @@ -0,0 +1,452 @@ +/* + Integrate.c + partition the integration region until each region + has approximately equal spread = 1/2 vol (max - min), + then do a main integration over all regions + this file is part of Divonne + last modified 15 Nov 11 th +*/ + + +#define INIDEPTH 3 +#define DEPTH 5 +#define POSTDEPTH 15 + +/*********************************************************************/ + +static int Integrate(This *t, real *integral, real *error, real *prob) +{ + TYPEDEFREGION; + + Totals totals[NCOMP]; + real nneed; + count dim, comp, iter, pass = 0, err, iregion; + number nwant, nmin = INT_MAX, neff; + int fail; + + if( VERBOSE > 1 ) { + char s[512]; + sprintf(s, "Divonne input parameters:\n" + " ndim " COUNT "\n ncomp " COUNT "\n" + " epsrel " REAL "\n epsabs " REAL "\n" + " flags %d\n seed %d\n" + " mineval " NUMBER "\n maxeval " NUMBER "\n" + " key1 %d\n key2 %d\n key3 %d\n maxpass " COUNT "\n" + " border " REAL "\n maxchisq " REAL "\n mindeviation " REAL "\n" + " ngiven " NUMBER "\n nextra " NUMBER, + t->ndim, t->ncomp, + t->epsrel, t->epsabs, + t->flags, t->seed, + t->mineval, t->maxeval, + t->key1, t->key2, t->key3, t->maxpass, + t->border.lower, t->maxchisq, t->mindeviation, + t->ngiven, t->nextra); + Print(s); + } + + if( BadComponent(t) ) return -2; + if( BadDimension(t, t->key1) || + BadDimension(t, t->key2) || + ((t->key3 & -2) && BadDimension(t, t->key3)) ) return -1; + + t->neval_opt = t->neval_cut = 0; + + AllocRegions(t); + for( dim = 0; dim < t->ndim; ++dim ) { + Bounds *b = &RegionPtr(0)->bounds[dim]; + b->lower = 0; + b->upper = 1; + } + t->nregions = 1; + + RuleIni(&t->rule7); + RuleIni(&t->rule9); + RuleIni(&t->rule11); + RuleIni(&t->rule13); + SamplesIni(&t->samples[0]); + SamplesIni(&t->samples[1]); + SamplesIni(&t->samples[2]); + + if( (fail = setjmp(t->abort)) ) goto abort; + + t->epsabs = Max(t->epsabs, NOTZERO); + + /* Step 1: partition the integration region */ + + if( VERBOSE ) Print("Partitioning phase:"); + + if( IsSobol(t->key1) || IsSobol(t->key2) || IsSobol(t->key3) ) + IniRandom(t); + + SamplesLookup(t, &t->samples[0], t->key1, + (number)47, (number)INT_MAX, (number)0); + SamplesAlloc(t, &t->samples[0]); + + t->totals = totals; + Zap(totals); + t->phase = 1; + + Iterate(t, 0, INIDEPTH, 0, NULL); + + for( iter = 1; ; ++iter ) { + Totals *maxtot; + count valid; + + for( comp = 0; comp < t->ncomp; ++comp ) { + Totals *tot = &totals[comp]; + tot->avg = tot->spreadsq = 0; + tot->spread = tot->secondspread = -INFTY; + } + + for( iregion = 0; iregion < t->nregions; ++iregion ) { + Region *region = RegionPtr(iregion); + for( comp = 0; comp < t->ncomp; ++comp ) { + cResult *r = ®ion->result[comp]; + Totals *tot = &totals[comp]; + tot->avg += r->avg; + tot->spreadsq += Sq(r->spread); + if( r->spread > tot->spread ) { + tot->secondspread = tot->spread; + tot->spread = r->spread; + tot->iregion = iregion; + } + else if( r->spread > tot->secondspread ) + tot->secondspread = r->spread; + } + } + + maxtot = totals; + valid = 0; + for( comp = 0; comp < t->ncomp; ++comp ) { + Totals *tot = &totals[comp]; + integral[comp] = tot->avg; + valid += tot->avg == tot->avg; + if( tot->spreadsq > maxtot->spreadsq ) maxtot = tot; + tot->spread = sqrt(tot->spreadsq); + error[comp] = tot->spread/t->samples[0].neff; + } + + if( VERBOSE ) { + char s[128 + 64*NCOMP], *p = s; + + p += sprintf(p, "\n" + "Iteration " COUNT " (pass " COUNT "): " COUNT " regions\n" + NUMBER7 " integrand evaluations so far,\n" + NUMBER7 " in optimizing regions,\n" + NUMBER7 " in finding cuts", + iter, pass, t->nregions, t->neval, t->neval_opt, t->neval_cut); + + for( comp = 0; comp < t->ncomp; ++comp ) + p += sprintf(p, "\n[" COUNT "] " + REAL " +- " REAL, + comp + 1, integral[comp], error[comp]); + + Print(s); + } + + if( valid == 0 ) goto abort; /* all NaNs */ + + if( t->neval > t->maxeval ) break; + + nneed = maxtot->spread/MaxErr(maxtot->avg); + if( nneed < MAXPRIME ) { + cnumber n = t->neval + t->nregions*(number)ceil(nneed); + if( n < nmin ) { + nmin = n; + pass = 0; + } + else if( ++pass > t->maxpass && n >= t->mineval ) break; + } + + Iterate(t, maxtot->iregion, DEPTH, -1, NULL); + } + + /* Step 2: do a "full" integration on each region */ + +/* nneed = t->samples[0].neff + 1; */ + nneed = 2*t->samples[0].neff; + for( comp = 0; comp < t->ncomp; ++comp ) { + Totals *tot = &totals[comp]; + creal maxerr = MaxErr(tot->avg); + tot->nneed = tot->spread/maxerr; + nneed = Max(nneed, tot->nneed); + tot->maxerrsq = Sq(maxerr); + tot->mindevsq = tot->maxerrsq*Sq(t->mindeviation); + } + nwant = (number)Min(ceil(nneed), MARKMASK/40.); + + err = SamplesLookup(t, &t->samples[1], t->key2, nwant, + (t->maxeval - t->neval)/t->nregions + 1, t->samples[0].n + 1); + + /* the number of points needed to reach the desired accuracy */ + fail = Unmark(err)*t->nregions; + + if( Marked(err) ) { + if( VERBOSE ) Print("\nNot enough samples left for main integration."); + for( comp = 0; comp < t->ncomp; ++comp ) + prob[comp] = -999; + neff = t->samples[0].neff; + } + else { + bool can_adjust = (t->key3 == 1 && t->samples[1].sampler != SampleRule && + (t->key2 < 0 || t->samples[1].neff < MAXPRIME)); + count df, nlimit; + + SamplesAlloc(t, &t->samples[1]); + + if( VERBOSE ) { + char s[128]; + sprintf(s, "\nMain integration on " COUNT + " regions with " NUMBER " samples per region.", + t->nregions, t->samples[1].neff); + Print(s); + } + + ResClear(integral); + ResClear(error); + ResClear(prob); + + nlimit = t->maxeval - t->nregions*t->samples[1].n; + df = 0; + +#define CopyPhaseResults(f) \ + for( comp = 0; comp < t->ncomp; ++comp ) { \ + PhaseResult *p = &totals[comp].phase[f]; \ + cResult *r = ®ion->result[comp]; \ + p->avg = r->avg; \ + p->err = r->err; \ + } + +#define Var2(f, res) Sq((res)->err ? (res)->err : r->spread/t->samples[f].neff) +#define Var(f) Var2(f, &tot->phase[f]) + + for( iregion = 0; iregion < t->nregions; ++iregion ) { + Region *region; + char s[64*NDIM + 256*NCOMP], *p = s; + int todo; + +refine: + region = RegionPtr(iregion); + CopyPhaseResults(0); + t->phase = 2; + region->isamples = 1; + t->samples[1].sampler(t, iregion); + CopyPhaseResults(1); + + if( can_adjust ) + for( comp = 0; comp < t->ncomp; ++comp ) + totals[comp].spreadsq -= Sq(region->result[comp].spread); + + nlimit += t->samples[1].n; + todo = 0; + + for( comp = 0; comp < t->ncomp; ++comp ) { + cResult *r = ®ion->result[comp]; + Totals *tot = &totals[comp]; + + if( t->neval < nlimit ) { + creal avg2 = tot->phase[1].avg; + creal diffsq = Sq(avg2 - tot->phase[0].avg); + + if( r->err*tot->nneed > r->spread || + diffsq > Max(t->maxchisq*(Var(0) + Var(1)), EPS*Sq(avg2)) ) { + if( t->key3 && diffsq > tot->mindevsq ) { + if( t->key3 == 1 ) { + if( VERBOSE > 2 ) Print("\nSplit"); + t->phase = 1; + Iterate(t, iregion, POSTDEPTH, 1, totals); + + if( can_adjust ) { + cnumber nnew = (tot->spreadsq/Sq(MARKMASK) > tot->maxerrsq) ? + MARKMASK : + (number)ceil(sqrt(tot->spreadsq/tot->maxerrsq)); + if( nnew > nwant + nwant/64 ) { + ccount err = SamplesLookup(t, &t->samples[1], t->key2, nnew, + (t->maxeval - t->neval)/t->nregions + 1, t->samples[1].n); + fail += Unmark(err)*t->nregions; + nwant = nnew; + SamplesFree(&t->samples[1]); + SamplesAlloc(t, &t->samples[1]); + + if( t->key2 > 0 && t->samples[1].neff >= MAXPRIME ) + can_adjust = false; + + if( VERBOSE > 2 ) { + char s[128]; + sprintf(s, "Sampling remaining " COUNT + " regions with " NUMBER " points per region.", + t->nregions, t->samples[1].neff); + Print(s); + } + } + } + goto refine; + } + todo |= 3; + } + todo |= 1; + } + } + } + + if( can_adjust ) { + for( comp = 0; comp < t->ncomp; ++comp ) + totals[comp].maxerrsq -= + Sq(region->result[comp].spread/t->samples[1].neff); + } + + switch( todo ) { + case 1: /* get spread right */ + region->isamples = 1; + Explore(t, iregion); + break; + + case 3: /* sample region again with more points */ + if( SamplesIniQ(&t->samples[2]) ) { + SamplesLookup(t, &t->samples[2], t->key3, + nwant, (number)INT_MAX, (number)0); + SamplesAlloc(t, &t->samples[2]); + } + t->phase = 3; + region->isamples = 2; + t->samples[2].sampler(t, iregion); + Explore(t, iregion); + ++region->depth; /* misused for df here */ + ++df; + } + + if( VERBOSE > 2 ) { + for( dim = 0; dim < t->ndim; ++dim ) { + cBounds *b = ®ion->bounds[dim]; + p += sprintf(p, + (dim == 0) ? "\nRegion (" REALF ") - (" REALF ")" : + "\n (" REALF ") - (" REALF ")", + b->lower, b->upper); + } + } + + for( comp = 0; comp < t->ncomp; ++comp ) { + Result *r = ®ion->result[comp]; + Totals *tot = &totals[comp]; + + creal x1 = tot->phase[0].avg; + creal v1 = Var(0); + creal x2 = tot->phase[1].avg; + creal v2 = Var(1); + creal r2 = v1 ? v2/v1 : + Sq(t->samples[1].neff/(real)t->samples[0].neff); + + real norm = 1 + r2; + real avg = x2 + r2*x1; + real sigsq = v2; + real chisq = Sq(x2 - x1); + real chiden = v1 + v2; + + if( todo == 3 ) { + creal x3 = r->avg; + creal v3 = Var2(2, r); + creal r3 = v2 ? v3/v2 : + Sq(t->samples[2].neff/(real)t->samples[1].neff); + + norm = 1 + r3*norm; + avg = x3 + r3*avg; + sigsq = v3; + chisq = v1*Sq(x3 - x2) + v2*Sq(x3 - x1) + v3*chisq; + chiden = v1*v2 + v3*chiden; + } + + avg = LAST ? r->avg : (sigsq *= norm = 1/norm, avg*norm); + if( chisq > EPS ) chisq /= Max(chiden, NOTZERO); + + if( VERBOSE > 2 ) { +#define Out2(f, res) (res)->avg, r->spread/t->samples[f].neff, (res)->err +#define Out(f) Out2(f, &tot->phase[f]) + p += sprintf(p, "\n[" COUNT "] " + REAL " +- " REAL "(" REAL ")\n " + REAL " +- " REAL "(" REAL ")", comp + 1, Out(0), Out(1)); + if( todo == 3 ) p += sprintf(p, "\n " + REAL " +- " REAL "(" REAL ")", Out2(2, r)); + p += sprintf(p, " \tchisq " REAL, chisq); + } + + integral[comp] += avg; + error[comp] += sigsq; + prob[comp] += chisq; + + r->avg = avg; + r->spread = sqrt(sigsq); + r->chisq = chisq; + } + + if( VERBOSE > 2 ) Print(s); + } + + for( comp = 0; comp < t->ncomp; ++comp ) + error[comp] = sqrt(error[comp]); + + df += t->nregions; + + if( VERBOSE > 2 ) { + char s[16 + 128*NCOMP], *p = s; + + p += sprintf(p, "\nTotals:"); + + for( comp = 0; comp < t->ncomp; ++comp ) + p += sprintf(p, "\n[" COUNT "] " + REAL " +- " REAL " \tchisq " REAL " (" COUNT " df)", + comp + 1, integral[comp], error[comp], prob[comp], df); + + Print(s); + } + + for( comp = 0; comp < t->ncomp; ++comp ) + prob[comp] = ChiSquare(prob[comp], df); + + neff = 1; + } + +#ifdef MLVERSION + if( REGIONS ) { + MLPutFunction(stdlink, "List", 2); + MLPutFunction(stdlink, "List", t->nregions); + for( iregion = 0; iregion < t->nregions; ++iregion ) { + Region *region = RegionPtr(iregion); + cBounds *b = region->bounds; + real lower[NDIM], upper[NDIM]; + + for( dim = 0; dim < t->ndim; ++dim ) { + lower[dim] = b[dim].lower; + upper[dim] = b[dim].upper; + } + + MLPutFunction(stdlink, "Cuba`Divonne`region", 4); + + MLPutRealList(stdlink, lower, t->ndim); + MLPutRealList(stdlink, upper, t->ndim); + + MLPutFunction(stdlink, "List", t->ncomp); + for( comp = 0; comp < t->ncomp; ++comp ) { + cResult *r = ®ion->result[comp]; + real res[] = {r->avg, r->spread/neff, r->chisq}; + MLPutRealList(stdlink, res, Elements(res)); + } + + MLPutInteger(stdlink, region->depth + 1); /* misused for df */ + } + } +#endif + +abort: + SamplesFree(&t->samples[2]); + SamplesFree(&t->samples[1]); + SamplesFree(&t->samples[0]); + RuleFree(&t->rule13); + RuleFree(&t->rule11); + RuleFree(&t->rule9); + RuleFree(&t->rule7); + + free(t->voidregion); + + return fail; +} + diff --git a/src/divonne/Iterate.c b/src/divonne/Iterate.c new file mode 100644 index 0000000..82c0e34 --- /dev/null +++ b/src/divonne/Iterate.c @@ -0,0 +1,138 @@ +/* + Iterate.c + recursion over regions + this file is part of Divonne + last modified 15 Nov 11 th +*/ + + +static void Iterate(This *t, count iregion, cint depth, cint isamples, + Totals *totals) +{ + TYPEDEFREGION; + Region *parent, *region; + typedef struct { + real avg, err, spread, spreadsq; + } Corr; + Corr corr[NCOMP]; + count ireg, mreg = iregion; + count comp, maxsplit; + int last, idest, isrc; + + region = RegionPtr(iregion); + region->depth = depth; + region->next = -iregion - 1; + if( isamples < 0 ) Split(t, iregion); + else { + region->isamples = isamples; + Explore(t, iregion); + } + + ireg = iregion + RegionPtr(iregion)->next; + + do { + region = RegionPtr(ireg); + if( region->depth > 0 ) { + --region->depth; +more: + ireg = ExploreParent(t, ireg); + if( ireg == -1 ) return; + region = RegionPtr(ireg); + } + if( region->depth < 0 ) mreg = IMax(mreg, ireg); + ireg += region->next; + } while( ireg > 0 ); + +#ifndef MLVERSION + if( t->running ) goto more; +#endif + + maxsplit = 1; + for( ireg = mreg; ireg >= iregion; --ireg ) { + parent = RegionPtr(ireg); + maxsplit -= NegQ(parent->depth); + if( parent->depth < 0 ) { + count xreg; + struct { + count from, to; + } todo[maxsplit], *tdmax = todo, *td; + count nsplit = 0; + real norm; + + memset(corr, 0, sizeof corr); + + tdmax->from = ireg + parent->next; + tdmax->to = tdmax->from - parent->depth; + ++tdmax; + for( td = todo; td < tdmax; ++td ) { + for( xreg = td->from; xreg < td->to; ++xreg ) { + Region *region = RegionPtr(xreg); + if( region->depth < 0 ) { + tdmax->from = xreg + region->next; + tdmax->to = tdmax->from - region->depth; + ++tdmax; + } + else { + ++nsplit; + for( comp = 0; comp < t->ncomp; ++comp ) { + cResult *r = ®ion->result[comp]; + Corr *c = &corr[comp]; + c->avg += r->avg; + c->err += r->err; + c->spread += Sq(r->spread); + } + } + } + } + + norm = 1./nsplit--; + for( comp = 0; comp < t->ncomp; ++comp ) { + Result *p = &parent->result[comp]; + Corr *c = &corr[comp]; + creal diff = fabs(p->avg - c->avg)*norm; + c->avg = diff*norm*nsplit; + c->err = (c->err == 0) ? 1 : 1 + diff/c->err; + c->spread = (c->spread == 0) ? 1 : 1 + diff/sqrt(c->spread); + } + + for( td = todo; td < tdmax; ++td ) + for( xreg = td->from; xreg < td->to; ++xreg ) { + Region *region = RegionPtr(xreg); + if( region->depth >= 0 ) { + cnumber neff = t->samples[region->isamples].neff; + for( comp = 0; comp < t->ncomp; ++comp ) { + Result *r = ®ion->result[comp]; + Corr *c = &corr[comp]; + if( r->err > 0 ) r->err = r->err*c->err + c->avg; + r->spread = r->spread*c->spread + c->avg*neff; + c->spreadsq += Sq(r->spread); + } + } + } + } + } + + if( totals ) + for( comp = 0; comp < t->ncomp; ++comp ) + totals[comp].spreadsq += corr[comp].spreadsq; + + for( last = -1, idest = isrc = iregion; iregion <= mreg; ++iregion ) { + Region *region = RegionPtr(iregion); + cint cur = NegQ(region->depth); + switch( cur - last ) { + case -1: + memmove(RegionPtr(idest), RegionPtr(isrc), + (iregion - isrc)*sizeof(Region)); + idest += iregion - isrc; + break; + case 1: + isrc = iregion; + } + last = cur; + } + + memmove(RegionPtr(idest), RegionPtr(iregion), + (t->nregions - iregion)*sizeof(Region)); + t->nregions += idest - iregion; +} + diff --git a/src/divonne/KorobovCoeff.c b/src/divonne/KorobovCoeff.c new file mode 100644 index 0000000..e8ace85 --- /dev/null +++ b/src/divonne/KorobovCoeff.c @@ -0,0 +1,881 @@ +#define KOROBOV_MINDIM 2 +#define KOROBOV_MAXDIM 33 +#define MAXPRIME 9689 + +#define Hash(x) ((19945 - x)*(-47 + x))/121634 + +static int prime[] = { + FIRST,47,53,59,67,71,79,83,89,97,103,109,113,127,131,137,139,149,151, + 157,163,173,179,181,191,193,199,211,223,227,229,233,239,241,251,257,263, + 269,277,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383, + 389,397,401,409,419,421,431,433,439,443,449,457,461,467,479,487,491,499, + 503,509,521,523,541,547,557,563,569,571,577,587,593,599,601,607,613,617, + 619,631,641,643,647,653,659,661,673,677,683,691,701,709,719,727,733,739, + 743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853,857,859, + 863,877,881,883,887,907,911,919,929,937,941,947,953,967,971,977,983,991, + 997,1009,1013,1019,1021,1031,1033,1039,1049,1051,1061,1063,1069,1087, + 1091,1093,1097,1103,1109,1117,1123,1129,1151,1153,1163,1171,1181,1187, + 1193,1201,1213,1217,1223,1229,1231,1237,1249,1259,1277,1279,1283,1289, + 1291,1297,1301,1303,1307,1319,1321,1327,1361,1367,1373,1381,1399,1409, + 1423,1427,1429,1433,1439,1447,1451,1453,1459,1471,1481,1483,1487,1489, + 1493,1499,1511,1523,1531,1543,1549,1553,1559,1567,1571,1579,1583,1597, + 1601,1607,1609,1613,1619,1621,1627,1637,1657,1663,1667,1669,1693,1697, + 1699,1709,1721,1723,1733,1741,1747,1753,1759,1777,1783,1787,1789,1801, + 1811,1823,1831,1847,1861,1867,1871,1873,1877,1879,1889,1901,1907,1913, + 1931,1933,1949,1951,1973,1979,1987,1993,1997,1999,2003,2011,2017,2027, + 2029,2039,2053,2063,2069,2081,2083,2087,2089,2099,2111,2113,2129,2131, + 2137,2141,2143,2153,2161,2179,2203,2207,2213,2221,2237,2239,2243,2251, + 2267,2269,2273,2281,2287,2293,2297,2309,2311,2333,2339,2341,2347,2351, + 2357,2371,2377,2381,2383,2389,2393,2399,2411,2417,2423,2437,2441,2447, + 2459,2467,2473,2477,2503,2521,2531,2539,2543,2549,2551,2557,2579,2591, + 2593,2609,2617,2621,2633,2647,2657,2659,2663,2671,2677,2683,2687,2689, + 2693,2699,2707,2711,2713,2719,2729,2731,2741,2749,2753,2767,2777,2789, + 2791,2797,2801,2803,2819,2833,2837,2843,2851,2857,2861,2879,2887,2897, + 2903,2909,2917,2927,2939,2953,2957,2963,2969,2971,2999,3001,3011,3019, + 3023,3037,3041,3049,3061,3067,3079,3083,3089,3109,3119,3121,3137,3163, + 3167,3169,3181,3187,3191,3203,3209,3217,3221,3229,3251,3253,3257,3259, + 3271,3299,3301,3307,3313,3319,3323,3329,3331,3343,3347,3359,3361,3371, + 3373,3389,3391,3407,3413,3433,3449,3457,3461,3463,3467,3469,3491,3499, + 3511,3517,3527,3529,3533,3539,3547,3557,3571,3581,3583,3593,3607,3613, + 3623,3631,3643,3659,3671,3673,3677,3691,3701,3709,3719,3733,3739,3761, + 3767,3769,3779,3793,3797,3803,3821,3833,3847,3851,3853,3863,3877,3889, + 3907,3911,3917,3929,3943,3947,3967,3989,4001,4003,4007,4013,4019,4027, + 4049,4051,4057,4073,4079,4091,4099,4111,4127,4133,4139,4153,4159,4177, + 4201,4211,4217,4219,4229,4231,4243,4259,4271,4283,4289,4297,4327,4337, + 4339,4349,4357,4363,4373,4391,4397,4409,4421,4423,4441,4451,4463,4481, + 4483,4493,4507,4517,4523,4547,4549,4561,4567,4583,4597,4603,4621,4637, + 4639,4651,4663,4673,4691,4703,4721,4723,4733,4751,4759,4783,4787,4789, + 4801,4813,4831,4861,4871,4877,4889,4903,4909,4919,4931,4933,4943,4957, + 4969,4987,4993,5003,5021,5023,5039,5051,5059,5077,5087,5101,5119,5147, + 5153,5167,5171,5179,5189,5209,5227,5231,5237,5261,5273,5281,5297,5309, + 5323,5333,5347,5351,5381,5387,5399,5413,5431,5437,5449,5471,5479,5501, + 5507,5519,5531,5557,5563,5573,5591,5623,5639,5641,5647,5657,5669,5683, + 5701,5711,5737,5743,5749,5779,5783,5801,5813,5827,5843,5857,5869,5881, + 5903,5923,5927,5953,5981,5987,6007,6011,6029,6037,6053,6067,6089,6101, + 6113,6131,6151,6163,6173,6197,6211,6229,6247,6257,6277,6287,6311,6323, + 6343,6359,6373,6389,6421,6427,6449,6469,6481,6491,6521,6529,6547,6563, + 6581,6599,6619,6637,6653,6673,6691,6709,6733,6737,6763,6781,6803,6823, + 6841,6863,6883,6899,6917,6947,6961,6983,7001,7019,7043,7057,7079,7103, + 7127,7151,7159,7187,7211,7229,7253,7283,7297,7321,7349,7369,7393,7417, + 7433,7459,7487,7507,7537,7561,7583,7607,7639,7669,7687,7717,7741,7759, + 7793,7823,7853,7883,7907,7937,7963,7993,8039,8059,8093,8123,8161,8191, + 8221,8263,8297,8329,8369,8419,8447,8501,8527,8563,8609,8663,8699,8747, + 8803,8849,8893,8963,9029,9091,9157,9239,9319,9413,9533,MarkLast(9689) +}; + +static short coeff[][32] = { + {13,11,10,3,9,2,2,2,2,9,2,2,7,2,2,2,2,2,2,6,2,2,2,13,11,10,3,9,2,2,2,2}, + {23,17,12,11,14,14,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,12,12,12,14,14,14}, + {18,14,5,14,2,2,19,19,25,25,18,18,18,2,2,2,2,2,2,2,2,2,2,2,25,6,2,2,2,18,14,5}, + {18,13,23,5,2,12,6,12,12,12,10,10,16,2,16,16,2,2,2,2,2,2,2,10,2,2,2,2,10,2,2,2}, + {21,22,7,21,2,20,20,2,2,2,2,22,2,2,2,2,2,2,2,6,6,21,2,2,2,2,2,2,2,2,6,6}, + {29,19,27,32,6,8,2,2,2,2,2,8,8,2,2,2,2,9,9,9,9,2,2,2,2,2,2,2,9,9,2,2}, + {30,19,24,16,22,8,2,2,22,5,9,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {34,28,13,28,27,27,2,4,2,2,2,16,16,4,20,20,36,20,36,5,5,5,36,36,5,5,5,7,5,7,7,2}, + {35,19,33,8,21,30,8,2,4,2,4,4,2,2,2,2,2,2,2,2,2,17,2,2,11,25,11,17,17,17,17,17}, + {37,21,35,29,27,19,19,2,2,2,5,15,2,2,15,15,19,19,19,19,19,2,2,2,2,2,19,2,2,2,2,2}, + {45,44,13,25,17,47,30,2,30,2,2,2,2,2,2,2,2,2,19,19,19,17,17,2,2,2,2,2,2,2,2,2}, + {35,22,37,9,35,12,35,8,2,2,50,50,2,2,32,32,32,31,13,8,8,8,2,22,50,9,9,9,22,22,22,10}, + {29,24,43,36,49,2,2,8,4,25,49,25,2,2,8,10,10,10,5,5,5,40,10,33,40,40,2,27,10,25,25,25}, + {50,18,32,39,21,2,2,2,4,4,36,36,14,14,14,14,2,2,2,17,17,17,16,16,2,14,14,14,14,2,2,2}, + {31,28,45,20,18,43,43,13,28,2,2,2,31,31,31,31,31,2,2,2,43,43,2,2,2,2,2,2,2,2,30,2}, + {39,15,41,7,24,2,2,30,40,2,2,25,25,25,25,2,2,2,2,2,2,6,6,2,25,2,5,2,2,25,2,2}, + {44,20,29,39,7,21,21,21,2,2,45,2,2,2,49,49,49,49,49,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {56,20,22,13,18,35,35,6,2,4,2,4,2,2,2,23,16,16,4,23,2,34,52,2,34,2,4,2,2,2,23,16}, + {46,32,17,18,29,27,31,31,31,2,2,4,15,2,2,2,2,2,2,2,2,2,2,2,2,2,23,32,32,32,15,15}, + {62,42,43,17,23,13,13,2,2,13,2,2,2,2,2,2,2,10,2,2,2,2,9,10,2,2,2,19,9,9,9,9}, + {64,34,16,28,16,51,47,2,2,2,6,18,39,39,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,12,12,2}, + {74,26,44,25,50,24,54,39,58,42,2,42,42,2,2,2,2,2,2,2,2,33,33,2,2,39,11,2,2,58,39,58}, + {70,22,50,22,16,9,25,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {74,21,17,25,35,33,10,2,10,20,20,57,57,57,2,2,57,2,2,2,2,2,2,2,13,2,2,2,2,2,2,2}, + {81,18,10,11,47,38,71,37,2,37,2,2,2,2,2,26,26,26,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {55,30,85,42,16,36,45,67,2,2,68,2,2,2,2,2,2,2,68,10,2,2,2,2,2,2,2,2,2,2,2,2}, + {64,17,24,26,49,12,10,39,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,59,2,2}, + {68,57,23,38,61,38,13,13,8,2,2,2,2,2,2,2,2,2,2,2,2,2,2,68,15,2,44,44,44,2,2,2}, + {94,28,58,29,13,5,15,8,66,2,2,2,39,39,15,66,2,2,6,6,2,2,66,66,66,66,2,2,2,2,2,66}, + {94,85,9,41,41,37,29,29,17,2,2,2,7,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,8,8}, + {89,32,75,77,77,13,2,30,30,2,2,2,2,2,2,2,2,2,2,67,67,2,2,2,2,2,2,2,2,8,19,32}, + {70,45,58,63,67,10,72,72,70,6,2,36,2,70,70,6,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {101,33,76,13,45,63,2,2,6,19,2,2,32,32,32,32,32,65,2,63,63,11,11,11,19,19,19,19,9,63,63,63}, + {70,89,44,37,19,45,2,2,2,8,10,8,54,54,80,80,80,80,80,2,116,2,116,2,2,80,40,51,100,100,8,2}, + {71,54,83,51,42,98,2,2,8,8,14,30,93,22,15,15,30,30,30,44,44,44,2,2,22,22,22,117,44,11,11,11}, + {109,37,51,113,17,10,2,2,17,17,55,2,55,55,55,55,55,55,2,2,2,57,48,48,55,55,2,2,55,2,2,55}, + {75,38,68,89,11,52,2,2,81,39,2,38,2,2,2,2,2,2,2,2,2,2,2,19,2,2,2,2,2,2,2,2}, + {81,84,35,34,20,93,2,12,12,12,2,96,2,96,96,2,96,2,2,2,2,2,2,2,2,2,2,2,2,56,56,56}, + {104,32,56,46,77,11,35,35,24,56,19,2,2,2,78,2,2,75,2,2,2,2,78,2,2,2,2,2,2,2,2,2}, + {81,103,25,35,28,15,20,20,20,2,2,2,2,20,20,20,107,107,2,2,2,2,2,2,2,2,2,2,2,2,13,13}, + {119,75,42,29,74,23,54,36,39,2,2,4,4,19,19,2,2,2,2,2,2,2,2,54,2,2,2,2,2,2,2,54}, + {115,73,22,102,75,138,16,73,50,16,2,50,2,2,2,133,2,2,2,2,2,2,2,2,2,2,2,2,2,33,33,33}, + {119,48,66,51,14,22,20,20,2,2,2,2,2,60,2,2,2,2,2,2,2,2,60,2,2,2,2,2,2,60,2,65}, + {121,94,80,29,51,69,42,36,14,14,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,17,2,2}, + {129,123,41,79,43,34,24,11,2,2,4,2,2,2,2,75,16,16,16,75,75,75,16,16,16,25,2,99,2,2,75,16}, + {128,33,35,68,22,8,62,94,2,2,2,62,62,2,98,2,2,4,98,2,2,32,81,32,32,32,98,98,98,98,98,98}, + {101,109,154,15,57,6,27,36,2,2,37,37,2,2,2,2,2,2,2,107,2,2,2,107,107,2,2,2,2,2,2,2}, + {106,40,24,38,61,118,106,106,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {149,111,58,79,127,13,41,33,27,16,30,2,61,2,72,2,2,2,2,2,2,2,2,2,2,2,2,75,75,2,2,2}, + {105,92,43,156,25,53,57,115,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {99,40,62,67,66,29,99,99,99,78,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,79}, + {109,42,96,95,66,41,103,84,13,103,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {111,72,16,89,25,86,117,29,14,14,2,2,2,2,2,60,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {106,72,49,94,140,44,97,157,75,2,2,4,123,123,2,2,123,123,123,123,2,2,2,2,2,2,2,2,2,2,2,2}, + {115,67,74,32,43,50,21,36,135,36,85,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {151,71,157,42,41,37,80,27,18,2,2,2,2,2,2,2,2,2,2,2,2,2,115,128,128,128,128,128,32,2,128,80}, + {119,91,38,30,92,44,32,76,22,2,34,2,2,2,2,2,2,2,2,2,2,2,2,2,2,129,2,2,129,2,2,2}, + {121,126,31,52,120,37,57,10,171,2,2,2,2,35,35,35,2,2,97,97,97,97,97,97,97,35,35,35,97,97,97,2}, + {155,86,49,104,87,94,64,45,61,91,91,91,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {164,121,44,166,47,33,7,15,13,2,2,122,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {128,120,133,17,71,52,25,107,42,21,21,2,2,2,2,4,4,96,2,9,9,2,9,94,94,94,94,94,94,94,94,96}, + {179,82,157,76,61,35,13,90,197,2,69,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,39,39}, + {136,136,148,63,66,10,169,95,95,163,30,28,28,2,41,130,2,2,2,21,2,2,2,2,2,2,2,2,2,2,2,36}, + {131,40,112,63,55,30,53,79,79,79,2,79,2,2,2,2,2,79,2,2,2,2,14,36,2,21,21,21,21,2,2,91}, + {165,81,92,48,9,110,12,40,40,34,2,2,2,107,107,107,2,107,2,2,2,2,2,2,2,2,2,2,2,15,41,41}, + {169,66,170,97,35,56,55,86,32,32,2,2,2,2,14,2,40,2,37,2,2,37,40,40,40,2,2,2,37,37,37,37}, + {135,63,126,156,70,18,49,143,6,117,2,109,109,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {193,59,51,68,68,15,170,170,170,143,143,12,2,2,2,63,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {145,101,56,65,23,76,110,2,4,4,4,146,146,146,2,146,2,2,2,2,2,2,2,2,2,2,2,2,2,2,146,146}, + {144,129,26,98,36,46,47,52,52,52,82,2,2,2,2,2,17,2,2,2,2,2,2,2,2,2,2,2,2,91,2,2}, + {145,78,166,171,56,20,63,2,2,33,33,33,33,2,78,47,47,47,47,47,2,2,2,2,2,78,78,78,2,2,2,2}, + {191,69,176,54,47,75,167,2,2,2,188,188,188,30,30,2,67,67,117,2,117,117,117,2,2,36,2,2,2,2,2,2}, + {186,96,29,122,47,96,170,157,157,157,157,108,159,2,195,195,26,26,26,26,26,2,2,2,2,132,132,132,2,2,2,2}, + {151,118,226,91,54,49,33,2,2,2,2,4,4,4,143,143,2,2,143,25,25,25,2,143,143,143,143,143,143,143,143,143}, + {144,91,237,82,81,75,138,163,163,163,117,117,44,2,44,136,136,136,136,2,2,2,2,2,122,122,122,122,2,2,2,136}, + {189,78,178,64,118,27,189,2,2,67,67,110,110,110,110,2,28,28,2,2,2,2,2,2,2,102,2,2,2,2,2,2}, + {165,202,83,76,125,65,42,2,44,44,23,2,23,23,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {209,204,92,75,85,146,104,2,7,18,8,2,2,2,204,95,95,95,2,2,2,95,95,95,95,95,95,95,2,2,2,95}, + {169,68,89,16,193,82,33,262,262,175,148,148,148,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {171,162,78,43,61,17,112,10,171,182,118,33,2,2,2,2,118,2,2,2,2,2,2,151,2,2,2,2,2,2,2,2}, + {211,121,119,55,90,211,96,89,225,25,178,36,36,36,2,2,108,2,2,2,2,2,2,2,2,2,2,2,2,184,2,2}, + {154,101,83,17,16,210,41,79,70,158,2,27,27,2,2,2,2,2,2,2,2,2,2,2,2,153,2,2,2,2,2,2}, + {169,179,130,79,148,180,136,17,47,119,2,119,119,169,169,2,169,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {241,171,148,31,172,34,66,60,156,140,2,2,2,75,75,2,2,2,2,2,2,2,190,190,2,2,2,30,2,2,2,2}, + {229,189,183,106,118,138,82,149,265,39,2,2,265,2,2,2,2,2,2,130,2,2,2,71,71,2,2,2,71,2,2,71}, + {165,157,127,21,64,15,80,130,130,130,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,74,2}, + {221,130,203,84,83,83,29,121,54,54,2,141,2,2,94,94,94,4,4,4,2,4,2,2,2,54,54,108,16,16,94,52}, + {230,166,20,160,121,102,153,94,16,67,2,2,2,2,2,2,97,97,97,2,2,97,97,2,97,97,97,97,97,97,97,97}, + {181,79,137,119,139,24,77,17,50,25,25,25,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {239,242,192,40,41,62,124,193,193,31,193,2,2,2,2,2,2,2,2,2,2,2,2,148,2,2,2,2,2,2,2,2}, + {239,178,73,122,239,51,95,48,78,88,78,2,2,2,2,2,2,2,2,2,2,2,144,144,2,2,144,144,144,2,144,144}, + {234,117,198,34,143,21,74,6,252,252,98,2,2,2,2,197,38,2,2,2,2,2,47,2,47,47,47,47,2,2,2,47}, + {179,110,38,28,58,39,16,29,42,125,202,8,8,129,4,4,2,2,2,67,67,2,2,2,2,2,2,8,67,67,2,2}, + {246,53,189,50,18,59,179,179,7,137,137,2,2,103,103,103,103,40,40,40,2,2,2,2,73,73,73,2,103,103,103,103}, + {239,133,87,92,193,12,206,238,238,238,31,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {191,244,60,193,18,32,193,104,74,125,125,66,2,2,2,2,2,2,2,2,2,2,125,125,2,125,125,125,2,2,2,2}, + {177,74,90,91,172,219,63,84,32,2,2,196,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {253,143,54,39,122,32,75,107,234,2,6,6,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {282,89,71,88,30,23,81,105,105,2,2,105,105,131,107,2,2,2,2,2,195,195,2,2,29,29,21,21,128,195,195,195}, + {259,115,171,40,156,71,67,24,24,2,2,2,24,4,4,4,2,234,2,2,2,2,2,2,2,2,2,74,74,2,2,2}, + {264,237,49,203,247,108,75,75,75,2,2,32,16,8,16,16,16,164,14,164,2,2,32,16,8,16,16,32,42,42,42,2}, + {264,106,89,51,29,226,23,286,286,151,151,151,151,151,2,2,2,2,2,2,31,31,31,2,2,2,2,2,2,2,2,284}, + {194,215,82,23,213,23,108,127,74,2,201,32,178,2,285,2,2,2,2,285,2,2,2,2,2,2,2,2,2,2,2,2}, + {196,267,251,111,231,14,30,52,95,2,154,53,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {266,67,22,101,102,157,53,95,130,2,42,76,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {281,205,107,178,236,122,122,316,76,215,215,2,60,2,2,2,2,2,2,227,2,2,2,2,2,2,2,2,27,2,2,2}, + {271,89,65,195,132,162,102,45,56,174,104,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {200,169,170,121,155,68,131,167,78,113,113,2,2,64,2,2,2,2,2,2,2,2,2,2,2,2,2,173,2,2,2,2}, + {288,143,265,264,71,19,231,169,27,27,27,2,2,2,2,2,2,2,2,2,2,2,2,2,51,2,2,2,2,2,2,2}, + {311,141,96,173,90,119,134,151,35,252,39,2,39,39,2,2,2,2,2,2,2,2,2,113,113,2,2,2,2,2,2,113}, + {311,230,52,138,225,346,162,216,216,91,160,182,91,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {275,167,128,244,184,184,44,210,237,139,139,139,139,2,2,2,2,2,2,2,2,2,2,73,2,2,2,2,2,2,2,2}, + {176,156,83,135,46,197,108,63,33,33,33,2,133,2,213,213,213,213,133,133,2,133,2,2,133,133,2,2,2,2,2,2}, + {283,125,141,192,89,181,106,208,124,124,2,112,112,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {289,191,171,152,191,173,54,13,21,56,56,56,2,2,2,2,2,2,2,2,2,220,2,2,2,2,2,2,2,2,2,2}, + {334,305,132,132,99,126,54,116,164,105,2,105,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,287,2,2,2,2}, + {240,166,44,193,153,333,15,99,246,99,2,2,99,99,2,2,2,2,195,195,195,2,195,195,2,263,263,2,195,195,195,263}, + {246,194,265,79,225,65,24,62,46,181,2,2,2,314,2,2,2,2,2,2,2,215,2,2,2,2,2,2,2,2,2,2}, + {229,334,285,302,21,26,24,97,64,40,2,2,2,231,231,231,231,65,2,148,2,2,2,2,2,2,2,2,2,2,2,2}, + {251,295,55,249,135,173,164,78,261,261,2,2,2,2,114,2,2,2,2,2,256,142,142,2,2,2,2,2,2,2,2,185}, + {232,153,55,60,181,79,107,70,29,35,2,2,58,58,2,58,2,2,2,2,61,61,2,61,61,2,2,61,61,90,2,90}, + {246,116,45,146,109,90,32,103,133,119,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {246,113,146,232,162,262,204,47,45,331,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {360,150,84,275,13,26,368,49,244,244,63,63,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {239,295,174,87,30,87,85,36,103,36,2,278,2,2,2,2,2,2,163,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {356,300,75,310,123,301,200,107,183,37,218,37,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {358,207,168,150,150,21,156,50,195,275,275,275,2,2,2,2,2,251,2,2,2,251,251,251,251,251,251,251,251,251,2,2}, + {322,194,234,62,236,147,239,400,255,255,80,4,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {326,276,134,100,143,113,115,221,13,339,194,194,194,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {337,132,27,45,14,81,110,84,238,224,211,2,29,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {192,213,113,174,403,117,342,342,311,35,35,2,2,2,2,2,2,2,2,101,2,2,2,2,2,2,2,2,2,101,101,101}, + {264,273,316,53,40,330,51,285,115,219,147,2,2,2,335,2,2,2,2,2,173,2,173,2,2,173,173,173,173,173,173,83}, + {254,293,407,118,54,296,160,231,4,4,93,2,2,2,2,2,60,61,2,2,120,127,127,127,88,88,88,88,88,88,88,88}, + {341,78,336,263,281,164,99,334,296,114,109,2,163,163,163,163,2,2,2,2,2,2,2,125,125,292,292,292,292,125,125,125}, + {355,87,212,100,89,210,133,344,120,45,45,138,138,138,138,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {274,141,46,219,158,284,38,79,73,185,35,6,81,2,2,2,2,53,2,2,81,81,2,81,2,2,2,53,53,53,53,53}, + {349,303,439,19,95,240,174,191,2,162,162,2,2,2,76,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {360,91,201,205,67,181,59,77,2,44,103,103,103,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,125}, + {283,154,261,91,77,147,227,105,116,311,256,256,2,116,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,32,2}, + {287,288,111,89,249,370,55,16,248,67,67,115,2,2,134,134,2,2,2,2,2,2,2,2,2,2,2,2,2,22,22,22}, + {284,270,282,37,29,181,160,49,285,285,374,250,2,374,374,2,2,2,179,179,35,2,179,179,2,179,179,2,2,285,285,285}, + {359,305,52,36,243,231,7,92,2,68,68,307,62,45,2,2,112,311,311,311,2,2,2,2,2,2,2,2,2,2,2,2}, + {288,119,218,137,364,38,27,380,2,2,211,23,33,2,2,2,2,2,225,225,225,2,2,225,225,225,2,2,2,2,2,2}, + {277,155,232,309,370,365,348,75,214,214,214,4,4,2,2,2,210,210,210,210,210,210,210,2,2,2,2,2,2,2,2,2}, + {292,204,91,41,124,190,107,322,125,125,125,125,125,25,25,62,2,2,146,146,2,2,62,146,2,146,114,146,114,2,2,2}, + {282,195,192,409,68,99,253,106,2,2,2,231,55,55,2,323,323,55,55,285,285,285,285,2,2,2,2,2,2,285,285,323}, + {299,122,174,403,113,77,63,275,2,2,2,138,276,227,38,227,2,237,2,2,2,2,2,2,2,2,2,2,352,352,352,2}, + {282,222,268,86,21,109,353,408,2,2,2,2,135,12,12,216,241,241,241,241,241,241,241,241,241,303,303,303,135,135,135,2}, + {374,94,89,257,137,246,186,196,2,2,2,2,2,454,122,122,122,122,2,2,2,28,28,94,94,94,94,94,122,122,122,122}, + {288,92,62,428,122,153,481,66,2,2,2,250,250,177,177,177,177,279,279,279,279,279,279,279,2,2,279,177,177,177,177,177}, + {288,370,141,284,207,192,450,67,2,2,2,183,217,217,217,183,183,167,202,202,202,202,167,167,2,2,2,164,164,80,167,167}, + {286,293,199,39,158,332,242,103,2,2,2,408,266,315,2,2,365,253,315,315,315,315,315,2,2,315,2,2,2,2,2,2}, + {407,83,435,187,40,16,52,65,2,2,244,39,77,119,119,2,2,2,119,342,342,2,2,2,2,2,342,2,2,58,58,119}, + {398,88,78,57,260,203,203,43,131,131,131,204,204,322,204,2,102,2,325,325,325,325,2,2,2,2,2,2,2,2,2,2}, + {390,174,70,155,163,67,225,49,2,34,34,151,151,2,2,111,2,2,111,111,2,2,2,2,2,2,2,2,2,2,2,2}, + {393,129,393,169,23,192,168,47,2,2,312,150,71,2,150,2,2,2,61,2,2,61,2,2,2,2,2,2,2,2,2,2}, + {408,136,71,63,63,159,222,68,181,181,124,227,14,14,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {294,169,79,242,160,123,178,290,186,186,56,399,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {415,228,69,68,193,122,21,362,33,22,362,57,2,2,2,2,46,46,196,196,196,2,196,196,196,2,196,2,2,2,2,2}, + {415,130,241,185,312,175,309,199,94,281,47,47,2,2,2,2,206,307,221,2,2,2,2,2,239,239,239,239,239,206,206,206}, + {417,238,147,165,346,19,92,164,266,291,291,43,2,2,2,345,2,2,2,345,345,2,2,2,2,2,345,2,2,2,2,2}, + {456,192,86,182,35,174,342,102,210,210,210,393,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,256,256,158}, + {307,255,92,38,325,61,103,246,176,319,80,89,2,241,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {432,168,63,154,166,46,479,145,144,288,288,288,288,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {341,256,113,85,188,233,161,29,110,167,91,91,253,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {311,360,312,158,73,16,106,209,472,48,24,203,203,2,2,2,2,234,234,234,2,234,234,203,2,2,2,234,234,234,234,234}, + {437,196,161,100,132,246,395,187,35,35,35,2,2,35,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {438,174,338,145,155,276,422,374,4,463,463,99,224,70,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {426,225,211,130,325,283,353,96,282,23,299,2,2,2,63,63,2,276,276,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {430,101,288,38,200,332,325,193,123,123,88,2,2,2,2,2,231,231,139,139,139,139,139,139,139,139,139,139,139,139,139,139}, + {434,143,308,389,365,363,174,63,121,125,260,2,2,260,260,2,2,2,2,2,2,2,2,2,2,258,2,2,2,258,2,2}, + {453,123,201,141,229,223,234,494,102,102,102,2,2,102,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,252}, + {438,168,65,264,304,74,168,88,114,132,187,2,127,127,2,2,2,2,2,81,81,56,2,2,2,307,2,2,2,2,81,81}, + {324,181,141,129,33,171,173,291,227,373,52,301,301,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {448,119,431,111,135,50,242,95,148,49,49,49,68,2,2,2,2,2,2,2,2,49,2,2,2,2,2,2,2,2,2,2}, + {335,114,55,47,33,173,287,345,198,198,136,238,238,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {468,377,243,237,332,512,27,167,22,169,14,14,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {456,162,188,223,408,209,28,164,299,299,258,186,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {445,391,115,226,96,456,239,214,556,158,158,282,2,2,2,2,2,2,2,2,2,2,2,2,2,331,2,2,2,2,2,2}, + {360,397,130,172,407,479,295,13,38,199,199,346,2,2,2,2,2,2,145,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {512,136,129,361,180,61,274,128,422,27,292,165,2,2,2,2,2,2,363,117,117,117,117,2,2,2,2,363,2,2,2,2}, + {478,433,483,302,200,227,273,27,171,171,371,102,2,2,2,2,2,20,2,2,2,2,2,2,2,2,403,403,2,2,2,2}, + {485,158,454,86,212,60,93,40,209,188,188,106,2,231,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {390,448,111,145,47,555,367,317,315,52,429,435,429,429,2,2,2,2,2,2,2,2,229,2,2,229,2,2,2,229,2,2}, + {490,331,187,398,407,373,497,219,423,423,378,378,2,419,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {378,406,112,198,539,550,516,59,240,240,23,316,2,122,2,2,2,2,2,2,2,2,2,2,111,111,2,2,2,95,2,2}, + {474,373,248,330,40,113,105,273,103,407,2,165,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {495,406,306,239,172,323,236,50,37,435,2,310,56,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {498,447,112,241,552,119,227,189,140,140,140,140,140,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {505,132,169,418,342,28,319,301,172,530,317,317,335,2,2,2,2,2,2,376,2,2,2,2,2,2,2,2,2,2,2,2}, + {397,393,191,269,462,151,264,134,307,307,2,163,163,2,2,2,2,2,2,2,2,2,2,2,2,2,159,2,2,2,2,2}, + {485,491,325,149,122,145,228,100,311,64,2,62,137,2,137,2,2,2,2,2,2,2,392,2,2,2,2,2,2,2,2,2}, + {364,462,360,383,182,187,123,69,129,146,2,156,149,2,149,2,2,2,2,2,2,2,303,303,303,2,2,2,2,2,149,266}, + {507,195,130,401,363,171,483,20,86,464,2,89,89,2,26,2,2,2,2,2,425,425,2,2,2,2,2,2,2,2,2,2}, + {380,220,87,122,242,78,207,371,95,305,2,2,2,2,440,440,445,358,358,331,331,358,445,445,445,445,445,445,445,445,445,445}, + {507,221,247,137,182,90,28,207,325,438,2,2,2,2,2,187,232,438,2,2,68,37,37,37,37,37,37,37,37,37,161,2}, + {509,265,101,126,203,86,152,416,352,85,2,2,2,284,391,368,2,2,152,2,2,2,325,2,2,2,2,2,2,2,2,2}, + {572,359,332,480,68,535,59,504,365,21,2,2,246,54,246,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {415,178,178,372,415,400,73,82,348,99,2,23,325,44,2,2,2,2,2,2,2,2,325,2,2,2,2,2,2,2,2,2}, + {430,275,236,361,42,552,368,236,653,74,65,458,288,307,307,2,2,2,2,2,2,2,65,65,2,2,2,2,2,2,2,2}, + {434,139,58,437,130,441,188,15,63,145,145,145,300,2,2,2,2,300,2,2,2,2,2,2,2,2,401,401,401,401,401,401}, + {542,138,266,514,552,202,103,197,574,48,2,96,96,2,2,96,96,217,2,2,2,2,2,2,2,2,2,2,2,2,2,217}, + {546,494,72,272,550,219,213,209,169,404,69,464,86,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {602,466,332,458,99,244,255,183,446,670,2,186,323,2,2,2,2,2,2,2,2,2,2,2,2,2,2,292,165,165,165,165}, + {422,413,561,110,242,62,436,478,18,150,606,88,643,2,249,2,2,2,2,456,2,2,2,2,2,2,2,2,2,2,2,456}, + {522,141,154,253,264,53,120,93,274,52,44,203,556,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {600,249,375,555,421,322,317,84,517,517,268,106,353,2,2,2,2,2,2,2,2,2,268,2,2,2,2,2,2,302,2,2}, + {555,516,310,438,290,559,52,265,248,193,285,441,285,285,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {555,300,232,386,470,300,355,177,57,407,450,279,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {544,177,79,306,256,402,205,496,398,115,115,43,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {534,274,194,220,575,81,206,544,341,85,137,429,429,429,429,344,2,2,2,2,2,315,315,315,315,315,315,72,72,72,2,2}, + {400,136,112,136,273,277,205,578,122,122,230,230,2,2,2,2,2,2,2,2,2,2,2,2,2,2,302,2,2,2,2,2}, + {576,421,115,52,253,373,17,657,43,178,178,58,485,485,485,485,485,485,2,2,2,159,159,159,159,2,619,2,2,2,2,2}, + {576,301,142,329,96,41,302,528,126,112,206,206,2,2,2,2,2,2,206,206,2,206,206,2,191,206,206,191,191,191,191,206}, + {548,538,508,250,539,102,73,285,119,433,480,480,2,2,2,480,480,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {622,526,294,56,498,176,237,351,25,26,474,55,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {446,163,469,481,240,278,51,373,491,13,22,419,2,2,2,2,2,2,2,2,2,176,176,2,2,2,2,2,2,2,2,2}, + {445,223,102,108,120,166,68,214,737,504,96,96,206,377,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,39,39,528}, + {453,121,489,84,434,505,78,575,468,372,468,468,83,468,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {425,355,128,58,194,82,438,117,10,34,34,35,112,107,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {432,479,328,443,253,634,271,429,406,543,406,543,543,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {433,294,192,205,152,70,99,68,392,169,309,390,390,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,199,2,2,2}, + {456,383,487,311,57,579,673,264,582,187,184,43,43,2,2,2,2,501,501,501,2,2,2,2,2,2,2,2,2,2,2,2}, + {437,561,384,619,363,420,614,117,217,247,405,142,142,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {434,372,239,508,478,26,375,255,151,151,650,112,251,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {437,133,516,423,305,90,135,25,266,487,6,286,286,2,2,2,2,2,2,2,2,2,2,2,2,510,510,2,2,2,2,2}, + {463,341,170,401,178,79,305,98,162,166,32,392,335,335,335,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {610,477,478,516,318,184,267,423,190,494,494,2,336,336,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {611,211,491,224,47,54,124,268,271,271,223,2,2,2,2,2,2,2,2,2,2,2,2,359,2,2,2,2,2,2,2,2}, + {590,463,461,162,162,622,167,254,29,377,377,75,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {478,388,612,404,491,561,180,80,262,58,94,2,2,275,2,2,2,2,2,151,2,2,2,2,2,312,312,312,2,2,2,275}, + {629,225,67,623,298,588,354,49,41,185,176,63,63,63,2,2,2,2,2,2,2,2,2,2,2,2,8,435,32,32,435,435}, + {671,275,392,298,612,328,337,215,58,58,124,2,2,490,392,2,2,2,125,457,457,2,2,2,2,2,2,2,2,2,2,457}, + {448,126,129,168,209,340,40,96,509,509,509,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {667,246,160,68,737,203,168,628,46,128,358,2,2,2,121,121,2,2,2,2,2,2,560,121,2,2,2,2,2,2,2,121}, + {635,212,284,356,187,591,275,361,194,317,488,2,2,2,2,2,2,97,6,2,6,247,2,2,2,2,2,2,2,2,2,6}, + {612,395,104,86,264,321,521,325,252,53,178,100,100,100,16,343,343,343,343,343,2,2,2,2,2,2,2,2,2,343,343,343}, + {486,428,287,472,292,141,504,178,585,98,282,2,2,2,2,2,2,2,2,2,2,2,2,284,284,284,78,284,2,2,2,2}, + {612,327,212,565,450,385,201,649,423,491,106,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {462,579,236,447,60,162,427,258,73,742,742,2,742,742,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {495,440,89,439,65,207,459,407,139,131,624,2,380,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {469,507,276,227,66,237,260,386,27,666,31,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {646,393,273,238,24,13,253,127,368,316,316,316,150,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {626,196,111,465,386,431,181,414,614,391,349,318,389,2,389,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {667,257,290,122,109,523,95,26,282,49,374,236,236,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,88,2,2}, + {653,169,261,533,488,282,213,443,337,480,503,174,534,2,2,2,2,2,534,2,2,2,2,534,2,2,2,2,534,2,2,2}, + {670,555,160,90,604,604,50,459,376,545,316,180,526,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {639,253,95,380,108,448,223,254,381,30,6,644,6,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {642,160,702,90,157,254,278,521,650,277,74,554,122,2,2,2,2,2,2,517,174,174,174,2,2,2,2,2,2,2,2,2}, + {678,254,190,197,637,49,130,25,374,357,357,411,643,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,537,2,2}, + {512,347,65,546,434,87,18,123,672,412,316,6,699,6,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {657,233,108,38,147,53,136,168,408,477,477,279,268,289,2,2,2,2,2,2,289,2,2,2,2,2,2,2,2,289,289,2}, + {498,431,217,101,78,143,111,113,181,825,458,140,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {660,624,376,472,165,66,158,308,492,779,305,305,2,576,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {521,249,388,155,467,245,134,311,72,312,312,623,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {408,348,216,299,302,668,347,63,172,141,272,168,678,2,2,2,512,2,2,2,2,4,2,2,2,494,64,64,64,128,16,512}, + {669,421,230,70,212,845,237,347,148,76,823,472,2,2,2,132,2,2,2,2,2,2,2,383,132,383,2,2,383,383,383,383}, + {693,530,139,82,780,416,270,278,330,484,484,200,2,2,2,2,137,94,2,2,2,2,2,2,2,2,484,2,2,2,2,2}, + {672,150,164,622,196,75,302,119,42,314,314,132,60,60,60,298,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {705,302,411,705,691,160,809,40,32,867,826,826,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {684,229,138,46,407,399,82,254,267,31,31,45,2,209,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {707,323,409,27,31,157,492,463,886,412,251,251,304,190,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {715,521,636,304,402,459,435,571,611,214,214,43,43,358,2,2,2,2,358,2,2,2,2,2,2,358,358,358,2,2,358,358}, + {768,224,219,425,467,147,151,643,316,263,263,263,263,263,2,2,2,2,2,272,139,2,2,2,2,2,2,2,2,2,272,53}, + {555,543,434,78,850,174,277,194,4,100,471,69,69,424,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {690,206,572,877,600,129,288,52,19,147,222,222,147,147,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {793,279,264,566,252,495,872,492,482,107,294,503,350,350,2,2,2,2,2,2,2,285,285,273,273,273,273,2,2,2,2,2}, + {703,427,225,320,136,47,103,547,239,217,73,68,68,204,204,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {791,275,60,137,352,839,67,476,356,216,216,563,563,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {703,312,472,588,228,512,386,668,477,617,389,389,389,2,296,2,2,2,2,343,343,2,2,343,343,2,2,617,617,617,617,2}, + {709,509,697,145,252,194,304,192,192,623,623,4,423,2,2,2,199,423,2,2,2,222,222,2,2,623,623,623,623,623,2,222}, + {587,453,117,107,672,86,248,568,568,294,294,513,78,2,2,164,82,2,2,2,2,22,2,2,2,2,2,2,2,2,2,2}, + {741,466,378,135,737,131,159,469,59,2,59,59,187,2,204,2,2,2,2,2,2,2,2,2,798,2,2,798,798,798,798,798}, + {539,310,463,103,553,45,609,326,197,2,62,113,272,2,62,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {750,703,182,242,92,335,272,466,594,2,701,569,474,129,140,140,2,507,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {547,210,113,361,584,121,65,307,98,2,2,552,514,514,2,514,207,514,514,514,2,2,2,2,2,2,2,2,2,2,2,2}, + {555,229,328,91,272,815,483,749,468,2,92,92,4,92,2,2,2,258,258,258,2,258,258,2,2,2,2,258,2,2,258,258}, + {580,145,358,434,630,73,604,366,366,2,2,398,398,207,2,207,487,2,2,487,207,2,2,207,207,207,2,2,2,2,207,207}, + {457,520,93,460,275,525,300,184,354,147,147,147,147,179,82,82,82,82,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {872,630,513,218,719,174,197,104,86,281,281,281,541,642,281,94,2,45,94,2,335,335,2,2,2,2,2,2,2,2,2,84}, + {765,421,129,298,867,365,222,476,401,142,90,22,22,88,226,657,2,2,477,2,2,2,2,2,226,226,2,226,2,2,2,226}, + {833,634,228,520,113,329,279,420,581,2,2,385,385,110,450,2,733,2,2,2,561,561,2,561,2,2,2,2,2,2,2,2}, + {587,553,360,539,227,800,312,143,536,2,2,2,64,64,64,2,2,2,179,179,493,2,2,184,184,184,58,2,2,2,493,493}, + {744,466,389,280,229,134,363,177,389,2,2,2,536,273,536,536,536,536,168,45,45,45,45,2,2,2,2,2,2,2,2,2}, + {841,222,158,469,253,91,347,241,766,2,2,2,88,88,88,439,439,439,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {462,653,478,67,269,150,474,711,220,669,669,669,669,669,390,352,325,2,229,545,545,545,545,545,545,545,545,2,545,352,309,352}, + {468,430,849,689,202,427,45,34,105,2,2,2,2,4,4,4,4,4,4,4,2,2,2,4,4,4,4,4,2,2,2,2}, + {610,289,503,744,775,512,605,454,484,2,2,2,444,466,145,631,2,631,631,631,631,631,631,631,631,631,2,2,631,631,631,858}, + {792,169,306,843,246,123,293,229,483,2,2,2,165,163,163,163,163,440,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {563,325,717,766,440,705,290,123,228,2,2,2,32,64,146,2,2,2,116,79,79,2,146,146,79,79,79,2,2,146,146,79}, + {795,185,350,211,82,537,106,680,62,2,2,537,423,423,423,2,2,501,501,2,501,2,501,2,2,2,2,2,2,2,2,2}, + {633,425,295,548,497,163,381,461,89,2,2,831,583,896,38,2,625,2,2,2,276,276,2,2,276,2,2,2,2,2,2,2}, + {767,318,84,97,208,387,423,196,417,2,396,396,396,396,396,128,128,2,2,2,328,328,4,4,4,4,101,2,2,328,82,16}, + {802,533,869,638,67,192,805,223,219,2,2,191,178,178,77,77,2,2,2,2,431,431,2,2,2,431,431,2,2,431,2,2}, + {781,638,410,399,336,465,856,426,28,2,4,4,6,6,2,2,2,449,372,372,449,449,449,2,2,449,449,449,449,449,449,2}, + {807,377,237,443,388,286,158,349,491,32,32,260,260,260,2,2,260,615,615,615,2,2,260,260,260,260,260,615,615,615,615,615}, + {780,359,766,618,41,596,86,636,287,707,707,96,49,373,613,373,2,2,2,2,2,2,2,613,613,613,2,2,2,2,2,2}, + {788,497,334,93,319,169,273,540,904,2,903,569,569,569,272,272,2,2,2,2,571,571,571,571,571,571,571,571,571,571,571,571}, + {622,309,913,550,994,90,257,588,29,526,526,526,496,496,576,2,2,2,2,2,182,182,182,2,2,447,447,447,447,447,447,182}, + {814,652,456,774,624,870,27,739,464,2,108,578,578,561,295,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {818,280,99,873,165,426,341,74,479,342,727,684,684,662,662,2,2,2,2,2,2,662,2,2,2,2,2,2,2,2,2,2}, + {593,411,953,203,89,57,785,354,349,424,424,707,707,707,829,2,2,2,2,2,670,670,670,2,2,424,424,424,2,2,670,424}, + {629,560,621,245,683,633,495,551,472,2,31,74,489,684,555,684,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {901,490,693,410,666,119,703,593,201,61,70,70,774,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,418,418}, + {669,321,391,548,189,157,337,42,796,871,276,622,30,2,2,2,2,2,2,2,580,580,107,2,2,2,2,2,434,434,434,434}, + {610,236,633,300,681,358,72,281,148,466,466,283,275,2,386,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {929,360,102,893,329,136,515,33,170,581,268,35,777,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {859,584,475,745,506,900,40,869,143,612,175,275,209,12,12,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {822,581,76,382,72,347,964,324,137,61,61,28,623,351,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {655,330,324,151,166,431,58,174,142,115,1003,66,724,778,2,2,2,503,503,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {867,820,301,252,61,331,105,309,562,218,365,326,768,672,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {623,330,182,489,212,223,741,490,40,412,801,681,681,801,2,2,71,2,2,2,2,2,2,427,2,2,2,2,2,2,2,2}, + {859,844,510,859,118,190,550,29,159,622,622,382,258,382,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {612,237,272,53,534,682,372,935,494,536,536,599,599,599,2,536,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {843,730,235,233,816,495,598,134,131,604,227,378,378,553,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {667,397,121,526,321,660,848,729,357,137,268,711,521,521,2,2,2,2,2,2,2,2,2,2,2,2,2,194,2,2,2,521}, + {939,783,796,676,259,643,103,289,15,471,80,80,2,239,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,888}, + {670,595,333,257,907,413,548,341,327,350,612,700,700,700,700,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {678,274,695,790,169,701,707,1084,470,123,846,846,217,121,317,2,2,2,83,83,83,83,83,83,83,83,83,2,2,2,2,2}, + {877,181,375,79,199,256,223,295,135,371,395,354,2,307,944,2,813,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {882,417,475,424,311,646,346,207,74,157,590,356,2,2,324,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {708,442,186,698,345,103,687,463,163,416,416,107,2,2,2,375,375,416,6,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {865,675,786,568,112,197,225,348,372,497,215,215,2,2,2,2,159,159,150,224,224,141,2,2,2,2,141,141,141,141,141,141}, + {844,244,672,489,839,263,14,233,422,392,8,392,2,2,2,2,2,2,815,815,815,815,257,257,105,105,2,2,2,815,815,815}, + {693,726,117,167,535,725,224,78,716,100,460,299,2,2,2,2,921,744,2,2,2,2,2,378,2,2,178,178,178,2,178,178}, + {898,559,396,742,51,143,411,221,116,756,756,756,2,2,2,701,701,2,2,2,2,240,225,256,322,322,240,240,240,240,240,322}, + {697,540,358,391,932,309,103,73,35,353,353,503,2,2,353,134,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {982,579,548,413,416,103,71,101,1039,526,684,684,2,2,656,2,2,2,2,2,2,2,2,2,2,2,656,656,656,2,656,656}, + {695,881,335,126,429,476,772,667,974,98,433,49,129,129,2,2,2,2,2,2,2,2,2,2,544,2,544,2,2,2,2,544}, + {859,361,215,569,255,378,543,436,220,34,105,105,816,816,816,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {705,770,134,178,940,944,654,600,46,797,797,591,2,145,616,2,2,2,2,2,2,389,389,2,122,2,2,2,389,389,909,389}, + {642,757,247,513,372,54,546,971,271,61,61,1018,2,143,332,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {988,271,675,163,379,108,48,472,870,485,485,18,2,485,528,528,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {865,827,614,74,725,685,724,190,178,272,835,722,2,35,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {923,397,722,186,203,575,24,144,36,526,206,787,12,100,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {718,359,103,558,684,560,67,35,120,342,680,265,265,265,2,2,265,2,2,2,2,2,2,2,2,2,430,2,2,2,2,2}, + {927,493,988,194,97,1006,377,578,105,248,707,784,98,784,2,2,2,2,2,2,2,2,2,370,370,2,370,2,2,2,2,2}, + {900,455,485,601,353,69,67,965,25,226,314,314,883,923,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {903,259,153,106,289,916,861,41,441,368,131,131,262,671,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {945,358,160,196,82,403,362,195,376,877,521,336,521,77,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {912,516,108,555,306,274,55,197,565,174,659,208,441,441,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {753,242,194,619,345,94,463,485,163,85,412,575,270,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {940,226,320,666,269,54,542,174,109,290,754,524,649,2,202,2,2,2,2,2,2,2,776,202,776,776,776,2,2,202,202,202}, + {915,210,456,377,303,237,225,521,621,175,569,20,124,2,601,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {672,652,792,253,796,404,171,90,406,433,43,159,72,2,2,372,2,540,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {733,439,537,37,149,650,916,443,743,621,921,664,664,2,2,2,2,2,682,523,523,523,2,2,523,523,523,523,523,523,523,523}, + {982,344,812,567,243,52,246,369,439,205,600,739,730,2,2,2,61,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {982,604,126,65,633,657,22,776,161,45,725,44,4,2,2,2,2,2,2,2,2,2,269,269,2,2,2,2,2,2,2,2}, + {745,600,284,1117,459,1135,300,52,845,331,334,334,334,2,334,334,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {936,409,217,57,574,395,481,245,548,268,447,598,375,2,192,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {986,241,233,45,721,325,350,222,35,1065,1065,1065,1065,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {755,796,877,981,259,194,1180,215,90,658,662,662,662,2,36,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {981,626,987,827,466,458,578,346,475,223,223,223,342,1058,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,728}, + {949,422,941,491,66,786,592,429,307,123,40,478,478,478,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {992,723,625,251,431,544,309,466,700,644,484,837,904,320,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1077,496,819,340,974,122,39,1209,819,18,461,648,648,394,2,2,2,2,2,2,61,2,2,2,2,2,394,2,2,2,2,394}, + {999,674,212,673,279,579,462,754,89,866,345,110,110,887,2,2,2,2,2,707,707,2,2,2,2,2,2,2,2,2,2,707}, + {1083,356,367,357,559,213,606,477,71,103,790,103,299,299,2,2,2,2,2,2,406,406,2,2,2,2,2,2,2,2,2,2}, + {1005,260,389,960,501,714,118,73,334,1019,704,204,504,205,822,822,2,2,2,2,2,2,2,2,2,2,684,2,2,2,2,2}, + {738,749,769,610,306,326,328,578,479,840,840,840,68,192,2,150,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1010,937,449,474,154,456,766,318,275,444,709,2,778,778,778,806,779,779,2,2,2,2,2,2,2,2,806,2,2,2,287,287}, + {1011,780,134,945,183,42,741,25,252,164,205,222,222,222,147,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1118,427,294,404,268,217,922,515,19,1045,1045,2,833,291,448,2,2,2,2,2,2,2,2,2,175,2,2,2,2,2,2,2}, + {1094,640,912,223,67,472,623,623,1244,65,1009,1209,1209,812,387,2,2,2,513,2,2,2,2,2,2,2,2,2,2,2,1209,234}, + {722,375,264,390,515,498,1161,391,884,551,238,2,2,825,549,2,2,2,551,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {792,250,299,210,496,682,94,207,220,227,227,2,2,227,73,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1022,409,93,359,983,345,280,280,104,940,940,2,2,382,1039,2,2,2,2,831,2,2,2,2,2,2,2,2,2,2,2,2}, + {1027,925,413,335,327,826,250,122,293,773,564,541,420,420,420,774,763,2,2,2,2,2,2,900,110,110,2,763,2,2,2,2}, + {1028,730,807,119,209,146,230,498,164,309,309,2,2,2,693,912,430,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {997,525,680,120,466,728,288,110,1082,544,572,2,2,663,290,290,2,2,754,2,2,2,2,2,582,582,582,582,582,2,2,2}, + {1055,395,795,561,222,85,294,433,377,89,89,2,2,2,456,821,2,2,821,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {997,614,240,638,755,575,874,321,600,235,665,2,2,2,154,154,767,767,2,767,2,2,2,2,2,2,2,2,2,2,2,2}, + {802,298,672,424,104,623,152,159,476,760,66,2,2,2,215,215,490,490,490,2,2,2,2,2,490,490,490,490,490,490,490,490}, + {1128,788,124,501,561,1015,419,787,48,620,705,2,2,2,2,88,18,2,215,215,215,2,2,215,215,2,2,2,215,2,2,2}, + {807,433,721,434,449,242,170,842,21,4,642,2,2,2,2,2,4,4,4,4,2,856,856,856,885,885,856,856,856,856,856,885}, + {755,612,235,265,369,855,414,362,478,518,518,2,2,64,16,8,32,4,16,8,8,1041,501,1041,2,2,64,16,8,8,16,270}, + {1004,719,1041,460,551,516,135,417,130,698,698,2,2,2,655,655,655,655,655,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1017,568,930,1113,556,1299,114,881,690,475,641,2,2,2,779,779,103,2,528,2,2,2,2,2,2,528,528,2,2,910,910,2}, + {814,473,286,752,476,779,420,569,742,164,490,2,2,2,793,812,812,812,2,812,812,2,2,526,526,812,526,2,2,2,526,526}, + {818,301,273,664,206,971,895,590,912,523,523,2,2,452,384,255,2,130,130,130,130,865,2,2,2,255,2,2,2,2,2,2}, + {820,249,292,1017,1017,143,403,37,433,456,515,2,2,69,640,2,2,2,2,2,2,2,2,2,2,2,2,824,824,824,2,2}, + {1078,527,589,244,170,892,827,606,1165,773,189,2,2,240,22,2,2,2,2,2,2,759,621,621,621,621,621,621,621,621,621,621}, + {865,1132,428,582,254,408,536,376,825,116,116,1266,1266,1266,705,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1062,268,389,1325,598,276,1270,48,572,439,302,2,544,609,544,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1065,517,247,1142,247,674,385,120,592,177,98,2,956,364,275,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {635,503,594,203,456,1246,221,396,1151,178,66,2,781,587,86,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1157,395,446,280,1130,695,668,271,111,882,477,615,615,615,2,2,2,2,2,2,2,2,615,615,615,615,615,615,615,2,305,2}, + {830,397,932,519,818,113,367,694,88,535,535,414,343,175,2,2,2,2,2,2,2,2,2,2,414,864,2,2,864,864,864,864}, + {793,463,329,730,390,551,968,92,511,470,424,563,672,563,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1050,749,809,479,87,757,288,172,597,722,4,418,418,390,2,2,2,2,2,390,390,2,2,2,2,2,2,2,2,2,2,2}, + {1084,402,130,1077,276,154,1068,779,511,853,83,757,757,38,2,2,2,2,2,202,2,2,2,2,2,2,2,2,2,2,757,2}, + {1090,255,271,110,159,235,158,236,271,815,1300,416,416,416,2,2,416,416,2,2,2,399,791,791,2,791,2,2,2,2,791,791}, + {1058,417,271,172,312,363,184,191,28,183,759,214,759,39,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1057,385,263,395,901,274,727,340,1117,263,813,870,858,429,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1102,846,985,1085,764,124,764,51,874,612,478,801,478,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1120,665,311,695,319,1033,511,297,602,1030,1030,714,240,240,2,2,2,2,2,2,2,2,2,2,2,2,2,953,2,2,2,2}, + {814,293,763,661,575,631,524,636,112,691,595,1103,405,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1106,662,258,190,1315,214,530,263,318,904,877,1317,318,2,510,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1206,469,299,1052,655,114,189,213,321,188,64,475,475,2,2,662,662,662,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1120,1159,358,347,838,207,357,167,476,52,672,38,822,2,2,2,2,2,2,213,2,2,2,2,2,2,2,2,2,2,2,2}, + {1076,596,553,545,79,727,881,121,298,169,639,368,695,115,115,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1130,177,84,673,350,543,543,95,128,954,430,884,884,2,884,884,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1222,412,430,707,691,746,131,607,311,607,112,217,912,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {820,461,681,382,273,273,358,274,274,91,887,676,386,2,676,676,2,2,2,2,2,2,2,200,2,2,2,2,200,2,2,2}, + {1096,1166,209,407,1127,400,974,322,428,906,631,134,171,2,2,2,2,664,664,664,2,2,2,2,2,2,2,2,2,2,2,2}, + {1091,946,437,51,527,802,597,639,587,645,510,586,586,2,2,2,2,2,2,2,2,2,2,2,2,2,2,168,168,168,168,168}, + {1148,585,868,1282,666,417,733,1231,515,332,1213,337,337,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1103,276,174,408,233,170,955,108,530,354,585,38,677,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,476}, + {1167,478,1169,1053,563,371,108,772,413,497,1338,991,660,2,2,2,2,2,2,2,2,2,2,2,2,27,2,2,2,2,2,2}, + {1108,437,1160,324,868,686,361,399,786,1161,1161,707,731,731,655,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1116,331,280,422,1109,341,570,243,849,241,566,61,608,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {898,782,478,1208,196,983,608,537,196,1141,141,296,715,715,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1099,1187,300,240,268,413,1366,634,184,768,773,365,783,224,783,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1158,945,300,1115,205,495,435,302,187,774,774,843,843,284,284,2,2,2,2,909,933,933,933,2,2,909,909,2,2,2,2,909}, + {904,660,1283,46,33,124,416,218,152,970,1241,305,307,307,307,260,894,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1127,553,287,58,739,99,514,739,766,42,580,241,598,598,936,936,936,629,629,629,629,2,2,2,2,2,2,2,2,2,2,2}, + {1142,370,287,925,307,1232,129,11,1284,1056,33,33,536,521,2,1286,2,2,2,2,2,2,2,2,2,2,2,2,847,847,847,847}, + {1140,814,528,677,84,1192,305,637,335,451,103,325,77,969,2,651,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1297,600,419,985,846,493,186,109,147,239,197,762,762,327,327,1004,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1181,615,482,653,238,130,313,506,98,1314,730,730,730,730,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {921,613,628,1288,111,150,191,233,633,83,387,602,105,394,2,2,2,2,2,2,2,351,2,2,351,351,351,2,2,2,351,351}, + {1192,555,586,516,1288,733,64,653,364,273,421,215,75,75,2,2,2,2,2,2,953,953,953,953,8,383,383,2,161,383,953,953}, + {1160,617,505,1205,374,906,23,408,194,91,91,91,585,984,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1203,1101,497,352,254,309,464,123,607,1080,265,1145,1145,1145,284,284,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1210,656,1026,782,802,442,1319,734,794,165,165,796,93,796,2,829,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {963,646,721,1161,219,667,1088,485,692,692,663,535,553,662,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,686,686,2}, + {966,590,140,297,189,844,633,12,847,742,742,244,281,34,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {941,231,1038,309,173,770,413,560,855,660,721,1103,721,721,721,2,2,2,2,2,2,2,2,2,2,2,2,174,2,2,2,2}, + {1213,305,656,983,1399,1196,692,986,9,339,754,308,2,308,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {972,768,1109,523,642,546,1452,29,1296,13,813,813,2,1496,2,2,2,2,2,2,2,165,165,165,165,165,165,2,2,2,2,544}, + {1330,671,528,831,1426,735,33,425,364,119,363,978,2,761,483,476,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1188,217,838,237,379,202,785,949,479,169,348,872,2,872,872,2,2,2,2,2,2,1028,2,2,2,2,2,2,2,2,2,2}, + {1190,286,513,881,390,215,387,130,749,554,1110,519,160,160,160,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1247,353,973,217,1044,1318,1115,319,203,390,1244,225,2,2,508,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {893,560,132,1420,721,191,568,799,412,22,322,93,2,4,4,4,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {987,774,678,175,145,264,588,97,1308,6,828,1129,2,2,2,45,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {955,980,753,96,574,259,1327,556,342,1415,1036,1036,478,478,478,478,507,2,2,2,2,151,2,2,317,2,2,2,2,2,595,595}, + {882,1038,211,110,942,337,1305,1225,661,183,381,381,2,2,2,2,347,2,2,2,2,2,2,2,600,431,431,431,431,431,431,431}, + {1208,486,343,725,677,1204,135,139,924,170,1111,317,2,2,2,2,202,706,202,107,107,107,2,2,706,706,107,107,2,2,2,706}, + {1259,1017,456,298,443,838,137,744,551,334,36,951,2,2,2,699,718,2,2,984,2,2,2,2,2,2,984,984,2,2,2,2}, + {1212,1186,641,284,565,636,895,82,690,117,184,184,2,2,2,397,902,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1263,370,132,635,381,47,537,179,1192,301,1282,33,2,2,2,1553,2,2,2,2,2,2,2,2,2,2,2,307,307,2,2,2}, + {1223,433,252,572,424,82,221,107,382,430,203,461,2,915,362,964,2,2,964,2,2,2,2,2,964,964,964,964,964,485,485,485}, + {1015,593,112,1408,51,104,199,221,931,1010,928,928,2,2,878,878,2,2,2,2,731,731,2,731,731,2,731,2,731,731,731,2}, + {1220,410,1193,352,260,434,469,41,1090,961,961,728,2,2,330,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {898,1043,391,1289,29,830,184,321,1136,85,1133,1082,864,864,2,2,2,2,2,2,2,2,789,789,2,789,789,2,2,789,789,2}, + {1223,434,851,152,140,1495,190,397,925,37,1080,430,2,2,204,2,759,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {989,1043,184,232,64,403,284,745,171,171,995,223,380,380,1400,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {939,1070,1288,254,973,901,321,109,568,713,336,988,2,946,262,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1276,636,569,258,325,675,342,85,88,579,833,833,833,833,520,2,2,2,2,520,520,2,2,2,2,2,2,2,2,2,2,2}, + {982,508,815,214,206,602,448,685,446,572,1549,8,1047,1047,1047,2,2,2,2,2,2,2,363,502,2,2,71,363,2,2,363,363}, + {1288,1398,789,514,151,600,1618,1194,1419,441,234,204,1191,438,828,2,857,857,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1002,342,1045,757,1008,979,322,240,1211,171,552,123,2,129,129,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1401,402,525,293,97,223,452,808,61,169,1023,1023,886,886,1023,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1001,644,263,164,136,939,624,95,489,1023,1107,331,331,10,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1011,475,845,532,567,951,663,295,877,1275,227,39,618,683,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1401,741,509,797,47,157,1256,482,1513,899,736,780,780,210,2,2,2,2,783,783,2,2,2,2,2,2,2,2,2,2,2,2}, + {1047,880,369,402,641,446,639,586,277,396,419,275,825,820,2,2,2,238,238,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1273,701,831,1294,1580,137,162,415,563,11,92,116,116,116,2,2,2,2,2,2,1029,1029,1029,504,504,877,877,877,877,877,1029,1029}, + {1335,400,315,412,172,125,568,1024,58,601,398,985,640,577,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1453,947,486,485,453,415,1164,684,504,605,422,998,727,727,2,2,2,1136,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1057,1198,146,529,284,1286,160,135,75,686,648,1425,821,586,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1052,442,936,64,132,1378,1323,161,161,161,230,131,12,12,2,2,2,2,2,2,2,2,2,2,2,998,998,998,998,2,2,2}, + {1422,838,234,554,736,243,344,526,1108,33,1303,699,249,305,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1279,681,423,813,806,269,412,420,985,485,761,1013,649,796,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {938,614,523,557,898,624,178,461,287,985,371,371,260,613,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1335,834,652,528,536,523,497,60,173,777,238,59,4,59,4,8,2,2,2,559,559,559,559,559,559,559,2,2,559,559,559,2}, + {1040,998,324,93,887,497,1326,443,152,1193,595,80,80,80,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1302,1116,283,1006,891,838,768,373,468,968,1178,1178,1269,1269,876,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1027,1128,114,395,357,417,848,22,389,1257,734,838,838,301,900,2,90,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1072,724,717,877,873,369,1031,698,917,1641,1641,1641,53,549,549,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {987,1243,424,240,53,1150,558,292,1107,574,814,1474,1474,1068,1186,2,2,2,2,2,2,2,2,2,2,2,2,2,2,859,2,2}, + {1040,420,960,882,64,661,292,146,976,427,689,248,248,248,638,2,2,2,2,2,2,2,2,2,2,2,2,861,861,861,2,861}, + {1040,522,666,398,78,208,293,818,134,867,147,147,482,2,4,629,629,629,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {987,1280,1245,1300,926,676,56,546,541,690,84,42,1000,1383,1383,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1324,588,1378,592,1445,1029,759,1296,739,931,363,704,312,704,704,704,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1058,454,1557,191,129,297,695,1390,1274,460,923,923,923,2,4,1059,2,2,2,2,2,2,2,2,2,2,2,2,2,1059,2,2}, + {1327,572,282,1022,907,1276,409,643,1050,633,187,187,187,2,228,45,2,2,2,2,2,2,320,2,2,2,2,2,2,2,2,2}, + {1395,958,237,101,559,891,560,47,524,747,197,589,589,917,887,887,887,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1395,529,461,402,194,392,122,781,111,162,780,593,593,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1378,541,848,848,347,856,102,104,183,156,395,130,1377,2,2,2,159,159,159,2,2,2,772,2,2,2,2,2,2,2,2,2}, + {1062,212,784,63,252,873,1302,1108,1380,84,1375,1375,1375,2,2,2,375,374,2,980,2,2,2,980,980,980,2,2,2,2,2,2}, + {1384,549,430,781,946,879,901,924,741,114,14,451,36,2,2,2,2,287,287,287,803,803,803,803,2,2,2,803,803,803,803,803}, + {1413,627,1329,1092,526,197,31,417,1149,981,964,1003,685,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,717,717,2,2}, + {1084,1174,1601,949,910,960,500,461,1290,23,1042,636,212,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1420,531,626,738,376,537,814,206,990,235,847,812,201,201,201,201,726,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1420,624,363,537,1436,278,292,377,263,820,376,382,382,2,654,655,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1091,793,1353,208,506,599,846,503,1011,247,289,61,1050,61,61,61,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1533,978,284,156,914,162,685,1184,252,1375,189,256,640,2,640,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1531,692,414,277,541,1371,1447,682,536,109,432,1240,1240,2,1022,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1128,398,791,1170,76,661,408,259,756,495,79,553,10,10,1532,1532,1532,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1040,704,618,854,374,1470,274,383,941,519,351,351,351,351,351,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1427,988,498,1529,99,678,1323,149,33,426,543,543,335,1507,772,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1102,349,490,266,144,220,599,437,743,764,647,1128,605,265,324,324,324,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1118,496,645,592,354,1133,935,428,72,532,182,182,1370,660,123,2,294,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1102,1042,315,745,1006,771,630,68,587,1187,295,295,295,408,408,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1118,724,1322,405,199,614,1087,885,1313,317,769,660,660,1158,535,2,2,2,373,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1120,772,743,488,346,126,784,584,943,153,311,133,133,969,605,605,2,2,2,2,2,2,2,605,2,2,2,2,2,2,2,2}, + {1404,284,176,590,1128,1371,322,543,1136,546,1315,174,174,777,777,891,2,2,2,2,2,2,2,579,579,579,579,2,2,2,2,2}, + {1441,791,233,141,141,316,89,296,462,1263,758,482,599,599,578,341,2,2,2,2,2,2,2,2,2,2,2,525,525,525,2,2}, + {1413,406,700,547,1166,250,518,543,104,331,205,205,691,691,2,2,118,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1485,400,1497,168,82,680,1103,554,249,702,493,101,296,236,2,236,236,944,944,2,2,394,2,2,2,2,2,2,2,2,2,2}, + {1127,869,558,533,1215,194,1762,784,593,777,1153,1079,1079,1079,2,2,2,330,1045,2,1045,1045,2,2,2,2,2,2,2,2,2,787}, + {1459,1243,467,533,266,1364,1031,890,1402,486,1678,1678,93,978,2,2,2,978,947,947,2,978,2,2,2,2,2,2,2,2,2,2}, + {1139,809,117,522,955,1096,1120,1470,116,184,1565,1565,557,557,2,2,2,2,2,829,1326,2,2,2,2,2,2,2,2,2,2,2}, + {1142,984,1044,590,340,241,662,357,366,1305,2,125,631,474,2,2,2,980,2,2,2,2,2,2,2,2,2,2,2,2,2,1273}, + {1469,1247,1277,616,209,486,106,552,219,217,471,272,272,1201,2,2,503,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1148,542,1478,496,950,464,1011,235,136,180,2,416,758,453,2,909,2,2,2,2,2,2,2,2,2,2,1019,1019,2,2,2,2}, + {1495,1178,874,415,1100,368,1057,1228,562,215,31,31,680,680,680,1208,2,2,2,2,2,2,2,2,1208,2,2,2,2,1208,1208,2}, + {1497,1166,1613,1403,107,803,993,539,1436,1289,2,240,334,634,532,1147,2,2,2,2,2,2,2,117,2,2,2,2,2,2,2,117}, + {1617,289,1033,169,355,260,30,45,721,906,88,44,44,418,417,218,2,2,846,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1628,721,400,239,728,1336,984,425,65,120,1232,463,463,640,349,616,616,2,2,2,2,2,2,2,2,147,147,2,2,2,2,2}, + {1628,286,541,530,1610,201,1220,1592,272,181,2,38,263,1586,1157,1157,1157,2,2,2,1157,1157,2,2,2,1157,2,2,2,2,1157,1157}, + {1531,621,210,755,482,82,1308,317,427,168,2,232,116,190,701,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,701}, + {1532,575,1245,360,249,630,133,1406,920,1539,63,63,76,82,82,2,2,2,770,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1102,785,118,93,1491,988,275,53,1328,26,2,2,240,647,240,761,761,761,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1478,722,113,1534,1751,115,1728,1234,777,282,508,508,1184,63,1184,855,855,2,2,2,738,738,578,578,2,2,2,2,2,2,2,2}, + {1480,536,1421,164,429,84,970,1673,548,497,2,2,530,156,156,128,245,2,2,2,2,260,2,2,2,2,2,2,2,2,2,2}, + {1533,1302,1286,538,619,526,1669,145,1034,125,2,1038,1038,388,388,387,729,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1482,961,1093,556,1746,628,427,689,510,751,684,37,37,1229,1256,882,1507,1507,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1472,852,533,433,924,57,53,1036,410,675,1212,1212,1212,600,600,1212,1259,1245,1245,1245,1245,2,2,2,2,2,2,2,2,2,2,2}, + {1228,425,1030,699,407,171,568,925,1104,97,2,2,1286,1286,1286,502,2,1219,1219,1219,1290,2,1219,1219,1219,1219,2,2,2,59,2,2}, + {1547,657,777,695,1254,224,933,367,212,385,2,2,2,1422,749,245,885,710,2,2,710,710,710,45,710,710,710,710,2,336,710,2}, + {1678,466,549,145,351,816,1041,334,192,192,2,2,348,1017,130,4,4,180,180,180,180,512,2,2,2,512,512,512,512,512,128,8}, + {1550,569,481,1041,1680,1114,1265,160,386,194,2,2,2,533,929,531,422,531,1355,1355,1355,1355,531,216,216,2,2,2,2,2,2,2}, + {1559,462,850,289,1570,71,512,858,810,835,2,2,2,2,1028,1205,1205,546,546,546,1205,1205,1205,1205,1205,2,1205,1205,1205,1205,1205,1205}, + {1192,888,701,164,131,613,282,237,525,366,2,2,2,2,1737,845,845,750,2,1062,1062,1062,1062,1062,1062,1062,1062,1062,1062,2,1261,1261}, + {1208,426,412,1072,274,248,1544,627,9,458,2,2,2,2,2,2,270,270,270,150,715,282,150,150,150,150,150,150,150,150,150,150}, + {1128,393,1522,96,160,581,540,120,441,176,2,2,2,2,2,2,1427,551,1102,1102,328,328,592,592,592,592,592,592,592,592,592,592}, + {1202,538,171,1177,1090,690,1566,746,1012,1012,2,2,2,2,313,781,808,313,1125,1117,930,1117,1117,1117,1117,1117,1117,1117,1117,2,2,2}, + {1567,1265,372,1633,613,484,243,1523,21,275,2,2,2,431,431,431,431,2,2,978,489,889,889,889,889,889,889,889,2,2,2,2}, + {1566,982,815,133,891,412,1179,831,651,268,2,2,2,367,366,367,367,63,63,767,2,2,2,2,2,2,2,2,2,2,2,2}, + {1522,1422,1017,124,499,451,731,1112,1355,1355,2,2,2,854,854,336,854,336,1297,2,2,2,193,193,193,193,193,2,2,2,2,2}, + {1160,1331,917,1696,401,547,122,592,863,863,2,2,703,703,703,703,495,495,495,2,2,495,495,495,495,495,269,2,2,2,269,269}, + {1538,814,1027,677,524,226,756,202,242,102,2,2,912,564,1289,682,2,1125,1125,1125,1125,2,1289,1289,1289,1125,1125,1125,2,1289,1289,1289}, + {1598,397,1471,1471,1162,866,236,948,1557,737,2,2,153,737,1408,765,765,608,2,2,2,171,608,608,608,608,2,608,608,2,2,2}, + {1598,434,107,270,148,1317,835,123,642,1236,2,2,67,633,771,878,771,878,878,2,2,2,771,2,2,2,2,2,2,2,2,2}, + {1628,1502,1042,822,80,403,1335,684,464,426,671,671,336,336,336,2,425,896,2,2,2,2,1337,1337,1337,1337,1337,1337,2,2,2,2}, + {1630,715,1368,1273,993,293,385,545,1267,896,1038,1038,270,1325,1325,2,2,961,961,961,961,961,961,2,2,961,961,2,2,961,2,961}, + {1612,723,409,641,796,1087,1228,1398,623,262,740,740,870,870,397,2,2,893,893,2,2,1367,328,2,328,2,2,2,2,2,2,2}, + {1614,588,652,105,441,844,734,912,532,878,1073,1073,62,1415,693,1431,1431,1431,1431,925,925,925,925,925,925,925,2,2,2,2,2,2}, + {1607,1503,1072,471,221,277,854,1236,263,752,2,694,1657,934,553,2,2,2,498,498,2,802,2,46,2,2,2,2,2,2,2,2}, + {1172,987,140,1964,584,600,852,1725,456,1199,718,718,791,981,791,2,2,2,2,2,1260,2,2,2,2,2,718,2,2,718,2,718}, + {1746,771,620,415,1057,437,613,1034,1662,837,2,1149,1466,1149,1149,1149,1466,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1263,835,1533,789,1259,174,1497,557,644,203,2,289,604,434,434,434,2,844,844,2,2,2,1111,1111,1111,2,2,2,2,2,2,2}, + {1272,884,388,1889,956,159,1172,595,219,645,2,629,107,107,1279,75,2,2,2,211,2,2,2,2,2,2,2,2,2,2,2,2}, + {1797,904,172,659,349,177,692,448,1141,990,640,99,1073,806,640,640,2,640,640,911,911,911,640,640,640,640,2,2,2,2,2,2}, + {1276,442,1008,1352,243,162,711,301,552,1002,668,668,384,71,384,384,2,2,2,2,2,727,727,727,777,777,777,777,777,777,2,777}, + {1600,1130,171,1113,813,722,117,990,37,24,969,94,825,1398,1398,1398,1398,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1198,496,714,609,644,1159,873,249,186,1539,136,239,379,1994,2,68,68,68,68,68,2,2,192,2,969,2,2,969,2,2,969,969}, + {1678,1316,460,1133,1003,150,1236,1316,1417,218,1763,1763,77,77,2,1491,771,771,771,771,771,2,771,2,2,2,2,2,2,2,2,2}, + {1682,449,1067,393,136,854,36,492,637,1053,247,1111,1111,1111,2,247,247,247,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1288,1690,702,760,420,333,1213,1911,805,351,67,67,1568,1568,2,2,604,142,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1288,1858,152,894,346,104,997,203,249,1006,1278,1489,1489,555,2,2,2,1074,1074,518,2,2,518,2,2,518,2,2,2,2,2,2}, + {1601,697,532,408,697,1140,1568,47,1499,780,1171,318,318,318,2,2,2,2,318,318,2,2,2,2,2,2,2,2,2,2,2,2}, + {1283,1078,791,873,655,412,389,835,292,958,1245,678,1611,1519,2,2,185,2,2,2,2,2,2,1245,1245,2,2,2,2,2,2,1245}, + {1685,1610,1447,1093,1255,937,703,431,522,1384,988,988,253,988,2,1892,1892,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1822,589,236,205,797,39,241,1048,181,386,102,102,102,111,1361,1361,1361,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1233,843,813,157,396,669,1531,439,640,733,996,996,996,1566,951,608,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1342,705,302,595,1200,52,83,647,519,139,103,103,103,513,2,513,2,2,2,2,2,513,2,2,2,2,2,2,2,2,2,2}, + {1630,1244,142,767,1299,719,629,1716,419,837,1145,1136,1148,1405,1405,1405,2,2,2,2,2,309,309,309,309,309,2,2,2,2,2,2}, + {1636,974,279,419,893,1608,1491,156,1486,115,730,730,863,509,924,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1318,1234,213,1089,1567,602,1330,404,467,718,249,215,354,177,59,332,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1732,1771,584,533,297,1056,669,293,146,311,1176,311,590,590,277,2,2,2,2,2,2,2,2,2,539,539,2,2,2,2,2,2}, + {1026,512,1196,394,1259,1313,762,549,311,1576,1576,465,465,140,465,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1360,383,1470,502,1010,153,1588,619,1246,396,1107,1107,112,423,423,2,2,2,2,2,202,2,2,2,2,2,2,2,2,2,2,2}, + {1320,1636,858,1210,509,194,1575,154,1424,455,1860,832,1075,581,262,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1670,1350,689,1074,437,956,587,642,1154,439,196,1108,1108,1108,990,2,2,2,2,2,1112,2,2,2,2,2,2,2,2,2,2,2}, + {1873,890,920,874,591,651,768,478,331,76,760,760,760,760,67,2,2,2,2,1241,1241,1241,1241,2,2,2,2,2,2,2,1241,1241}, + {1682,867,333,102,628,891,654,506,995,684,961,563,1313,1313,1313,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1672,1248,429,813,262,92,809,1248,560,1365,1392,753,753,1259,1261,2,2,2,2,2,2,2,2,177,177,2,2,2,2,2,2,2}, + {1391,1598,1112,590,797,584,1354,47,1473,1291,1874,48,491,463,990,463,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1875,1576,924,677,461,134,1525,1619,44,701,299,743,728,791,791,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,791,791}, + {1267,904,1187,1595,765,1451,494,1573,950,909,87,1265,757,1371,1005,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1360,1091,1478,1237,97,578,1616,494,1422,223,865,1092,359,2,1080,4,2,2,2,688,1965,2,1965,2,2,2,2,2,2,2,2,2}, + {1750,386,393,840,723,791,1707,1319,1525,83,1302,571,280,2,280,73,2,2,2,1207,2,2,2,2,2,2,2,2,2,2,2,2}, + {1763,1018,1859,432,717,723,874,1294,1050,1800,1237,619,1074,2,10,1237,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1376,652,461,225,361,936,1073,1279,149,619,983,511,1994,2,2,1076,1076,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1947,393,495,946,1375,391,2128,582,1143,695,1872,760,760,2,2,1456,974,974,435,974,974,435,974,2,974,974,2,2,2,2,2,2}, + {1768,1463,531,1008,95,1677,362,1105,985,177,1682,1682,244,2,2,1234,1041,1041,1041,2,2,2,1041,1041,2,2,2,2,2,2,1894,2}, + {1780,1739,1357,1684,1586,736,208,966,1691,339,339,128,128,2,2,128,128,128,2,2,128,2,2,2,2,1929,2,2,338,2,2,338}, + {1387,1459,358,1409,1919,917,777,223,313,1847,1012,1024,1024,2,2,2,2,1420,1420,1428,1420,2,1420,1420,2,2,2,1420,1117,1117,1117,1117}, + {1289,907,228,665,1695,1735,489,214,762,1777,321,1674,932,2,2,2,2,1358,709,2,1959,1959,372,2,2,372,372,2,2,372,372,372}, + {1378,680,1117,1367,759,62,319,563,505,1138,1093,345,693,2,2,2,780,780,2,2,2,729,729,729,2,2,2,2,2,2,729,729}, + {1802,1645,453,1079,604,618,334,855,541,167,37,88,849,2,2,518,518,2,2,530,2,2,2,2,2,2,2,119,119,2,2,2}, + {1275,1612,143,1586,502,987,555,436,2236,1826,494,494,358,2,2,213,2,2,2,2,2,2,1585,1585,1585,1585,1585,1585,1585,1585,1585,1585}, + {1322,512,560,432,365,87,1835,1137,515,1271,1739,309,309,1229,1229,1229,2,2,2,2,2,2,2,2,416,416,416,416,2,2,2,2}, + {1758,835,287,888,391,875,1834,516,1432,1171,98,408,302,976,976,1963,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1444,394,1613,796,645,1406,186,158,402,1364,314,588,606,2,577,117,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1337,1391,137,371,165,87,1026,20,419,99,572,572,918,854,918,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1323,589,526,1555,1636,1172,86,42,1545,57,627,1769,1769,2,867,343,2,2,2,2,2,2,2,724,2,2,2,2,724,724,2,2}, + {1323,1647,384,301,270,549,1098,1144,1066,55,88,1805,683,2,945,120,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1327,1075,539,1017,926,350,1102,236,494,1268,286,286,1293,267,227,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1472,661,1538,487,94,2209,563,138,881,1735,718,203,1382,1473,1473,1473,1473,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1413,766,349,1471,45,625,733,1082,170,58,1268,207,1081,1081,1081,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1861,1487,419,97,799,1791,458,1029,370,627,57,414,414,1540,247,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1881,716,268,387,2138,1212,999,408,1363,434,1429,1429,1648,1648,1007,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1480,1131,1089,1688,340,962,505,1816,139,44,1350,403,1385,1996,173,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1868,650,1146,1690,948,43,497,692,1628,1302,1302,108,462,731,731,2,2,2,2,2,2,2,185,185,185,2,2,2,2,2,2,2}, + {2023,1204,531,733,1054,618,668,363,783,218,1302,2055,559,2055,2055,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1328,601,601,617,554,467,391,1545,162,1361,807,1565,1565,243,1344,2,725,510,510,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1808,1525,1129,652,1195,329,1410,558,1322,911,161,536,737,94,306,2,2,2,2,2,2,2,2,2,541,541,541,2,2,2,2,2}, + {1911,1338,639,1106,854,128,19,1353,847,253,618,517,2054,2054,93,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1825,850,180,1483,864,953,50,81,106,432,1372,1372,1212,10,10,10,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1892,441,977,228,1252,604,735,136,889,878,1319,1319,2127,2127,1963,367,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1820,1553,536,1351,425,1268,227,1742,429,348,1397,552,1151,1151,2,180,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1453,1044,556,833,305,1493,989,1158,726,1790,532,1229,1229,1229,2,2,2,2,2,2,2,2,2,2,259,2,2,2,2,2,420,2}, + {2059,592,492,973,137,1331,392,334,635,1480,2254,1796,1796,284,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1460,986,709,268,755,824,83,893,115,656,2071,1323,1001,144,2,2,2,2,2,2,2,2,1527,1527,1527,1527,1527,1527,1527,2,801,801}, + {1850,1476,792,840,2037,229,1578,526,431,1485,1450,1001,1001,1001,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1922,1383,813,346,1247,666,1931,1111,2042,79,682,501,1349,1930,2,2,681,681,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1922,542,1739,625,88,1376,259,49,338,318,505,788,1314,657,2,2,2,1314,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1948,1530,576,582,1069,119,2131,41,1178,1677,1677,1677,325,346,2,2,2,2,2,2,1401,2,33,2,2,2,2,2,2,71,71,71}, + {1928,1111,168,1252,1467,1083,1927,603,1278,714,1027,50,751,1970,2,2,2,2,621,2,100,2,2,10,10,2,2,2,2,793,793,793}, + {1394,896,674,2350,1375,1599,1858,135,762,722,628,685,705,28,2,2,2,2,2,2,2,2,2,2,2,855,2,2,2,2,2,2}, + {1540,791,518,419,1130,1068,299,1386,1378,134,859,859,71,162,2,71,71,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2107,709,828,154,542,184,1094,1665,307,1549,177,2007,85,773,2,2,2,2,2,2,2,2,2,2,2,697,2,2,2,2,697,2}, + {1977,1218,244,365,576,666,761,238,629,913,1907,986,1351,986,704,1257,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1496,1912,1291,1053,510,2322,1048,1530,2223,673,894,594,628,332,2,2,2,2,2,295,295,295,2,2,2,2,2,2,2,2,2,2}, + {1520,1107,1082,687,484,1732,676,1595,467,653,1091,428,2113,332,332,332,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1905,612,920,848,562,2032,230,1305,1073,851,731,798,798,357,516,2,2,2,2,2,2,2,2,1465,1465,373,2,2,2,2,2,2}, + {1428,1062,1016,75,297,1130,533,768,464,753,48,1510,1510,418,375,1626,2,221,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1396,729,1710,337,371,489,1341,2117,132,1870,853,853,408,1079,328,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1978,1051,977,588,1423,1001,508,409,825,497,659,1063,384,463,463,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1534,854,2007,1207,947,1773,1571,1505,909,1471,1655,1655,2334,1327,409,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2157,2106,679,238,378,49,1101,588,811,1313,1556,2301,475,812,812,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2187,1515,549,1416,1073,1613,47,1046,390,252,1214,1404,1404,933,1013,2,2,2,1025,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2145,1069,662,709,737,1141,1737,827,1384,1628,107,107,1032,277,277,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2173,1379,155,393,1578,610,1911,899,697,58,185,597,597,1249,1369,2,2,2,2,1369,2,2,2,2,2,2,2,2,2,2,2,2}, + {1413,1589,1603,2268,520,333,1416,859,1619,867,1154,512,1291,413,413,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1600,1823,1698,1268,623,583,1932,1674,522,529,1862,1281,246,989,246,2,246,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1559,992,174,1313,612,1487,1487,461,702,37,1660,839,2,95,1628,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2066,1719,710,1294,2041,377,1005,690,132,756,1618,187,187,726,187,615,615,2,2,2,2,851,2,2,2,2,2,2,2,2,2,744}, + {2192,1029,310,1609,592,1542,265,117,2006,82,162,205,2,2009,2009,1201,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1570,1504,1414,1143,1999,1932,1015,1015,556,514,626,79,2,79,1795,1461,1461,2,2,2,2,2,2,2,1461,1461,1461,1461,1461,2,2,2}, + {1562,937,1964,934,1349,378,459,109,1676,1655,1339,1809,2,768,768,188,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1965,949,1057,1043,2256,1571,970,348,69,1324,1174,485,105,105,105,2172,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2044,1869,838,1424,1097,155,1142,230,1335,420,235,1510,2,431,425,622,2,2,2,2,2,625,2,2,2,625,625,2,2,2,2,2}, + {1976,1433,820,504,421,1007,388,1083,635,82,1524,750,2,2,870,106,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1990,1948,1138,1787,253,115,312,1912,341,1624,260,1783,1315,1315,790,790,790,790,790,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1993,585,327,1393,1013,1671,1758,1436,1989,1217,1109,1476,2,2,1042,756,1042,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2056,1062,1605,1943,680,445,113,857,650,1388,2016,1231,2,2,1292,1292,1292,2,1039,1039,1039,1039,1039,2,2,2,2,2,2,2,2,2}, + {2008,1773,416,1954,1314,742,1694,505,202,1747,785,375,2,2,2,477,1538,477,2,2,2,2,2,1309,1309,1309,1309,2,2,2,1309,2}, + {1658,1008,258,749,427,1071,2052,263,1047,2152,1602,1602,2,2,2,1311,669,669,2,1897,1897,1897,669,669,669,669,669,669,669,669,669,669}, + {2258,1887,1875,1021,863,604,543,1115,509,1243,312,213,2,2,2,2,335,770,770,2,1143,567,2,2,567,567,567,411,2,2,2,411}, + {2266,1872,991,1468,1168,939,907,833,624,701,386,1713,2,2,2,2,2,931,861,381,1299,2,861,2,2,2,861,861,861,861,861,2}, + {2273,1510,803,2278,842,1245,1389,230,822,1564,113,1276,2,2,2,2,1350,273,273,2,2,2,2,2,1281,1281,1281,2,2,1281,1281,1281}, + {2278,1028,548,373,190,1443,614,2386,1940,930,557,2069,2,2,2,558,112,103,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2108,776,1568,342,2215,1882,681,1292,1601,586,1481,618,1930,1930,1930,1930,2146,89,89,2,2,2,2,2,2,2,1171,2,2,2,2,2}, + {2139,2177,1652,392,715,605,778,632,472,1619,64,64,2,2,2,1747,859,2,2,2,2,2,216,216,216,216,1747,1747,1747,1747,1747,1747}, + {1492,448,271,135,1288,417,130,83,235,2313,482,746,2,2,746,609,611,611,611,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1628,846,1504,138,464,401,501,506,967,1027,1540,1035,2,1921,1539,1539,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1630,1677,1624,301,1038,909,887,374,411,143,1021,174,2,1393,19,634,2,2,2,2,2,2,2,873,2,2,873,873,2,2,2,2}, + {1654,1131,2054,994,2170,548,801,252,87,219,488,2239,2,1232,1839,1822,2,2,2,968,2,2,2,2,2,2,2,2,2,2,2,2}, + {2065,1520,1423,1797,899,1425,1801,776,2365,58,646,695,2,998,998,1342,2,2,2,2,2,2,2,2,2,2,2,2,1150,1150,2,2}, + {2304,1948,316,1063,237,607,1143,2575,1388,1022,127,251,2,438,1570,1570,1570,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2177,710,1912,617,809,1078,199,905,673,519,457,52,2,1348,1348,410,2,2,2,2,2,340,2,2,2,2,2,2,2,2,2,2}, + {2073,1543,1586,1296,2466,753,455,46,119,1694,2035,1592,206,206,206,2,2,2,2,2,2,2,2,1172,2,2,2,2,2,2,2,2}, + {2075,1056,874,2101,566,1790,1333,386,538,1560,2254,331,717,717,717,454,454,2,2,2,2,2,2,2,2,2,2,2,454,454,2,2}, + {1670,977,1540,553,855,1729,239,757,191,62,732,549,1092,1092,199,199,199,199,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2096,1155,2257,125,1986,245,1208,2146,2287,680,1413,73,467,1410,1410,2,2,2,2,2,133,133,133,2,2,2,2,2,2,2,2,2}, + {1538,1026,2157,1457,1784,2559,184,29,614,273,697,697,1922,697,697,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2106,856,1025,382,389,272,425,672,1021,216,601,292,510,510,876,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1662,608,2478,266,1330,505,40,2058,964,724,596,1221,1221,310,42,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1600,1338,196,1510,1371,1138,957,169,545,1176,1131,2460,1708,541,541,2,363,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1682,1008,737,444,822,999,2066,283,646,1860,1008,778,1178,1178,458,1743,2,2,2,2,2,2,2,2,2,2,2,2,2,1743,1743,1743}, + {2132,756,1097,166,202,411,640,717,514,1389,633,633,633,633,633,633,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2386,748,620,478,647,898,320,53,1115,190,60,1860,1860,802,802,2,2,2,2,1264,1346,1346,2,2,2,2,2,2,2,2,2,2}, + {2125,996,1081,124,1140,628,1668,1913,151,2495,523,430,260,708,2190,2190,2190,2,2,2,2,2,1660,2,2,497,497,497,497,497,497,2}, + {1602,1489,895,383,56,698,2081,1728,794,789,16,16,797,302,52,2,2,2,2,2,2,2,2,2,2,797,797,797,797,797,797,1808}, + {2210,606,901,547,131,1924,1852,1271,194,766,390,390,520,795,1429,1429,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1731,599,817,724,718,1038,1082,2503,1341,936,421,1802,1304,1304,1491,186,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1614,1058,847,689,749,1028,1047,1474,117,1369,1442,1442,1540,700,104,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1733,679,2041,2420,326,934,1172,1431,193,370,1073,1073,1073,260,2,2,2,2,2,2,2,2,2,2,2,2,1193,2,2,2,2,2}, + {2168,1532,769,2570,1303,357,1793,1633,1226,1025,205,1218,1984,764,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2234,1706,356,581,532,933,1704,387,1345,1345,34,135,350,307,614,614,307,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1656,2093,354,310,306,1553,106,459,175,55,1482,958,254,254,2,356,356,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1654,1035,330,533,1446,953,499,142,1527,1748,265,1437,265,510,2,2,2,2,2,2,2,1835,1835,1835,1835,2,2,2,2,2,2,2}, + {1600,479,1457,246,2025,618,1612,2139,169,1492,1097,1327,2007,2007,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1744,447,845,2145,748,1555,1193,1312,916,1770,1294,546,794,323,2,2,2,2,2,1733,1733,2,2,1730,2,1733,1733,2,2,1733,551,551}, + {1766,1558,1901,1393,987,1859,815,1165,50,2065,88,88,1453,1453,2,2,2,995,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1615,1267,1242,1494,399,663,68,1209,1573,528,640,1200,248,640,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1678,592,1351,509,312,721,163,1597,1262,199,2643,1330,1661,992,2,2,719,2,2,2,2,2,2,2,2,2,2,2,2,2,1704,2}, + {2207,970,838,2043,1016,561,267,329,584,608,679,303,832,1613,959,959,959,1409,1409,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2217,352,447,914,1200,561,614,1616,509,2292,1114,1114,1229,52,1053,1053,1053,2,2,2,2,2,2,2,2,2,2,2,2,2,795,795}, + {2313,595,1593,1951,133,282,372,2396,1117,226,2104,267,374,267,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2312,1231,1604,997,652,1096,1070,320,481,662,911,1610,342,2527,606,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2245,1541,1828,783,615,428,1282,1892,848,1219,2465,314,314,314,2,2,2,2,2,2,2,2,2,2,2,2,1323,2,2,2,2,1323}, + {2522,1030,324,1264,628,1339,480,234,2351,1085,1979,2333,1339,1356,1356,2286,2,2,2,2,2,2,2,2,2,2,2,2,2530,2,2,2}, + {2519,1136,612,209,994,1179,1060,2621,130,485,661,1444,2122,124,258,1114,2,2,806,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2524,1894,253,2072,1242,355,888,1362,28,480,452,1216,595,545,354,1145,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2569,1356,1053,410,437,58,1508,831,2272,383,1725,615,1191,1191,1191,2493,186,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2386,1106,709,251,784,929,1551,2481,304,2148,1546,955,2453,866,866,2,2,2264,2264,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2228,1163,1995,649,1000,680,325,1591,774,767,711,711,1418,524,711,401,976,2,2,2,2,2005,2005,2,2,2,2,2,2,1390,1390,2}, + {2362,1706,564,1088,1296,1267,70,1015,496,1298,758,154,240,240,154,154,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1715,2260,357,557,783,1195,2288,1997,1120,144,247,175,1277,203,203,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2294,2360,1353,748,1439,226,940,2316,1112,1527,214,1406,1429,712,1124,2,595,595,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2303,1018,316,280,1616,909,97,1126,1295,736,216,54,2045,726,1673,2,2,2,2,2,779,779,2,2,2,2,2,2,2,2,2,2}, + {2390,491,1217,1148,2314,2250,2180,308,613,662,1346,1346,1346,1280,778,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1732,527,1303,664,71,294,404,917,1074,180,2618,2412,441,1987,1750,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1695,1287,1346,1181,1412,1653,830,2025,957,1720,1614,887,964,964,964,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1756,2308,1986,101,957,633,1940,1002,390,1237,95,1441,95,95,705,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2461,1412,540,1183,229,300,47,585,518,402,1863,1863,560,1326,1326,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1707,717,366,287,1883,50,599,1371,474,1551,947,2142,1885,947,2008,1004,1004,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2632,567,1149,1227,1156,2052,643,1585,1197,581,63,718,699,149,149,1940,2,2,2,2,2,2,2,2,2,2,2,2146,2,2,2,2}, + {1773,2024,377,340,1938,103,1180,600,199,848,2449,2449,506,506,762,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2435,1920,394,1482,266,1637,911,1697,1689,1249,1085,1085,397,2292,1355,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2345,662,270,324,1061,1080,1952,593,1480,2111,2667,2093,2059,2120,955,1447,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1777,455,1487,1190,455,1542,977,2308,437,1129,410,856,1420,412,412,766,2,2034,2034,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2662,2224,1142,656,59,598,730,458,226,1151,741,1286,1015,2,688,2017,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2666,768,529,990,2329,130,1678,2466,318,1083,387,1524,511,2,731,731,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2501,1216,246,1278,718,704,2019,88,273,1203,67,1488,1828,2,2,1489,1489,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2474,2292,1818,2061,2833,751,2172,1708,1210,1675,370,131,163,2,2,163,163,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1751,1575,889,828,82,1956,712,499,1420,1686,339,2326,2035,2,2,558,558,2,1234,2,2,2,2,2,2,2,2,2,2,2,1239,1239}, + {2522,1148,1943,168,218,252,543,1535,2004,130,353,353,42,2,2,2,1173,1173,2,1547,2,2,2,2,2,2,2,1547,1173,1547,1547,2}, + {2695,432,1213,579,865,1637,1857,84,447,155,2492,347,1980,2,2,2,1155,1155,1155,2,1933,1933,1933,2,2,2,2,2,2,2,1901,1901}, + {1808,1683,474,1761,106,602,1416,217,1351,1602,366,393,1966,2,2,2,2,2,378,378,606,606,606,2,2,2,2,2,919,919,919,919}, + {2428,1576,1692,449,2012,240,1167,418,272,1557,2197,645,645,2,2,2,2,2,2150,2150,2,2,562,715,2,2,2,81,81,2,2,2}, + {2727,781,1689,1709,997,2563,1032,468,44,992,1214,725,75,2,2,2,2,360,360,380,2,2,2,2,2,2,2,2,2,2,2,2}, + {1948,1085,1344,2090,1435,2389,3193,1007,1003,244,667,1838,2062,2,2,2,1802,299,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2433,932,689,818,2014,1498,749,1645,867,1627,47,1766,2193,2,2,2030,2030,2,430,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2463,712,1525,2092,2942,352,761,242,2178,2339,483,1905,1347,2,2,65,529,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2547,920,386,925,74,579,323,2319,520,2332,1535,751,1591,2,770,770,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2452,2588,2055,665,818,2622,413,1260,965,211,989,1219,166,2,1251,1251,2,1256,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1957,2311,993,276,293,2826,1087,880,927,1811,1122,2974,2974,2,2,590,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2552,998,533,827,1619,831,1861,918,750,1955,241,1899,448,2151,2151,449,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1860,579,1000,1575,898,170,185,1032,293,2754,438,459,459,2,1199,1199,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2602,2417,1888,2528,1410,669,1543,233,814,2478,225,1449,1449,224,1671,1671,2,2,931,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1889,2527,1366,1371,387,925,1751,162,250,1064,292,467,467,546,1244,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2492,1186,1350,1616,2749,1962,33,708,279,813,1390,489,1203,268,173,2410,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2500,1575,423,541,561,380,262,1564,1923,1242,2084,1758,1283,2213,924,924,2,2,2,2,2,2,2,2,2,2,2,1827,1827,2,2,2}, + {1842,1736,489,743,1539,1681,683,1412,1418,312,2778,2778,1975,1975,803,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2824,1183,2201,278,241,2230,1591,1648,1036,818,1321,1312,754,813,813,813,813,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1900,2506,952,1059,163,870,681,1235,1271,1188,2071,1705,1183,648,404,2,2,2,2,2,2,2,2,2236,2236,2,2,2,2,2,2,2}, + {2662,1443,2327,132,490,1149,1572,744,429,621,1763,2383,1903,1246,964,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2673,2182,1307,1776,1233,1828,1828,340,249,216,503,160,160,582,926,2129,2129,2129,2129,2129,2129,2129,2129,2129,2129,2129,2129,2,1018,1018,1103,1103}, + {2042,620,1074,2057,2758,859,815,1127,766,1693,252,808,981,416,416,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2102,881,2170,1673,705,101,58,1712,1568,214,758,488,1007,269,243,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2692,2665,961,1478,324,429,1311,376,1648,130,2083,1047,409,343,343,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2122,1087,563,1669,647,2996,151,2458,250,310,71,1348,355,965,2815,1333,1333,2,2,2,2,2218,2,2,2,2,2,2,2,2,2,2}, + {1952,1968,2260,2945,2464,1055,2626,570,1316,1828,1828,970,970,221,220,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2072,1947,1779,254,2822,1552,855,804,3452,202,695,82,684,208,1270,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1947,1699,1341,486,1765,1960,264,899,1082,1674,987,1878,930,1008,930,930,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1953,1527,1643,591,1517,2427,1232,1555,2542,495,675,2534,2534,3106,83,3106,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2744,1728,2213,792,761,1667,1908,31,447,442,815,2865,762,762,762,762,2,2,2,2,2,2,2,2,2,2,2,2,649,649,649,2}, + {2722,1406,1257,807,2191,3017,1330,1023,602,2124,794,530,733,733,1083,2528,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {1963,1525,437,398,609,393,2420,3059,435,1251,1977,1672,450,1960,1954,1960,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2626,2468,2838,845,2060,218,1080,912,911,1973,1365,920,1316,1316,2,1316,1316,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2734,1727,1743,1026,809,1154,779,244,1238,1616,812,784,825,1810,1810,1810,1810,1559,1559,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2180,2262,1651,204,3193,2121,2725,1016,629,1834,603,2848,26,26,728,728,728,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2648,1328,2578,133,1377,105,2485,2139,323,1045,145,761,1201,1848,2,814,814,814,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2999,358,250,1379,102,2349,1491,2074,42,376,2811,1220,296,296,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2810,1274,499,742,1724,425,190,1561,1302,2603,2255,917,661,661,2,2,2,495,2,2,2,2,2,2,2,2,2,2575,2,2,2,2}, + {2150,589,876,1616,2655,432,902,1028,433,1375,574,1400,1400,1400,2,2,2,2,2,1529,1529,1529,1529,1529,1529,1529,1529,1529,1529,1529,1529,1529}, + {1665,1856,201,824,796,249,1217,590,1375,1175,1599,824,824,3319,2,2,2,601,1961,1961,2,2,2,1961,2,2,2,2,2,2,1961,2}, + {2704,2239,1260,140,2161,2781,1840,574,2353,343,3218,61,2108,2038,1873,2,1833,1408,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2173,876,802,2197,3338,176,1783,224,1763,1160,1264,1264,2864,554,2,552,552,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2708,1663,2279,824,836,1598,2101,1620,1202,1606,1368,1079,1167,1999,2848,2848,2848,1101,1101,1101,2,2,2,662,2,2,2,2272,2,2,2,2}, + {1987,1463,2328,1890,1443,2086,283,2895,522,1577,1514,1657,2605,891,2,1181,1181,2,2,2121,2,2,2,2,2,2,2,2,2,2,2,2}, + {2173,1637,1139,905,1802,1378,296,439,1507,1017,1427,209,708,462,1508,1508,1508,2,2,2,2,2,2,2,2,2,2,2240,2240,1459,1459,1459}, + {2206,1526,628,2877,802,2587,1253,1258,1044,2195,3246,40,2898,2898,1704,598,2,2145,2,2,2,2145,2,2,2,2,2,2,2,2,2,2}, + {2182,618,1022,1433,1138,1580,2590,149,796,2090,743,294,294,1117,720,3003,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2025,1805,1466,1213,2006,1903,568,1700,1355,865,1783,1006,1006,1070,1070,268,2,2,2,2,2388,2388,845,845,845,2,2,2,2,2,2,2}, + {2185,1038,3050,1461,2270,2159,958,1637,233,2483,525,987,437,437,437,3065,2,2,2160,2160,2,2,2,2,2,2,2,2,2,2,2160,2160}, + {2083,1465,847,1450,502,447,2168,794,1761,1324,162,188,2853,2853,636,973,2,563,2,2,2,2,2,2,2,2089,2089,2089,2089,2089,2,2}, + {2923,2303,203,508,472,648,3169,269,515,3147,2415,1700,1700,1700,1461,1461,1461,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2300,1116,1555,2794,1095,998,1999,894,963,753,324,2130,2675,2675,554,2045,2,2,2,2,2,2,2,2130,2130,2130,2,2,2,2,2,2}, + {2103,768,702,1548,1486,2228,2846,861,665,1497,1046,1046,2252,394,394,1901,1155,2,2,2,2,2,2,2,2,2,2,2,192,192,192,192}, + {2923,640,661,2179,1207,182,872,171,738,269,1372,222,908,2069,2069,2,1550,516,2,2,2,2,2,2,2,2,2,2,1109,2,2,2}, + {2833,2005,387,733,562,468,317,224,94,478,1606,2522,1606,2001,1087,2,2,1087,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2125,2479,1749,1226,1169,1681,459,652,1087,2211,1613,686,2213,1689,2446,2,2,2925,2925,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2953,1059,205,3093,138,132,2148,1345,1499,216,151,1296,2446,1610,1632,2,2,2,2,4,4,2,2,2,987,987,2,2,2,2,2,2}, + {3199,1431,593,2050,2785,507,1540,1103,1740,459,62,1766,1781,1121,1600,2,1600,1600,125,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2258,1714,415,373,1919,2605,693,827,1918,496,1479,1903,86,1083,415,2,2,38,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {3289,2032,329,2169,2323,1599,517,1704,1847,804,632,40,40,40,40,40,40,40,40,2,2,1600,2,2,1600,2,2,2,2,2,2,2}, + {2165,2725,2293,368,705,3063,494,103,12,1332,175,2331,3144,2165,1709,1709,2090,2,2,2,2,1363,1363,2,566,2,2,2,2,2,2,2}, + {2300,1070,2169,2540,734,1002,912,1386,2215,224,1285,880,2052,2052,1301,959,563,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {3267,1852,1037,648,611,1250,432,853,1467,179,715,2,2033,841,2607,2607,2607,2607,2,2,2,2,1874,1874,2,1874,899,2,2,2,2,2}, + {2348,2565,794,859,1740,1596,532,462,457,1014,1227,2,2761,954,249,249,2,458,458,2,2,2,2,2,909,909,2,2,2,2,2,2}, + {3038,2399,1450,1276,1222,727,552,646,1055,2351,686,63,252,504,3166,1802,2,2,1165,1165,1165,2,2,2,2,1165,1165,1165,2,2,1165,2}, + {3038,2519,1494,107,2597,802,535,1669,1695,1928,1940,1580,1580,85,2274,1551,2,2,2431,560,560,560,2,2,1098,2,2,2,2,2,2,2}, + {3040,1044,1927,1952,1479,3124,1373,1990,588,2550,1277,2,629,2671,1842,2712,840,1702,2,1669,2,1347,2,2,2,2,1669,1669,1669,1669,2,1669}, + {3056,1567,691,1243,653,751,248,842,1954,480,458,2,2,2451,934,3172,3556,2259,2312,2,2562,2562,2,2,2562,2562,2562,2562,2562,2,2,2}, + {2959,2553,1333,877,2492,3169,2498,686,2030,2820,3233,1313,1313,1471,1471,1471,1471,2,2,1471,1471,2,2,1481,2,1887,2,2,2,2,2,2}, + {3398,964,862,301,1705,2002,310,644,144,1091,1507,2,2,2460,496,496,2517,2517,1842,2,2,1964,2,2,2,2,2,1676,2,2,2,2}, + {2379,3034,166,302,2108,1078,2976,68,158,134,1567,2,2,1514,1514,1514,1883,1883,2,2,1883,1883,1883,1883,1883,1883,1883,1883,1883,2,2,2}, + {2386,1270,1204,1032,1474,224,496,2296,1536,1219,311,2,2,2,2,1238,2108,2108,2108,2108,2108,2108,2108,2108,1444,1444,1444,1444,1444,1444,1444,2}, + {2431,739,2488,1386,1632,2107,2602,2139,1751,349,3147,2,2,64,16,8,32,4,4,32,728,728,728,728,2,2,64,16,8,180,180,180}, + {3405,2142,1621,110,2112,2097,807,740,747,282,372,2,2,2,2,2493,2493,2493,1299,2,132,1872,2,1843,2,2,2,2,2,2,2,2}, + {3157,1230,685,1513,663,1335,2100,1441,1826,1670,1539,2,2,2,2899,2899,1378,54,2,46,46,2,2,1362,1362,2,2,2,2,2,2,2}, + {2415,822,3658,449,1980,891,129,823,1787,621,514,2668,2668,2668,2668,2668,666,269,2830,2,2,2,2,241,370,370,370,370,2,2,2,2}, + {2463,2664,2825,1208,882,629,428,428,356,343,1730,2,769,769,769,1714,769,2,2,955,769,2,2,955,955,955,2,2,2,955,955,955}, + {2447,1588,1077,831,1413,2362,1499,1812,1112,815,129,1034,1034,1867,194,518,1454,723,723,1251,2,160,2,2,1251,1251,2,2,2,2,2,2}, + {3094,1638,1514,843,1503,1884,1481,727,723,1319,226,2,676,2401,1699,562,639,639,1176,2,2,2,2,824,2,2,2,2,2,2,2,2}, + {3125,2004,547,2986,2919,471,948,1747,201,1862,802,2,1238,1277,1277,1277,2,2,1245,1245,1245,2,2,2743,1245,1245,2,2,2,2,2,2}, + {2582,2469,533,1726,1575,1505,2448,2031,1257,427,588,1633,202,3553,1938,672,195,195,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2378,636,1958,1628,1255,2285,2208,1626,719,2944,1086,1436,1436,1719,2111,655,2637,2637,2,2,2,2637,2637,2,2637,2637,2637,2637,2637,2,2637,2637}, + {2372,3079,2161,515,368,847,955,1257,1937,315,2666,1938,1723,1252,1252,362,362,2,2205,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2518,2060,1055,362,1455,1899,1105,1560,2237,2451,2080,181,2346,181,1829,1829,1829,2,2,1509,1509,1509,2,1509,2,2,2,2,2,2,2,2}, + {3580,1671,674,1838,814,1409,323,3021,1047,2579,2579,2968,2968,102,2656,2638,2638,4006,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {3194,1576,1084,859,2879,1600,953,1429,471,867,1105,1490,293,293,293,2,2,198,2619,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {3215,2004,3333,2271,3283,1660,2135,1696,1413,1362,834,253,253,253,3802,2,2,2,1881,690,690,2,2,2,1881,1881,1881,1881,1881,2,1881,1881}, + {3719,2441,2094,1665,1707,1827,1310,230,1635,143,386,1029,1070,1062,1062,2,1062,1062,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {3249,1309,1232,472,711,2557,1479,1027,145,489,1377,2928,2928,3522,3522,3522,968,415,415,2,2,2,2,1332,1332,1332,2,1332,2891,2,1332,2891}, + {2462,1962,257,2244,1966,1905,204,262,799,319,752,1696,971,971,3781,1426,1426,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {3434,3131,1399,3413,1533,281,3288,1242,810,135,2506,2506,1742,946,1015,1044,1044,1044,2,2,2,2,1044,1837,1837,1837,1837,1837,2,2,2,2}, + {2518,1200,631,596,1946,365,2960,413,592,3878,242,2714,2364,1402,1402,2322,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {3362,2012,1759,2002,1365,150,3120,471,1590,3246,1296,196,196,196,2984,2323,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {3382,899,3140,2860,1155,1840,2822,355,1753,1856,1018,822,52,52,52,1102,1102,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2728,1334,274,1330,2674,2614,931,2250,883,1506,2193,1345,1089,500,2,219,390,2,2,2,2,2,2,2,390,2,2,2,2,2,2,2}, + {3911,3343,202,675,1733,71,166,176,1323,2864,899,2155,1108,2172,2,2,1829,2172,1107,2,2,2,2,1107,1107,1107,2,2,2,1107,2,2}, + {2757,3466,1411,1168,340,2760,1053,524,53,2090,1227,26,260,830,2,2,2,1139,2,2,2,2,2,2,2,2,2,2,2,2,2,2}, + {2662,902,2371,1920,1097,1476,1008,1012,3556,468,3374,2560,591,1446,2,298,298,149,149,149,149,149,3135,3135,3135,3135,3135,2,2,2,2,2}, + {2861,1407,1848,245,2186,1209,164,2577,625,132,657,2333,2333,2213,2213,2213,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2} +}; diff --git a/src/divonne/Rule.c b/src/divonne/Rule.c new file mode 100644 index 0000000..f5d9337 --- /dev/null +++ b/src/divonne/Rule.c @@ -0,0 +1,685 @@ +/* + Rule.c + integration with cubature rules + code lifted with minor modifications from DCUHRE + by J. Berntsen, T. Espelid, and A. Genz + this file is part of Divonne + last modified 6 Nov 11 th +*/ + + +enum { nrules = 5 }; + +#define TYPEDEFSET \ + typedef struct { \ + count n; \ + real weight[nrules], scale[nrules], norm[nrules]; \ + real gen[NDIM]; \ + } Set + +/*********************************************************************/ + +static void Rule13Alloc(This *t) +{ + static creal w[][nrules] = { + { .00844923090033615, .3213775489050763, .3372900883288987, + -.8264123822525677, .6539094339575232 }, + { .023771474018994404, -.1767341636743844, -.1644903060344491, + .306583861409436, -.2041614154424632}, + { .02940016170142405, .07347600537466073, .07707849911634623, + .002389292538329435, -.174698151579499 }, + { .006644436465817374, -.03638022004364754, -.03804478358506311, + -.1343024157997222, .03937939671417803 }, + { .0042536044255016, .021252979220987123, .02223559940380806, + .08833366840533902, .006974520545933992 }, + { 0, .1460984204026913, .1480693879765931, + 0, 0 }, + { .0040664827465935255, .017476132861520992, 4.467143702185815e-6, + .0009786283074168292, .0066677021717782585 }, + { .03362231646315497, .1444954045641582, .150894476707413, + -.1319227889147519, .05512960621544304 }, + { .033200804136503725, .0001307687976001325, 3.6472001075162155e-5, + .00799001220015063, .05443846381278608 }, + { .014093686924979677, .0005380992313941161, .000577719899901388, + .0033917470797606257, .02310903863953934 }, + { .000977069770327625, .0001042259576889814, .0001041757313688177, + .0022949157182832643, .01506937747477189 }, + { .007531996943580376, -.001401152865045733, -.001452822267047819, + -.01358584986119197, -.060570216489018905 }, + { .02577183086722915, .008041788181514763, .008338339968783704, + .04025866859057809, .04225737654686337}, + { .015625, -.1420416552759383, -.147279632923196, + .003760268580063992, .02561989142123099 } + }; + + static creal g[] = { + .12585646717265545, .3506966822267133, + .4795480315809981, .4978005239276064, + .25, .07972723291487795, + .1904495567970094, .3291384627633596, + .43807365825146577, .499121592026599, + .4895111329084231, .32461421628226944, + .43637106005656195, .1791307322940614, + .2833333333333333, .1038888888888889 }; + + enum { nsets = 14, ndim = 2 }; + + TYPEDEFSET; + + count n, r; + Set *first, *last, *s, *x; + + Alloc(first, nsets); + Clear(first, nsets); + + last = first; + n = last->n = 1; + Copy(last->weight, w[0], nrules); + + ++last; + n += last->n = 2*ndim; + Copy(last->weight, w[1], nrules); + last->gen[0] = g[0]; + + ++last; + n += last->n = 2*ndim; + Copy(last->weight, w[2], nrules); + last->gen[0] = g[1]; + + ++last; + n += last->n = 2*ndim; + Copy(last->weight, w[3], nrules); + last->gen[0] = g[2]; + + ++last; + n += last->n = 2*ndim; + Copy(last->weight, w[4], nrules); + last->gen[0] = g[3]; + + ++last; + n += last->n = 2*ndim; + Copy(last->weight, w[5], nrules); + last->gen[0] = g[4]; + + ++last; + n += last->n = 2*ndim*(ndim - 1); + Copy(last->weight, w[6], nrules); + last->gen[0] = g[5]; + last->gen[1] = g[5]; + + ++last; + n += last->n = 2*ndim*(ndim - 1); + Copy(last->weight, w[7], nrules); + last->gen[0] = g[6]; + last->gen[1] = g[6]; + + ++last; + n += last->n = 2*ndim*(ndim - 1); + Copy(last->weight, w[8], nrules); + last->gen[0] = g[7]; + last->gen[1] = g[7]; + + ++last; + n += last->n = 2*ndim*(ndim - 1); + Copy(last->weight, w[9], nrules); + last->gen[0] = g[8]; + last->gen[1] = g[8]; + + ++last; + n += last->n = 2*ndim*(ndim - 1); + Copy(last->weight, w[10], nrules); + last->gen[0] = g[9]; + last->gen[1] = g[9]; + + ++last; + n += last->n = 4*ndim*(ndim - 1); + Copy(last->weight, w[11], nrules); + last->gen[0] = g[10]; + last->gen[1] = g[11]; + + ++last; + n += last->n = 4*ndim*(ndim - 1); + Copy(last->weight, w[12], nrules); + last->gen[0] = g[12]; + last->gen[1] = g[13]; + + ++last; + n += last->n = 4*ndim*(ndim - 1); + Copy(last->weight, w[13], nrules); + last->gen[0] = g[14]; + last->gen[1] = g[15]; + + t->rule13.first = first; + t->rule13.last = last; + t->rule13.errcoeff[0] = 10; + t->rule13.errcoeff[1] = 1; + t->rule13.errcoeff[2] = 5; + t->rule13.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( x = first; x <= last; ++x ) + sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]); + s->scale[r] = scale; + s->norm[r] = 1/sum; + } +} + +/*********************************************************************/ + +static void Rule11Alloc(This *t) +{ + static creal w[][nrules] = { + { .0009903847688882167, 1.715006248224684, 1.936014978949526, + .517082819560576, 2.05440450381852 }, + { .0084964717409851, -.3755893815889209, -.3673449403754268, + .01445269144914044, .013777599884901202 }, + { .00013587331735072814, .1488632145140549, .02929778657898176, + -.3601489663995932, -.576806291790441 }, + { .022982920777660364, -.2497046640620823, -.1151883520260315, + .3628307003418485, .03726835047700328 }, + { .004202649722286289, .1792501419135204, .05086658220872218, + .007148802650872729, .0068148789397772195 }, + { .0012671889041675774, .0034461267589738897, .04453911087786469, + -.09222852896022966, .057231697338518496 }, + { .0002109560854981544, -.005140483185555825, -.022878282571259, + .01719339732471725, -.044930187438112855 }, + { .016830857056410086, .006536017839876424, .02908926216345833, + -.102141653746035, .027292365738663484 }, + { .00021876823557504823, -.00065134549392297, -.002898884350669207, + -.007504397861080493, .000354747395055699 }, + { .009690420479796819, -.006304672433547204, -.028059634133074954, + .01648362537726711, .01571366799739551 }, + { .030773311284628138, .01266959399788263, .05638741361145884, + .05234610158469334, .049900992192785674 }, + { .0084974310856038, -.005454241018647931, -.02427469611942451, + .014454323316130661, .0137791555266677 }, + { .0017749535291258914, .004826995274768427, .021483070341828822, + .003019236275367777, .0028782064230998723 } + }; + + static creal g[] = { + .095, .25, + .375, .4, + .4975, .49936724991757, + .38968518428362114, .49998494965443835, + .3951318612385894, .22016983438253684, + .4774686911397297, .2189239229503431, + .4830546566815374, .2288552938881567 }; + + enum { nsets = 13, ndim = 3 }; + + TYPEDEFSET; + + count n, r; + Set *first, *last, *s, *x; + + Alloc(first, nsets); + Clear(first, nsets); + + last = first; + n = last->n = 1; + Copy(last->weight, w[0], nrules); + + ++last; + n += last->n = 2*ndim; + Copy(last->weight, w[1], nrules); + last->gen[0] = g[0]; + + ++last; + n += last->n = 2*ndim; + Copy(last->weight, w[2], nrules); + last->gen[0] = g[1]; + + ++last; + n += last->n = 2*ndim; + Copy(last->weight, w[3], nrules); + last->gen[0] = g[2]; + + ++last; + n += last->n = 2*ndim; + Copy(last->weight, w[4], nrules); + last->gen[0] = g[3]; + + ++last; + n += last->n = 2*ndim; + Copy(last->weight, w[5], nrules); + last->gen[0] = g[4]; + + ++last; + n += last->n = 2*ndim*(ndim - 1); + Copy(last->weight, w[6], nrules); + last->gen[0] = g[5]; + last->gen[1] = g[5]; + + ++last; + n += last->n = 2*ndim*(ndim - 1); + Copy(last->weight, w[7], nrules); + last->gen[0] = g[6]; + last->gen[1] = g[6]; + + ++last; + n += last->n = 4*ndim*(ndim - 1)*(ndim - 2)/3; + Copy(last->weight, w[8], nrules); + last->gen[0] = g[7]; + last->gen[1] = g[7]; + last->gen[2] = g[7]; + + ++last; + n += last->n = 4*ndim*(ndim - 1)*(ndim - 2)/3; + Copy(last->weight, w[9], nrules); + last->gen[0] = g[8]; + last->gen[1] = g[8]; + last->gen[2] = g[8]; + + ++last; + n += last->n = 4*ndim*(ndim - 1)*(ndim - 2)/3; + Copy(last->weight, w[10], nrules); + last->gen[0] = g[9]; + last->gen[1] = g[9]; + last->gen[2] = g[9]; + + ++last; + n += last->n = 4*ndim*(ndim - 1)*(ndim - 2); + Copy(last->weight, w[11], nrules); + last->gen[0] = g[10]; + last->gen[1] = g[11]; + last->gen[2] = g[11]; + + ++last; + n += last->n = 4*ndim*(ndim - 1)*(ndim - 2); + Copy(last->weight, w[12], nrules); + last->gen[0] = g[12]; + last->gen[1] = g[12]; + last->gen[2] = g[13]; + + t->rule11.first = first; + t->rule11.last = last; + t->rule11.errcoeff[0] = 4; + t->rule11.errcoeff[1] = .5; + t->rule11.errcoeff[2] = 3; + t->rule11.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( x = first; x <= last; ++x ) + sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]); + s->scale[r] = scale; + s->norm[r] = 1/sum; + } +} + +/*********************************************************************/ + +static void Rule9Alloc(This *t) +{ + static creal w[] = { + -.0023611709677855117884, .11415390023857325268, + -.63833920076702389094, .74849988504685208004, + -.0014324017033399125142, .057471507864489725949, + -.14225104571434243234, -.062875028738286979989, + .254591133248959089, -1.207328566678236261, + .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 = t->ndim; + ccount twondim = 1 << ndim; + count dim, n, r; + Set *first, *last, *s, *x; + + 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]; + + t->rule9.first = first; + t->rule9.last = last; + t->rule9.errcoeff[0] = 5; + t->rule9.errcoeff[1] = 1; + t->rule9.errcoeff[2] = 5; + t->rule9.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( x = first; x <= last; ++x ) + sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]); + s->scale[r] = scale; + s->norm[r] = 1/sum; + } +} + +/*********************************************************************/ + +static void Rule7Alloc(This *t) +{ + 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 = t->ndim; + ccount twondim = 1 << ndim; + count dim, n, r; + Set *first, *last, *s, *x; + + 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]; + + t->rule7.first = first; + t->rule7.last = last; + t->rule7.errcoeff[0] = 5; + t->rule7.errcoeff[1] = 1; + t->rule7.errcoeff[2] = 5; + t->rule7.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( x = first; x <= last; ++x ) + sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]); + s->scale[r] = scale; + s->norm[r] = 1/sum; + } +} + +/*********************************************************************/ + +static inline void RuleIni(Rule *rule) +{ + rule->first = NULL; +} + +/*********************************************************************/ + +static inline bool RuleIniQ(Rule *rule) +{ + return rule->first == NULL; +} + +/*********************************************************************/ + +static inline void RuleFree(Rule *rule) +{ + free(rule->first); +} + +/*********************************************************************/ + +static real *ExpandFS(cThis *t, cBounds *b, real *g, real *x) +{ + count dim, ndim = t->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, j = dim, ix = 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 SampleRule(This *t, ccount iregion) +{ + SAMPLERDEFS; + TYPEDEFSET; + Set *first = (Set *)samples->rule->first; + Set *last = (Set *)samples->rule->last; + Set *s; + creal *errcoeff = samples->rule->errcoeff; + count comp, rul, sn; + + for( s = first; s <= last; ++s ) + if( s->n ) x = ExpandFS(t, b, s->gen, x); + + DoSample(t, n, samples->x, f, t->ndim); + + for( comp = 0; comp < t->ncomp; ++comp ) { + real sum[nrules]; + creal *f1 = f++; + + Zap(sum); + for( s = first; s <= last; ++s ) + for( sn = s->n; sn; --sn ) { + creal fun = *f1; + f1 += t->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[comp].avg = region->vol*sum[0]; + r[comp].err = region->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]) ); + } +} + diff --git a/src/divonne/Sample.c b/src/divonne/Sample.c new file mode 100644 index 0000000..2d0123a --- /dev/null +++ b/src/divonne/Sample.c @@ -0,0 +1,261 @@ +/* + Sample.c + most of what is related to sampling + this file is part of Divonne + last modified 7 Nov 11 th +*/ + + +#define MARKMASK 0xfffffff +#define Marked(x) ((x) & ~MARKMASK) +#define Unmark(x) ((x) & MARKMASK) + +#define EXTRAPOLATE_EPS (.25*t->border.lower) +/*#define EXTRAPOLATE_EPS 0x1p-26*/ + +/*********************************************************************/ + +static inline void SamplesIni(Samples *samples) +{ + samples->x = NULL; +} + +/*********************************************************************/ + +static inline bool SamplesIniQ(cSamples *samples) +{ + return samples->x == NULL; +} + +/*********************************************************************/ + +static inline void SamplesFree(cSamples *samples) +{ + free(samples->x); +} + +/*********************************************************************/ + +static void SampleSobol(This *t, ccount iregion) +{ + SAMPLERDEFS; + real avg[NCOMP], norm; + number i; + count dim, comp; + + for( i = 0; i < n; ++i ) { + t->rng.getrandom(t, x); + for( dim = 0; dim < t->ndim; ++x, ++dim ) + *x = b[dim].lower + *x*(b[dim].upper - b[dim].lower); + } + + DoSample(t, n, samples->x, f, t->ndim); + + ResCopy(avg, f); + f += t->ncomp; + for( i = 2; i < n; ++i ) + for( comp = 0; comp < t->ncomp; ++comp ) + avg[comp] += *f++; + + norm = region->vol/samples->neff; + for( comp = 0; comp < t->ncomp; ++comp ) { + r[comp].avg = norm*avg[comp]; + r[comp].err = 0; + } +} + +/*********************************************************************/ + +static void SampleKorobov(This *t, ccount iregion) +{ + SAMPLERDEFS; + real *xlast = x + t->ndim, *flast = f + t->ncomp; + real avg[NCOMP], norm; + cnumber neff = samples->neff; + number nextra = 0, i; + real dist = 0; + count dim, comp; + + for( i = 1; i < n; ++i ) { + number c = i; + for( dim = 0; dim < t->ndim; ++dim ) { + creal dx = abs(2*c - neff)/(real)neff; + *xlast++ = b[dim].lower + dx*(b[dim].upper - b[dim].lower); + c = c*samples->coeff % neff; + } + } + + for( dim = 0; dim < t->ndim; ++dim ) { + creal dx = (x[dim] = b[dim].upper) - t->border.upper; + if( dx > 0 ) dist += Sq(dx); + } + + if( dist > 0 ) { + dist = sqrt(dist)/EXTRAPOLATE_EPS; + for( dim = 0; dim < t->ndim; ++dim ) { + real x2 = x[dim], dx = x2 - t->border.upper; + if( dx > 0 ) { + x[dim] = t->border.upper; + x2 = t->border.upper - dx/dist; + } + xlast[dim] = x2; + } + nextra = 1; + } + + DoSample(t, n + nextra, x, f, t->ndim); + + ResCopy(avg, flast); + flast += t->ncomp; + for( i = 2; i < n; ++i ) + for( comp = 0; comp < t->ncomp; ++comp ) + avg[comp] += *flast++; + + if( nextra ) { + for( comp = 0; comp < t->ncomp; ++comp ) + f[comp] += dist*(f[comp] - flast[comp]); + for( dim = 0; dim < t->ndim; ++dim ) + x[dim] = b[dim].upper; + } + + norm = region->vol/samples->neff; + for( comp = 0; comp < t->ncomp; ++comp ) { + r[comp].avg = norm*(avg[comp] + avg[comp] + f[comp]); + r[comp].err = 0; + } +} + +/*********************************************************************/ + +#define IsSobol(k) NegQ(k) +#define IsRule(k, d) (k == 9 || k == 7 || (k == 11 && d == 3) || (k == 13 && d == 2)) + +/* The following coding is used for key1, key2, key3: + 0 = for key1, key2: use default, + for key3: do nothing, + 1 = for key3: split region again, + 7 = degree-7 cubature rule, + 9 = degree-9 cubature rule, + 11 = degree-11 cubature rule (only in 3 dims), + 13 = degree-13 cubature rule (only in 2 dims), + -inf..-40 = absolute # of points, Sobol numbers, + -39..-1 = multiplicator, Sobol numbers, + 1..39 = multiplicator, Korobov numbers, + 40..inf = absolute # of points, Korobov numbers. */ + +static count SamplesLookup(This *t, Samples *samples, cint key, + cnumber nwant, cnumber nmax, number nmin) +{ + number n; + + if( key == 13 && t->ndim == 2 ) { + if( RuleIniQ(&t->rule13) ) Rule13Alloc(t); + samples->rule = &t->rule13; + samples->n = n = nmin = t->rule13.n; + samples->sampler = SampleRule; + } + else if( key == 11 && t->ndim == 3 ) { + if( RuleIniQ(&t->rule11) ) Rule11Alloc(t); + samples->rule = &t->rule11; + samples->n = n = nmin = t->rule11.n; + samples->sampler = SampleRule; + } + else if( key == 9 ) { + if( RuleIniQ(&t->rule9) ) Rule9Alloc(t); + samples->rule = &t->rule9; + samples->n = n = nmin = t->rule9.n; + samples->sampler = SampleRule; + } + else if( key == 7 ) { + if( RuleIniQ(&t->rule7) ) Rule7Alloc(t); + samples->rule = &t->rule7; + samples->n = n = nmin = t->rule7.n; + samples->sampler = SampleRule; + } + else { + n = Abs1(key); + if( n < 40 ) n *= nwant; + samples->sampler = (key < 0) ? SampleSobol : + (n = n/2 + 1, SampleKorobov); + samples->n = IMin(n, nmax); + } + + samples->neff = samples->n; + + return IDim(n - nmax) | Marked(nmax - nmin); +} + +/*********************************************************************/ + +static void SamplesAlloc(cThis *t, Samples *samples) +{ +#define FIRST -INT_MAX +#define MarkLast(x) (x | Marked(INT_MAX)) + +#include "KorobovCoeff.c" + + number nx, nf; + + if( samples->sampler == SampleKorobov ) { + enum { max = Elements(prime) - 2 }; + cint n = IMin(2*samples->n - 1, MAXPRIME); + int i = Hash(n), p; + count shift = 2 + NegQ(n - 1000); + + while( i = IMin(IDim(i), max), + n > (p = prime[i + 1]) || n <= prime[i] ) { + cint d = (n - Unmark(p)) >> ++shift; + i += Min1(d); + } + + samples->coeff = coeff[i][t->ndim - KOROBOV_MINDIM]; + samples->neff = p = Unmark(p); + samples->n = p/2 + 1; + } + + nx = t->ndim*(samples->n + 1); /* need 1 for extrapolation */ + nf = t->ncomp*(samples->n + 1); + + Alloc(samples->x, nx + nf + t->ncomp + t->ncomp); + samples->f = samples->x + nx; +} + +/*********************************************************************/ + +static real Sample(This *t, creal *x0) +{ + real xtmp[2*NDIM], ftmp[2*NCOMP], *xlast = xtmp, f; + real dist = 0; + count dim, comp; + number n = 1; + + for( dim = 0; dim < t->ndim; ++dim ) { + creal x1 = *xlast++ = Min(Max(*x0++, 0.), 1.); + real dx; + if( (dx = x1 - t->border.lower) < 0 || + (dx = x1 - t->border.upper) > 0 ) dist += Sq(dx); + } + + if( dist > 0 ) { + dist = sqrt(dist)/EXTRAPOLATE_EPS; + for( dim = 0; dim < t->ndim; ++dim ) { + real x2 = xtmp[dim], dx, b; + if( (dx = x2 - (b = t->border.lower)) < 0 || + (dx = x2 - (b = t->border.upper)) > 0 ) { + xtmp[dim] = b; + x2 = b - dx/dist; + } + *xlast++ = x2; + } + n = 2; + } + + DoSample(t, n, xtmp, ftmp, t->ndim); + + comp = Untag(t->selectedcomp); + f = ftmp[comp]; + if( n > 1 ) f += dist*(f - ftmp[comp + t->ncomp]); + + return Sign(t->selectedcomp)*f; +} + diff --git a/src/divonne/Split.c b/src/divonne/Split.c new file mode 100644 index 0000000..e5883e6 --- /dev/null +++ b/src/divonne/Split.c @@ -0,0 +1,311 @@ +/* + Split.c + determine optimal cuts for splitting a region + this file is part of Divonne + last modified 15 Nov 11 th +*/ + + +#define BNDTOL .05 +#define FRACT .5 +#define BIG 1e10 +#define SINGTOL 1e-4 + +#define LHSTOL .1 +#define GAMMATOL .1 + +/* the next four macros must be in sync with the typedef of Bounds! */ +#define Lower(d) (2*d) +#define Upper(d) (2*d + 1) +#define Dim(i) ((i) >> 1) +#define SignedDelta(i) (2*(i & 1) - 1)*delta[i] + +typedef struct { + count i; + real save, delta; + real f, df, fold; + real lhs, row, sol; +} Cut; + +/*********************************************************************/ + +static inline real Div(creal a, creal b) +{ + return (b != 0 && fabs(b) < BIG*fabs(a)) ? a/b : a; +} + +/*********************************************************************/ + +static void SomeCut(This *t, Cut *cut, Bounds *b) +{ + count dim, maxdim; + static count nextdim = 0; + real xmid[NDIM], ymid, maxdev; + + for( dim = 0; dim < t->ndim; ++dim ) + xmid[dim] = .5*(b[dim].upper + b[dim].lower); + ymid = Sample(t, xmid); + + maxdev = 0; + maxdim = 0; + for( dim = 0; dim < t->ndim; ++dim ) { + real ylower, yupper, dev; + creal x = xmid[dim]; + xmid[dim] = b[dim].lower; + ylower = Sample(t, xmid); + xmid[dim] = b[dim].upper; + yupper = Sample(t, xmid); + xmid[dim] = x; + + dev = fabs(ymid - .5*(ylower + yupper)); + if( dev >= maxdev ) { + maxdev = dev; + maxdim = dim; + } + } + + if( maxdev > 0 ) nextdim = 0; + else maxdim = nextdim++ % t->ndim; + + cut->i = Upper(maxdim); + cut->save = b[maxdim].upper; + b[maxdim].upper = xmid[maxdim]; +} + +/*********************************************************************/ + +static inline real Volume(cThis *t, creal *delta) +{ + real vol = 1; + count dim; + for( dim = 0; dim < t->ndim; ++dim ) + vol *= delta[Lower(dim)] + delta[Upper(dim)]; + return vol; +} + +/*********************************************************************/ + +static inline real SetupEqs(Cut *cut, ccount ncuts, real f) +{ + real sqsum = 0; + Cut *c = &cut[ncuts]; + while( --c >= cut ) { + sqsum += Sq(c->lhs = f - c->f); + f = c->f; + } + return sqsum; +} + +/*********************************************************************/ + +static inline void SolveEqs(Cut *cut, count ncuts, + creal *delta, creal diff) +{ + real last = 0; + real r = 1; + Cut *c; + + for( c = cut; ; ++c ) { + ccount dim = Dim(c->i); + c->row = r -= + Div(diff, (delta[Lower(dim)] + delta[Upper(dim)])*c->df); + if( --ncuts == 0 ) break; + last += r*c->lhs; + } + + last = Div(c->lhs - last, r); + + for( ; c >= cut; last += (--c)->lhs ) { + creal delmin = -(c->delta = delta[c->i]); + creal delmax = FRACT*(delmin + c->save); + c->sol = Div(last, c->df); + if( c->sol > delmax ) c->sol = .75*delmax; + if( c->sol < delmin ) c->sol = .75*delmin; + } +} + +/*********************************************************************/ + +static count FindCuts(This *t, Cut *cut, Bounds *bounds, creal vol, + real *xmajor, creal fmajor, creal fdiff) +{ + cint sign = (fdiff < 0) ? -1 : 1; + + count ncuts = 0, icut; + real delta[2*NDIM]; + real gamma, fgamma, lhssq; + count dim, div; + + for( dim = 0; dim < t->ndim; ++dim ) { + cBounds *b = &bounds[dim]; + creal xsave = xmajor[dim]; + real dist = b->upper - xsave; + if( dist >= BNDTOL*(b->upper - b->lower) ) { + Cut *c = &cut[ncuts++]; + c->i = Upper(dim); + c->save = dist; + xmajor[dim] += dist *= FRACT; + c->f = Sample(t, xmajor); + xmajor[dim] = xsave; + } + delta[Upper(dim)] = dist; + } + + for( dim = 0; dim < t->ndim; ++dim ) { + cBounds *b = &bounds[dim]; + creal xsave = xmajor[dim]; + real dist = xsave - b->lower; + if( dist >= BNDTOL*(b->upper - b->lower) ) { + Cut *c = &cut[ncuts++]; + c->i = Lower(dim); + c->save = dist; + xmajor[dim] -= dist *= FRACT; + c->f = Sample(t, xmajor); + xmajor[dim] = xsave; + } + delta[Lower(dim)] = dist; + } + + if( ncuts == 0 ) { + SomeCut(t, cut, bounds); + return 1; + } + + for( ; ; ) { + real mindiff = INFTY; + Cut *mincut = cut; + + for( icut = 0; icut < ncuts; ++icut ) { + Cut *c = &cut[icut]; + creal diff = fabs(fmajor - c->f); + if( diff <= mindiff ) { + mindiff = diff; + mincut = c; + } + } + + gamma = Volume(t, delta)/vol; + fgamma = fmajor + (gamma - 1)*fdiff; + + if( sign*(mincut->f - fgamma) < 0 ) break; + + if( --ncuts == 0 ) { + SomeCut(t, cut, bounds); + return 1; + } + + delta[mincut->i] = mincut->save; + memcpy(mincut, mincut + 1, (char *)&cut[ncuts] - (char *)mincut); + } + + for( icut = 0; icut < ncuts; ++icut ) { + Cut *c = &cut[icut]; + c->fold = c->f; + c->df = (c->f - fmajor)/delta[c->i]; + } + + lhssq = SetupEqs(cut, ncuts, fgamma); + +repeat: + SolveEqs(cut, ncuts, delta, gamma*fdiff); + + for( div = 1; div <= 16; div *= 4 ) { + real gammanew, lhssqnew; + + for( icut = 0; icut < ncuts; ++icut ) { + Cut *c = &cut[icut]; + real *x = &xmajor[Dim(c->i)]; + creal xsave = *x; + delta[c->i] = c->delta + c->sol/div; + *x += SignedDelta(c->i); + c->f = Sample(t, xmajor); + *x = xsave; + } + + gammanew = Volume(t, delta)/vol; + fgamma = fmajor + (gammanew - 1)*fdiff; + lhssqnew = SetupEqs(cut, ncuts, fgamma); + + if( lhssqnew <= lhssq ) { + real fmax; + + if( fabs(gammanew - gamma) < GAMMATOL*gamma ) break; + gamma = gammanew; + + fmax = fabs(fgamma); + for( icut = 0; icut < ncuts; ++icut ) { + Cut *c = &cut[icut]; + creal dfmin = SINGTOL*c->df; + creal sol = c->sol/div; + real df = c->f - c->fold; + df = (fabs(sol) < BIG*fabs(df)) ? df/sol : 1; + c->df = (fabs(df) < fabs(dfmin)) ? dfmin : df; + fmax = Max(fmax, fabs(c->f)); + c->fold = c->f; + } + + if( lhssqnew < Sq((1 + fmax)*LHSTOL) ) break; + lhssq = lhssqnew; + goto repeat; + } + } + + for( icut = 0; icut < ncuts; ++icut ) { + Cut *c = &cut[icut]; + real *b = (real *)bounds + c->i; + c->save = *b; + *b = xmajor[Dim(c->i)] + SignedDelta(c->i); + } + + return ncuts; +} + +/*********************************************************************/ + +static void Split(This *t, ccount iregion) +{ + TYPEDEFREGION; + Region *region = RegionPtr(iregion); + Cut cut[2*NDIM], *c; + Errors errors[NCOMP]; + count ncuts, succ; + int depth; + real *b, tmp; + + t->selectedcomp = region->cutcomp; + t->neval_cut -= t->neval; + ncuts = FindCuts(t, cut, region->bounds, region->vol, + (real *)region->result + region->xmajor, region->fmajor, + region->fmajor - region->fminor); + t->neval_cut += t->neval; + + depth = region->depth - ncuts; + + EnlargeRegions(t, ++ncuts); + region = RegionPtr(iregion); + region->depth = -ncuts; + succ = iregion + region->next; + region->next = t->nregions - iregion; + b = (real *)region->bounds; + + region = RegionPtr(t->nregions); + VecCopy(region->bounds, b); + region->depth = IDim(depth) + 1; + region->next = 1; + region->isamples = 0; + for( c = cut; --ncuts; ++c ) { + ccount ci = c->i; + creal tmp = b[ci ^ 1]; + b[ci ^ 1] = b[ci]; + b[ci] = c->save; + region = RegionPtr(++t->nregions); + VecCopy(region->bounds, b); + region->depth = IDim(depth) + 1; + region->next = 1; + region->isamples = 0; + ++depth; + b[ci ^ 1] = tmp; + } + region->next = succ - t->nregions++; +} + diff --git a/src/divonne/common.c b/src/divonne/common.c new file mode 100644 index 0000000..8da4d00 --- /dev/null +++ b/src/divonne/common.c @@ -0,0 +1,34 @@ +/* + common.c + includes most of the modules + this file is part of Divonne + last modified 7 Nov 11 th +*/ + + +#include "Random.c" +#include "ChiSquare.c" +#include "Rule.c" +#include "Sample.c" +#include "FindMinimum.c" +#include "Split.c" +#include "Explore.c" +#include "Iterate.c" + +static inline bool BadDimension(cThis *t, ccount key) +{ + if( t->ndim > NDIM ) return true; + if( IsSobol(key) ) return + t->ndim < SOBOL_MINDIM || (t->seed == 0 && t->ndim > SOBOL_MAXDIM); + if( IsRule(key, t->ndim) ) return t->ndim < 1; + return t->ndim < KOROBOV_MINDIM || t->ndim > KOROBOV_MAXDIM; +} + +static inline bool BadComponent(cThis *t) +{ + if( t->ncomp > NCOMP ) return true; + return t->ncomp < 1; +} + +#include "Integrate.c" + diff --git a/src/divonne/decl.h b/src/divonne/decl.h new file mode 100644 index 0000000..99fac6d --- /dev/null +++ b/src/divonne/decl.h @@ -0,0 +1,127 @@ +/* + decl.h + Type declarations + this file is part of Divonne + last modified 14 Nov 11 th +*/ + + +#include "stddecl.h" + +#define Tag(x) ((x) | INT_MIN) +#define Untag(x) ((x) & INT_MAX) + +typedef struct { + real lower, upper; +} Bounds; + +typedef const Bounds cBounds; + +typedef struct { + real avg, err; +} PhaseResult; + +typedef struct { + real avg, spreadsq; + real spread, secondspread; + real nneed, maxerrsq, mindevsq; + PhaseResult phase[2]; + int iregion; +} Totals; + +typedef struct { + void *first, *last; + real errcoeff[3]; + count n; +} Rule; + +typedef const Rule cRule; + +typedef struct samples { + real *x, *f; + void (*sampler)(struct _this *t, ccount); + cRule *rule; + number n, neff; + count coeff; +} Samples; + +typedef const Samples cSamples; + +typedef struct { + real diff, err, spread; +} Errors; + +typedef const Errors cErrors; + +typedef int (*Integrand)(ccount *, creal *, ccount *, real *, void *, cint *); + +typedef void (*PeakFinder)(ccount *, cBounds *, number *, real *); + +typedef struct _this { + count ndim, ncomp; +#ifndef MLVERSION + Integrand integrand; + void *userdata; + PeakFinder peakfinder; + int ncores, *child; + int running, nchildren; + fd_set children; +#endif + real epsrel, epsabs; + int flags, seed; + number mineval, maxeval; + int key1, key2, key3; + count maxpass; + Bounds border; + real maxchisq, mindeviation; + number ngiven, nextra; + real *xgiven, *xextra, *fgiven, *fextra; + count ldxgiven; + count nregions; + number neval, neval_opt, neval_cut; + count phase; + count selectedcomp, size; + Samples samples[3]; + Totals *totals; + Rule rule7, rule9, rule11, rule13; + RNGState rng; + void *voidregion; + jmp_buf abort; +} This; + +typedef const This cThis; + +#define TYPEDEFREGION \ + typedef struct { \ + real avg, err, spread, chisq; \ + real fmin, fmax; \ + real xmin[NDIM], xmax[NDIM]; \ + } Result; \ + typedef const Result cResult; \ + typedef struct region { \ + int depth, next; \ + count isamples, cutcomp, xmajor; \ + real fmajor, fminor, vol; \ + Bounds bounds[NDIM]; \ + Result result[NCOMP]; \ + } Region + +#define RegionPtr(n) (&((Region *)t->voidregion)[n]) + +#define CHUNKSIZE 4096 + +#define AllocRegions(t) \ + MemAlloc((t)->voidregion, ((t)->size = CHUNKSIZE)*sizeof(Region)) + +#define EnlargeRegions(t, n) if( (t)->nregions + n > (t)->size ) \ + ReAlloc((t)->voidregion, ((t)->size += CHUNKSIZE)*sizeof(Region)) + +#define SAMPLERDEFS \ + TYPEDEFREGION; \ + Region *region = RegionPtr(iregion); \ + cBounds *b = region->bounds; \ + Result *r = region->result; \ + cSamples *samples = &t->samples[region->isamples]; \ + real *x = samples->x, *f = samples->f; \ + cnumber n = samples->n + |