Correct endpoints when index is before the epoch

The endpoints C code casts the double index to long, which truncates
it toward zero. This behavior is desired when the index is positive,
because it moves the endpoint *back* in time. But when the index is
negative, truncating toward zero moves the endpoint *forward* in time.

This is also an issue if the index value is stored as integer, since the
C99 specification states that integer division truncates toward zero.

If the first index value is less than zero, branch into a special case
to handle pre-epoch index values. This avoids performance degradation
if all index values are after the epoch.

If the index value is less than zero, simply add 1 to offset the
truncation toward zero. We also need to furthre adjust the potential
endpoint value if the index is exactly equal to zero.

Fixes #144.

Author: joshuaulrich

src/endpoints.c

@@ -36,26 +36,70 @@ SEXP endpoints (SEXP _x, SEXP _on, SEXP _k, SEXP _addlast /* TRUE */)
       /*int_index = INTEGER(getAttrib(_x, install("index")));*/
       int_index = INTEGER(_x);
       ep[0] = 0;
-      for(i=1,j=1; i<nr; i++) {
-        int_tmp[0] = int_index[i] / on / k +1;
-        int_tmp[1] = int_index[i-1] / on / k +1;
-        if( (int_tmp[0] - int_tmp[1]) != 0 ) {
-          ep[j] = i;
-          j++;
-        }
+      /* special handling if index values < 1970-01-01 00:00:00 UTC */
+      if(int_index[0] < 0) {
+          for(i=1,j=1; i<nr; i++) {
+            int idx_i0 = int_index[i];
+            int idx_i1 = int_index[i-1];
+
+            // if index equals epoch
+            int epoch_adj = (idx_i0 == 0);
+            // cast truncates toward zero
+            // (i.e. *forward* in time if index < 0)
+            if(idx_i0 < 0) idx_i0++;
+            if(idx_i1 < 0) idx_i1++;
+
+            int_tmp[0] = idx_i0 / on / k + epoch_adj;
+            int_tmp[1] = idx_i1 / on / k;
+            if( (int_tmp[0] - int_tmp[1]) != 0 ) {
+              ep[j] = i;
+              j++;
+            }
+          }
+      } else {
+          for(i=1,j=1; i<nr; i++) {
+            int_tmp[0] = int_index[i] / on / k +1;
+            int_tmp[1] = int_index[i-1] / on / k +1;
+            if( (int_tmp[0] - int_tmp[1]) != 0 ) {
+              ep[j] = i;
+              j++;
+            }
+          }
       }
       break;
     case REALSXP:
       /*real_index = REAL(getAttrib(_x, install("index")));*/
       real_index = REAL(_x);
       ep[0] = 0;
-      for(i=1,j=1; i<nr; i++) {
-        real_tmp[0] = (long)real_index[i] / on / k +1;
-        real_tmp[1] = (long)real_index[i-1] / on / k +1;
-        if( (real_tmp[0] - real_tmp[1]) != 0 ) {
-          ep[j] = i;
-          j++;
-        }
+      /* special handling if index values < 1970-01-01 00:00:00 UTC */
+      if(real_index[0] < 0) {
+          for(i=1,j=1; i<nr; i++) {
+            double idx_i0 = real_index[i];
+            double idx_i1 = real_index[i-1];
+
+            // if index equals epoch
+            int epoch_adj = (idx_i0 == 0);
+            // cast truncates toward zero
+            // (i.e. *forward* in time if index < 0)
+            if(idx_i0 < 0) idx_i0 += 1.0;
+            if(idx_i1 < 0) idx_i1 += 1.0;
+
+            real_tmp[0] = (long)(idx_i0 / on / k + epoch_adj);
+            real_tmp[1] = (long)(idx_i1 / on / k);
+            if( (real_tmp[0] - real_tmp[1]) != 0) {
+              ep[j] = i;
+              j++;
+            }
+          }
+      } else {
+          for(i=1,j=1; i<nr; i++) {
+            real_tmp[0] = (long)real_index[i] / on / k +1;
+            real_tmp[1] = (long)real_index[i-1] / on / k +1;
+            if( (real_tmp[0] - real_tmp[1]) != 0 ) {
+              ep[j] = i;
+              j++;
+            }
+          }
       }
       break;
     default: