1 /*
2 * rangefunc.c: Functions on ranges (data sets).
3 *
4 * Authors:
5 * Morten Welinder <terra@gnome.org>
6 * Andreas J. Guelzow <aguelzow@taliesin.ca>
7 *
8 * This library is free software; you can redistribute it and/or
9 * modify it under the terms of the GNU General Public License as
10 * published by the Free Software Foundation; either version 2 of the
11 * License, or (at your option) version 3.
12 *
13 * This library is distributed in the hope that it will be useful, but
14 * WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 * Library General Public License for more details.
17 *
18 * You should have received a copy of the GNU Library General Public
19 * License along with this library; if not, write to the Free Software
20 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
21 * USA.
22 */
23
24 #include <goffice/goffice-config.h>
25 #include <goffice/goffice.h>
26
27 #include <math.h>
28 #include <stdlib.h>
29 #include <string.h>
30
31 #ifndef DOUBLE
32
33 #define DOUBLE double
34 #ifdef HAVE_LONG_DOUBLE
35 #define LDOUBLE long double
36 #endif
37 #define SUFFIX(_n) _n
38
39 #ifdef GOFFICE_WITH_LONG_DOUBLE
40 #include "go-rangefunc.c"
41 #undef DOUBLE
42 #undef LDOUBLE
43 #undef SUFFIX
44
45 #ifdef HAVE_SUNMATH_H
46 #include <sunmath.h>
47 #endif
48 #define DOUBLE long double
49 #define SUFFIX(_n) _n ## l
50 #endif
51
52 #endif
53
54 #ifndef LDOUBLE
55 #define LDOUBLE DOUBLE
56 #endif
57
58 /* ------------------------------------------------------------------------- */
59
SUFFIX(GOAccumulator)60 static SUFFIX(GOAccumulator)*
61 SUFFIX(sum_helper) (DOUBLE const *xs, int n)
62 {
63 SUFFIX(GOAccumulator) *acc = SUFFIX(go_accumulator_new) ();
64 while (n > 0) {
65 n--;
66 SUFFIX(go_accumulator_add) (acc, xs[n]);
67 }
68 return acc;
69 }
70
71
72 /* Arithmetic sum. */
73 int
SUFFIX(go_range_sum)74 SUFFIX(go_range_sum) (DOUBLE const *xs, int n, DOUBLE *res)
75 {
76 void *state = SUFFIX(go_accumulator_start) ();
77 SUFFIX(GOAccumulator) *acc = SUFFIX(sum_helper) (xs, n);
78 *res = SUFFIX(go_accumulator_value) (acc);
79 SUFFIX(go_accumulator_free) (acc);
80 SUFFIX(go_accumulator_end) (state);
81 return 0;
82 }
83
84 /* Arithmetic sum of squares. */
85 int
SUFFIX(go_range_sumsq)86 SUFFIX(go_range_sumsq) (DOUBLE const *xs, int n, DOUBLE *res)
87 {
88 void *state = SUFFIX(go_accumulator_start) ();
89 SUFFIX(GOAccumulator) *acc = SUFFIX(go_accumulator_new) ();
90 while (n > 0) {
91 DOUBLE x = xs[--n];
92 SUFFIX(GOQuad) q;
93 SUFFIX(go_quad_mul12) (&q, x, x);
94 SUFFIX(go_accumulator_add_quad) (acc, &q);
95 }
96 *res = SUFFIX(go_accumulator_value) (acc);
97 SUFFIX(go_accumulator_free) (acc);
98 SUFFIX(go_accumulator_end) (state);
99 return 0;
100 }
101
102 /* Arithmetic average. */
103 int
SUFFIX(go_range_average)104 SUFFIX(go_range_average) (DOUBLE const *xs, int n, DOUBLE *res)
105 {
106 if (n <= 0)
107 return 1;
108
109 if (SUFFIX(go_range_constant) (xs, n)) {
110 *res = xs[0];
111 return 0;
112 }
113
114 SUFFIX(go_range_sum) (xs, n, res);
115 *res /= n;
116 return 0;
117 }
118
119 /* Minimum element. */
120 int
SUFFIX(go_range_min)121 SUFFIX(go_range_min) (DOUBLE const *xs, int n, DOUBLE *res)
122 {
123 if (n > 0) {
124 DOUBLE min = xs[0];
125 int i;
126
127 for (i = 1; i < n; i++)
128 if (xs[i] < min)
129 min = xs[i];
130 *res = min;
131 return 0;
132 } else
133 return 1;
134 }
135
136 /* Maximum element. */
137 int
SUFFIX(go_range_max)138 SUFFIX(go_range_max) (DOUBLE const *xs, int n, DOUBLE *res)
139 {
140 if (n > 0) {
141 DOUBLE max = xs[0];
142 int i;
143
144 for (i = 1; i < n; i++)
145 if (xs[i] > max)
146 max = xs[i];
147 *res = max;
148 return 0;
149 } else
150 return 1;
151 }
152
153 /* Maximum absolute element. */
154 int
SUFFIX(go_range_maxabs)155 SUFFIX(go_range_maxabs) (DOUBLE const *xs, int n, DOUBLE *res)
156 {
157 if (n > 0) {
158 DOUBLE max = SUFFIX(fabs) (xs[0]);
159 int i;
160
161 for (i = 1; i < n; i++) {
162 DOUBLE xabs = SUFFIX(fabs) (xs[i]);
163 if (xabs > max)
164 max = xabs;
165 }
166 *res = max;
167 return 0;
168 } else
169 return 1;
170 }
171
172 /* Sum of square deviations from mean. */
173 int
SUFFIX(go_range_devsq)174 SUFFIX(go_range_devsq) (DOUBLE const *xs, int n, DOUBLE *res)
175 {
176 if (SUFFIX(go_range_constant) (xs, n))
177 *res = 0;
178 else {
179 void *state;
180 SUFFIX(GOAccumulator) *acc;
181 DOUBLE sum, sumh, suml;
182 SUFFIX(GOQuad) qavg, qtmp, qn;
183
184 state = SUFFIX(go_accumulator_start) ();
185 acc = SUFFIX(sum_helper) (xs, n);
186 sumh = SUFFIX(go_accumulator_value) (acc);
187 SUFFIX(go_accumulator_add) (acc, -sumh);
188 suml = SUFFIX(go_accumulator_value) (acc);
189
190 SUFFIX(go_quad_init) (&qavg, sumh);
191 SUFFIX(go_quad_init) (&qtmp, suml);
192 SUFFIX(go_quad_add) (&qavg, &qavg, &qtmp);
193 SUFFIX(go_range_sum) (xs, n, &sum);
194 SUFFIX(go_quad_init) (&qn, n);
195 SUFFIX(go_quad_div) (&qavg, &qavg, &qn);
196 /*
197 * The just-generated qavg isn't exact, but should be
198 * good to near-GOQuad precision. We hope that's good
199 * enough.
200 */
201
202 SUFFIX(go_accumulator_clear) (acc);
203 while (n > 0) {
204 SUFFIX(GOQuad) q;
205 n--;
206 SUFFIX(go_quad_init) (&q, xs[n]);
207 SUFFIX(go_quad_sub) (&q, &q, &qavg);
208 SUFFIX(go_quad_mul) (&q, &q, &q);
209 SUFFIX(go_accumulator_add_quad) (acc, &q);
210 }
211 *res = SUFFIX(go_accumulator_value) (acc);
212 SUFFIX(go_accumulator_free) (acc);
213 SUFFIX(go_accumulator_end) (state);
214 }
215 return 0;
216 }
217
218 static int
SUFFIX(float_compare)219 SUFFIX(float_compare) (DOUBLE const *a, DOUBLE const *b)
220 {
221 if (*a < *b)
222 return -1;
223 else if (*a == *b)
224 return 0;
225 else
226 return 1;
227 }
228
229 DOUBLE *
SUFFIX(go_range_sort)230 SUFFIX(go_range_sort) (DOUBLE const *xs, int n)
231 {
232 if (n <= 0)
233 return NULL;
234 else {
235 DOUBLE *ys = g_new (DOUBLE, n);
236 memcpy (ys, xs, n * sizeof (DOUBLE));
237 qsort (ys, n, sizeof (ys[0]), (int (*) (const void *, const void *))&SUFFIX(float_compare));
238 return ys;
239 }
240 }
241 /* This requires sorted data. */
242 int
SUFFIX(go_range_fractile_inter_sorted)243 SUFFIX(go_range_fractile_inter_sorted) (DOUBLE const *xs, int n, DOUBLE *res, DOUBLE f)
244 {
245 DOUBLE fpos, residual;
246 int pos;
247
248 if (n <= 0 || f < 0.0 || f > 1.0)
249 return 1;
250
251 fpos = (n - 1) * f;
252 pos = (int)fpos;
253 residual = fpos - pos;
254
255 if (residual == 0.0 || pos + 1 >= n)
256 *res = xs[pos];
257 else
258 *res = (1 - residual) * xs[pos] + residual * xs[pos + 1];
259
260 return 0;
261 }
262
263 int
SUFFIX(go_range_fractile_inter)264 SUFFIX(go_range_fractile_inter) (DOUBLE const *xs, int n, DOUBLE *res, DOUBLE f)
265 {
266 DOUBLE *ys = SUFFIX(go_range_sort) (xs, n);
267 int error = SUFFIX(go_range_fractile_inter_sorted) (ys, n, res, f);
268 g_free (ys);
269 return error;
270 }
271
272 int
SUFFIX(go_range_fractile_inter_nonconst)273 SUFFIX(go_range_fractile_inter_nonconst) (DOUBLE *xs, int n, DOUBLE *res, DOUBLE f)
274 {
275 qsort (xs, n, sizeof (xs[0]), (int (*) (const void *, const void *))&SUFFIX(float_compare));
276 return SUFFIX(go_range_fractile_inter_sorted) (xs, n, res, f);
277 }
278
279 int
SUFFIX(go_range_median_inter)280 SUFFIX(go_range_median_inter) (DOUBLE const *xs, int n, DOUBLE *res)
281 {
282 return SUFFIX(go_range_fractile_inter) (xs, n, res, 0.5);
283 }
284
285 int
SUFFIX(go_range_median_inter_sorted)286 SUFFIX(go_range_median_inter_sorted) (DOUBLE const *xs, int n, DOUBLE *res)
287 {
288 return SUFFIX(go_range_fractile_inter_sorted) (xs, n, res, 0.5);
289 }
290
291 int
SUFFIX(go_range_median_inter_nonconst)292 SUFFIX(go_range_median_inter_nonconst) (DOUBLE *xs, int n, DOUBLE *res)
293 {
294 return SUFFIX(go_range_fractile_inter_nonconst) (xs, n, res, 0.5);
295 }
296
297 int
SUFFIX(go_range_constant)298 SUFFIX(go_range_constant) (DOUBLE const *xs, int n)
299 {
300 int i;
301 for (i = 1; i < n; i++)
302 if (xs[0] != xs[i])
303 return 0;
304 return 1;
305 }
306
307 int
SUFFIX(go_range_increasing)308 SUFFIX(go_range_increasing) (DOUBLE const *xs, int n)
309 {
310 int i = 0;
311 DOUBLE last;
312 g_return_val_if_fail (n == 0 || xs != NULL, 0);
313 while ( i < n && isnan (xs[i]))
314 i++;
315 if (i == n)
316 return 0;
317 last = xs[i];
318 for (i = i + 1; i < n; i++) {
319 if (isnan (xs[i]))
320 continue;
321 if (last >= xs[i])
322 return 0;
323 last = xs[i];
324 }
325 return 1;
326 }
327
328 int
SUFFIX(go_range_decreasing)329 SUFFIX(go_range_decreasing) (DOUBLE const *xs, int n)
330 {
331 int i = 0;
332 DOUBLE last;
333 g_return_val_if_fail (n == 0 || xs != NULL, 0);
334 while ( i < n && isnan (xs[i]))
335 i++;
336 if (i == n)
337 return 0;
338 last = xs[i];
339 for (i = i + 1; i < n; i++) {
340 if (isnan (xs[i]))
341 continue;
342 if (last <= xs[i])
343 return 0;
344 last = xs[i];
345 }
346 return 1;
347 }
348
349 int
SUFFIX(go_range_vary_uniformly)350 SUFFIX(go_range_vary_uniformly) (DOUBLE const *xs, int n)
351 {
352 return SUFFIX(go_range_increasing) (xs, n) || SUFFIX(go_range_decreasing) (xs, n);
353 }
354