Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F8308886
Rule.c
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
20 KB
Subscribers
None
Rule.c
View Options
/*
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 8 Jun 10 th
*/
enum
{
nrules
=
5
};
#define TYPEDEFSET \
typedef struct { \
count n; \
real weight[nrules], scale[nrules], norm[nrules]; \
real gen[1000]; \
} Set
//#TT used to be real gen[NDIM];
/*********************************************************************/
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
->
rule
.
first
=
first
;
t
->
rule
.
last
=
last
;
t
->
rule
.
errcoeff
[
0
]
=
10
;
t
->
rule
.
errcoeff
[
1
]
=
1
;
t
->
rule
.
errcoeff
[
2
]
=
5
;
t
->
rule
.
n
=
n
;
for
(
s
=
first
;
s
<=
last
;
++
s
)
for
(
r
=
1
;
r
<
nrules
-
1
;
++
r
)
{
creal
scale
=
(
s
->
weight
[
r
]
==
0
)
?
100
:
-
s
->
weight
[
r
+
1
]
/
s
->
weight
[
r
];
real
sum
=
0
;
for
(
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
->
rule
.
first
=
first
;
t
->
rule
.
last
=
last
;
t
->
rule
.
errcoeff
[
0
]
=
4
;
t
->
rule
.
errcoeff
[
1
]
=
.5
;
t
->
rule
.
errcoeff
[
2
]
=
3
;
t
->
rule
.
n
=
n
;
for
(
s
=
first
;
s
<=
last
;
++
s
)
for
(
r
=
1
;
r
<
nrules
-
1
;
++
r
)
{
creal
scale
=
(
s
->
weight
[
r
]
==
0
)
?
100
:
-
s
->
weight
[
r
+
1
]
/
s
->
weight
[
r
];
real
sum
=
0
;
for
(
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
->
rule
.
first
=
first
;
t
->
rule
.
last
=
last
;
t
->
rule
.
errcoeff
[
0
]
=
5
;
t
->
rule
.
errcoeff
[
1
]
=
1
;
t
->
rule
.
errcoeff
[
2
]
=
5
;
t
->
rule
.
n
=
n
;
for
(
s
=
first
;
s
<=
last
;
++
s
)
for
(
r
=
1
;
r
<
nrules
-
1
;
++
r
)
{
creal
scale
=
(
s
->
weight
[
r
]
==
0
)
?
100
:
-
s
->
weight
[
r
+
1
]
/
s
->
weight
[
r
];
real
sum
=
0
;
for
(
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
->
rule
.
first
=
first
;
t
->
rule
.
last
=
last
;
t
->
rule
.
errcoeff
[
0
]
=
5
;
t
->
rule
.
errcoeff
[
1
]
=
1
;
t
->
rule
.
errcoeff
[
2
]
=
5
;
t
->
rule
.
n
=
n
;
for
(
s
=
first
;
s
<=
last
;
++
s
)
for
(
r
=
1
;
r
<
nrules
-
1
;
++
r
)
{
creal
scale
=
(
s
->
weight
[
r
]
==
0
)
?
100
:
-
s
->
weight
[
r
+
1
]
/
s
->
weight
[
r
];
real
sum
=
0
;
for
(
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
RuleAlloc
(
This
*
t
)
{
if
(
t
->
key
==
13
&&
t
->
ndim
==
2
)
Rule13Alloc
(
t
);
else
if
(
t
->
key
==
11
&&
t
->
ndim
==
3
)
Rule11Alloc
(
t
);
else
if
(
t
->
key
==
9
)
Rule9Alloc
(
t
);
else
if
(
t
->
key
==
7
)
Rule7Alloc
(
t
);
else
{
if
(
t
->
ndim
==
2
)
Rule13Alloc
(
t
);
else
if
(
t
->
ndim
==
3
)
Rule11Alloc
(
t
);
else
Rule9Alloc
(
t
);
}
Alloc
(
t
->
rule
.
x
,
t
->
rule
.
n
*
(
t
->
ndim
+
t
->
ncomp
));
t
->
rule
.
f
=
t
->
rule
.
x
+
t
->
rule
.
n
*
t
->
ndim
;
}
/*********************************************************************/
static
inline
void
RuleFree
(
cThis
*
t
)
{
free
(
t
->
rule
.
x
);
free
(
t
->
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
Sample
(
This
*
t
,
void
*
voidregion
)
{
TYPEDEFREGION
;
TYPEDEFSET
;
Region
*
const
region
=
(
Region
*
)
voidregion
;
creal
vol
=
ldexp
(
1.
,
-
region
->
div
);
real
*
x
=
t
->
rule
.
x
,
*
f
=
t
->
rule
.
f
;
Set
*
first
=
(
Set
*
)
t
->
rule
.
first
,
*
last
=
(
Set
*
)
t
->
rule
.
last
,
*
s
;
creal
*
errcoeff
=
t
->
rule
.
errcoeff
;
creal
ratio
=
Sq
(
first
[
2
].
gen
[
0
]
/
first
[
1
].
gen
[
0
]);
ccount
offset
=
2
*
t
->
ndim
*
t
->
ncomp
;
count
dim
,
comp
,
rul
,
n
,
maxdim
=
0
;
real
maxrange
=
0
;
for
(
dim
=
0
;
dim
<
t
->
ndim
;
++
dim
)
{
cBounds
*
b
=
&
region
->
bounds
[
dim
];
creal
range
=
b
->
upper
-
b
->
lower
;
if
(
range
>
maxrange
)
{
maxrange
=
range
;
maxdim
=
dim
;
}
}
for
(
s
=
first
;
s
<=
last
;
++
s
)
if
(
s
->
n
)
x
=
ExpandFS
(
t
,
region
->
bounds
,
s
->
gen
,
x
);
DoSample
(
t
,
t
->
rule
.
n
,
t
->
rule
.
x
,
f
);
for
(
comp
=
0
;
comp
<
t
->
ncomp
;
++
comp
)
{
Result
*
r
=
&
region
->
result
[
comp
];
real
sum
[
nrules
];
creal
*
f1
=
f
;
creal
base
=
*
f1
*
2
*
(
1
-
ratio
);
real
maxdiff
=
0
;
count
bisectdim
=
maxdim
;
for
(
dim
=
0
;
dim
<
t
->
ndim
;
++
dim
)
{
creal
*
fp
=
f1
+
t
->
ncomp
;
creal
*
fm
=
fp
+
t
->
ncomp
;
creal
fourthdiff
=
fabs
(
base
+
ratio
*
(
fp
[
0
]
+
fm
[
0
])
-
(
fp
[
offset
]
+
fm
[
offset
]));
f1
=
fm
;
if
(
fourthdiff
>
maxdiff
)
{
maxdiff
=
fourthdiff
;
bisectdim
=
dim
;
}
}
r
->
bisectdim
=
bisectdim
;
f1
=
f
++
;
Zap
(
sum
);
for
(
s
=
first
;
s
<=
last
;
++
s
)
for
(
n
=
s
->
n
;
n
;
--
n
)
{
creal
fun
=
*
f1
;
f1
+=
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
->
avg
=
vol
*
sum
[
0
];
r
->
err
=
vol
*
(
(
errcoeff
[
0
]
*
sum
[
1
]
<=
sum
[
2
]
&&
errcoeff
[
0
]
*
sum
[
2
]
<=
sum
[
3
])
?
errcoeff
[
1
]
*
sum
[
1
]
:
errcoeff
[
2
]
*
Max
(
Max
(
sum
[
1
],
sum
[
2
]),
sum
[
3
])
);
}
if
(
VERBOSE
>
2
)
{
char
s
[
64
*
NDIM
+
128
*
NCOMP
],
*
p
=
s
;
for
(
dim
=
0
;
dim
<
t
->
ndim
;
++
dim
)
{
cBounds
*
b
=
&
region
->
bounds
[
dim
];
p
+=
sprintf
(
p
,
(
dim
==
0
)
?
"
\n
Region ("
REALF
") - ("
REALF
")"
:
"
\n
("
REALF
") - ("
REALF
")"
,
b
->
lower
,
b
->
upper
);
}
for
(
comp
=
0
;
comp
<
t
->
ncomp
;
++
comp
)
{
cResult
*
r
=
&
region
->
result
[
comp
];
p
+=
sprintf
(
p
,
"
\n
["
COUNT
"] "
REAL
" +- "
REAL
,
comp
+
1
,
r
->
avg
,
r
->
err
);
}
Print
(
s
);
}
}
File Metadata
Details
Attached
Mime Type
text/plain
Expires
Sat, Dec 21, 1:29 PM (19 h, 27 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4022961
Default Alt Text
Rule.c (20 KB)
Attached To
rSARTRESVN sartresvn
Event Timeline
Log In to Comment