Skip to content

Commit c78488f

Browse files
SebKrantzSebKrantz
committed
Refactor prangeR function to enhance NA handling and ensure double precision for range calculations, improving type conversion and memory management.
1 parent db64cec commit c78488f

File tree

1 file changed

+85
-161
lines changed

1 file changed

+85
-161
lines changed

‎src/psum.c‎

Lines changed: 85 additions & 161 deletions
Original file line numberDiff line numberDiff line change
@@ -1359,39 +1359,42 @@ SEXP prangeR(SEXP na, SEXP args) {
13591359
if (hasFactor) {
13601360
error("Function 'prange' is not meaningful for factors.");
13611361
}
1362+
// For range, we need to compute max - min, so we need double precision to avoid overflow
1363+
if(anstype != REALSXP) anstype = REALSXP;
13621364
int nprotect=1;
1363-
const bool init_convert = (anstype != type0);
1364-
SEXP ans = init_convert ? PROTECT(allocVector(anstype, len0)) : PROTECT(duplicate(args0));
1365-
if (init_convert) {
1366-
copyMostAttrib(args0, ans);
1367-
}
1365+
SEXP ans = PROTECT(allocVector(REALSXP, len0));
1366+
copyMostAttrib(args0, ans);
13681367
const bool narm = asLogical(na);
1369-
switch(anstype) {
1370-
case LGLSXP:
1371-
case INTSXP: {
1372-
int *restrict pans = (anstype == LGLSXP) ? LOGICAL(ans) : INTEGER(ans);
1373-
if (init_convert) {
1374-
const int *restrict p0 = LOGICAL(args0);
1375-
for (ssize_t j = 0; j < len0; ++j) {
1376-
pans[j] = p0[j];
1377-
}
1368+
double *restrict pans = REAL(ans);
1369+
// Initialize with first argument
1370+
if (type0 != REALSXP) {
1371+
const int *restrict p0 = (type0 == LGLSXP) ? LOGICAL(args0) : INTEGER(args0);
1372+
for (ssize_t j = 0; j < len0; ++j) {
1373+
pans[j] = p0[j] == NA_INTEGER ? R_NaN : (double)p0[j];
13781374
}
1379-
int *restrict pmax = (int*)R_alloc(len0, sizeof(int));
1380-
if (narm) {
1381-
int *restrict found = (int*)R_alloc(len0, sizeof(int));
1382-
for (ssize_t j = 0; j < len0; ++j) {
1383-
if (pans[j] == NA_INTEGER) {
1384-
found[j] = 0;
1385-
} else {
1386-
found[j] = 1;
1387-
pmax[j] = pans[j];
1388-
}
1375+
} else {
1376+
const double *restrict p0 = REAL(args0);
1377+
for (ssize_t j = 0; j < len0; ++j) {
1378+
pans[j] = p0[j];
1379+
}
1380+
}
1381+
double *restrict pmax = (double*)R_alloc(len0, sizeof(double));
1382+
if (narm) {
1383+
int *restrict found = (int*)R_alloc(len0, sizeof(int));
1384+
for (ssize_t j = 0; j < len0; ++j) {
1385+
if (ISNAN(pans[j])) {
1386+
found[j] = 0;
1387+
} else {
1388+
found[j] = 1;
1389+
pmax[j] = pans[j];
13891390
}
1390-
for (int i = 1; i < n; ++i) {
1391-
const SEXP argi = PTR_ETL(args, i);
1392-
const int *restrict pa = (UTYPEOF(argi) == LGLSXP) ? LOGICAL(argi) : INTEGER(argi);
1391+
}
1392+
for (int i = 1; i < n; ++i) {
1393+
const SEXP argi = PTR_ETL(args, i);
1394+
if (UTYPEOF(argi) == REALSXP) {
1395+
const double *restrict pa = REAL(argi);
13931396
for (ssize_t j = 0; j < len0; ++j) {
1394-
if (pa[j] != NA_INTEGER) {
1397+
if (!ISNAN(pa[j])) {
13951398
if (!found[j]) {
13961399
pans[j] = pa[j];
13971400
pmax[j] = pa[j];
@@ -1406,159 +1409,80 @@ SEXP prangeR(SEXP na, SEXP args) {
14061409
}
14071410
}
14081411
}
1409-
}
1410-
for (ssize_t j = 0; j < len0; ++j) {
1411-
if (!found[j]) {
1412-
pans[j] = NA_INTEGER;
1413-
} else {
1414-
pans[j] = pmax[j] - pans[j];
1415-
}
1416-
}
1417-
} else {
1418-
for (ssize_t j = 0; j < len0; ++j) {
1419-
pmax[j] = pans[j];
1420-
}
1421-
for (int i = 1; i < n; ++i) {
1422-
const SEXP argi = PTR_ETL(args, i);
1412+
} else {
14231413
const int *restrict pa = (UTYPEOF(argi) == LGLSXP) ? LOGICAL(argi) : INTEGER(argi);
14241414
for (ssize_t j = 0; j < len0; ++j) {
1425-
if (pans[j] == NA_INTEGER || pa[j] == NA_INTEGER) {
1426-
pans[j] = NA_INTEGER;
1427-
pmax[j] = NA_INTEGER;
1428-
} else {
1429-
if (pa[j] < pans[j]) {
1430-
pans[j] = pa[j];
1431-
}
1432-
if (pa[j] > pmax[j]) {
1433-
pmax[j] = pa[j];
1415+
if (pa[j] != NA_INTEGER) {
1416+
const double v = (double)pa[j];
1417+
if (!found[j]) {
1418+
pans[j] = v;
1419+
pmax[j] = v;
1420+
found[j] = 1;
1421+
} else {
1422+
if (v < pans[j]) {
1423+
pans[j] = v;
1424+
}
1425+
if (v > pmax[j]) {
1426+
pmax[j] = v;
1427+
}
14341428
}
14351429
}
14361430
}
14371431
}
1438-
for (ssize_t j = 0; j < len0; ++j) {
1439-
if (pans[j] == NA_INTEGER || pmax[j] == NA_INTEGER) {
1440-
pans[j] = NA_INTEGER;
1441-
} else {
1442-
pans[j] = pmax[j] - pans[j];
1443-
}
1444-
}
14451432
}
1446-
} break;
1447-
case REALSXP: {
1448-
double *restrict pans = REAL(ans);
1449-
if (init_convert) {
1450-
const int *restrict p0 = (type0 == LGLSXP) ? LOGICAL(args0) : INTEGER(args0);
1451-
for (ssize_t j = 0; j < len0; ++j) {
1452-
pans[j] = p0[j] == NA_INTEGER ? R_NaN : (double)p0[j];
1433+
for (ssize_t j = 0; j < len0; ++j) {
1434+
if (!found[j]) {
1435+
pans[j] = R_NaN;
1436+
} else {
1437+
pans[j] = pmax[j] - pans[j];
14531438
}
14541439
}
1455-
double *restrict pmax = (double*)R_alloc(len0, sizeof(double));
1456-
if (narm) {
1457-
int *restrict found = (int*)R_alloc(len0, sizeof(int));
1458-
for (ssize_t j = 0; j < len0; ++j) {
1459-
if (ISNAN(pans[j])) {
1460-
found[j] = 0;
1461-
} else {
1462-
found[j] = 1;
1463-
pmax[j] = pans[j];
1464-
}
1465-
}
1466-
for (int i = 1; i < n; ++i) {
1467-
const SEXP argi = PTR_ETL(args, i);
1468-
if (UTYPEOF(argi) == REALSXP) {
1469-
const double *restrict pa = REAL(argi);
1470-
for (ssize_t j = 0; j < len0; ++j) {
1471-
if (!ISNAN(pa[j])) {
1472-
if (!found[j]) {
1473-
pans[j] = pa[j];
1474-
pmax[j] = pa[j];
1475-
found[j] = 1;
1476-
} else {
1477-
if (pa[j] < pans[j]) {
1478-
pans[j] = pa[j];
1479-
}
1480-
if (pa[j] > pmax[j]) {
1481-
pmax[j] = pa[j];
1482-
}
1483-
}
1440+
} else {
1441+
for (ssize_t j = 0; j < len0; ++j) {
1442+
pmax[j] = pans[j];
1443+
}
1444+
for (int i = 1; i < n; ++i) {
1445+
const SEXP argi = PTR_ETL(args, i);
1446+
if (UTYPEOF(argi) == REALSXP) {
1447+
const double *restrict pa = REAL(argi);
1448+
for (ssize_t j = 0; j < len0; ++j) {
1449+
if (ISNAN(pans[j]) || ISNAN(pa[j])) {
1450+
pans[j] = R_NaN;
1451+
pmax[j] = R_NaN;
1452+
} else {
1453+
if (pa[j] < pans[j]) {
1454+
pans[j] = pa[j];
14841455
}
1485-
}
1486-
} else {
1487-
const int *restrict pa = (UTYPEOF(argi) == LGLSXP) ? LOGICAL(argi) : INTEGER(argi);
1488-
for (ssize_t j = 0; j < len0; ++j) {
1489-
if (pa[j] != NA_INTEGER) {
1490-
const double v = (double)pa[j];
1491-
if (!found[j]) {
1492-
pans[j] = v;
1493-
pmax[j] = v;
1494-
found[j] = 1;
1495-
} else {
1496-
if (v < pans[j]) {
1497-
pans[j] = v;
1498-
}
1499-
if (v > pmax[j]) {
1500-
pmax[j] = v;
1501-
}
1502-
}
1456+
if (pa[j] > pmax[j]) {
1457+
pmax[j] = pa[j];
15031458
}
15041459
}
15051460
}
1506-
}
1507-
for (ssize_t j = 0; j < len0; ++j) {
1508-
if (!found[j]) {
1509-
pans[j] = R_NaN;
1510-
} else {
1511-
pans[j] = pmax[j] - pans[j];
1512-
}
1513-
}
1514-
} else {
1515-
for (ssize_t j = 0; j < len0; ++j) {
1516-
pmax[j] = pans[j];
1517-
}
1518-
for (int i = 1; i < n; ++i) {
1519-
const SEXP argi = PTR_ETL(args, i);
1520-
if (UTYPEOF(argi) == REALSXP) {
1521-
const double *restrict pa = REAL(argi);
1522-
for (ssize_t j = 0; j < len0; ++j) {
1523-
if (ISNAN(pans[j]) || ISNAN(pa[j])) {
1524-
pans[j] = R_NaN;
1525-
pmax[j] = R_NaN;
1526-
} else {
1527-
if (pa[j] < pans[j]) {
1528-
pans[j] = pa[j];
1529-
}
1530-
if (pa[j] > pmax[j]) {
1531-
pmax[j] = pa[j];
1532-
}
1461+
} else {
1462+
const int *restrict pa = (UTYPEOF(argi) == LGLSXP) ? LOGICAL(argi) : INTEGER(argi);
1463+
for (ssize_t j = 0; j < len0; ++j) {
1464+
if (ISNAN(pans[j]) || pa[j] == NA_INTEGER) {
1465+
pans[j] = R_NaN;
1466+
pmax[j] = R_NaN;
1467+
} else {
1468+
const double v = (double)pa[j];
1469+
if (v < pans[j]) {
1470+
pans[j] = v;
15331471
}
1534-
}
1535-
} else {
1536-
const int *restrict pa = (UTYPEOF(argi) == LGLSXP) ? LOGICAL(argi) : INTEGER(argi);
1537-
for (ssize_t j = 0; j < len0; ++j) {
1538-
if (ISNAN(pans[j]) || pa[j] == NA_INTEGER) {
1539-
pans[j] = R_NaN;
1540-
pmax[j] = R_NaN;
1541-
} else {
1542-
const double v = (double)pa[j];
1543-
if (v < pans[j]) {
1544-
pans[j] = v;
1545-
}
1546-
if (v > pmax[j]) {
1547-
pmax[j] = v;
1548-
}
1472+
if (v > pmax[j]) {
1473+
pmax[j] = v;
15491474
}
15501475
}
15511476
}
15521477
}
1553-
for (ssize_t j = 0; j < len0; ++j) {
1554-
if (ISNAN(pans[j]) || ISNAN(pmax[j])) {
1555-
pans[j] = R_NaN;
1556-
} else {
1557-
pans[j] = pmax[j] - pans[j];
1558-
}
1478+
}
1479+
for (ssize_t j = 0; j < len0; ++j) {
1480+
if (ISNAN(pans[j]) || ISNAN(pmax[j])) {
1481+
pans[j] = R_NaN;
1482+
} else {
1483+
pans[j] = pmax[j] - pans[j];
15591484
}
15601485
}
1561-
} break;
15621486
}
15631487
UNPROTECT(nprotect);
15641488
return ans;

0 commit comments

Comments
 (0)