Author Archives: Jonathan Carroll

Counting Digits Quickly

By: Jonathan Carroll

Re-posted from: https://jcarroll.com.au/2025/06/29/counting-digits-quickly/

When things run slower than we’d like in R we tend to reach for another, usually
compiled, language, and move our code there. What if it “just happened”? What
started out as a silly exploration of how to count digits ended up with a race
to see which language does it fastest. Maybe some surprises here for some, maybe
some bad implementations on my part – let’s find out.

I saw some recent activity on {quickr};
Tomasz Kalinowski’s R to Fortran transpiler – I had starred the repo a long time
ago (and in haste, accidentally unstarred it, then re-starred it) but never
really played with it. I’m familiar with slightly older Fortran; nowadays it’s
called “modern Fortran”, but I did my PhD using Fortran95 in the late 2000’s.
I’ve even pushed some of my postdoc code to GitHub
after getting it working again for a recent student.

I figured now was a great chance to have a proper play with the package.

{quickr} “transpiles” R code which means it takes R code converts the syntax
into Fortran syntax using the same variables and equivalent functions where
available. The idea being that when R isn’t working fast enough for you, instead
of re-writing your function in something like C++ (via {Rcpp}) it can
automatically write a Fortran version of your code and compile that into a
highly performant function which can be called with the same arguments. Faster
running code with no additional effort – sounds great!

The README for {quickr} has some examples highlighting how it can improve the
performance of some functions beyond what {Rcpp} can offer, in some cases
approaching C speeds. That’s not surprising to those who know Fortran – it’s
still very much used in theoretical physics partly because of the performance,
partly due to the existing support in that field, but also partly because despite
being an ‘old’ language, it’s actually pretty nice to use.

One of the big advantages of Fortran I found when learning other languages
after learning Fortran was that there’s no manual memory management. If you
want a vector or an array/tensor with many dimensions, you just ask for it
(specifying a size along each dimension or dynamically sizing, but never
manually freeing memory). R is known for its statistics chops, but under the
hood of some of these functions still call out to Fortran code.

I wanted some example code to try out myself and see if I even recognise the
Fortran it produces. I didn’t just want to use the example code from the package,
so what could I use?

In this post I
celebrated the fact that Julia has a ndigits() function, while in R I cheated
and used nchar() which works fine provided you’re dealing with non-negative
integers up to 99999, outside of which it doesn’t do what you want

nchar(99)       # 99    = 2 characters
## [1] 2
nchar(99999)    # 99999 = 5 characters
## [1] 5
nchar(99999+1)  # 1e+05 = 5 characters
## [1] 5
nchar(-99)      # -99   = 3 characters
## [1] 3

I had some interesting discussions on Mastodon
about different ways to implement ndigits() properly for R and in the end,
re-implementing the Julia solution seemed to work great for all edge cases. I
decided to use this for my Fortran testing with {quickr}.

I got the package installed and the compiler hooked up correctly so that I could
run the example code, then tried adapting it to the ndigits() problem.

R

The R code I started with was

nd_R <- function(x) {
  out <- double(length(x))
  x <- abs(x / 10)
  for (v in seq_along(x)) {
    d <- 1
    m <- 1
    while (m <= x[v]) {
      m <- m * 10
      d <- d + 1
    }
    out[v] <- d
  }
  out
}

nd_R(c(123456, 234, -72))
## [1] 6 3 2

and while this looks like a moderate amount of code, in essence it’s taking the
absolute value of the input (since we want to ignore negatives, and which is
nicely vectorised in R), dividing by 10, checking if we’ve exceeded the input
yet, and if not, stepping through successive multiples of 10 until we do, which
finds the first power of 10 that is greater than our value, indicating the
number of digits. For what it’s worth, this is why in that post I noted an
alternative route to achieving this; ceil(log10(x)).

Fortran

Hoping to immediately transpile this to Fortran, I immediately hit my first
snag; {quickr} hasn’t yet implemented while() so I can’t transpile this
exactly as I have it. There’s no early return() or break either, so I can’t
just exit an oversized loop early. Without an alternative, I’m going to cheat a
bit and just run a loop 12 times – this puts an upper limit on the input to a 12
digit number, but I can live with that.

Update while writing the post: I suppose good things come to those who wait –
on digging through some source code for this post I saw that while has been
implemented in the last week, so I’m going to pretend that was always the case.

The other piece this transpiler needs is a type declaration for the input; R
is fully dynamic in that a function can take any type of object and it’s up to
the function to decide what to do with it. Fortran is a bit stricter, and
requires types to be annotated, so I need to add a declare(type()) to the code.

nd2f <- function(x) {
  declare(
    type(x = double(NA))
  )
  out <- double(length(x))
  x <- abs(x / 10)
  for (v in seq_along(x)) {
    d <- 1
    m <- 1
    while (m <= x[v]) {
      m <- m * 10
      d <- d + 1
    }
    out[v] <- d
  }
  out
}

Note that this is still very much R code at this point – I can even run it
in R and get the same answers as before

nd2f(c(123456, 234, -72))
## [1] 6 3 2

What surprised me here is that declare() is a base R function (not from
{quickr}) intended for “specifying information about R code for use by the
interpreter, compiler, and code analysis tools”. I was originally thinking it
would be neat to be able to leverage that for some type-checking on the R side
as well as being informative to the Fortran code, but it “ignores the arguments
and returns NULL invisibly”, so no go on this throwing an error from R

int_id <- function(x) {
  declare(type(x = integer(NA)))
  x
}

int_id(3L)
## [1] 3
int_id(1.5)
## [1] 1.5

The magic happens when we ask {quickr} to do the transpilation.

The type information is used in the Fortran code, so compiling the id() example
produces something that is more restrictive on types

int_id_F <- quickr::quick(int_id)
int_id_F(3L)
## [1] 3
int_id_F(1.5)
## Error in int_id_F(1.5): typeof(x) must be 'integer', not 'double'

I can inspect the generated code with r2f(), though one wouldn’t normally need
to – it’s interesting to see what the Fortran code looks like

quickr:::r2f(int_id)
## subroutine int_id(x, x__len_) bind(c)
##   use iso_c_binding, only: c_int, c_ptrdiff_t
##   implicit none
## 
##   ! manifest start
##   ! sizes
##   integer(c_ptrdiff_t), intent(in), value :: x__len_
## 
##   ! args
##   integer(c_int), intent(in out) :: x(x__len_)
##   ! manifest end
## 
## 
## end subroutine
## 
## @r: function (x)
##   {
##       declare(type(x = integer(NA)))
##       x
##   }
## @closure: function (x)
##   {
##       declare(type(x = integer(NA)))
##       x
##   }

But of course, this just returns the value and that’s not particularly
enlightening. Doing the same for the ndigits code

quickr:::r2f(nd2f)
## subroutine nd2f(x, out, x__len_) bind(c)
##   use iso_c_binding, only: c_double, c_int, c_ptrdiff_t
##   implicit none
## 
##   ! manifest start
##   ! sizes
##   integer(c_ptrdiff_t), intent(in), value :: x__len_
## 
##   ! args
##   real(c_double), intent(in out) :: x(x__len_)
##   real(c_double), intent(out) :: out(x__len_)
## 
##   ! locals
##   integer(c_int) :: v
##   real(c_double) :: d
##   real(c_double) :: m
##   ! manifest end
## 
## 
##   out = 0
##   x = abs((x / 10.0_c_double))
##   do v = 1, size(x)
##     d = 1.0_c_double
##     m = 1.0_c_double
##     do while ((m <= x(v)))
##       m = (m * 10.0_c_double)
##       d = (d + 1.0_c_double)
##     end do
##     out(v) = d
##   end do
## end subroutine
## 
## @r: function (x)
##   {
##       declare(type(x = double(NA)))
##       out <- double(length(x))
##       x <- abs(x/10)
##       for (v in seq_along(x)) {
##           d <- 1
##           m <- 1
##           while (m <= x[v]) {
##               m <- m * 10
##               d <- d + 1
##           }
##           out[v] <- d
##       }
##       out
##   }
## @closure: function (x)
##   {
##       declare(type(x = double(NA)))
##       out <- double(length(x))
##       x <- abs(x/10)
##       for (v in seq_along(x)) {
##           d <- 1
##           m <- 1
##           while (m <= x[v]) {
##               m <- m * 10
##               d <- d + 1
##           }
##           out[v] <- d
##       }
##       out
##   }

