csd <- function(x) {
  xsq <- x^2                             # Calculate the square of each element in x
  n <- length(x)                         # Determine the length of x
  cumsumx <- cumsum(x)                   # Calculate the cumulative sum of x
  cumxbar <- cumsumx / 1:length(x)       # Calculate the cumulative mean of x
  p1 <- cumsum(xsq)                      # Calculate the cumulative sum of x squared
  p2 <- -2 * cumxbar * cumsumx           # Calculate the second part of the variance
  p3 <- 1:n * (cumxbar^2)                # Calculate the third part of the variance
  cumvar <- (p1 + p2 + p3) / (1:n - 1)   # Calculate the cumulative variance
  cumsd <- sqrt(cumvar)                  # Calculate the cumulative standard deviation
  
  return(cumsd)                          # Return the cumulative standard deviation
}

> microbenchmark::microbenchmark(
+   csd(x),
+   csd2(x),
+   times = 100L
+ )
Unit: microseconds
    expr     min      lq      mean   median       uq     max neval cld
  csd(x) 18470.6 19005.5 20086.684 19541.55 20052.30 31079.3   100  a 
 csd2(x)    18.0    20.8    28.965    29.00    32.45    97.5   100   b
> 
> rbenchmark::benchmark(
+   current_func = csd(x),
+   optimized_func = csd2(x),
+   replications = 100L,
+   columns = c("test","replications","elapsed",
+               "relative","user.self","sys.self")
+ )
            test replications elapsed relative user.self sys.self
1   current_func          100    2.04       NA      1.96     0.06
2 optimized_func          100    0.00       NA      0.00     0.00
> 19541/29
[1] 673.8276
> 20086/28.965
[1] 693.4576