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
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
|
/* integration/util.c
*
* Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3 of the License, or (at
* your option) any later version.
*
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
static inline
void update (gsl_integration_workspace * workspace,
double a1, double b1, double area1, double error1,
double a2, double b2, double area2, double error2);
static inline void
retrieve (const gsl_integration_workspace * workspace,
double * a, double * b, double * r, double * e);
static inline
void update (gsl_integration_workspace * workspace,
double a1, double b1, double area1, double error1,
double a2, double b2, double area2, double error2)
{
double * alist = workspace->alist ;
double * blist = workspace->blist ;
double * rlist = workspace->rlist ;
double * elist = workspace->elist ;
size_t * level = workspace->level ;
const size_t i_max = workspace->i ;
const size_t i_new = workspace->size ;
const size_t new_level = workspace->level[i_max] + 1;
/* append the newly-created intervals to the list */
if (error2 > error1)
{
alist[i_max] = a2; /* blist[maxerr] is already == b2 */
rlist[i_max] = area2;
elist[i_max] = error2;
level[i_max] = new_level;
alist[i_new] = a1;
blist[i_new] = b1;
rlist[i_new] = area1;
elist[i_new] = error1;
level[i_new] = new_level;
}
else
{
blist[i_max] = b1; /* alist[maxerr] is already == a1 */
rlist[i_max] = area1;
elist[i_max] = error1;
level[i_max] = new_level;
alist[i_new] = a2;
blist[i_new] = b2;
rlist[i_new] = area2;
elist[i_new] = error2;
level[i_new] = new_level;
}
workspace->size++;
if (new_level > workspace->maximum_level)
{
workspace->maximum_level = new_level;
}
qpsrt (workspace) ;
}
static inline void
retrieve (const gsl_integration_workspace * workspace,
double * a, double * b, double * r, double * e)
{
const size_t i = workspace->i;
double * alist = workspace->alist;
double * blist = workspace->blist;
double * rlist = workspace->rlist;
double * elist = workspace->elist;
*a = alist[i] ;
*b = blist[i] ;
*r = rlist[i] ;
*e = elist[i] ;
}
static inline double
sum_results (const gsl_integration_workspace * workspace);
static inline double
sum_results (const gsl_integration_workspace * workspace)
{
const double * const rlist = workspace->rlist ;
const size_t n = workspace->size;
size_t k;
double result_sum = 0;
for (k = 0; k < n; k++)
{
result_sum += rlist[k];
}
return result_sum;
}
static inline int
subinterval_too_small (double a1, double a2, double b2);
static inline int
subinterval_too_small (double a1, double a2, double b2)
{
const double e = GSL_DBL_EPSILON;
const double u = GSL_DBL_MIN;
double tmp = (1 + 100 * e) * (fabs (a2) + 1000 * u);
int status = fabs (a1) <= tmp && fabs (b2) <= tmp;
return status;
}
|
