r — округление удваивается до n десятичных знаков, используя «округление до четного»; в переполнении стека

Я хочу реализовать функцию в C ++:

double round_to_even(double num, int decimal_places);

/* rounds the first argument to three decimal places
using round-to-even (unbiased rounding) */

который ведет себя так, как работает округление в R, используя метод, известный как округление до четного, непредвзятое округление или округление статистики:
https://stat.ethz.ch/R-manual/R-devel/library/base/html/Round.html
Лучший способ проиллюстрировать это на нескольких примерах из R, где decimal_place = 3:

> round(0.1247,3) # rounding the first argument to three decimal places
[1] 0.125

> round(0.1244,3)
[1] 0.124

Обратите внимание на то, как он ведет себя так, как и следовало ожидать, округляя в большую сторону, когда цифра справа от указанного десятичного_значения больше 5 (здесь равно 7), и обрезая, когда она меньше 5 (здесь равна 4). Однако, когда цифра справа от указанного decimal_place равна 5 (или в «средней точке»), она ведет себя следующим образом:

> round(0.1255,3) # will round up
[1] 0.126

> round(.1275,3) # will round up
[1] 0.128

> round(0.1265,3) # will truncate
[1] 0.126

> round(0.1245,3) # will truncate
[1] 0.124

путем округления, когда цифра слева от decimal_place нечетная, и усечения, когда четная.

Есть ли краткий способ реализовать это, и в частности без преобразования цифр в символы?

Редактировать:
Это лучшее, что я мог придумать:

double round_to_even(double number, int decimal_points)
/* rounds the first argument to three decimal places
using round-to-even (unbiased rounding) */
{
double num_left = number, num_right = number;
int digit_left, digit_right;

num_left *= pow(10, decimal_points + 1);
digit_left = fmod(num_left, 10);

if (digit_left == 5)
{
num_right *= pow(10, decimal_points);
digit_right = fmod(num_right, 10);

if (digit_right % 2 == 0) // if even
return floor(number * pow(10, decimal_points)) / pow(10, decimal_points);
else // otherwise it's odd
return ceil(number * pow(10, decimal_points)) / pow(10, decimal_points);
}
else { // standard round-to-nearest
return round(number * pow(10, decimal_points)) / pow(10, decimal_points);
}
}

И я проверил это:

std::vector<double> test_vector({ 0.1247, 0.1244, 0.1255, 0.1275, 0.1265, 0.1245 });
std::vector<double> expected_values({ 0.125, 0.124, 0.126, 0.128, 0.126, 0.124 });

for (std::vector<double>::size_type i = 0; i < test_vector.size(); i++)
std::cout << "expected: " << expected_values[i] << "\t got: " << round_to_even(test_vector[i], 3) << std::endl;

Что дает следующий вывод:

expected: 0.125  got: 0.125
expected: 0.124  got: 0.124
expected: 0.126  got: 0.126
expected: 0.128  got: 0.128
expected: 0.126  got: 0.126
expected: 0.124  got: 0.124