1 /*
2 * R : A Computer Language for Statistical Data Analysis
3 * Copyright (C) 1995, 1996 Robert Gentleman and Ross Ihaka
4 * Copyright (C) 1997--2015 The R Core Team
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, a copy is available at
18 * https://www.R-project.org/Licenses/
19 */
20
21 #ifdef HAVE_CONFIG_H
22 #include <config.h>
23 #endif
24
25 #include <Defn.h>
26 #include <Internal.h>
27
cumsum(SEXP x,SEXP s)28 static SEXP cumsum(SEXP x, SEXP s)
29 {
30 LDOUBLE sum = 0.;
31 double *rx = REAL(x), *rs = REAL(s);
32 for (R_xlen_t i = 0 ; i < XLENGTH(x) ; i++) {
33 sum += rx[i]; /* NA and NaN propagated */
34 rs[i] = (double) sum;
35 }
36 return s;
37 }
38
39 /* We need to ensure that overflow gives NA here */
icumsum(SEXP x,SEXP s)40 static SEXP icumsum(SEXP x, SEXP s)
41 {
42 int *ix = INTEGER(x), *is = INTEGER(s);
43 double sum = 0.0;
44 for (R_xlen_t i = 0 ; i < XLENGTH(x) ; i++) {
45 if (ix[i] == NA_INTEGER) break;
46 sum += ix[i];
47 if(sum > INT_MAX || sum < 1 + INT_MIN) { /* INT_MIN is NA_INTEGER */
48 warning(_("integer overflow in 'cumsum'; use 'cumsum(as.numeric(.))'"));
49 break;
50 }
51 is[i] = (int) sum;
52 }
53 return s;
54 }
55
ccumsum(SEXP x,SEXP s)56 static SEXP ccumsum(SEXP x, SEXP s)
57 {
58 Rcomplex sum;
59 sum.r = 0;
60 sum.i = 0;
61 for (R_xlen_t i = 0 ; i < XLENGTH(x) ; i++) {
62 sum.r += COMPLEX(x)[i].r;
63 sum.i += COMPLEX(x)[i].i;
64 COMPLEX(s)[i].r = sum.r;
65 COMPLEX(s)[i].i = sum.i;
66 }
67 return s;
68 }
69
cumprod(SEXP x,SEXP s)70 static SEXP cumprod(SEXP x, SEXP s)
71 {
72 LDOUBLE prod;
73 double *rx = REAL(x), *rs = REAL(s);
74 prod = 1.0;
75 for (R_xlen_t i = 0 ; i < XLENGTH(x) ; i++) {
76 prod *= rx[i]; /* NA and NaN propagated */
77 rs[i] = (double) prod;
78 }
79 return s;
80 }
81
ccumprod(SEXP x,SEXP s)82 static SEXP ccumprod(SEXP x, SEXP s)
83 {
84 Rcomplex prod, tmp;
85 prod.r = 1;
86 prod.i = 0;
87 for (R_xlen_t i = 0 ; i < XLENGTH(x) ; i++) {
88 tmp.r = prod.r;
89 tmp.i = prod.i;
90 prod.r = COMPLEX(x)[i].r * tmp.r - COMPLEX(x)[i].i * tmp.i;
91 prod.i = COMPLEX(x)[i].r * tmp.i + COMPLEX(x)[i].i * tmp.r;
92 COMPLEX(s)[i].r = prod.r;
93 COMPLEX(s)[i].i = prod.i;
94 }
95 return s;
96 }
97
cummax(SEXP x,SEXP s)98 static SEXP cummax(SEXP x, SEXP s)
99 {
100 double max, *rx = REAL(x), *rs = REAL(s);
101 max = R_NegInf;
102 for (R_xlen_t i = 0 ; i < XLENGTH(x) ; i++) {
103 if(ISNAN(rx[i]) || ISNAN(max))
104 max = max + rx[i]; /* propagate NA and NaN */
105 else
106 max = (max > rx[i]) ? max : rx[i];
107 rs[i] = max;
108 }
109 return s;
110 }
111
cummin(SEXP x,SEXP s)112 static SEXP cummin(SEXP x, SEXP s)
113 {
114 double min, *rx = REAL(x), *rs = REAL(s);
115 min = R_PosInf; /* always positive, not NA */
116 for (R_xlen_t i = 0 ; i < XLENGTH(x) ; i++ ) {
117 if (ISNAN(rx[i]) || ISNAN(min))
118 min = min + rx[i]; /* propagate NA and NaN */
119 else
120 min = (min < rx[i]) ? min : rx[i];
121 rs[i] = min;
122 }
123 return s;
124 }
125
icummax(SEXP x,SEXP s)126 static SEXP icummax(SEXP x, SEXP s)
127 {
128 int *ix = INTEGER(x);
129 if(ix[0] == NA_INTEGER)
130 return s; // all NA
131 int *is = INTEGER(s), max = ix[0];
132 is[0] = max;
133 for (R_xlen_t i = 1 ; i < XLENGTH(x) ; i++) {
134 if(ix[i] == NA_INTEGER) break;
135 is[i] = max = (max > ix[i]) ? max : ix[i];
136 }
137 return s;
138 }
139
icummin(SEXP x,SEXP s)140 static SEXP icummin(SEXP x, SEXP s)
141 {
142 int *ix = INTEGER(x), *is = INTEGER(s);
143 int min = ix[0];
144 is[0] = min;
145 for (R_xlen_t i = 1 ; i < XLENGTH(x) ; i++ ) {
146 if(ix[i] == NA_INTEGER) break;
147 is[i] = min = (min < ix[i]) ? min : ix[i];
148 }
149 return s;
150 }
151
do_cum(SEXP call,SEXP op,SEXP args,SEXP env)152 SEXP attribute_hidden do_cum(SEXP call, SEXP op, SEXP args, SEXP env)
153 {
154 SEXP s, t, ans;
155 R_xlen_t i, n;
156 checkArity(op, args);
157 if (DispatchGroup("Math", call, op, args, env, &ans))
158 return ans;
159 if (isComplex(CAR(args))) {
160 t = CAR(args);
161 n = XLENGTH(t);
162 PROTECT(s = allocVector(CPLXSXP, n));
163 setAttrib(s, R_NamesSymbol, getAttrib(t, R_NamesSymbol));
164 UNPROTECT(1);
165 if(n == 0) return s;
166 /* no need to initialize s, ccum* set all elements */
167 switch (PRIMVAL(op) ) {
168 case 1: /* cumsum */
169 return ccumsum(t, s);
170 break;
171 case 2: /* cumprod */
172 return ccumprod(t, s);
173 break;
174 case 3: /* cummax */
175 errorcall(call, _("'cummax' not defined for complex numbers"));
176 break;
177 case 4: /* cummin */
178 errorcall(call, _("'cummin' not defined for complex numbers"));
179 break;
180 default:
181 errorcall(call, "unknown cumxxx function");
182 }
183 } else if( ( isInteger(CAR(args)) || isLogical(CAR(args)) ) &&
184 PRIMVAL(op) != 2) {
185 PROTECT(t = coerceVector(CAR(args), INTSXP));
186 n = XLENGTH(t);
187 PROTECT(s = allocVector(INTSXP, n));
188 setAttrib(s, R_NamesSymbol, getAttrib(t, R_NamesSymbol));
189 if(n == 0) {
190 UNPROTECT(2); /* t, s */
191 return s;
192 }
193 for(i = 0 ; i < n ; i++) INTEGER(s)[i] = NA_INTEGER;
194 switch (PRIMVAL(op) ) {
195 case 1: /* cumsum */
196 ans = icumsum(t,s);
197 break;
198 case 3: /* cummax */
199 ans = icummax(t,s);
200 break;
201 case 4: /* cummin */
202 ans = icummin(t,s);
203 break;
204 default:
205 errorcall(call, _("unknown cumxxx function"));
206 ans = R_NilValue;
207 }
208 UNPROTECT(2); /* t, s */
209 return ans;
210 } else {
211 PROTECT(t = coerceVector(CAR(args), REALSXP));
212 n = XLENGTH(t);
213 PROTECT(s = allocVector(REALSXP, n));
214 setAttrib(s, R_NamesSymbol, getAttrib(t, R_NamesSymbol));
215 UNPROTECT(2);
216 if(n == 0) return s;
217 /* no need to initialize s, cum* set all elements */
218 switch (PRIMVAL(op) ) {
219 case 1: /* cumsum */
220 return cumsum(t,s);
221 break;
222 case 2: /* cumprod */
223 return cumprod(t,s);
224 break;
225 case 3: /* cummax */
226 return cummax(t,s);
227 break;
228 case 4: /* cummin */
229 return cummin(t,s);
230 break;
231 default:
232 errorcall(call, _("unknown cumxxx function"));
233 }
234 }
235 return R_NilValue; /* for -Wall */
236 }
237