The subroutine itself looks a lot like the R code; sure, some type annotations
are sprinkled around, do v = 1, size(x) replaces for v in seq_along(x) and
do while replaces while, but I don’t think it’s entirely alien.

What might surprise some is the line

x = abs((x / 10.0_c_double))

Notice there’s no loop around this? Fortran is an array language…
Rank-polymorphism, baby! I covered this in
another post of mine
but thanks to this, abs() is vectorised wherever needed

program test_abs
  implicit none
  integer, dimension(5) :: i = [-1, 2, -3, 4, -5]
  write(*,*) abs(i)
end program test_abs
#           1           2           3           4           5

Generating the compiled Fortran code from nd2f is as easy as

nd_F <- quickr::quick(nd2f)
nd_F
## function (x) 
## .External(<pointer: 0x11fddd73c>, x)

and we see that it’s referencing some external code. This can be called

nd_F(c(123456, 234, -72))
## [1] 6 3 2

with the big benefit that now it’s a LOT faster!

Generating a million random values and excluding any zero values, we can see the
40x performance increase (!!!)

set.seed(1)
nums <- round(runif(1e6, -1, 1) * 1e6)
nums <- nums[nums != 0]

b0 <- bench::mark(
  R = nd_R(nums),
  Fortran = nd_F(nums),
  min_iterations = 10
)

dplyr::arrange(b0[, 1:8], median)
## # A tibble: 2 × 6
##   expression      min   median `itr/sec` mem_alloc `gc/sec`
##   <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
## 1 Fortran      3.56ms   3.88ms    255.      15.3MB    266. 
## 2 R          154.27ms 154.42ms      6.48    15.3MB     25.9
plot(b0)

For those not familiar, this benchmark plot shows the individual times taken for
repeated executions of the code in each ‘expression’, grouped vertically by the
‘expression’ itself (annotated as the language here) with some random scatter to show
the spread of execution times. Points to the left are faster. It’s also worth
noting that bench::mark() defaults to check = TRUE so we can rest assured that
the results from each of the different languages we’re about to explore are
consistent and it’s not some artifact of one language doing less work.

If you run these yourself you’ll get slightly different results. I’m running them
on a newish M3 Macbook Pro.

All that performance increase from just adding one line to the R code and
wrapping it with one other function (resulting in an entirely different program
being written and compiled, producing the correct results).

I should note that in the first iteration of this post (in which while was not yet
supported) I used an excessive for loop which resulted in a
not-as-impressive-but-still-very-impressive 15x performance boost.

R (compiled)

If compiled code is so great, what about just compiling the R code with, e.g.
compiler::cmpfun()?

nd_comp = compiler::cmpfun(nd_R)

nd_comp(c(123456, 234, -72))
## [1] 6 3 2
b1 <- bench::mark(
  compiled = nd_comp(nums),
  R = nd_R(nums),
  Fortran = nd_F(nums),
  min_iterations = 10
)

dplyr::arrange(b1[, 1:8], median)
## # A tibble: 3 × 6
##   expression      min   median `itr/sec` mem_alloc `gc/sec`
##   <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
## 1 Fortran      3.58ms   3.93ms    253.      15.3MB   123.  
## 2 compiled   147.86ms    149ms      6.66    15.3MB     6.66
## 3 R          154.69ms 155.25ms      6.30    15.3MB     2.70
plot(b1)

That doesn’t help; by the time the benchmark was running the nd_R function had
been called enough times for it to be JIT compiled, anyway.

This did get me thinking, though – what about other compiled alternatives?

C

Since I’m going through Harvard’s CS50 ‘Introduction to Computer Science’ course
with R Contributors
to learn a bit more structured C I figured I’d add that via coolbutuseless’
{callme} package. This surely isn’t
the world’s greatest C code, but it compiles and runs…

callme::compile(
  "
#include <R.h>
#include <Rinternals.h>
#include <stdlib.h>
#include <math.h>

SEXP nd_C(SEXP vec) {
  double *vec_ptr = REAL(vec);
  SEXP res = PROTECT(allocVector(REALSXP, length(vec)));
  double *res_ptr = REAL(res);
  for (int i = 0; i < length(vec); i++) {
    double abs_x = fabs(vec_ptr[i] / 10.0);
        int d = 1;
        double m = 1.0;
        while (m <= abs_x) {
            m *= 10.0;
            d++;
        }
        res_ptr[i] = d;
  }

  UNPROTECT(1);
  return res;
}
"
)

nd_C(c(123456, 234, -72))
## [1] 6 3 2

So, how does it compare?

b2 <- bench::mark(
  C = nd_C(nums),
  R = nd_R(nums),
  Fortran = nd_F(nums),
  min_iterations = 10
)

dplyr::arrange(b2[, 1:8], median)
## # A tibble: 3 × 6
##   expression      min   median `itr/sec` mem_alloc `gc/sec`
##   <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
## 1 Fortran       3.3ms   3.58ms    277.     15.26MB   107.  
## 2 C            3.92ms   4.14ms    240.      7.63MB    50.6 
## 3 R          147.99ms 149.25ms      6.71   15.26MB     4.47
plot(b2)

Whoa – automatically transpiled Fortran runs faster than (my) C… That’s fast.

Impressively fast

Impressively fast

C++

What about C++ via {Rcpp}? Dealing with vectors is made easier by {Rcpp} having
pre-built types compatible with R, and this otherwise looks very similar to the
R code

nd_Rcpp <- Rcpp::cppFunction(
  "
IntegerVector nd(NumericVector x) {
    int n = x.size();
    IntegerVector out(n);

    for (int v = 0; v < n; v++) {
        double abs_x = std::abs(x[v] / 10.0);
        int d = 1;
        int m = 1;
        while (m <= abs_x) {
            m *= 10;
            d++;
        }
        out[v] = d;
    }

    return out;
}
"
)

nd_Rcpp(c(123456, 234, -72))
## [1] 6 3 2
b3 <- bench::mark(
  `C++` = nd_Rcpp(nums),
  C = nd_C(nums),
  R = nd_R(nums),
  Fortran = nd_F(nums),
  min_iterations = 10
)

dplyr::arrange(b3[, 1:8], median)
## # A tibble: 4 × 6
##   expression      min   median `itr/sec` mem_alloc `gc/sec`
##   <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
## 1 C++          3.02ms   3.12ms    320.      3.82MB    16.6 
## 2 Fortran      3.29ms   3.58ms    279.     15.26MB    92.9 
## 3 C            3.75ms   3.93ms    255.      7.63MB    24.9 
## 4 R          148.52ms 148.85ms      6.71   15.26MB     1.68
plot(b3)

This one seems to wander around a bit; on different runs I’ve seen performance
equal or better to the C code and on others, about 3x as long, but generally
pretty fast.

Julia

After all of this, I remembered that I was comparing the Julia implementation –
how does that perform? Julia is a JIT/AOT compiled language, so maybe it’s not
too bad… I can still call that directly from R

JuliaCall::julia_eval("ndigits.([123456, 234, -72])")
## [1] 6 3 2

keeping in mind that the Julia function ndigits (the implementation for which
I’ve borrowed for all of the examples, so we are dealing with the same
algorithm in each case) is in fact compiled, but available as ndigits(). As
long as I make the vector available in a Julia session (as integers; the
function is only defined for integers) I can run this

JuliaCall::julia_assign("nums", as.integer(nums))

b4 <- bench::mark(
  Julia = JuliaCall::julia_eval("ndigits.(nums)"),
  `C++` = nd_Rcpp(nums),
  C = nd_C(nums),
  R = nd_R(nums),
  Fortran = nd_F(nums),
  min_iterations = 10
)

