Skip to content

Commit 4074e12

Browse files
jangoreckijangorecki
committed
zeros and num stability
1 parent cf898ab commit 4074e12

File tree

1 file changed

+99
-16
lines changed

1 file changed

+99
-16
lines changed

‎src/froll.c‎

Lines changed: 99 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -888,7 +888,13 @@ void frollminExact(const double *x, uint64_t nx, ans_t *ans, int k, double fill,
888888
#undef PROD_WINDOW_STEP_FRONT
889889
#define PROD_WINDOW_STEP_FRONT \
890890
if (R_FINITE(x[i])) { \
891-
w *= x[i]; \
891+
if (x[i] == 0.0) { \
892+
zerc++; \
893+
} else { \
894+
w += log(fabsl(x[i])); \
895+
if (x[i] < 0.0) \
896+
negc++; \
897+
} \
892898
} else if (ISNAN(x[i])) { \
893899
nc++; \
894900
} else if (x[i]==R_PosInf) { \
@@ -899,7 +905,13 @@ void frollminExact(const double *x, uint64_t nx, ans_t *ans, int k, double fill,
899905
#undef PROD_WINDOW_STEP_BACK
900906
#define PROD_WINDOW_STEP_BACK \
901907
if (R_FINITE(x[i-k])) { \
902-
w /= x[i-k]; \
908+
if (x[i-k] == 0.0) { \
909+
zerc--; \
910+
} else { \
911+
w -= log(fabsl(x[i-k])); \
912+
if (x[i-k] < 0.0) \
913+
negc--; \
914+
} \
903915
} else if (ISNAN(x[i-k])) { \
904916
nc--; \
905917
} else if (x[i-k]==R_PosInf) { \
@@ -911,26 +923,57 @@ void frollminExact(const double *x, uint64_t nx, ans_t *ans, int k, double fill,
911923
#define PROD_WINDOW_STEP_VALUE \
912924
if (nc == 0) { \
913925
if (pinf == 0 && ninf == 0) { \
914-
ans->dbl_v[i] = (double) w; \
926+
if (zerc) { \
927+
ans->dbl_v[i] = 0.0; \
928+
} else { \
929+
long double res = expl(w); \
930+
if (negc % 2) \
931+
res = -res; \
932+
ans->dbl_v[i] = (double) res; \
933+
} \
915934
} else { \
916-
ans->dbl_v[i] = (ninf+(w<0))%2 ? R_NegInf : R_PosInf; \
935+
if (zerc) { \
936+
ans->dbl_v[i] = R_NaN; \
937+
} else { \
938+
if ((ninf + (negc%2)) % 2) { \
939+
ans->dbl_v[i] = R_NegInf; \
940+
} else { \
941+
ans->dbl_v[i] = R_PosInf; \
942+
} \
943+
} \
917944
} \
918945
} else if (nc == k) { \
919946
ans->dbl_v[i] = narm ? 1.0 : NA_REAL; \
920947
} else { \
921948
if (narm) { \
922949
if (pinf == 0 && ninf == 0) { \
923-
ans->dbl_v[i] = (double) w; \
950+
if (zerc) { \
951+
ans->dbl_v[i] = 0.0; \
952+
} else { \
953+
long double res = expl(w); \
954+
if (negc % 2) \
955+
res = -res; \
956+
ans->dbl_v[i] = (double) res; \
957+
} \
924958
} else { \
925-
ans->dbl_v[i] = (ninf+(w<0))%2 ? R_NegInf : R_PosInf; \
959+
if (zerc) { \
960+
ans->dbl_v[i] = R_NaN; \
961+
} else { \
962+
if ((ninf + (negc%2)) % 2) { \
963+
ans->dbl_v[i] = R_NegInf; \
964+
} else { \
965+
ans->dbl_v[i] = R_PosInf; \
966+
} \
967+
} \
926968
} \
927969
} else { \
928970
ans->dbl_v[i] = NA_REAL; \
929971
} \
930972
}
931973

932974
/* fast rolling prod - fast
933-
* same as mean fast
975+
* work in log space: sum rather than prod
976+
* track zero and negative
934977
*/
935978
void frollprodFast(const double *x, uint64_t nx, ans_t *ans, int k, double fill, bool narm, int hasnf, bool verbose) {
936979
if (verbose)
@@ -943,35 +986,75 @@ void frollprodFast(const double *x, uint64_t nx, ans_t *ans, int k, double fill,
943986
}
944987
return;
945988
}
946-
long double w = 1.0;
989+
long double w = 0.0;
990+
int zerc = 0;
991+
int negc = 0;
947992
bool truehasnf = hasnf>0;
948993
if (!truehasnf) {
949994
int i;
950995
for (i=0; i<k-1; i++) { // #loop_counter_not_local_scope_ok
951-
w *= x[i];
996+
if (x[i] == 0.0) {
997+
zerc++;
998+
} else {
999+
w += log(fabsl(x[i]));
1000+
if (x[i] < 0.0)
1001+
negc++;
1002+
}
9521003
ans->dbl_v[i] = fill;
9531004
}
954-
w *= x[i];
955-
ans->dbl_v[i] = (double) w;
1005+
if (x[i] == 0.0) {
1006+
zerc++;
1007+
} else {
1008+
w += log(fabsl(x[i]));
1009+
if (x[i] < 0.0)
1010+
negc++;
1011+
}
1012+
if (zerc) {
1013+
ans->dbl_v[i] = 0.0;
1014+
} else {
1015+
long double res = expl(w);
1016+
if (negc % 2)
1017+
res = -res;
1018+
ans->dbl_v[i] = (double) res;
1019+
}
9561020
if (R_FINITE((double) w)) {
9571021
for (uint64_t i=k; i<nx; i++) {
958-
w /= x[i-k];
959-
w *= x[i];
960-
ans->dbl_v[i] = (double) w;
1022+
if (x[i-k] == 0.0) {
1023+
zerc--;
1024+
} else {
1025+
w -= log(fabsl(x[i-k]));
1026+
if (x[i-k] < 0.0)
1027+
negc--;
1028+
}
1029+
if (x[i] == 0.0) {
1030+
zerc++;
1031+
} else {
1032+
w += log(fabsl(x[i]));
1033+
if (x[i] < 0.0)
1034+
negc++;
1035+
}
1036+
if (zerc) {
1037+
ans->dbl_v[i] = 0.0;
1038+
} else {
1039+
long double res = expl(w);
1040+
if (negc % 2)
1041+
res = -res;
1042+
ans->dbl_v[i] = (double) res;
1043+
}
9611044
}
9621045
if (!R_FINITE((double) w)) {
9631046
if (hasnf==-1)
9641047
ansSetMsg(ans, 2, "%s: has.nf=FALSE used but non-finite values are present in input, use default has.nf=NA to avoid this warning", __func__);
9651048
if (verbose)
9661049
ansSetMsg(ans, 0, "%s: non-finite values are present in input, re-running with extra care for NFs\n", __func__);
967-
w = 1.0; truehasnf = true;
1050+
w = 0.0; zerc = 0; negc = 0; truehasnf = true;
9681051
}
9691052
} else {
9701053
if (hasnf==-1)
9711054
ansSetMsg(ans, 2, "%s: has.nf=FALSE used but non-finite values are present in input, use default has.nf=NA to avoid this warning", __func__);
9721055
if (verbose)
9731056
ansSetMsg(ans, 0, "%s: non-finite values are present in input, skip non-finite inaware attempt and run with extra care for NFs straighaway\n", __func__);
974-
w = 1.0; truehasnf = true;
1057+
w = 0.0; zerc = 0; negc = 0; truehasnf = true;
9751058
}
9761059
}
9771060
if (truehasnf) {

0 commit comments

Comments
 (0)