1 /*
2 * R : A Computer Language for Statistical Data Analysis
3 * Copyright (C) 2001-2016 the R Core Team
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, a copy is available at
17 * https://www.R-project.org/Licenses/
18 */
19
20 #ifdef HAVE_CONFIG_H
21 #include <config.h>
22 #endif
23
24 #include <math.h> // for isfinite
25 #include <Rinternals.h>
26 #include <R_ext/Applic.h>
27
28 #ifdef ENABLE_NLS
29 #include <libintl.h>
30 #define _(String) dgettext ("stats", String)
31 #else
32 #define _(String) (String)
33 #endif
34
35 /* called via .External(.) :*/
36 SEXP call_dqags(SEXP args);
37 SEXP call_dqagi(SEXP args);
38
39 typedef struct int_struct
40 {
41 SEXP f; /* function */
42 SEXP env; /* where to evaluate the calls */
43 } int_struct, *IntStruct;
44
45
46 /* This is *the* ``integr_fn f'' used when called from R : */
Rintfn(double * x,int n,void * ex)47 static void Rintfn(double *x, int n, void *ex)
48 {
49 SEXP args, resultsxp, tmp;
50 int i;
51 IntStruct IS = (IntStruct) ex;
52
53 PROTECT(args = allocVector(REALSXP, n));
54 for(i = 0; i < n; i++) REAL(args)[i] = x[i];
55
56 PROTECT(tmp = lang2(IS->f , args));
57 PROTECT(resultsxp = eval(tmp, IS->env));
58
59 if(length(resultsxp) != n)
60 error("evaluation of function gave a result of wrong length");
61 if(TYPEOF(resultsxp) == INTSXP) {
62 resultsxp = coerceVector(resultsxp, REALSXP);
63 } else if(TYPEOF(resultsxp) != REALSXP)
64 error("evaluation of function gave a result of wrong type");
65 for(i = 0; i < n; i++) {
66 x[i] = REAL(resultsxp)[i];
67 if(!R_FINITE(x[i]))
68 error("non-finite function value");
69 }
70 UNPROTECT(3);
71 return;
72 }
73
call_dqags(SEXP args)74 SEXP call_dqags(SEXP args)
75 {
76 int_struct is;
77 SEXP ans, ansnames;
78 double lower, upper, epsabs, epsrel, result, abserr, *work;
79 int neval, ier, limit, lenw, last, *iwork;
80
81 args = CDR(args);
82 is.f = CAR(args); args = CDR(args);
83 is.env = CAR(args); args = CDR(args);
84 if(length(CAR(args)) > 1) error(_("'%s' must be of length one"), "lower");
85 lower = asReal(CAR(args)); args = CDR(args);
86 if(length(CAR(args)) > 1) error(_("'%s' must be of length one"), "upper");
87 upper = asReal(CAR(args)); args = CDR(args);
88 epsabs = asReal(CAR(args)); args = CDR(args);
89 epsrel = asReal(CAR(args)); args = CDR(args);
90 limit = asInteger(CAR(args)); args = CDR(args);
91 lenw = 4 * limit;
92 iwork = (int *) R_alloc((size_t) limit, sizeof(int));
93 work = (double *) R_alloc((size_t) lenw, sizeof(double));
94
95 Rdqags(Rintfn, (void*)&is,
96 &lower, &upper, &epsabs, &epsrel, &result,
97 &abserr, &neval, &ier, &limit, &lenw, &last, iwork, work);
98
99 PROTECT(ans = allocVector(VECSXP, 4));
100 PROTECT(ansnames = allocVector(STRSXP, 4));
101 SET_STRING_ELT(ansnames, 0, mkChar("value"));
102 SET_VECTOR_ELT(ans, 0, allocVector(REALSXP, 1));
103 REAL(VECTOR_ELT(ans, 0))[0] = result;
104 SET_STRING_ELT(ansnames, 1, mkChar("abs.error"));
105 SET_VECTOR_ELT(ans, 1, allocVector(REALSXP, 1));
106 REAL(VECTOR_ELT(ans, 1))[0] = abserr;
107 SET_STRING_ELT(ansnames, 2, mkChar("subdivisions"));
108 SET_VECTOR_ELT(ans, 2, allocVector(INTSXP, 1));
109 INTEGER(VECTOR_ELT(ans, 2))[0] = last;
110 SET_STRING_ELT(ansnames, 3, mkChar("ierr"));
111 SET_VECTOR_ELT(ans, 3, allocVector(INTSXP, 1));
112 INTEGER(VECTOR_ELT(ans, 3))[0] = ier;
113 setAttrib(ans, R_NamesSymbol, ansnames);
114 UNPROTECT(2);
115 return ans;
116 }
117
call_dqagi(SEXP args)118 SEXP call_dqagi(SEXP args)
119 {
120 int_struct is;
121 SEXP ans, ansnames;
122 double bound, epsabs, epsrel, result, abserr, *work;
123 int inf, neval, ier, limit, lenw, last, *iwork;
124
125 args = CDR(args);
126 is.f = CAR(args); args = CDR(args);
127 is.env = CAR(args); args = CDR(args);
128 if(length(CAR(args)) > 1) error(_("'%s' must be of length one"), "bound");
129 bound = asReal(CAR(args)); args = CDR(args);
130 inf = asInteger(CAR(args)); args = CDR(args);
131 epsabs = asReal(CAR(args)); args = CDR(args);
132 epsrel = asReal(CAR(args)); args = CDR(args);
133 limit = asInteger(CAR(args)); args = CDR(args);
134 lenw = 4 * limit;
135 iwork = (int *) R_alloc((size_t) limit, sizeof(int));
136 work = (double *) R_alloc((size_t) lenw, sizeof(double));
137
138 Rdqagi(Rintfn, (void*)&is, &bound,&inf,&epsabs,&epsrel,&result,
139 &abserr,&neval,&ier,&limit,&lenw,&last,iwork,work);
140
141 PROTECT(ans = allocVector(VECSXP, 4));
142 PROTECT(ansnames = allocVector(STRSXP, 4));
143 SET_STRING_ELT(ansnames, 0, mkChar("value"));
144 SET_VECTOR_ELT(ans, 0, allocVector(REALSXP, 1));
145 REAL(VECTOR_ELT(ans, 0))[0] = result;
146 SET_STRING_ELT(ansnames, 1, mkChar("abs.error"));
147 SET_VECTOR_ELT(ans, 1, allocVector(REALSXP, 1));
148 REAL(VECTOR_ELT(ans, 1))[0] = abserr;
149 SET_STRING_ELT(ansnames, 2, mkChar("subdivisions"));
150 SET_VECTOR_ELT(ans, 2, allocVector(INTSXP, 1));
151 INTEGER(VECTOR_ELT(ans, 2))[0] = last;
152 SET_STRING_ELT(ansnames, 3, mkChar("ierr"));
153 SET_VECTOR_ELT(ans, 3, allocVector(INTSXP, 1));
154 INTEGER(VECTOR_ELT(ans, 3))[0] = ier;
155 setAttrib(ans, R_NamesSymbol, ansnames);
156 UNPROTECT(2);
157 return ans;
158 }
159