Skip to content

Commit 0c443e2

Browse files
mns-nordicalsmns-nordicals
authored
[2/15] Fix wrong-kind literals and kind-less casts in logical division, integer %/%, and Re() (#138)
* Emit real literals in the logical-to-double cast maybe_cast_double() wrote merge(1_c_double, 0_c_double, x), but a kind-suffixed 1_c_double is an integer(8) literal, not a double, so the "cast" produced integers and TRUE / FALSE was an integer division by zero (undefined behavior; gfortran -O2 happens to return 1) instead of R's Inf. Spell the literals 1.0_c_double/0.0_c_double, matching the correct spelling already used by floor()/ceiling(). Every consumer of the helper (/, as.double(), trunc(), BLAS operand casts) inherits the fix. * Compute integer %/% in double precision The integer branch of %/% emitted int(floor(real(a) / real(b))); kind-less real() of an integer is default (single-precision) real, which cannot represent odd integers above 2^24, so e.g. 16777219L %/% 2L returned 8388610 instead of R's 8388609 — a silent wrong answer. Cast both operands with kind=c_double and the result with kind=c_int. * Emit Re() with an explicit c_double kind Re() lowered to kind-less real(arg). For a complex argument that is correct (F2018: the result takes the argument's kind), but the handler accepts any mode, and for a double or integer argument kind-less real() is single precision — so Re(0.1) returned 0.100000001490116, silently rounded through a float while declared real(c_double). Emit real(arg, kind=c_double), which is exact for double/integer arguments and a no-op change of meaning for complex ones.
1 parent 7a36ea3 commit 0c443e2

9 files changed

Lines changed: 277 additions & 5 deletions

File tree

‎R/r2f-arithmetic.R‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -83,7 +83,9 @@ r2f_handlers[["%/%"]] <- function(args, scope, ...) {
8383

8484
expr <- switch(
8585
out_val@mode,
86-
integer = glue("int(floor(real({left}) / real({right})))"),
86+
integer = glue(
87+
"int(floor(real({left}, kind=c_double) / real({right}, kind=c_double)), kind=c_int)"
88+
),
8789
double = glue("floor({left} / {right})"),
8890
stop("%/% only implemented for numeric types")
8991
)

‎R/r2f-math.R‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -131,7 +131,9 @@ r2f_handlers[["abs"]] <- function(args, scope, ...) {
131131
register_unary_intrinsic(
132132
"Re",
133133
mode_fun = function(arg) "double",
134-
expr_fun = function(arg, intrinsic) glue("real({arg})")
134+
# kind-less real() of an integer or double argument is single precision;
135+
# kind=c_double is a no-op change of meaning for complex arguments (F2018).
136+
expr_fun = function(arg, intrinsic) glue("real({arg}, kind=c_double)")
135137
)
136138

137139
register_unary_intrinsic(

‎R/r2f-operators-helpers.R‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -30,7 +30,7 @@ booleanize_logical_as_int <- function(x) {
3030
maybe_cast_double <- function(x) {
3131
if (x@value@mode == "logical") {
3232
Fortran(
33-
glue("merge(1_c_double, 0_c_double, {x})"),
33+
glue("merge(1.0_c_double, 0.0_c_double, {x})"),
3434
Variable("double", x@value@dims)
3535
)
3636
} else if (x@value@mode == "integer") {

‎tests/testthat/_snaps/div-cast.md‎

Lines changed: 76 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -188,7 +188,7 @@
188188
! manifest end
189189
190190
191-
out_ = (a / merge(1_c_double, 0_c_double, (b/=0)))
191+
out_ = (a / merge(1.0_c_double, 0.0_c_double, (b/=0)))
192192
end subroutine
193193
Code
194194
cat(cwrapper)
@@ -241,6 +241,81 @@
241241
return out_;
242242
}
243243

244+
# division casts logical by logical
245+
246+
Code
247+
fn
248+
Output
249+
function(a, b) {
250+
declare(type(a = logical(1)), type(b = logical(1)))
251+
a / b
252+
}
253+
<environment: 0x0>
254+
Code
255+
cat(fsub)
256+
Output
257+
subroutine fn(a, b, out_) bind(c)
258+
use iso_c_binding, only: c_double, c_int
259+
implicit none
260+
261+
! manifest start
262+
! args
263+
integer(c_int), intent(in) :: a ! logical
264+
integer(c_int), intent(in) :: b ! logical
265+
real(c_double), intent(out) :: out_
266+
! manifest end
267+
268+
269+
out_ = (merge(1.0_c_double, 0.0_c_double, (a/=0)) / merge(1.0_c_double, 0.0_c_double, (b/=0)))
270+
end subroutine
271+
Code
272+
cat(cwrapper)
273+
Output
274+
#define R_NO_REMAP
275+
#include <R.h>
276+
#include <Rinternals.h>
277+
278+
279+
extern void fn(
280+
const int* const a__,
281+
const int* const b__,
282+
double* const out___);
283+
284+
SEXP fn_(SEXP _args) {
285+
// a
286+
_args = CDR(_args);
287+
SEXP a = CAR(_args);
288+
if (TYPEOF(a) != LGLSXP) {
289+
Rf_error("typeof(a) must be 'logical', not '%s'", Rf_type2char(TYPEOF(a)));
290+
}
291+
const int* const a__ = LOGICAL(a);
292+
const R_xlen_t a__len_ = Rf_xlength(a);
293+
294+
// b
295+
_args = CDR(_args);
296+
SEXP b = CAR(_args);
297+
if (TYPEOF(b) != LGLSXP) {
298+
Rf_error("typeof(b) must be 'logical', not '%s'", Rf_type2char(TYPEOF(b)));
299+
}
300+
const int* const b__ = LOGICAL(b);
301+
const R_xlen_t b__len_ = Rf_xlength(b);
302+
303+
if (a__len_ != 1)
304+
Rf_error("length(a) must be 1, not %0.f",
305+
(double)a__len_);
306+
if (b__len_ != 1)
307+
Rf_error("length(b) must be 1, not %0.f",
308+
(double)b__len_);
309+
const R_xlen_t out___len_ = (1);
310+
SEXP out_ = PROTECT(Rf_allocVector(REALSXP, out___len_));
311+
double* out___ = REAL(out_);
312+
313+
fn(a__, b__, out___);
314+
315+
UNPROTECT(1);
316+
return out_;
317+
}
318+
244319
# division casts complex
245320

246321
Code

‎tests/testthat/_snaps/div-mod.md‎

Lines changed: 75 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -160,3 +160,78 @@
160160
return out_;
161161
}
162162

163+
# integer %/% is exact beyond 2^24
164+
165+
Code
166+
fn
167+
Output
168+
function(a, b) {
169+
declare(type(a = integer(1)), type(b = integer(1)))
170+
a %/% b
171+
}
172+
<environment: 0x0>
173+
Code
174+
cat(fsub)
175+
Output
176+
subroutine fn(a, b, out_) bind(c)
177+
use iso_c_binding, only: c_double, c_int
178+
implicit none
179+
180+
! manifest start
181+
! args
182+
integer(c_int), intent(in) :: a
183+
integer(c_int), intent(in) :: b
184+
integer(c_int), intent(out) :: out_
185+
! manifest end
186+
187+
188+
out_ = int(floor(real(a, kind=c_double) / real(b, kind=c_double)), kind=c_int)
189+
end subroutine
190+
Code
191+
cat(cwrapper)
192+
Output
193+
#define R_NO_REMAP
194+
#include <R.h>
195+
#include <Rinternals.h>
196+
197+
198+
extern void fn(
199+
const int* const a__,
200+
const int* const b__,
201+
int* const out___);
202+
203+
SEXP fn_(SEXP _args) {
204+
// a
205+
_args = CDR(_args);
206+
SEXP a = CAR(_args);
207+
if (TYPEOF(a) != INTSXP) {
208+
Rf_error("typeof(a) must be 'integer', not '%s'", Rf_type2char(TYPEOF(a)));
209+
}
210+
const int* const a__ = INTEGER(a);
211+
const R_xlen_t a__len_ = Rf_xlength(a);
212+
213+
// b
214+
_args = CDR(_args);
215+
SEXP b = CAR(_args);
216+
if (TYPEOF(b) != INTSXP) {
217+
Rf_error("typeof(b) must be 'integer', not '%s'", Rf_type2char(TYPEOF(b)));
218+
}
219+
const int* const b__ = INTEGER(b);
220+
const R_xlen_t b__len_ = Rf_xlength(b);
221+
222+
if (a__len_ != 1)
223+
Rf_error("length(a) must be 1, not %0.f",
224+
(double)a__len_);
225+
if (b__len_ != 1)
226+
Rf_error("length(b) must be 1, not %0.f",
227+
(double)b__len_);
228+
const R_xlen_t out___len_ = (1);
229+
SEXP out_ = PROTECT(Rf_allocVector(INTSXP, out___len_));
230+
int* out___ = INTEGER(out_);
231+
232+
fn(a__, b__, out___);
233+
234+
UNPROTECT(1);
235+
return out_;
236+
}
237+

‎tests/testthat/_snaps/unary-intrinsics.md‎

Lines changed: 64 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1691,7 +1691,7 @@
16911691
! manifest end
16921692
16931693
1694-
out = real(z)
1694+
out = real(z, kind=c_double)
16951695
end subroutine
16961696
Code
16971697
cat(cwrapper)
@@ -1982,6 +1982,69 @@
19821982
return out;
19831983
}
19841984

1985+
# Re() on a double argument keeps double precision
1986+
1987+
Code
1988+
fn
1989+
Output
1990+
function(x) {
1991+
declare(type(x = double(n)))
1992+
out <- Re(x)
1993+
out
1994+
}
1995+
<environment: 0x0>
1996+
Code
1997+
cat(fsub)
1998+
Output
1999+
subroutine fn(x, out, x__len_) bind(c)
2000+
use iso_c_binding, only: c_double, c_ptrdiff_t
2001+
implicit none
2002+
2003+
! manifest start
2004+
! sizes
2005+
integer(c_ptrdiff_t), intent(in), value :: x__len_
2006+
2007+
! args
2008+
real(c_double), intent(in) :: x(x__len_)
2009+
real(c_double), intent(out) :: out(x__len_)
2010+
! manifest end
2011+
2012+
2013+
out = real(x, kind=c_double)
2014+
end subroutine
2015+
Code
2016+
cat(cwrapper)
2017+
Output
2018+
#define R_NO_REMAP
2019+
#include <R.h>
2020+
#include <Rinternals.h>
2021+
2022+
2023+
extern void fn(
2024+
const double* const x__,
2025+
double* const out__,
2026+
const R_xlen_t x__len_);
2027+
2028+
SEXP fn_(SEXP _args) {
2029+
// x
2030+
_args = CDR(_args);
2031+
SEXP x = CAR(_args);
2032+
if (TYPEOF(x) != REALSXP) {
2033+
Rf_error("typeof(x) must be 'double', not '%s'", Rf_type2char(TYPEOF(x)));
2034+
}
2035+
const double* const x__ = REAL(x);
2036+
const R_xlen_t x__len_ = Rf_xlength(x);
2037+
2038+
const R_xlen_t out__len_ = x__len_;
2039+
SEXP out = PROTECT(Rf_allocVector(REALSXP, out__len_));
2040+
double* out__ = REAL(out);
2041+
2042+
fn(x__, out__, x__len_);
2043+
2044+
UNPROTECT(1);
2045+
return out;
2046+
}
2047+
19852048
# logical not local used as ifelse mask compiles and runs
19862049

19872050
Code

‎tests/testthat/test-div-cast.R‎

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -42,6 +42,27 @@ test_that("division casts logical", {
4242
})
4343

4444

45+
test_that("division casts logical by logical", {
46+
div_lgl_lgl <- function(a, b) {
47+
declare(type(a = logical(1)), type(b = logical(1)))
48+
a / b
49+
}
50+
51+
# locks the 1.0_c_double/0.0_c_double literals: with integer-kind
52+
# 1_c_double literals this was an integer division, and TRUE/FALSE
53+
# divided by zero instead of yielding Inf
54+
expect_translation_snapshots(div_lgl_lgl)
55+
56+
expect_quick_equal(
57+
div_lgl_lgl,
58+
list(TRUE, TRUE),
59+
list(TRUE, FALSE), # Inf
60+
list(FALSE, TRUE),
61+
list(FALSE, FALSE) # NaN
62+
)
63+
})
64+
65+
4566
test_that("division casts complex", {
4667
div_cplx <- function(a, b) {
4768
declare(type(a = complex(n)), type(b = complex(n)))

‎tests/testthat/test-div-mod.R‎

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -36,3 +36,22 @@ test_that("%% and %/%", {
3636
}
3737
expect_quick_identical(div_int, x)
3838
})
39+
40+
test_that("integer %/% is exact beyond 2^24", {
41+
div_int_big <- function(a, b) {
42+
declare(type(a = integer(1)), type(b = integer(1)))
43+
a %/% b
44+
}
45+
46+
# locks the kind=c_double casts: kind-less real() is single precision,
47+
# which cannot represent odd integers above 2^24, so e.g.
48+
# 16777219L %/% 2L came back as 8388610 instead of 8388609
49+
expect_translation_snapshots(div_int_big)
50+
51+
expect_quick_identical(
52+
div_int_big,
53+
list(16777219L, 2L),
54+
list(-16777219L, 2L),
55+
list(.Machine$integer.max, 7L)
56+
)
57+
})

‎tests/testthat/test-unary-intrinsics.R‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -106,6 +106,21 @@ test_that("complex unary intrinsics", {
106106
})
107107

108108

109+
test_that("Re() on a double argument keeps double precision", {
110+
fn <- function(x) {
111+
declare(type(x = double(n)))
112+
out <- Re(x)
113+
out
114+
}
115+
116+
# locks real(x, kind=c_double): kind-less real() of a double argument
117+
# is single precision, so Re() silently rounded its input
118+
expect_translation_snapshots(fn)
119+
120+
expect_quick_identical(fn, list(c(0.1, 1 / 3, -2.5)))
121+
})
122+
123+
109124
test_that("unary logical 'not'-operator on vector", {
110125
fn <- function(x) {
111126
declare(type(x = integer(n)))

0 commit comments

Comments
 (0)