rrd_graph_rpn: add a MEDIAN operator. -- Aaron Gallagher <_@habnab.it>
[rrdtool.git] / src / rrd_rpncalc.c
1 /****************************************************************************
2  * RRDtool 1.4.3  Copyright by Tobi Oetiker, 1997-2010
3  ****************************************************************************
4  * rrd_rpncalc.c  RPN calculator functions
5  ****************************************************************************/
6
7 #include <limits.h>
8 #include <locale.h>
9 #include <stdlib.h>
10
11 #include "rrd_tool.h"
12 #include "rrd_rpncalc.h"
13 // #include "rrd_graph.h"
14
15 short     addop2str(
16     enum op_en op,
17     enum op_en op_type,
18     char *op_str,
19     char **result_str,
20     unsigned short *offset);
21 int       tzoffset(
22     time_t);            /* used to implement LTIME */
23
24 short rpn_compact(
25     rpnp_t *rpnp,
26     rpn_cdefds_t **rpnc,
27     short *count)
28 {
29     short     i;
30
31     *count = 0;
32     /* count the number of rpn nodes */
33     while (rpnp[*count].op != OP_END)
34         (*count)++;
35     if (++(*count) > DS_CDEF_MAX_RPN_NODES) {
36         rrd_set_error("Maximum %d RPN nodes permitted. Got %d RPN nodes at present.",
37                       DS_CDEF_MAX_RPN_NODES-1,(*count)-1);
38         return -1;
39     }
40
41     /* allocate memory */
42     *rpnc = (rpn_cdefds_t *) calloc(*count, sizeof(rpn_cdefds_t));
43     for (i = 0; rpnp[i].op != OP_END; i++) {
44         (*rpnc)[i].op = (char) rpnp[i].op;
45         if (rpnp[i].op == OP_NUMBER) {
46             /* rpnp.val is a double, rpnc.val is a short */
47             double    temp = floor(rpnp[i].val);
48
49             if (temp < SHRT_MIN || temp > SHRT_MAX) {
50                 rrd_set_error
51                     ("constants must be integers in the interval (%d, %d)",
52                      SHRT_MIN, SHRT_MAX);
53                 free(*rpnc);
54                 return -1;
55             }
56             (*rpnc)[i].val = (short) temp;
57         } else if (rpnp[i].op == OP_VARIABLE || rpnp[i].op == OP_PREV_OTHER) {
58             (*rpnc)[i].val = (short) rpnp[i].ptr;
59         }
60     }
61     /* terminate the sequence */
62     (*rpnc)[(*count) - 1].op = OP_END;
63     return 0;
64 }
65
66 rpnp_t   *rpn_expand(
67     rpn_cdefds_t *rpnc)
68 {
69     short     i;
70     rpnp_t   *rpnp;
71
72     /* DS_CDEF_MAX_RPN_NODES is small, so at the expense of some wasted
73      * memory we avoid any reallocs */
74     rpnp = (rpnp_t *) calloc(DS_CDEF_MAX_RPN_NODES, sizeof(rpnp_t));
75     if (rpnp == NULL) {
76         rrd_set_error("failed allocating rpnp array");
77         return NULL;
78     }
79     for (i = 0; rpnc[i].op != OP_END; ++i) {
80         rpnp[i].op = (enum op_en)rpnc[i].op;
81         if (rpnp[i].op == OP_NUMBER) {
82             rpnp[i].val = (double) rpnc[i].val;
83         } else if (rpnp[i].op == OP_VARIABLE || rpnp[i].op == OP_PREV_OTHER) {
84             rpnp[i].ptr = (long) rpnc[i].val;
85         }
86     }
87     /* terminate the sequence */
88     rpnp[i].op = OP_END;
89     return rpnp;
90 }
91
92 /* rpn_compact2str: convert a compact sequence of RPN operator nodes back
93  * into a CDEF string. This function is used by rrd_dump.
94  * arguments:
95  *  rpnc: an array of compact RPN operator nodes
96  *  ds_def: a pointer to the data source definition section of an RRD header
97  *   for lookup of data source names by index
98  *  str: out string, memory is allocated by the function, must be freed by the
99  *   the caller */
100 void rpn_compact2str(
101     rpn_cdefds_t *rpnc,
102     ds_def_t *ds_def,
103     char **str)
104 {
105     unsigned short i, offset = 0;
106     char      buffer[7];    /* short as a string */
107
108     for (i = 0; rpnc[i].op != OP_END; i++) {
109         if (i > 0)
110             (*str)[offset++] = ',';
111
112 #define add_op(VV,VVV) \
113           if (addop2str((enum op_en)(rpnc[i].op), VV, VVV, str, &offset) == 1) continue;
114
115         if (rpnc[i].op == OP_NUMBER) {
116             /* convert a short into a string */
117 #if defined(_WIN32) && !defined(__CYGWIN__) && !defined(__CYGWIN32__)
118             _itoa(rpnc[i].val, buffer, 10);
119 #else
120             sprintf(buffer, "%d", rpnc[i].val);
121 #endif
122             add_op(OP_NUMBER, buffer)
123         }
124
125         if (rpnc[i].op == OP_VARIABLE) {
126             char     *ds_name = ds_def[rpnc[i].val].ds_nam;
127
128             add_op(OP_VARIABLE, ds_name)
129         }
130
131         if (rpnc[i].op == OP_PREV_OTHER) {
132             char     *ds_name = ds_def[rpnc[i].val].ds_nam;
133
134             add_op(OP_VARIABLE, ds_name)
135         }
136 #undef add_op
137
138 #define add_op(VV,VVV) \
139           if (addop2str((enum op_en)rpnc[i].op, VV, #VVV, str, &offset) == 1) continue;
140
141         add_op(OP_ADD, +)
142             add_op(OP_SUB, -)
143             add_op(OP_MUL, *)
144             add_op(OP_DIV, /)
145             add_op(OP_MOD, %)
146             add_op(OP_SIN, SIN)
147             add_op(OP_COS, COS)
148             add_op(OP_LOG, LOG)
149             add_op(OP_FLOOR, FLOOR)
150             add_op(OP_CEIL, CEIL)
151             add_op(OP_EXP, EXP)
152             add_op(OP_DUP, DUP)
153             add_op(OP_EXC, EXC)
154             add_op(OP_POP, POP)
155             add_op(OP_LT, LT)
156             add_op(OP_LE, LE)
157             add_op(OP_GT, GT)
158             add_op(OP_GE, GE)
159             add_op(OP_EQ, EQ)
160             add_op(OP_IF, IF)
161             add_op(OP_MIN, MIN)
162             add_op(OP_MAX, MAX)
163             add_op(OP_LIMIT, LIMIT)
164             add_op(OP_UNKN, UNKN)
165             add_op(OP_UN, UN)
166             add_op(OP_NEGINF, NEGINF)
167             add_op(OP_NE, NE)
168             add_op(OP_PREV, PREV)
169             add_op(OP_INF, INF)
170             add_op(OP_ISINF, ISINF)
171             add_op(OP_NOW, NOW)
172             add_op(OP_LTIME, LTIME)
173             add_op(OP_TIME, TIME)
174             add_op(OP_ATAN2, ATAN2)
175             add_op(OP_ATAN, ATAN)
176             add_op(OP_SQRT, SQRT)
177             add_op(OP_SORT, SORT)
178             add_op(OP_REV, REV)
179             add_op(OP_TREND, TREND)
180             add_op(OP_TRENDNAN, TRENDNAN)
181             add_op(OP_PREDICT, PREDICT)
182             add_op(OP_PREDICTSIGMA, PREDICTSIGMA)
183             add_op(OP_RAD2DEG, RAD2DEG)
184             add_op(OP_DEG2RAD, DEG2RAD)
185             add_op(OP_AVG, AVG)
186             add_op(OP_ABS, ABS)
187             add_op(OP_ADDNAN, ADDNAN)
188 #undef add_op
189     }
190     (*str)[offset] = '\0';
191
192 }
193
194 short addop2str(
195     enum op_en op,
196     enum op_en op_type,
197     char *op_str,
198     char **result_str,
199     unsigned short *offset)
200 {
201     if (op == op_type) {
202         short     op_len;
203
204         op_len = strlen(op_str);
205         *result_str = (char *) rrd_realloc(*result_str,
206                                            (op_len + 1 +
207                                             *offset) * sizeof(char));
208         if (*result_str == NULL) {
209             rrd_set_error("failed to alloc memory in addop2str");
210             return -1;
211         }
212         strncpy(&((*result_str)[*offset]), op_str, op_len);
213         *offset += op_len;
214         return 1;
215     }
216     return 0;
217 }
218
219 void parseCDEF_DS(
220     const char *def,
221     rrd_t *rrd,
222     int ds_idx)
223 {
224     rpnp_t   *rpnp = NULL;
225     rpn_cdefds_t *rpnc = NULL;
226     short     count, i;
227
228     rpnp = rpn_parse((void *) rrd, def, &lookup_DS);
229     if (rpnp == NULL) {
230         rrd_set_error("failed to parse computed data source");
231         return;
232     }
233     /* Check for OP nodes not permitted in COMPUTE DS.
234      * Moved this check from within rpn_compact() because it really is
235      * COMPUTE DS specific. This is less efficient, but creation doesn't
236      * occur too often. */
237     for (i = 0; rpnp[i].op != OP_END; i++) {
238         if (rpnp[i].op == OP_TIME || rpnp[i].op == OP_LTIME ||
239             rpnp[i].op == OP_PREV || rpnp[i].op == OP_COUNT ||
240             rpnp[i].op == OP_TREND || rpnp[i].op == OP_TRENDNAN ||
241             rpnp[i].op == OP_PREDICT || rpnp[i].op ==  OP_PREDICTSIGMA ) {
242             rrd_set_error
243                 ("operators TIME, LTIME, PREV COUNT TREND TRENDNAN PREDICT PREDICTSIGMA are not supported with DS COMPUTE");
244             free(rpnp);
245             return;
246         }
247     }
248     if (rpn_compact(rpnp, &rpnc, &count) == -1) {
249         free(rpnp);
250         return;
251     }
252     /* copy the compact rpn representation over the ds_def par array */
253     memcpy((void *) &(rrd->ds_def[ds_idx].par[DS_cdef]),
254            (void *) rpnc, count * sizeof(rpn_cdefds_t));
255     free(rpnp);
256     free(rpnc);
257 }
258
259 /* lookup a data source name in the rrd struct and return the index,
260  * should use ds_match() here except:
261  * (1) need a void * pointer to the rrd
262  * (2) error handling is left to the caller
263  */
264 long lookup_DS(
265     void *rrd_vptr,
266     char *ds_name)
267 {
268     unsigned int i;
269     rrd_t    *rrd;
270
271     rrd = (rrd_t *) rrd_vptr;
272
273     for (i = 0; i < rrd->stat_head->ds_cnt; ++i) {
274         if (strcmp(ds_name, rrd->ds_def[i].ds_nam) == 0)
275             return i;
276     }
277     /* the caller handles a bad data source name in the rpn string */
278     return -1;
279 }
280
281 /* rpn_parse : parse a string and generate a rpnp array; modified
282  * str2rpn() originally included in rrd_graph.c
283  * arguments:
284  * key_hash: a transparent argument passed to lookup(); conceptually this
285  *    is a hash object for lookup of a numeric key given a variable name
286  * expr: the string RPN expression, including variable names
287  * lookup(): a function that retrieves a numeric key given a variable name
288  */
289 rpnp_t   *rpn_parse(
290     void *key_hash,
291     const char *const expr_const,
292     long      (*lookup) (void *,
293                          char *))
294 {
295     int       pos = 0;
296     char     *expr;
297     long      steps = -1;
298     rpnp_t   *rpnp;
299     char      vname[MAX_VNAME_LEN + 10];
300     char     *old_locale;
301
302     old_locale = setlocale(LC_NUMERIC, NULL);
303     setlocale(LC_NUMERIC, "C");
304
305     rpnp = NULL;
306     expr = (char *) expr_const;
307
308     while (*expr) {
309         if ((rpnp = (rpnp_t *) rrd_realloc(rpnp, (++steps + 2) *
310                                            sizeof(rpnp_t))) == NULL) {
311             setlocale(LC_NUMERIC, old_locale);
312             return NULL;
313         }
314
315         else if ((sscanf(expr, "%lf%n", &rpnp[steps].val, &pos) == 1)
316                  && (expr[pos] == ',')) {
317             rpnp[steps].op = OP_NUMBER;
318             expr += pos;
319         }
320 #define match_op(VV,VVV) \
321         else if (strncmp(expr, #VVV, strlen(#VVV))==0 && ( expr[strlen(#VVV)] == ',' || expr[strlen(#VVV)] == '\0' )){ \
322             rpnp[steps].op = VV; \
323             expr+=strlen(#VVV); \
324         }
325
326 #define match_op_param(VV,VVV) \
327         else if (sscanf(expr, #VVV "(" DEF_NAM_FMT ")",vname) == 1) { \
328           int length = 0; \
329           if ((length = strlen(#VVV)+strlen(vname)+2, \
330               expr[length] == ',' || expr[length] == '\0') ) { \
331              rpnp[steps].op = VV; \
332              rpnp[steps].ptr = (*lookup)(key_hash,vname); \
333              if (rpnp[steps].ptr < 0) { \
334                            free(rpnp); \
335                            return NULL; \
336                          } else expr+=length; \
337           } \
338         }
339
340         match_op(OP_ADD, +)
341             match_op(OP_SUB, -)
342             match_op(OP_MUL, *)
343             match_op(OP_DIV, /)
344             match_op(OP_MOD, %)
345             match_op(OP_SIN, SIN)
346             match_op(OP_COS, COS)
347             match_op(OP_LOG, LOG)
348             match_op(OP_FLOOR, FLOOR)
349             match_op(OP_CEIL, CEIL)
350             match_op(OP_EXP, EXP)
351             match_op(OP_DUP, DUP)
352             match_op(OP_EXC, EXC)
353             match_op(OP_POP, POP)
354             match_op(OP_LTIME, LTIME)
355             match_op(OP_LT, LT)
356             match_op(OP_LE, LE)
357             match_op(OP_GT, GT)
358             match_op(OP_GE, GE)
359             match_op(OP_EQ, EQ)
360             match_op(OP_IF, IF)
361             match_op(OP_MIN, MIN)
362             match_op(OP_MAX, MAX)
363             match_op(OP_LIMIT, LIMIT)
364             /* order is important here ! .. match longest first */
365             match_op(OP_UNKN, UNKN)
366             match_op(OP_UN, UN)
367             match_op(OP_NEGINF, NEGINF)
368             match_op(OP_NE, NE)
369             match_op(OP_COUNT, COUNT)
370             match_op_param(OP_PREV_OTHER, PREV)
371             match_op(OP_PREV, PREV)
372             match_op(OP_INF, INF)
373             match_op(OP_ISINF, ISINF)
374             match_op(OP_NOW, NOW)
375             match_op(OP_TIME, TIME)
376             match_op(OP_ATAN2, ATAN2)
377             match_op(OP_ATAN, ATAN)
378             match_op(OP_SQRT, SQRT)
379             match_op(OP_SORT, SORT)
380             match_op(OP_REV, REV)
381             match_op(OP_TREND, TREND)
382             match_op(OP_TRENDNAN, TRENDNAN)
383             match_op(OP_PREDICT, PREDICT)
384             match_op(OP_PREDICTSIGMA, PREDICTSIGMA)
385             match_op(OP_RAD2DEG, RAD2DEG)
386             match_op(OP_DEG2RAD, DEG2RAD)
387             match_op(OP_AVG, AVG)
388             match_op(OP_ABS, ABS)
389             match_op(OP_ADDNAN, ADDNAN)
390             match_op(OP_MEDIAN, MEDIAN)
391 #undef match_op
392             else if ((sscanf(expr, DEF_NAM_FMT "%n", vname, &pos) == 1)
393                      && ((rpnp[steps].ptr = (*lookup) (key_hash, vname)) !=
394                          -1)) {
395             rpnp[steps].op = OP_VARIABLE;
396             expr += pos;
397         }
398
399         else {
400             setlocale(LC_NUMERIC, old_locale);
401             free(rpnp);
402             return NULL;
403         }
404
405         if (*expr == 0)
406             break;
407         if (*expr == ',')
408             expr++;
409         else {
410             setlocale(LC_NUMERIC, old_locale);
411             free(rpnp);
412             return NULL;
413         }
414     }
415     rpnp[steps + 1].op = OP_END;
416     setlocale(LC_NUMERIC, old_locale);
417     return rpnp;
418 }
419
420 void rpnstack_init(
421     rpnstack_t *rpnstack)
422 {
423     rpnstack->s = NULL;
424     rpnstack->dc_stacksize = 0;
425     rpnstack->dc_stackblock = 100;
426 }
427
428 void rpnstack_free(
429     rpnstack_t *rpnstack)
430 {
431     if (rpnstack->s != NULL)
432         free(rpnstack->s);
433     rpnstack->dc_stacksize = 0;
434 }
435
436 static int rpn_compare_double(
437     const void *x,
438     const void *y)
439 {
440     double    diff = *((const double *) x) - *((const double *) y);
441
442     return (diff < 0) ? -1 : (diff > 0) ? 1 : 0;
443 }
444
445 /* rpn_calc: run the RPN calculator; also performs variable substitution;
446  * moved and modified from data_calc() originally included in rrd_graph.c 
447  * arguments:
448  * rpnp : an array of RPN operators (including variable references)
449  * rpnstack : the initialized stack
450  * data_idx : when data_idx is a multiple of rpnp.step, the rpnp.data pointer
451  *   is advanced by rpnp.ds_cnt; used only for variable substitution
452  * output : an array of output values; OP_PREV assumes this array contains
453  *   the "previous" value at index position output_idx-1; the definition of
454  *   "previous" depends on the calling environment
455  * output_idx : an index into the output array in which to store the output
456  *   of the RPN calculator
457  * returns: -1 if the computation failed (also calls rrd_set_error)
458  *           0 on success
459  */
460 short rpn_calc(
461     rpnp_t *rpnp,
462     rpnstack_t *rpnstack,
463     long data_idx,
464     rrd_value_t *output,
465     int output_idx)
466 {
467     int       rpi;
468     long      stptr = -1;
469
470     /* process each op from the rpn in turn */
471     for (rpi = 0; rpnp[rpi].op != OP_END; rpi++) {
472         /* allocate or grow the stack */
473         if (stptr + 5 > rpnstack->dc_stacksize) {
474             /* could move this to a separate function */
475             rpnstack->dc_stacksize += rpnstack->dc_stackblock;
476             rpnstack->s = (double*)rrd_realloc(rpnstack->s,
477                                       (rpnstack->dc_stacksize) *
478                                       sizeof(*(rpnstack->s)));
479             if (rpnstack->s == NULL) {
480                 rrd_set_error("RPN stack overflow");
481                 return -1;
482             }
483         }
484 #define stackunderflow(MINSIZE)                         \
485         if(stptr<MINSIZE){                              \
486             rrd_set_error("RPN stack underflow");       \
487             return -1;                                  \
488         }
489
490         switch (rpnp[rpi].op) {
491         case OP_NUMBER:
492             rpnstack->s[++stptr] = rpnp[rpi].val;
493             break;
494         case OP_VARIABLE:
495         case OP_PREV_OTHER:
496             /* Sanity check: VDEFs shouldn't make it here */
497             if (rpnp[rpi].ds_cnt == 0) {
498                 rrd_set_error("VDEF made it into rpn_calc... aborting");
499                 return -1;
500             } else {
501                 /* make sure we pull the correct value from
502                  * the *.data array. Adjust the pointer into
503                  * the array acordingly. Advance the ptr one
504                  * row in the rra (skip over non-relevant
505                  * data sources)
506                  */
507                 if (rpnp[rpi].op == OP_VARIABLE) {
508                     rpnstack->s[++stptr] = *(rpnp[rpi].data);
509                 } else {
510                     if ((output_idx) <= 0) {
511                         rpnstack->s[++stptr] = DNAN;
512                     } else {
513                         rpnstack->s[++stptr] =
514                             *(rpnp[rpi].data - rpnp[rpi].ds_cnt);
515                     }
516
517                 }
518                 if (data_idx % rpnp[rpi].step == 0) {
519                     rpnp[rpi].data += rpnp[rpi].ds_cnt;
520                 }
521             }
522             break;
523         case OP_COUNT:
524             rpnstack->s[++stptr] = (output_idx + 1);    /* Note: Counter starts at 1 */
525             break;
526         case OP_PREV:
527             if ((output_idx) <= 0) {
528                 rpnstack->s[++stptr] = DNAN;
529             } else {
530                 rpnstack->s[++stptr] = output[output_idx - 1];
531             }
532             break;
533         case OP_UNKN:
534             rpnstack->s[++stptr] = DNAN;
535             break;
536         case OP_INF:
537             rpnstack->s[++stptr] = DINF;
538             break;
539         case OP_NEGINF:
540             rpnstack->s[++stptr] = -DINF;
541             break;
542         case OP_NOW:
543             rpnstack->s[++stptr] = (double) time(NULL);
544             break;
545         case OP_TIME:
546             /* HACK: this relies on the data_idx being the time,
547              ** which the within-function scope is unaware of */
548             rpnstack->s[++stptr] = (double) data_idx;
549             break;
550         case OP_LTIME:
551             rpnstack->s[++stptr] =
552                 (double) tzoffset(data_idx) + (double) data_idx;
553             break;
554         case OP_ADD:
555             stackunderflow(1);
556             rpnstack->s[stptr - 1] = rpnstack->s[stptr - 1]
557                 + rpnstack->s[stptr];
558             stptr--;
559             break;
560         case OP_ADDNAN:
561             stackunderflow(1);
562             if (isnan(rpnstack->s[stptr - 1])) {
563                 rpnstack->s[stptr - 1] = rpnstack->s[stptr];
564             } else if (isnan(rpnstack->s[stptr])) {
565                 /* NOOP */
566                 /* rpnstack->s[stptr - 1] = rpnstack->s[stptr - 1]; */
567             } else {
568                 rpnstack->s[stptr - 1] = rpnstack->s[stptr - 1]
569                     + rpnstack->s[stptr];
570             }
571
572             stptr--;
573             break;
574         case OP_SUB:
575             stackunderflow(1);
576             rpnstack->s[stptr - 1] = rpnstack->s[stptr - 1]
577                 - rpnstack->s[stptr];
578             stptr--;
579             break;
580         case OP_MUL:
581             stackunderflow(1);
582             rpnstack->s[stptr - 1] = (rpnstack->s[stptr - 1])
583                 * (rpnstack->s[stptr]);
584             stptr--;
585             break;
586         case OP_DIV:
587             stackunderflow(1);
588             rpnstack->s[stptr - 1] = rpnstack->s[stptr - 1]
589                 / rpnstack->s[stptr];
590             stptr--;
591             break;
592         case OP_MOD:
593             stackunderflow(1);
594             rpnstack->s[stptr - 1] = fmod(rpnstack->s[stptr - 1]
595                                           , rpnstack->s[stptr]);
596             stptr--;
597             break;
598         case OP_SIN:
599             stackunderflow(0);
600             rpnstack->s[stptr] = sin(rpnstack->s[stptr]);
601             break;
602         case OP_ATAN:
603             stackunderflow(0);
604             rpnstack->s[stptr] = atan(rpnstack->s[stptr]);
605             break;
606         case OP_RAD2DEG:
607             stackunderflow(0);
608             rpnstack->s[stptr] = 57.29577951 * rpnstack->s[stptr];
609             break;
610         case OP_DEG2RAD:
611             stackunderflow(0);
612             rpnstack->s[stptr] = 0.0174532952 * rpnstack->s[stptr];
613             break;
614         case OP_ATAN2:
615             stackunderflow(1);
616             rpnstack->s[stptr - 1] = atan2(rpnstack->s[stptr - 1],
617                                            rpnstack->s[stptr]);
618             stptr--;
619             break;
620         case OP_COS:
621             stackunderflow(0);
622             rpnstack->s[stptr] = cos(rpnstack->s[stptr]);
623             break;
624         case OP_CEIL:
625             stackunderflow(0);
626             rpnstack->s[stptr] = ceil(rpnstack->s[stptr]);
627             break;
628         case OP_FLOOR:
629             stackunderflow(0);
630             rpnstack->s[stptr] = floor(rpnstack->s[stptr]);
631             break;
632         case OP_LOG:
633             stackunderflow(0);
634             rpnstack->s[stptr] = log(rpnstack->s[stptr]);
635             break;
636         case OP_DUP:
637             stackunderflow(0);
638             rpnstack->s[stptr + 1] = rpnstack->s[stptr];
639             stptr++;
640             break;
641         case OP_POP:
642             stackunderflow(0);
643             stptr--;
644             break;
645         case OP_EXC:
646             stackunderflow(1);
647             {
648                 double    dummy;
649
650                 dummy = rpnstack->s[stptr];
651                 rpnstack->s[stptr] = rpnstack->s[stptr - 1];
652                 rpnstack->s[stptr - 1] = dummy;
653             }
654             break;
655         case OP_EXP:
656             stackunderflow(0);
657             rpnstack->s[stptr] = exp(rpnstack->s[stptr]);
658             break;
659         case OP_LT:
660             stackunderflow(1);
661             if (isnan(rpnstack->s[stptr - 1]));
662             else if (isnan(rpnstack->s[stptr]))
663                 rpnstack->s[stptr - 1] = rpnstack->s[stptr];
664             else
665                 rpnstack->s[stptr - 1] = rpnstack->s[stptr - 1] <
666                     rpnstack->s[stptr] ? 1.0 : 0.0;
667             stptr--;
668             break;
669         case OP_LE:
670             stackunderflow(1);
671             if (isnan(rpnstack->s[stptr - 1]));
672             else if (isnan(rpnstack->s[stptr]))
673                 rpnstack->s[stptr - 1] = rpnstack->s[stptr];
674             else
675                 rpnstack->s[stptr - 1] = rpnstack->s[stptr - 1] <=
676                     rpnstack->s[stptr] ? 1.0 : 0.0;
677             stptr--;
678             break;
679         case OP_GT:
680             stackunderflow(1);
681             if (isnan(rpnstack->s[stptr - 1]));
682             else if (isnan(rpnstack->s[stptr]))
683                 rpnstack->s[stptr - 1] = rpnstack->s[stptr];
684             else
685                 rpnstack->s[stptr - 1] = rpnstack->s[stptr - 1] >
686                     rpnstack->s[stptr] ? 1.0 : 0.0;
687             stptr--;
688             break;
689         case OP_GE:
690             stackunderflow(1);
691             if (isnan(rpnstack->s[stptr - 1]));
692             else if (isnan(rpnstack->s[stptr]))
693                 rpnstack->s[stptr - 1] = rpnstack->s[stptr];
694             else
695                 rpnstack->s[stptr - 1] = rpnstack->s[stptr - 1] >=
696                     rpnstack->s[stptr] ? 1.0 : 0.0;
697             stptr--;
698             break;
699         case OP_NE:
700             stackunderflow(1);
701             if (isnan(rpnstack->s[stptr - 1]));
702             else if (isnan(rpnstack->s[stptr]))
703                 rpnstack->s[stptr - 1] = rpnstack->s[stptr];
704             else
705                 rpnstack->s[stptr - 1] = rpnstack->s[stptr - 1] ==
706                     rpnstack->s[stptr] ? 0.0 : 1.0;
707             stptr--;
708             break;
709         case OP_EQ:
710             stackunderflow(1);
711             if (isnan(rpnstack->s[stptr - 1]));
712             else if (isnan(rpnstack->s[stptr]))
713                 rpnstack->s[stptr - 1] = rpnstack->s[stptr];
714             else
715                 rpnstack->s[stptr - 1] = rpnstack->s[stptr - 1] ==
716                     rpnstack->s[stptr] ? 1.0 : 0.0;
717             stptr--;
718             break;
719         case OP_IF:
720             stackunderflow(2);
721             rpnstack->s[stptr - 2] = (isnan(rpnstack->s[stptr - 2])
722                                       || rpnstack->s[stptr - 2] ==
723                                       0.0) ? rpnstack->s[stptr] : rpnstack->
724                 s[stptr - 1];
725             stptr--;
726             stptr--;
727             break;
728         case OP_MIN:
729             stackunderflow(1);
730             if (isnan(rpnstack->s[stptr - 1]));
731             else if (isnan(rpnstack->s[stptr]))
732                 rpnstack->s[stptr - 1] = rpnstack->s[stptr];
733             else if (rpnstack->s[stptr - 1] > rpnstack->s[stptr])
734                 rpnstack->s[stptr - 1] = rpnstack->s[stptr];
735             stptr--;
736             break;
737         case OP_MAX:
738             stackunderflow(1);
739             if (isnan(rpnstack->s[stptr - 1]));
740             else if (isnan(rpnstack->s[stptr]))
741                 rpnstack->s[stptr - 1] = rpnstack->s[stptr];
742             else if (rpnstack->s[stptr - 1] < rpnstack->s[stptr])
743                 rpnstack->s[stptr - 1] = rpnstack->s[stptr];
744             stptr--;
745             break;
746         case OP_LIMIT:
747             stackunderflow(2);
748             if (isnan(rpnstack->s[stptr - 2]));
749             else if (isnan(rpnstack->s[stptr - 1]))
750                 rpnstack->s[stptr - 2] = rpnstack->s[stptr - 1];
751             else if (isnan(rpnstack->s[stptr]))
752                 rpnstack->s[stptr - 2] = rpnstack->s[stptr];
753             else if (rpnstack->s[stptr - 2] < rpnstack->s[stptr - 1])
754                 rpnstack->s[stptr - 2] = DNAN;
755             else if (rpnstack->s[stptr - 2] > rpnstack->s[stptr])
756                 rpnstack->s[stptr - 2] = DNAN;
757             stptr -= 2;
758             break;
759         case OP_UN:
760             stackunderflow(0);
761             rpnstack->s[stptr] = isnan(rpnstack->s[stptr]) ? 1.0 : 0.0;
762             break;
763         case OP_ISINF:
764             stackunderflow(0);
765             rpnstack->s[stptr] = isinf(rpnstack->s[stptr]) ? 1.0 : 0.0;
766             break;
767         case OP_SQRT:
768             stackunderflow(0);
769             rpnstack->s[stptr] = sqrt(rpnstack->s[stptr]);
770             break;
771         case OP_SORT:
772             stackunderflow(0);
773             {
774                 int       spn = (int) rpnstack->s[stptr--];
775
776                 stackunderflow(spn - 1);
777                 qsort(rpnstack->s + stptr - spn + 1, spn, sizeof(double),
778                       rpn_compare_double);
779             }
780             break;
781         case OP_REV:
782             stackunderflow(0);
783             {
784                 int       spn = (int) rpnstack->s[stptr--];
785                 double   *p, *q;
786
787                 stackunderflow(spn - 1);
788
789                 p = rpnstack->s + stptr - spn + 1;
790                 q = rpnstack->s + stptr;
791                 while (p < q) {
792                     double    x = *q;
793
794                     *q-- = *p;
795                     *p++ = x;
796                 }
797             }
798             break;
799         case OP_PREDICT:
800         case OP_PREDICTSIGMA:
801             stackunderflow(2);
802             {
803                 /* the local averaging window (similar to trend, but better here, as we get better statistics thru numbers)*/
804                 int   locstepsize = rpnstack->s[--stptr];
805                 /* the number of shifts and range-checking*/
806                 int     shifts = rpnstack->s[--stptr];
807                 stackunderflow(shifts);
808                 // handle negative shifts special
809                 if (shifts<0) {
810                     stptr--;
811                 } else {
812                     stptr-=shifts;
813                 }
814                 /* the real calculation */
815                 double val=DNAN;
816                 /* the info on the datasource */
817                 time_t  dsstep = (time_t) rpnp[rpi - 1].step;
818                 int    dscount = rpnp[rpi - 1].ds_cnt;
819                 int   locstep = (int)ceil((float)locstepsize/(float)dsstep);
820
821                 /* the sums */
822                 double    sum = 0;
823                 double    sum2 = 0;
824                 int       count = 0;
825                 /* now loop for each position */
826                 int doshifts=shifts;
827                 if (shifts<0) { doshifts=-shifts; }
828                 for(int loop=0;loop<doshifts;loop++) {
829                     /* calculate shift step */
830                     int shiftstep=1;
831                     if (shifts<0) {
832                         shiftstep = loop*rpnstack->s[stptr];
833                     } else { 
834                         shiftstep = rpnstack->s[stptr+loop]; 
835                     }
836                     if(shiftstep <0) {
837                         rrd_set_error("negative shift step not allowed: %i",shiftstep);
838                         return -1;
839                     }
840                     shiftstep=(int)ceil((float)shiftstep/(float)dsstep);
841                     /* loop all local shifts */
842                     for(int i=0;i<=locstep;i++) {
843                         /* now calculate offset into data-array - relative to output_idx*/
844                         int offset=shiftstep+i;
845                         /* and process if we have index 0 of above */
846                         if ((offset>=0)&&(offset<output_idx)) {
847                             /* get the value */
848                             val =rpnp[rpi - 1].data[-dscount * offset];
849                             /* and handle the non NAN case only*/
850                             if (! isnan(val)) {
851                                 sum+=val;
852                                 sum2+=val*val;
853                                 count++;
854                             }
855                         }
856                     }
857                 }
858                 /* do the final calculations */
859                 val=DNAN;
860                 if (rpnp[rpi].op == OP_PREDICT) {  /* the average */
861                     if (count>0) {
862                         val = sum/(double)count;
863                     } 
864                 } else {
865                     if (count>1) { /* the sigma case */
866                         val=count*sum2-sum*sum;
867                         if (val<0) {
868                             val=DNAN;
869                         } else {
870                             val=sqrt(val/((float)count*((float)count-1.0)));
871                         }
872                     }
873                 }
874                 rpnstack->s[stptr] = val;
875             }
876             break;
877         case OP_TREND:
878         case OP_TRENDNAN:
879             stackunderflow(1);
880             if ((rpi < 2) || (rpnp[rpi - 2].op != OP_VARIABLE)) {
881                 rrd_set_error("malformed trend arguments");
882                 return -1;
883             } else {
884                 time_t    dur = (time_t) rpnstack->s[stptr];
885                 time_t    step = (time_t) rpnp[rpi - 2].step;
886
887                 if (output_idx + 1 >= (int) ceil((float) dur / (float) step)) {
888                     int       ignorenan = (rpnp[rpi].op == OP_TREND);
889                     double    accum = 0.0;
890                     int       i = 0;
891                     int       count = 0;
892
893                     do {
894                         double    val =
895                             rpnp[rpi - 2].data[rpnp[rpi - 2].ds_cnt * i--];
896                         if (ignorenan || !isnan(val)) {
897                             accum += val;
898                             ++count;
899                         }
900
901                         dur -= step;
902                     } while (dur > 0);
903
904                     rpnstack->s[--stptr] =
905                         (count == 0) ? DNAN : (accum / count);
906                 } else
907                     rpnstack->s[--stptr] = DNAN;
908             }
909             break;
910         case OP_AVG:
911             stackunderflow(0);
912             {
913                 int       i = (int) rpnstack->s[stptr--];
914                 double    sum = 0;
915                 int       count = 0;
916
917                 stackunderflow(i - 1);
918                 while (i > 0) {
919                     double    val = rpnstack->s[stptr--];
920
921                     i--;
922                     if (isnan(val)) {
923                         continue;
924                     }
925                     count++;
926                     sum += val;
927                 }
928                 /* now push the result back on stack */
929                 if (count > 0) {
930                     rpnstack->s[++stptr] = sum / count;
931                 } else {
932                     rpnstack->s[++stptr] = DNAN;
933                 }
934             }
935             break;
936         case OP_ABS:
937             stackunderflow(0);
938             rpnstack->s[stptr] = fabs(rpnstack->s[stptr]);
939             break;
940         case OP_MEDIAN:
941             stackunderflow(0);
942             {
943                 int elements = (int) rpnstack->s[stptr--];
944                 int final_elements = elements;
945                 double *element_ptr = rpnstack->s + stptr - elements + 1;
946                 double *goodvals = element_ptr;
947                 double *badvals = element_ptr + elements - 1;
948
949                 stackunderflow(elements - 1);
950
951                 /* move values to consolidate the non-NANs for sorting, keeping
952                  * track of how many NANs we encounter. */
953                 while (goodvals < badvals) {
954                     if (isnan(*goodvals)) {
955                         *goodvals = *badvals--;
956                         --final_elements;
957                     } else {
958                         ++goodvals;
959                     }
960                 }
961
962                 stptr -= elements;
963                 if (!final_elements) {
964                     /* no non-NAN elements; push NAN */
965                     rpnstack->s[++stptr] = DNAN;
966                 } else {                   
967                     /* when goodvals and badvals meet, they might have met on a
968                      * NAN, which wouldn't decrease final_elements. so, check
969                      * that now. */
970                     if (isnan(*goodvals)) --final_elements;
971                     /* and finally, take the median of the remaining non-NAN
972                      * elements. */
973                     qsort(element_ptr, final_elements, sizeof(double),
974                           rpn_compare_double);
975                     if (final_elements % 2 == 1){
976                        rpnstack->s[++stptr] = element_ptr[ final_elements / 2 ];
977                     }
978                     else {
979                        rpnstack->s[++stptr] = 0.5 * ( element_ptr[ final_elements / 2 ] + element_ptr[ final_elements / 2 - 1 ] );
980                     }
981                 }
982             }
983             break;
984        
985         case OP_END:
986             break;
987         }
988 #undef stackunderflow
989     }
990     if (stptr != 0) {
991         rrd_set_error("RPN final stack size != 1");
992         return -1;
993     }
994
995     output[output_idx] = rpnstack->s[0];
996     return 0;
997 }
998
999 /* figure out what the local timezone offset for any point in
1000    time was. Return it in seconds */
1001 int tzoffset(
1002     time_t now)
1003 {
1004     int       gm_sec, gm_min, gm_hour, gm_yday, gm_year,
1005         l_sec, l_min, l_hour, l_yday, l_year;
1006     struct tm t;
1007     int       off;
1008
1009     gmtime_r(&now, &t);
1010     gm_sec = t.tm_sec;
1011     gm_min = t.tm_min;
1012     gm_hour = t.tm_hour;
1013     gm_yday = t.tm_yday;
1014     gm_year = t.tm_year;
1015     localtime_r(&now, &t);
1016     l_sec = t.tm_sec;
1017     l_min = t.tm_min;
1018     l_hour = t.tm_hour;
1019     l_yday = t.tm_yday;
1020     l_year = t.tm_year;
1021     off =
1022         (l_sec - gm_sec) + (l_min - gm_min) * 60 + (l_hour - gm_hour) * 3600;
1023     if (l_yday > gm_yday || l_year > gm_year) {
1024         off += 24 * 3600;
1025     } else if (l_yday < gm_yday || l_year < gm_year) {
1026         off -= 24 * 3600;
1027     }
1028     return off;
1029 }