typedef double (__cdecl *func)(const double);
double rombf(const double x);
double romb(double upper_limit, double lower_limit, double precision, func f);
int main(int argc, char *argv[])
{
double upper_limit = 0.0; // 積分上限
double lower_limit = 1.0; // 積分下限
double precision = 0.000001; // 積分精度
double result = 0.0; // 積分結果
result = romb(upper_limit, lower_limit, precision, &rombf);
printf("result = %e \n\n", result);
return 0;
}
double romb(double upper_limit, double lower_limit, double precision, func f)
{
int m, n, i, k;
double y[10], h, ep, p, x, s, q;
h = lower_limit - upper_limit;
y[0] = h * ((*f)(upper_limit) + (*f)(lower_limit)) / 2.0;
m = 1; n = 1; ep = precision + 1.0;
while ((ep >= precision) && (m <= 9))
{
p = 0.0;
for (i = 0; i <= n - 1; i++)
{
x = upper_limit + (i + 0.5) * h;
p = p + (*f)(x);
}
p = (y[0] + h*p) / 2.0;
s = 1.0;
for (k = 1; k <= m; k++)
{
s = 4.0 * s;
q = (s * p - y[k-1]) / (s - 1.0);
y[k-1] = p; p = q;
}
ep = fabs(q - y[m-1]);
m = m + 1; y[m-1] = q; n = n + n; h = h / 2.0;
}
return(q);
}
double rombf(double x)
{
double y;
y = x / (4.0 + x * x);
return(y);
}