1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
|
/*
* John C. Bowman and Supakorn "Jamie" Rassameemasmuang <jamievlin@outlook.com>
* University of Alberta
* CUDA Adaptive Simpson integration
*/
#include <cuda.h>
#include "utils.cuh"
// Compute a numerical approximation to an integral via adaptive Simpson's Rule
// This routine ignores underflow.
__device__ constexpr float sixth=1.0/6.0;
__device__ constexpr int depth=10;
__device__ constexpr int acc = 1.0/256;
template<typename T>
struct TABLE {
bool left; // left interval?
float dat;
T psum, f1t, f2t, f3t, estr;
};
template<typename T=glm::vec3, typename TUtil=Vec3Utility, typename Tf>
__device__ inline T
simpson(Tf f, // Function to be integrated.
float a, float b, // Lower, upper limits of integration.
float acc) // Desired relative accuracy of integral.
// Try to make |error| <= acc*abs(integral).
{
T integral,diff,area,estl,estr,est,fv0,fv1,fv2,fv3,fv4;
float dx;
TABLE<T> table[depth],*p,*pstop;
p=table;
pstop=table+depth-1;
p->left=true;
p->psum=TUtil::init();
float alpha=a;
float da=b-a;
fv0=f(alpha);
fv2=f(alpha+0.5f*da);
fv4=f(alpha+da);
float wt=sixth*da;
est=wt*(fv0+4.0f*fv2+fv4);
area=est;
float acc2=acc*acc;
// Have estimate est of integral on (alpha, alpha+da).
// Bisect and compute estimates on left and right half intervals.
// integral is the best value for the integral.
for(;;) {
dx=0.5f*da;
float arg=alpha+0.5f*dx;
fv1=f(arg);
fv3=f(arg+dx);
wt=sixth*dx;
estl=wt*(fv0+4.0f*fv1+fv2);
estr=wt*(fv2+4.0f*fv3+fv4);
integral=estl+estr;
diff=est-integral;
area -= diff;
if(p >= pstop || (TUtil::abs2(diff) <= acc2*TUtil::abs2(area))) {
// Accept approximate integral.
// If it was a right interval, add results to finish at this level.
// If it was a left interval, process right interval.
for(;;) {
if(p->left == false) { // process right-half interval
alpha += da;
p->left=true;
p->psum=integral;
fv0=p->f1t;
fv2=p->f2t;
fv4=p->f3t;
da=p->dat;
est=p->estr;
break;
}
integral += p->psum;
if(--p <= table) return integral;
}
} else {
// Raise level and store information for processing right-half interval.
++p;
da=dx;
est=estl;
p->left=false;
p->f1t=fv2;
p->f2t=fv3;
p->f3t=fv4;
p->dat=dx;
p->estr=estr;
fv4=fv2;
fv2=fv1;
}
}
}
|