dplyr::arrange(b4[, 1:8], median)
## # A tibble: 5 × 6
##   expression      min   median `itr/sec` mem_alloc `gc/sec`
##   <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
## 1 C++          3.02ms   3.06ms    324.      3.81MB    26.6 
## 2 Julia         3.1ms   3.32ms    266.      3.81MB    22.0 
## 3 Fortran      3.29ms    3.6ms    270.     15.26MB    84.2 
## 4 C            3.74ms   3.94ms    251.      7.63MB    50.7 
## 5 R          157.75ms 157.91ms      6.33   15.26MB     2.71
plot(b4)

Ten points to Julia – remember, this is an interpreted language.

That’s really fast!

That’s really fast!

I should note there’s work being done towards making Julia binaries out of scripts, but this still has a startup time
of a few dozen milliseconds for even a Hello, World example.

Rust

One more? What about Rust? We can use {rextendr} to call Rust code inline,
making sure to target the release profile for maximum performance

rextendr::rust_function(
  r"(
  fn nd_Rust(x: &[f64]) -> Vec<i32> {
    let mut out = vec![0; x.len()];
    for v in 0..x.len() {
        let abs_x = (x[v].abs() / 10.0);
        let mut d = 1;
        let mut m = 1.0;
        while m <= abs_x {
            m *= 10.0;
            d += 1;
        }
        out[v] = d;
    }
    out
  }
)",
  profile = "release"
)
## ℹ build directory: '/private/var/folders/1h/k6c5hb4d2qx07m8kfqb54f9c0000gn/T/RtmppiKq7x/file8b2f28646e2'
## ✔ Writing '/private/var/folders/1h/k6c5hb4d2qx07m8kfqb54f9c0000gn/T/RtmppiKq7x/file8b2f28646e2/target/extendr_wrappers.R'
nd_Rust(c(123456, 234, -72))
## [1] 6 3 2
b5 <- bench::mark(
  Rust = nd_Rust(nums),
  Julia = JuliaCall::julia_eval("ndigits.(nums)"),
  `C++` = nd_Rcpp(nums),
  C = nd_C(nums),
  R = nd_R(nums),
  Fortran = nd_F(nums),
  min_iterations = 10
)

dplyr::arrange(b5[, 1:8], median)
## # A tibble: 6 × 6
##   expression      min   median `itr/sec` mem_alloc `gc/sec`
##   <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
## 1 Rust         2.87ms   3.07ms    322.      3.82MB   24.9  
## 2 C++          3.03ms   3.13ms    318.      3.81MB   22.1  
## 3 Fortran      3.28ms   3.56ms    283.     15.26MB   65.5  
## 4 Julia        3.24ms   3.66ms    265.      3.81MB   17.8  
## 5 C            3.75ms   3.94ms    255.      7.63MB   22.5  
## 6 R          150.69ms 151.09ms      6.61   15.26MB    0.734
plot(b5)

Ridiculous speeds!

Ridiculous speeds!

We are truly spoiled for choice these days – not only do we have a plethora of
languages we can call directly from R, but several languages which run faster than
even (at least my implementation of) C and count number of digits of
a million values in under 4ms.

Python

Just for funsies, what about Python? It’s not a compiled language, but maybe if
I use numpy it will be fast … ? It’s at least another language I can call from
R that is generally considered ‘faster’. Is it?

