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
|
#include <stdio.h>
#include <stdarg.h>
// define MPFR_USE_VA_LIST to enable the use of mpfr_inits and mpfr_clears
#define MPFR_USE_VA_LIST
#include <mpfr.h>
#include <libinum.h>
int func(mpfr_t* out, mpfr_t in, void* extra_args);
int main(int argc, const char* argv[]){
int ret;
mpfr_t tolerance, val, lower, upper;
array_mpfr abcissa, weights;
// precision of MPFR floats
mpfr_set_default_prec(53);
// tolerance for computing weights
mpfr_init(tolerance);
mpfr_set_d(tolerance, 1.e-11, MPFR_RNDN);
// compute weights
ret=gauss_legendre_weights_mpfr(10, tolerance, 10000, &abcissa, &weights);
// return codes
if(ret==LIBINUM_ERROR_MAXITER){
printf("error: maximum number of iterations reached when computing the integration abcissa\n");
mpfr_clear(tolerance);
return(ret);
}
else if(ret==LIBINUM_ERROR_NAN){
printf("error: infinity encountered when computing the integration abcissa\n");
mpfr_clear(tolerance);
return(ret);
}
else{
printf("abcissa:\n");
array_mpfr_print(abcissa);
printf("\nweights:\n");
array_mpfr_print(weights);
printf("\n");
}
mpfr_clear(tolerance);
mpfr_inits(val, lower, upper, NULL);
mpfr_set_ui(lower, 0, MPFR_RNDN);
mpfr_set_ui(upper, 2, MPFR_RNDN);
ret=integrate_gauss_mpfr(&val, &func, lower, upper, abcissa, weights, NULL);
// return codes
if(ret==LIBINUM_ERROR_SIZE_MISMATCH){
printf("error: the mpfr_arrays 'abcissa' and 'weights' must have the same length\n");
mpfr_clears(val, upper, lower, NULL);
array_mpfr_free(abcissa);
array_mpfr_free(weights);
return(ret);
}
else if(ret==LIBINUM_ERROR_NAN){
printf("error: the integrand is singular\n");
array_mpfr_free(abcissa);
array_mpfr_free(weights);
mpfr_clears(val, upper, lower, NULL);
return(ret);
}
else{
printf("\\int_0^2 x/sin(x) dx = ");
fprint_mpfr(stdout, val);
printf("\n");
}
mpfr_clears(val, upper, lower, NULL);
array_mpfr_free(abcissa);
array_mpfr_free(weights);
return(0);
}
/*
// e^x
int func(mpfr_t* out, mpfr_t in, void* extra_args){
mpfr_exp(*out, in, MPFR_RNDN);
return(0);
}
*/
// x/sin(x)
int func(mpfr_t* out, mpfr_t in, void* extra_args){
mpfr_sin(*out, in, MPFR_RNDN);
mpfr_div(*out, in, *out, MPFR_RNDN);
return(0);
}
|