library(reticulate)
reticulate::py_run_string('
import numpy as np
def nd_python(x):
    x = np.asarray(x)
    out = np.zeros(len(x), dtype=int)

    for v in range(len(x)):
        abs_x = abs(x[v] / 10.0)
        d = 1
        m = 1
        while m <= abs_x:
            m *= 10
            d += 1
        out[v] = d

    return out.tolist()
')

py$nd_python(c(123456, 234, -72))
## [1] 6 3 2
b6 <- bench::mark(
  Python = py$nd_python(nums),
  Rust = nd_Rust(nums),
  Julia = JuliaCall::julia_eval("ndigits.(nums)"),
  `C++` = nd_Rcpp(nums),
  C = nd_C(nums),
  R = nd_R(nums),
  Fortran = nd_F(nums),
  min_iterations = 10
)

dplyr::arrange(b6[, 1:8], median)
## # A tibble: 7 × 6
##   expression      min   median `itr/sec` mem_alloc `gc/sec`
##   <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
## 1 Rust         2.87ms   3.08ms    322.      3.81MB   13.0  
## 2 C++          3.02ms   3.13ms    321.      3.81MB    6.21 
## 3 Fortran      3.27ms   3.44ms    285.     15.26MB   32.7  
## 4 Julia        3.33ms   3.53ms    280.      3.81MB    7.49 
## 5 C            3.74ms   3.94ms    255.      7.63MB   10.6  
## 6 R           157.9ms 158.51ms      6.31   15.26MB    0.701
## 7 Python     269.08ms 271.27ms      3.64    3.81MB    0
plot(b6)

In fairness, there’s overhead here involved with calling it from R, but I think
that’s apples-to-apples considering I’m doing the same with all the compiled
languages.

Does it scale?

I’ve been running these benchmarks for a million numbers, but how do the results
scale with that size? What if it’s just a handful of numbers? What about in
between these extremes? Running the benchmarks at various scales should show this.

n_vals <- 10^(1:7)
scales <- purrr::map_df(n_vals, ~{
  set.seed(1)
  nums <- round(runif(.x, -1, 1) * .x)
  nums <- nums[nums != 0]
  JuliaCall::julia_assign("nums", as.integer(nums))
  b <- bench::mark(
    Python = py$nd_python(nums),
    Rust = nd_Rust(nums),
    Julia = JuliaCall::julia_eval("ndigits.(nums)"),
    `C++` = nd_Rcpp(nums),
    C = nd_C(nums),
    R = nd_R(nums),
    Fortran = nd_F(nums),
    min_iterations = 10,
    check = TRUE
  )
  dplyr::bind_cols(vec_len = .x, b[, 1:8])
})

library(ggplot2)
ggplot(scales,
       aes(x = vec_len,
           y = 1e6*as.numeric(median),
           col = as.character(expression)
       )) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  scale_x_log10() +
  scale_y_log10() +
  scale_color_discrete(palette = "Set2") +
  labs(
    title = "Scaling of Counting ndigits Benchmarks",
    x = "Vector Length",
    y = "Microseconds",
    color = "Language"
  ) +
  theme_bw()

What a nice, log-log linear result with that one exception – Julia is pretty
constant up until 1000, after which it starts to follow the same trajectory as
the other languages – presumably that’s just the overhead of starting up the
Julia runtime, which is a known bottleneck.

There’s definitely a clear divide between the interpreted languages (R and Python)
and the compiled ones.

At lower vector lengths there’s a little bit of a spread with Fortran really showing
off at the lowest lengths

dplyr::arrange(scales[scales$vec_len == 10, ], median)
## # A tibble: 7 × 7
##   vec_len expression      min   median `itr/sec` mem_alloc `gc/sec`
##     <dbl> <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
## 1      10 Fortran     122.9ns 205.12ns  4081331.        0B     0   
## 2      10 C++         287.1ns 369.04ns  2585126.        0B     0   
## 3      10 C             410ns 491.97ns  1914470.        0B     0   
## 4      10 Rust        532.9ns  697.1ns  1367702.        0B     0   
## 5      10 R             943ns   1.11µs   878760.        0B     0   
## 6      10 Python       20.7µs  22.92µs    43281.        0B     4.33
## 7      10 Julia        69.8µs  71.59µs    13636.        0B     2.03

but we’re looking at sub microsecond differences – what will you do with all
that free time?

By the time we’re looking at 1000 values, the compiled languages are all about
the same

dplyr::arrange(scales[scales$vec_len == 1000, ], median)
## # A tibble: 7 × 7
##   vec_len expression      min   median `itr/sec` mem_alloc `gc/sec`
##     <dbl> <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
## 1    1000 C++          1.68µs   1.84µs   487745.    3.95KB     0   
## 2    1000 Rust         2.13µs   2.46µs   397120.    3.95KB     0   
## 3    1000 Fortran      1.97µs   2.79µs   347513.   15.72KB    34.8 
## 4    1000 C            2.79µs   3.08µs   299955.    7.86KB    30.0 
## 5    1000 Julia       73.06µs  77.49µs    12381.    3.95KB     0   
## 6    1000 R           77.53µs  81.22µs    12306.   15.72KB     0   
## 7    1000 Python      213.9µs 232.35µs     4282.    3.95KB     2.06

At ten million values it’s a complete wash the compiled languages with maybe a
slight drop for C

dplyr::arrange(scales[scales$vec_len == 1e7, ], median)
## # A tibble: 7 × 7
##    vec_len expression      min   median `itr/sec` mem_alloc `gc/sec`
##      <dbl> <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
## 1 10000000 Rust        33.78ms  34.58ms    28.9      38.1MB   2.07  
## 2 10000000 C++         35.05ms  35.38ms    28.3      38.1MB   2.02  
## 3 10000000 Julia       37.63ms  38.56ms    18.1      38.1MB   2.01  
## 4 10000000 Fortran     39.29ms   41.9ms    23.8     152.6MB  17.0   
## 5 10000000 C           43.99ms  44.19ms    22.5      76.3MB   2.05  
## 6 10000000 R              1.8s    1.81s     0.550   152.6MB   0.236 
## 7 10000000 Python        3.06s     3.1s     0.322    38.1MB   0.0357

All very interesting!

It would probably be worthwhile digging into the memory usage of all of these
since there’s a big difference that likely indicates something different is
happening, but that’s beyond my understanding – feel free to let me know!

So, what might be the reason for Rust and Julia to be so fast, even compared to
C? These are newer languages with a lot of focus on their compilers, and it’s
entirely possible that they’re able to make some better optimisations compared
to a very general C compiler, but more likely that’s the upper limit of what a
computer can do in that much time and my C code is non-optimal.

Conclusions

Back to the original point, though – the transpilation does an amazing job
of improving the code without having to write more code in a different
language
. Sure, Julia solves this ‘two language problem’ by just being
ridiculously fast to begin with, but if I am writing R code, it’s fantastic to
see there’s an option for just “making it go brrr” without actually doing
anything extra.

Not all of R has been translated to Fortran so there’s a lot of code that won’t
transpile just yet, but it’s a truly inspiring project that I’ll surely be
keeping a close eye on.

I’d love to hear what people think about these comparisons – are there points I’ve
overlooked? Better ways to do it? Improvements to my implementations which change
the results? Other considerations I’ve missed? As always, I can be found on
Mastodon and the comment section below.

devtools::session_info()
## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value
##  version  R version 4.4.1 (2024-06-14)
##  os       macOS 15.5
##  system   aarch64, darwin20
##  ui       X11
##  language (EN)
##  collate  en_US.UTF-8
##  ctype    en_US.UTF-8
##  tz       Australia/Adelaide
##  date     2025-06-29
##  pandoc   3.4 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/aarch64/ (via rmarkdown)
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package      * version    date (UTC) lib source
##  beeswarm       0.4.0      2021-06-01 [1] CRAN (R 4.4.1)
##  bench          1.1.4      2025-01-16 [1] CRAN (R 4.4.1)
##  blogdown       1.21.1     2025-06-28 [1] Github (rstudio/blogdown@33313a5)
##  bookdown       0.41       2024-10-16 [1] CRAN (R 4.4.1)
##  brio           1.1.5      2024-04-24 [1] CRAN (R 4.4.0)
##  bslib          0.8.0      2024-07-29 [1] CRAN (R 4.4.0)
##  cachem         1.1.0      2024-05-16 [1] CRAN (R 4.4.0)
##  callme         0.1.10     2024-07-27 [1] CRAN (R 4.4.0)
##  cli            3.6.4      2025-02-13 [1] CRAN (R 4.4.1)
##  codetools      0.2-20     2024-03-31 [1] CRAN (R 4.4.1)
##  devtools       2.4.5      2022-10-11 [1] CRAN (R 4.4.0)
##  dichromat      2.0-0.1    2022-05-02 [1] CRAN (R 4.4.1)
##  digest         0.6.37     2024-08-19 [1] CRAN (R 4.4.1)
##  dotty          0.1.0      2024-08-30 [1] CRAN (R 4.4.1)
##  dplyr          1.1.4      2023-11-17 [1] CRAN (R 4.4.0)
##  ellipsis       0.3.2      2021-04-29 [1] CRAN (R 4.4.0)
##  evaluate       1.0.3      2025-01-10 [1] CRAN (R 4.4.1)
##  farver         2.1.2      2024-05-13 [1] CRAN (R 4.4.0)
##  fastmap        1.2.0      2024-05-15 [1] CRAN (R 4.4.0)
##  fs             1.6.5      2024-10-30 [1] CRAN (R 4.4.1)
##  generics       0.1.3      2022-07-05 [1] CRAN (R 4.4.0)
##  ggbeeswarm     0.7.2      2023-04-29 [1] CRAN (R 4.4.0)
##  ggplot2      * 3.5.2.9001 2025-06-15 [1] Github (tidyverse/ggplot2@9f80c8c)
##  glue           1.8.0      2024-09-30 [1] CRAN (R 4.4.1)
##  gtable         0.3.6      2024-10-25 [1] CRAN (R 4.4.1)
##  here           1.0.1      2020-12-13 [1] CRAN (R 4.4.0)
##  htmltools      0.5.8.1    2024-04-04 [1] CRAN (R 4.4.0)
##  htmlwidgets    1.6.4      2023-12-06 [1] CRAN (R 4.4.0)
##  httpuv         1.6.15     2024-03-26 [1] CRAN (R 4.4.0)
##  jquerylib      0.1.4      2021-04-26 [1] CRAN (R 4.4.0)
##  jsonlite       2.0.0      2025-03-27 [1] CRAN (R 4.4.1)
##  JuliaCall      0.17.6     2024-12-07 [1] CRAN (R 4.4.1)
##  knitr          1.50       2025-03-16 [1] CRAN (R 4.4.1)
##  later          1.4.1      2024-11-27 [1] CRAN (R 4.4.1)
##  lattice        0.22-6     2024-03-20 [1] CRAN (R 4.4.1)
##  lifecycle      1.0.4      2023-11-07 [1] CRAN (R 4.4.0)
##  magrittr       2.0.3      2022-03-30 [1] CRAN (R 4.4.0)
##  Matrix         1.7-1      2024-10-18 [1] CRAN (R 4.4.1)
##  memoise        2.0.1      2021-11-26 [1] CRAN (R 4.4.0)
##  mime           0.12       2021-09-28 [1] CRAN (R 4.4.0)
##  miniUI         0.1.1.1    2018-05-18 [1] CRAN (R 4.4.0)
##  pillar         1.10.1     2025-01-07 [1] CRAN (R 4.4.1)
##  pkgbuild       1.4.7      2025-03-24 [1] CRAN (R 4.4.1)
##  pkgconfig      2.0.3      2019-09-22 [1] CRAN (R 4.4.0)
##  pkgload        1.4.0      2024-06-28 [1] CRAN (R 4.4.0)
##  png            0.1-8      2022-11-29 [1] CRAN (R 4.4.0)
##  processx       3.8.6      2025-02-21 [1] CRAN (R 4.4.1)
##  profmem        0.7.0      2025-05-02 [1] CRAN (R 4.4.1)
##  profvis        0.4.0      2024-09-20 [1] CRAN (R 4.4.1)
##  promises       1.3.2      2024-11-28 [1] CRAN (R 4.4.1)
##  ps             1.9.0      2025-02-18 [1] CRAN (R 4.4.1)
##  purrr          1.0.4      2025-02-05 [1] CRAN (R 4.4.1)
##  quickr         0.1.0.9000 2025-06-29 [1] Github (t-kalinowski/quickr@254b4d0)
##  R6             2.6.1      2025-02-15 [1] CRAN (R 4.4.1)
##  RColorBrewer   1.1-3      2022-04-03 [1] CRAN (R 4.4.0)
##  Rcpp           1.0.14     2025-01-12 [1] CRAN (R 4.4.1)
##  remotes        2.5.0      2024-03-17 [1] CRAN (R 4.4.1)
##  reticulate   * 1.42.0     2025-03-25 [1] CRAN (R 4.4.1)
##  rextendr       0.3.1      2023-06-20 [1] CRAN (R 4.4.0)
##  rlang          1.1.5      2025-01-17 [1] CRAN (R 4.4.1)
##  rmarkdown      2.28       2024-08-17 [1] CRAN (R 4.4.0)
##  rprojroot      2.0.4      2023-11-05 [1] CRAN (R 4.4.0)
##  rstudioapi     0.17.1     2024-10-22 [1] CRAN (R 4.4.1)
##  S7             0.2.0      2024-11-07 [1] CRAN (R 4.4.1)
##  sass           0.4.9      2024-03-15 [1] CRAN (R 4.4.0)
##  scales         1.4.0      2025-04-24 [1] CRAN (R 4.4.1)
##  sessioninfo    1.2.2      2021-12-06 [1] CRAN (R 4.4.0)
##  shiny          1.9.1      2024-08-01 [1] CRAN (R 4.4.0)
##  stringi        1.8.4      2024-05-06 [1] CRAN (R 4.4.0)
##  tibble         3.2.1      2023-03-20 [1] CRAN (R 4.4.0)
##  tidyr          1.3.1      2024-01-24 [1] CRAN (R 4.4.0)
##  tidyselect     1.2.1      2024-03-11 [1] CRAN (R 4.4.0)
##  urlchecker     1.0.1      2021-11-30 [1] CRAN (R 4.4.0)
##  usethis        3.1.0.9000 2025-03-31 [1] Github (r-lib/usethis@a653d6e)
##  utf8           1.2.4      2023-10-22 [1] CRAN (R 4.4.0)
##  vctrs          0.6.5      2023-12-01 [1] CRAN (R 4.4.0)
##  vipor          0.4.7      2023-12-18 [1] CRAN (R 4.4.1)
##  withr          3.0.2      2024-10-28 [1] CRAN (R 4.4.1)
##  xfun           0.51       2025-02-19 [1] CRAN (R 4.4.1)
##  xtable         1.8-4      2019-04-21 [1] CRAN (R 4.4.0)
##  yaml           2.3.10     2024-07-26 [1] CRAN (R 4.4.0)
## 
##  [1] /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library
## 
## ─ Python configuration ───────────────────────────────────────────────────────
##  python:         /Users/jono/Library/Caches/org.R-project.R/R/reticulate/uv/cache/archive-v0/2tc-cviHm3ODucI_hIfUb/bin/python3
##  libpython:      /Users/jono/Library/Caches/org.R-project.R/R/reticulate/uv/python/cpython-3.11.12-macos-aarch64-none/lib/libpython3.11.dylib
##  pythonhome:     /Users/jono/Library/Caches/org.R-project.R/R/reticulate/uv/cache/archive-v0/2tc-cviHm3ODucI_hIfUb:/Users/jono/Library/Caches/org.R-project.R/R/reticulate/uv/cache/archive-v0/2tc-cviHm3ODucI_hIfUb
##  virtualenv:     /Users/jono/Library/Caches/org.R-project.R/R/reticulate/uv/cache/archive-v0/2tc-cviHm3ODucI_hIfUb/bin/activate_this.py
##  version:        3.11.12 (main, Apr  9 2025, 03:49:53) [Clang 20.1.0 ]
##  numpy:          /Users/jono/Library/Caches/org.R-project.R/R/reticulate/uv/cache/archive-v0/2tc-cviHm3ODucI_hIfUb/lib/python3.11/site-packages/numpy
##  numpy_version:  2.3.1
##  
##  NOTE: Python version was forced by py_require()
## 
## ──────────────────────────────────────────────────────────────────────────────

Rotation with Modulo

By: Jonathan Carroll

Re-posted from: https://jcarroll.com.au/2025/05/03/rotation-with-modulo/

How well do you know your fundamental operators in different languages? ‘Easy’
examples help to fortify that knowledge, and comparing across languages makes
for some neat implementation detail discoveries.

I saw this toot from @gregeganSF on Mastodon

Post by @gregeganSF@mathstodon.xyz
View on Mastodon

which says

To rotate a j-digit number n by k digits, if n≠10^j-1 there is a simple formula:

rot(n,j,k) = (n*10^k) mod (10^j-1)

e.g. 1234 * 100 mod 9999 = 3412

Why? The mod subtracts (10^j-1) times the leftmost k digits of n*10^k, removing them from the left and adding them on the right.

and my first thought was “ooh, that’s cool”, but my second thought was “I’m
going to implement this in a bunch of languages!”. Sure, it’s a very small bit
of math to implement without any particular sharp edges of iteration/recursion,
but that means I’ll be working with some basic functionality and I believe it’s
very important to have that locked in comfortably. Let’s see how it goes!

I’m not aware of a “name” for this as such. Writing it out a little more styled

\[
rot(n, j, k) = (n\times10^k)\ \%\ (10^j-1)
\]

it looks like I’ll just need a ‘power’ and a ‘modulo’. The \(j\) there is the
number of digits in \(n\) and sure, we could count that ourselves and pass it as
an argument, but even better might be to calculate it as well. That means
figuring out how many digits are in a number.

As always, my go-to is R, so let’s start there.

R

In R the power operator is ^. Also **, but that’s almost never used –
there’s even a note about that in the documentation

** is translated in the parser to ^, but this was undocumented for many years.

Modulo is %% which I can never remember because it’s similar to integer
division which is %/%.

To get the number of digits we can use the fact that nchar() will first convert
its input into a character vector, so 12345 becomes "12345" and thus the
number of characters of that is the number of digits. If that wasn’t the case
I could still do

ceiling(log10(314159))
## [1] 6

but the nchar() approach seems fine. Putting those pieces together into a
function which takes the number and how many places to move it, I get

rot <- function(n, k) {
  (n*10^k) %% (10^nchar(n)-1)
}
rot(12345, 3)
## [1] 45123

You can see that the values are ‘rotated’ cyclically by 3 places (to the left).

R doesn’t have a built-in way to achieve this even for a vector of values
(rotating ‘digits’ of an integer is a toy problem that is unlikely to actually
come up in real situations). One solution is my
{vec} package which does implement a
ring-buffer sort of effect for vectors

# remotes::install_github("jonocarroll/vec")
library(vec)
v <- as_vec(1:5)
rotate(v, 3)
## [1] 4 5 1 2 3

Under the hood this uses modulo on the indices

x <- 1:5; n <- 3
x[(((n - 1 + seq_along(x)) %% length(x))) + 1]
## [1] 4 5 1 2 3

One extra feature of this is it also takes negative values to shift the other way

rotate(v, -3)
## [1] 3 4 5 1 2

whereas the rot() implementation above can’t.

Julia

As is almost always the case, the Julia functionality looks closer to what one
might “expect” from translating maths or stats; the power operator remains ^,
but the modulo operator is a more familiar %. There is an actual ndigits()
function to get the number of digits which, as far as I can tell from the source,
doesn’t first convert to character. The examples for that function do highlight
a failing of the R approach – if a value is negative then as.character() will
produce the wrong number of digits. In R:

as.character(-1234)
## [1] "-1234"
nchar(-1234)
## [1] 5

while in Julia

ndigits(-1234)
## 4

We’re not dealing with negatives here, but that’s certainly a gotcha.

Implementing the Julia function can be done in a single line so no need for a
function keyword

rot(n,k) = (n*10^k) % (10^ndigits(n)-1)
## rot (generic function with 1 method)
rot(12345, 3)
## 45123

If we were working with vectors, Julia also has a built-in way to do the
cyclic rotation, though it seems to shift in the opposite direction

x = [1, 2, 3, 4, 5]
## 5-element Vector{Int64}:
##  1
##  2
##  3
##  4
##  5
circshift(x, -3)
## 5-element Vector{Int64}:
##  4
##  5
##  1
##  2
##  3
circshift(x, 3)
## 5-element Vector{Int64}:
##  3
##  4
##  5
##  1
##  2

Here ends the Algol-language portion of the post, and I’ll move on to some
languages where these operations are even more fundamental…

APL

When I see ‘straightforward’ math operations on numbers I now think APL because
that’s what it was originally built for and it works so very well; each math
operation you could perform on an array of values has a ‘glyph’ representing it
so should be a better ‘translation’ of the math on a blackboard directly to code.

One thing R users will immediately recognise is the assignment glyph ; yes
that’s a single glyph, not R’s <-, but it works the same as in R.

You’re probably familiar with the multiplication glyph × and addition +.
The power/exponentiation glyph is *. Nothing too surprising there, I hope.

Because - is the ‘subtract’ operation, there’s a distinct glyph for ‘negative’
¯ (a raised hyphen) so it isn’t confused with subtraction. Modulo takes some
more inspiration from math and is |. The only other potentially confusing
glyphs are length and format which, when combined, do something very
similar to R’s “how many characters does this use?”. Again the ‘negative number’
problem is here, but we’re not worried about that in this case.

Putting those pieces together requires knowing that APL evaluates right-to-left,
with the argument to the right of an operator in a “function” being denoted by
omega () and the one to the left being alpha (). The function I came up
with looks like this

    rot←{(¯1+10*≢⍕⍵)|⍵×10*⍺}

It’s entirely possible that it can be shortened or improved; I have a tendency
to overlook where parentheses are really required and opportunities for
simplification. Nonetheless, if you read from right-to-left it spells out the
calculation we want. Applying it to some value means placing it between its two
arguments

    3 rot 12345
45123

Try it out yourself!

Most of the difficulty I faced when building this was dealing with order-of-operations
which need to be right-to-left. There are ways to ‘swap’ the order of arguments
to a function (such as modulo) to make it read more similarly to the hand-written
expression, but I both couldn’t get that to produce the right answer and didn’t
feel it was necessary.

In terms of working with vectors, that’s where APL shines. There is a rotate
glyph which when given just one argument reverses a vector, but with a second
argument does exactly what we want; rotates by that many places

    x←1 2 3 4 5
    3 ⌽ x
4 5 1 2 3

If you don’t look too closely at the type of the data, we can use this to rotate
a string made of character digits by again using format to make a vector of
characters from the number

    3 ⌽ ⍕12345
45123

(the whitespace here is purely for demonstration; 3⌽⍕12345 works just the same).

Uiua

Uiua is a much newer language that has a lot of support for operating on data,
but it behaves differently to all of the above; it’s a stack-based language so
you work with data on a stack, not as variables. I’ve played around with it and
really enjoy working with it – there’s even an
Exercism track now – but in trying
to write this solution I realised that I’d only ever worked with one ‘thing’ on
the stack, even if that was an entire vector of values. This problem invites us
to work with a value to be rotated and how many places to rotate it; that
meant I learned a lot about figuring out which value from the stack I want.

Entering operators into Uiua is greatly eased by having translation and
auto-complete; you don’t have to figure out how to type you can start typing
mod and as long as it’s a unique completion, Uiua will convert it to the
appropriate glyph. Additionally, there are some formatting tricks such as taking
a double-underscore suffix to a function to make a combined glyph with a preset
value; log__10 translates to ₙ₁₀.

The operators I need here are modulo , multiply ×, power , log10 ₙ₁₀,
ceiling , and subtract -, with the additional dip to use the second
value on the stack, and backward ˜ to swap the arguments of modulo. I couldn’t
immediately think of a way to cleanly get the number of digits of a value (I did
later, which I’ll come back to) so I went the log10 route and my solution is
(again, right-to-left)

3 12345
˜◿⊙(-1˜ⁿ10⌈ₙ₁₀)⊸טⁿ10
45123

Try it out yourself!

Working with vectors is much cleaner here, too, and there’s a simple rotate
that does the same as APL

↻ 3 [1 2 3 4 5]
[4 5 1 2 3]

Reading this, the vector is placed on the stack, then the value 3 is put on the
top of the stack, then rotate takes two values from the stack and performs that
operation, leaving one value (the result) on the stack.

Uiua also has a really cool feature of “undoing” an operation, where the inverse
can be calculated. If I wanted to turn a string into a number I would use parse
and I can do the opposite by prefixing it with un ° to make “unparse”
which converts a number to a string. Since a string is just a vector of
characters (in this case digits) the rotate works naturally, albeit returning a
string

↻ 3 °⋕ 12345
"45123"

Uiua goes one step further with an under which takes some value, performs
some transformation, applies a function, then undoes that transformation. In my
case, I want to “unparse then re-parse” which seems like a great fit for this

⍜°⋕(↻ 3) 12345
45123

and returns a number again, because under applies the back-transformation of parse.

The unparse °⋕ is recognised as a compound and is passed as the first argument
to under, while I need the rotate and 3 to go together with parentheses. If I
always wanted to shift by 3 places I could use the double underscore to ‘attach’
the 3 to the rotate producing rotate__3 ↻₃ but what I have above allows for
changing that number.

Looking into it this way, it’s more obvious that I could get the number of digits
with lengthunparse as ⧻°⋕, but exploring how to go the long way around was
entirely worthwhile, and not necessarily longer with ceilinglog__10 as ⌈ₙ₁₀.


I knew I’d find lots of interesting things when I saw this and I was right – just
spending the time going through the documentation of these ‘basic’ functions
reminded me of things I’ve forgotten and some new things I don’t think I knew
before.

I’d love to hear what people think about these comparisons – are there points I’ve
overlooked? Better ways to do it? Different functions in some other languages?
Considerations I’ve missed? As always, I can be found on
Mastodon and the comment section below.

devtools::session_info()
## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value
##  version  R version 4.4.1 (2024-06-14)
##  os       macOS 15.4.1
##  system   aarch64, darwin20
##  ui       X11
##  language (EN)
##  collate  en_US.UTF-8
##  ctype    en_US.UTF-8
##  tz       Australia/Adelaide
##  date     2025-05-03
##  pandoc   3.4 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/aarch64/ (via rmarkdown)
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package     * version    date (UTC) lib source
##  blogdown      1.19       2024-02-01 [1] CRAN (R 4.4.0)
##  bookdown      0.41       2024-10-16 [1] CRAN (R 4.4.1)
##  bslib         0.8.0      2024-07-29 [1] CRAN (R 4.4.0)
##  cachem        1.1.0      2024-05-16 [1] CRAN (R 4.4.0)
##  cli           3.6.4      2025-02-13 [1] CRAN (R 4.4.1)
##  devtools      2.4.5      2022-10-11 [1] CRAN (R 4.4.0)
##  digest        0.6.37     2024-08-19 [1] CRAN (R 4.4.1)
##  ellipsis      0.3.2      2021-04-29 [1] CRAN (R 4.4.0)
##  evaluate      1.0.3      2025-01-10 [1] CRAN (R 4.4.1)
##  fastmap       1.2.0      2024-05-15 [1] CRAN (R 4.4.0)
##  fs            1.6.5      2024-10-30 [1] CRAN (R 4.4.1)
##  glue          1.8.0      2024-09-30 [1] CRAN (R 4.4.1)
##  htmltools     0.5.8.1    2024-04-04 [1] CRAN (R 4.4.0)
##  htmlwidgets   1.6.4      2023-12-06 [1] CRAN (R 4.4.0)
##  httpuv        1.6.15     2024-03-26 [1] CRAN (R 4.4.0)
##  jquerylib     0.1.4      2021-04-26 [1] CRAN (R 4.4.0)
##  jsonlite      2.0.0      2025-03-27 [1] CRAN (R 4.4.1)
##  JuliaCall     0.17.6     2024-12-07 [1] CRAN (R 4.4.1)
##  knitr         1.48       2024-07-07 [1] CRAN (R 4.4.0)
##  later         1.4.1      2024-11-27 [1] CRAN (R 4.4.1)
##  lifecycle     1.0.4      2023-11-07 [1] CRAN (R 4.4.0)
##  magrittr      2.0.3      2022-03-30 [1] CRAN (R 4.4.0)
##  memoise       2.0.1      2021-11-26 [1] CRAN (R 4.4.0)
##  mime          0.12       2021-09-28 [1] CRAN (R 4.4.0)
##  miniUI        0.1.1.1    2018-05-18 [1] CRAN (R 4.4.0)
##  pkgbuild      1.4.7      2025-03-24 [1] CRAN (R 4.4.1)
##  pkgload       1.4.0      2024-06-28 [1] CRAN (R 4.4.0)
##  profvis       0.4.0      2024-09-20 [1] CRAN (R 4.4.1)
##  promises      1.3.2      2024-11-28 [1] CRAN (R 4.4.1)
##  purrr         1.0.4      2025-02-05 [1] CRAN (R 4.4.1)
##  R6            2.6.1      2025-02-15 [1] CRAN (R 4.4.1)
##  Rcpp          1.0.14     2025-01-12 [1] CRAN (R 4.4.1)
##  remotes       2.5.0      2024-03-17 [1] CRAN (R 4.4.1)
##  rlang         1.1.5      2025-01-17 [1] CRAN (R 4.4.1)
##  rmarkdown     2.28       2024-08-17 [1] CRAN (R 4.4.0)
##  rstudioapi    0.17.1     2024-10-22 [1] CRAN (R 4.4.1)
##  sass          0.4.9      2024-03-15 [1] CRAN (R 4.4.0)
##  sessioninfo   1.2.2      2021-12-06 [1] CRAN (R 4.4.0)
##  shiny         1.9.1      2024-08-01 [1] CRAN (R 4.4.0)
##  urlchecker    1.0.1      2021-11-30 [1] CRAN (R 4.4.0)
##  usethis       3.1.0.9000 2025-03-31 [1] Github (r-lib/usethis@a653d6e)
##  vctrs         0.6.5      2023-12-01 [1] CRAN (R 4.4.0)
##  vec         * 0.1.0      2025-01-07 [1] Github (jonocarroll/vec@595b07e)
##  xfun          0.51       2025-02-19 [1] CRAN (R 4.4.1)
##  xtable        1.8-4      2019-04-21 [1] CRAN (R 4.4.0)
##  yaml          2.3.10     2024-07-26 [1] CRAN (R 4.4.0)
## 
##  [1] /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library
## 
## ──────────────────────────────────────────────────────────────────────────────

These Languages are Accumulating

By: Jonathan Carroll

Re-posted from: https://jcarroll.com.au/2024/11/28/these-languages-are-accumulating/

I keep saying that the more programming languages you know, the more you will
understand all the others you know – I’m now at the point where I want to solve
every problem I see in a handful of different languages. They all offer
different functionality, and some are certainly more suited to particular
problems than others, but there’s a world of difference between two characters
and importing from two libraries.

A newsletter I follow (and can’t find online copies of) that demonstrates neat
things in python (gotta learn it, despite not loving it) recently covered
accumulate, showing that sum and accumulate were sort of related

>>> my_list = [42, 73, 0, 16, 10]
>>> sum(my_list)
141
>>> from itertools import accumulate
>>> list(accumulate(my_list))
[42, 115, 115, 131, 141]

sum adds up all the elements of the list, while accumulate does the same but
keeps each successive partial sum.

It rounds out the demo with an alternative function being used in accumulate

>>> from itertools import accumulate
>>> from operator import mul  # mul(2, 3) == 6
>>> initial_investment = 1000
>>> rates = [1.01, 1.01, 1.02, 1.025, 1.035, 1.035, 1.06]
>>> list(
...     accumulate(rates, mul, initial=initial_investment)
... )
[1000, 1010.0, 1020.1, 1040.502, 1066.515, 1103.843, 1142.477, 1211.026]

Now, firstly… from operator import mul??? It looks like there’s no way to
pass * as an argument to a function. I could define a function that performs
the same on known arguments, e.g. lambda x, y: x * y

>>> list(accumulate(rates, lambda x, y: x*y, initial=initial_investment))
[1000, 1010.0, 1020.1, 1040.502, 1066.5145499999999, 1103.8425592499998, 1142.4770488237498, 1211.0256717531747]

but… ew.

It’s possible that there’s a different way to approach this. A list
comprehension comes to mind, e.g. something like

>>> [sum(my_list[0:i]) for i in range(1, len(my_list)+1)]
[42, 115, 115, 131, 141]

but that requires performing a sum for each sub-interval, so performance would
not scale well (admittedly, that was not a consideration here at all). I also
don’t believe there’s a built-in prod so one must import math in order to do
similar

>>> import math
>>> x = [initial_investment] + rates
>>> [math.prod(x[0:i]) for i in range(1, len(x)+1)]
[1000, 1010.0, 1020.1, 1040.502, 1066.5145499999999, 1103.8425592499998, 1142.4770488237498, 1211.0256717531747]

In R that could use the built-in cumprod for the cumulative product

initial_investment <- 1000
rates = c(1.01, 1.01, 1.02, 1.025, 1.035, 1.035, 1.06)

cumprod(c(initial_investment, rates))
## [1] 1000.000 1010.000 1020.100 1040.502 1066.515 1103.843 1142.477 1211.026

but that has the ‘multiply’ operation hardcoded. cumsum uses + as the
function… hmmm. Maybe R doesn’t have a generalised accumulate?

I’ve been playing around with Haskell lately, so recursive functions to the
rescue! One feature of recursive functions in R that I really like is Recall
which calls the function in which it is defined with a new set of arguments –
perfect for recursion!

accumulate_recall <- function(x, f, i=x[1]) {
  if (!length(x)) return(NULL)
  c(i, Recall(tail(x, -1), f, f(i, x[2])))
}

It’s also robust against renaming the function; the body doesn’t actually call
accumulate_recall by name at all.

This might be inefficient, though – it’s not uncommon to blow out the stack, so
a new Tailcall function (which doesn’t have the same elegance of being robust
against renaming) helps with flagging this as something that can be optimised

accumulate <- function(x, f, i=x[1]) {
  if (!length(x)) return(NULL)
  c(i, Tailcall(accumulate, tail(x, -1), f, f(i, x[2])))
}

With this, I can emulate the cumsum() and cumprod() functions

cumprod(1:6)
## [1]   1   2   6  24 120 720
accumulate(1:6, `*`)
## [1]   1   2   6  24 120 720
cumsum(2:6)
## [1]  2  5  9 14 20
accumulate(2:6, `+`)
## [1]  2  5  9 14 20

unless I try to calculate something too big…

cumprod(5:15)
##  [1]           5          30         210        1680       15120      151200
##  [7]     1663200    19958400   259459200  3632428800 54486432000
accumulate(5:15, `*`)
## Warning in f(i, x[2]): NAs produced by integer overflow
##  [1]         5        30       210      1680     15120    151200   1663200
##  [8]  19958400 259459200        NA        NA

It appears that the built-in functions convert to numeric. That’s easily fixed
on input

accumulate(as.numeric(5:15), `*`)
##  [1]           5          30         210        1680       15120      151200
##  [7]     1663200    19958400   259459200  3632428800 54486432000

In any case, there’s a generalised accumulate that takes the bare functions as
arguments.

But it can be so much cleaner than this!

In APL you won’t find any function named “sum” because it is just a reduction
(Reduce in R) with the function +

      sum←+/
      
      sum ⍳6 ⍝ sum the values 1:6
21

      sum 1↓⍳6 ⍝ sum the values 2:6
20

which in R is

sum(1:6)
## [1] 21
sum(2:6)
## [1] 20

Why would you write sum if you can just use +/? It’s fewer
characters to write out the implementation than the name!

For accumulate the terminology in APL is scan which uses a very similar
glyph because the operation itself is very similar; a reduce (/) is just the
last value of a scan (\) which keeps the progressive values. In both cases,
the operator (either slash) takes a binary function as the left argument and
produces a modified function – in these examples, effectively sum and prod
which is then applied to values on the right. The scan version does the same

      +\⍳6
1 3 6 10 15 21

      ×\⍳6
1 2 6 24 120 720
accumulate(1:6, `+`)
## [1]  1  3  6 10 15 21
accumulate(1:6, `*`)
## [1]   1   2   6  24 120 720

As for the rates example above, we concatenate the initial value with catenate
(,) just like the R example, but otherwise this works fine

      rates ← 1.01 1.01 1.02 1.025 1.035 1.035 1.06
      inv ← 1000
      
      ×/inv, rates
1211.025672

      ×\inv, rates
1000 1010 1020.1 1040.502 1066.51455 1103.842559 1142.477049 1211.025672

So all of that recursive R code made to generalise the cumulative application of
a function provided as an argument is boiled down to just the single glyph \.
Outstanding!

What’s more, there are lots of binary functions one would use this with, all
of which have spelled-out names in other languages

      +/ ⍝ sum (add)
      ×/ ⍝ prod (multiply)
      ∧/ ⍝ all (and)
      ∨/ ⍝ any (or)
      ⌈/ ⍝ maximum (max)
      ⌊/ ⍝ minimum (min)

In summary, it seems that looking across these languages, the available options
range from a single glyph for scan along with the bare binary operator, e.g.
×/; a cumprod() function which isn’t well-generalised but works out of the
box; and then there’s whatever mess this is (once you’ve installed these)

>>> from itertools import accumulate
>>> from operator import mul
>>> list(accumulate(rates, mul, initial=initial_investment))

Where did we go so wrong?

For what it’s worth, Julia has a reduce and an accumulate that behave very
nicely; generalised for the binary function as an argument

julia> reduce(+, 1:6)
21

julia> reduce(*, 1:6)
720

julia> accumulate(+, 1:6)
6-element Vector{Int64}:
  1
  3
  6
 10
 15
 21

julia> accumulate(*, 1:6)
6-element Vector{Int64}:
   1
   2
   6
  24
 120
 720

This is extremely close to the APL approach, but with longer worded names for
the reduce and scan operators. It also defines the more convenient sum,
prod, cumsum, and cumprod; no shortage of ways to do this in Julia!

In Haskell, foldl and scanl are the (left-associative) version of reduce
and accumulate, and passing an infix as an argument necessitates wrapping it
in parentheses

ghci> foldl (+) 0 [1..6]
21

ghci> scanl (+) 0 [1..6]
[0,1,3,6,10,15,21]

ghci> foldl (*) 1 [1..6]
720

ghci> scanl (*) 1 [1..6]
[1,1,2,6,24,120,720]

This requires an explicit starting value, unless one uses the specialised
versions which use the first value as an initial value

ghci> foldl1 (+) [1..6]
21

ghci> scanl1 (+) [1..6]
[1,3,6,10,15,21]

ghci> foldl1 (*) [1..6]
720

ghci> scanl1 (*) [1..6]
[1,2,6,24,120,720]

I started this post hoping to demonstrate how nice the APL syntax was for this,
but the detour through generalising the R function was a lot of unexpected fun
as well.

Comments, improvements, or your own solutions are most welcome. I can be found
on Mastodon or use the comments below.

Addendums

It should probably be noted that R does have a function scan but it’s for
reading data into a vector – if you ever spot someone using it for that… run.
I have war stories about that function.

I’d love to hear how this is accomplished in some other languages, too – does it
have a built-in accumulate that takes a binary function?

devtools::session_info()
## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value
##  version  R version 4.4.1 (2024-06-14)
##  os       macOS Sonoma 14.6
##  system   aarch64, darwin20
##  ui       X11
##  language (EN)
##  collate  en_US.UTF-8
##  ctype    en_US.UTF-8
##  tz       Australia/Adelaide
##  date     2024-11-28
##  pandoc   3.5 @ /opt/homebrew/bin/ (via rmarkdown)
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package     * version    date (UTC) lib source
##  blogdown      1.19       2024-02-01 [1] CRAN (R 4.4.0)
##  bookdown      0.41       2024-10-16 [1] CRAN (R 4.4.1)
##  bslib         0.8.0      2024-07-29 [1] CRAN (R 4.4.0)
##  cachem        1.1.0      2024-05-16 [1] CRAN (R 4.4.0)
##  cli           3.6.3      2024-06-21 [1] CRAN (R 4.4.0)
##  devtools      2.4.5      2022-10-11 [1] CRAN (R 4.4.0)
##  digest        0.6.37     2024-08-19 [1] CRAN (R 4.4.1)
##  ellipsis      0.3.2      2021-04-29 [1] CRAN (R 4.4.0)
##  evaluate      1.0.1      2024-10-10 [1] CRAN (R 4.4.1)
##  fastmap       1.2.0      2024-05-15 [1] CRAN (R 4.4.0)
##  fs            1.6.5      2024-10-30 [1] CRAN (R 4.4.1)
##  glue          1.8.0      2024-09-30 [1] CRAN (R 4.4.1)
##  htmltools     0.5.8.1    2024-04-04 [1] CRAN (R 4.4.0)
##  htmlwidgets   1.6.4      2023-12-06 [1] CRAN (R 4.4.0)
##  httpuv        1.6.15     2024-03-26 [1] CRAN (R 4.4.0)
##  jquerylib     0.1.4      2021-04-26 [1] CRAN (R 4.4.0)
##  jsonlite      1.8.9      2024-09-20 [1] CRAN (R 4.4.1)
##  knitr         1.48       2024-07-07 [1] CRAN (R 4.4.0)
##  later         1.3.2      2023-12-06 [1] CRAN (R 4.4.0)
##  lifecycle     1.0.4      2023-11-07 [1] CRAN (R 4.4.0)
##  magrittr      2.0.3      2022-03-30 [1] CRAN (R 4.4.0)
##  memoise       2.0.1      2021-11-26 [1] CRAN (R 4.4.0)
##  mime          0.12       2021-09-28 [1] CRAN (R 4.4.0)
##  miniUI        0.1.1.1    2018-05-18 [1] CRAN (R 4.4.0)
##  pkgbuild      1.4.5      2024-10-28 [1] CRAN (R 4.4.1)
##  pkgload       1.4.0      2024-06-28 [1] CRAN (R 4.4.0)
##  profvis       0.4.0      2024-09-20 [1] CRAN (R 4.4.1)
##  promises      1.3.0      2024-04-05 [1] CRAN (R 4.4.0)
##  purrr         1.0.2      2023-08-10 [1] CRAN (R 4.4.0)
##  R6            2.5.1      2021-08-19 [1] CRAN (R 4.4.0)
##  Rcpp          1.0.13-1   2024-11-02 [1] CRAN (R 4.4.1)
##  remotes       2.5.0.9000 2024-11-03 [1] Github (r-lib/remotes@5b7eb08)
##  rlang         1.1.4      2024-06-04 [1] CRAN (R 4.4.0)
##  rmarkdown     2.28       2024-08-17 [1] CRAN (R 4.4.0)
##  rstudioapi    0.17.1     2024-10-22 [1] CRAN (R 4.4.1)
##  sass          0.4.9      2024-03-15 [1] CRAN (R 4.4.0)
##  sessioninfo   1.2.2      2021-12-06 [1] CRAN (R 4.4.0)
##  shiny         1.9.1      2024-08-01 [1] CRAN (R 4.4.0)
##  urlchecker    1.0.1      2021-11-30 [1] CRAN (R 4.4.0)
##  usethis       3.0.0      2024-07-29 [1] CRAN (R 4.4.0)
##  vctrs         0.6.5      2023-12-01 [1] CRAN (R 4.4.0)
##  xfun          0.49       2024-10-31 [1] CRAN (R 4.4.1)
##  xtable        1.8-4      2019-04-21 [1] CRAN (R 4.4.0)
##  yaml          2.3.10     2024-07-26 [1] CRAN (R 4.4.0)
## 
##  [1] /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library
## 
## ──────────────────────────────────────────────────────────────────────────────