From a4d7776da747c53efc415bcfdaa5101eb088fb63 Mon Sep 17 00:00:00 2001 From: Peter Ullrich Date: Thu, 30 Apr 2026 21:39:23 +0200 Subject: [PATCH 1/7] Replace Float.round/2 algorithm with Cox 2026 uscale algorithm --- lib/elixir/lib/float.ex | 197 ++++++++++++++++++++-------------------- 1 file changed, 98 insertions(+), 99 deletions(-) diff --git a/lib/elixir/lib/float.ex b/lib/elixir/lib/float.ex index bb7b194575e..aa6da2688f6 100644 --- a/lib/elixir/lib/float.ex +++ b/lib/elixir/lib/float.ex @@ -344,9 +344,6 @@ defmodule Float do """ @spec round(float, precision_range) :: float - # This implementation is slow since it relies on big integers. - # Faster implementations are available on more recent papers - # and could be implemented in the future. def round(float, precision \\ 0) def round(float, 0) when float == 0.0, do: float @@ -366,135 +363,137 @@ defmodule Float do raise ArgumentError, invalid_precision_message(precision) end + # Float rounding via Russ Cox's "unrounded scaling" algorithm specialised for + # fixed-precision rounding. Reference: https://research.swtch.com/fp + # + # 1. Decompose `f = (-1)^s * m * 2^e` (m has 53 bits incl. implicit leading 1). + # 2. Bounded multiply `prod = m * 10^p` (≤103 bits, exact for p in 0..15). + # 3. Round `prod / 2^-e` to integer `n` per the requested rounding mode. + # 4. Emit float closest to `n / 10^p`: + # - fast path: when n < 2^53, IEEE float division is correctly rounded. + # - slow path: bignum scaling + manual mantissa extraction. defp round(num, _precision, _rounding) when is_float(num) and num == 0.0, do: num - defp round(float, precision, rounding) do - <> = <> - {num, count} = decompose(significant, 1) - count = count - exp + 1023 + defp round(float, precision, mode) do + <> = <> cond do - # Precision beyond 15 digits - count >= 104 -> - case rounding do - :ceil when sign === 0 -> 1 / power_of_10(precision) - :floor when sign === 1 -> -1 / power_of_10(precision) - :ceil when sign === 1 -> -0.0 - :half_up when sign === 1 -> -0.0 - _ -> 0.0 - end + # Subnormal — tiny but non-zero; treat per-mode (ceil(+) and floor(-) bump + # to 10^-precision; everything else rounds to signed zero). + exp_field == 0 -> + tiny_round(sign, precision, mode) - # We are asking more precision than we have - count <= precision -> + # |f| >= 2^52 — already integer-valued at any precision >= 1. + exp_field - 1075 >= 0 -> float true -> - # Difference in precision between float and asked precision - # We subtract 1 because we need to calculate the remainder too - diff = count - precision - 1 - - # Get up to latest so we calculate the remainder - power_of_10 = power_of_10(diff) - - # Convert the numerand to decimal base - num = num * power_of_5(count) - - # Move to the given precision - 1 - num = div(num, power_of_10) - div = div(num, 10) - num = rounding(rounding, sign, num, div) - - # Convert back to float without loss - # https://www.exploringbinary.com/correct-decimal-to-floating-point-using-big-integers/ - den = power_of_10(precision) - boundary = den <<< 52 - - cond do - num == 0 and sign == 1 -> - -0.0 - - num == 0 -> - 0.0 - - num >= boundary -> - {den, exp} = scale_down(num, boundary, 52) - decimal_to_float(sign, num, den, exp) - - true -> - {num, exp} = scale_up(num, boundary, 52) - decimal_to_float(sign, num, den, exp) - end + do_round(sign, @power_of_2_to_52 ||| mant, 1075 - exp_field, precision, mode) end end - defp decompose(significant, initial) do - decompose(significant, 1, 0, initial) - end - - defp decompose(<<1::1, bits::bitstring>>, count, last_count, acc) do - decompose(bits, count + 1, count, (acc <<< (count - last_count)) + 1) - end - - defp decompose(<<0::1, bits::bitstring>>, count, last_count, acc) do - decompose(bits, count + 1, last_count, acc) + # |f * 10^p| < 0.5 — integer round is 0; ceil/floor still bump per sign. + defp do_round(sign, _m, shift, precision, mode) when shift >= 104 do + tiny_round(sign, precision, mode) end - defp decompose(<<>>, _count, last_count, acc) do - {acc, last_count} - end + defp do_round(sign, m, shift, precision, mode) do + pow = power_of_10(precision) + prod = m * pow + half = 1 <<< (shift - 1) + q = prod >>> shift + rem = prod - (q <<< shift) + n = round_step(mode, sign, q, rem, half) - defp scale_up(num, boundary, exp) when num >= boundary, do: {num, exp} - defp scale_up(num, boundary, exp), do: scale_up(num <<< 1, boundary, exp - 1) + cond do + n == 0 -> + signed_zero(sign) - defp scale_down(num, den, exp) do - new_den = den <<< 1 + n < @power_of_2_to_52 <<< 1 -> + # Both n and pow fit in 53 bits, so IEEE float division is correctly rounded. + r = :erlang.float(n) / :erlang.float(pow) + if sign == 1, do: -r, else: r - if num < new_den do - {den >>> 52, exp} - else - scale_down(num, new_den, exp + 1) + true -> + bignum_to_float(sign, n, pow) end end - defp decimal_to_float(sign, num, den, exp) do + defp round_step(:half_up, _sign, q, rem, half) do + if rem >= half, do: q + 1, else: q + end + + defp round_step(:floor, 0, q, _rem, _half), do: q + defp round_step(:floor, 1, q, rem, _half) when rem > 0, do: q + 1 + defp round_step(:floor, 1, q, _rem, _half), do: q + defp round_step(:ceil, 0, q, rem, _half) when rem > 0, do: q + 1 + defp round_step(:ceil, 0, q, _rem, _half), do: q + defp round_step(:ceil, 1, q, _rem, _half), do: q + + defp signed_zero(0), do: 0.0 + defp signed_zero(1), do: -0.0 + + # Result of rounding a non-zero float whose |f * 10^precision| < 0.5. + # ceil(+) → +10^-precision, floor(-) → -10^-precision, others → signed 0. + defp tiny_round(0, precision, :ceil), do: 1.0 / :erlang.float(power_of_10(precision)) + defp tiny_round(1, precision, :floor), do: -1.0 / :erlang.float(power_of_10(precision)) + defp tiny_round(sign, _precision, _mode), do: signed_zero(sign) + + # Slow path: emit float closest to (sign * n / pow) when n >= 2^53. + # The binary emission step is always IEEE round-to-nearest-even, regardless + # of the integer-rounding mode used for `n`. + defp bignum_to_float(sign, n, pow) do + s = bit_length(n) - bit_length(pow) - 53 + {num, den, exp} = align(n, pow, s) + quo = div(num, den) rem = num - quo * den - - tmp = - case den >>> 1 do - den when rem > den -> quo + 1 - den when rem < den -> quo - _ when (quo &&& 1) === 1 -> quo + 1 - _ -> quo + half = den >>> 1 + + mant = + cond do + rem > half -> quo + 1 + rem < half -> quo + (quo &&& 1) === 1 -> quo + 1 + true -> quo end - tmp = tmp - @power_of_2_to_52 - <> = <> - tmp + <> = <> + f end - defp rounding(:floor, 1, _num, div), do: div + 1 - defp rounding(:ceil, 0, _num, div), do: div + 1 + # Pick (num, den, exp) so that num/den ∈ [2^52, 2^53) and the resulting + # float = num/den * 2^(exp-52). + defp align(n, pow, s) when s >= 0 do + new_pow = pow <<< s - defp rounding(:half_up, _sign, num, div) do - case rem(num, 10) do - rem when rem < 5 -> div - rem when rem >= 5 -> div + 1 - end + if n < new_pow <<< 53, + do: {n, new_pow, 52 + s}, + else: {n, new_pow <<< 1, 53 + s} end - defp rounding(_, _, _, div), do: div + defp align(n, pow, s) do + new_n = n <<< -s + boundary = pow <<< 52 + + if new_n >= boundary, + do: {new_n, pow, 52 + s}, + else: {new_n <<< 1, pow, 51 + s} + end - Enum.reduce(0..104, 1, fn x, acc -> + defp bit_length(0), do: 0 + defp bit_length(n) when n > 0, do: bit_length(n, 0) + defp bit_length(n, acc) when n >= 1 <<< 64, do: bit_length(n >>> 64, acc + 64) + defp bit_length(n, acc) when n >= 1 <<< 16, do: bit_length(n >>> 16, acc + 16) + defp bit_length(n, acc) when n >= 1 <<< 4, do: bit_length(n >>> 4, acc + 4) + defp bit_length(n, acc) when n >= 1, do: bit_length(n >>> 1, acc + 1) + defp bit_length(_, acc), do: acc + + Enum.reduce(0..15, 1, fn x, acc -> defp power_of_10(unquote(x)), do: unquote(acc) acc * 10 end) - Enum.reduce(0..104, 1, fn x, acc -> - defp power_of_5(unquote(x)), do: unquote(acc) - acc * 5 - end) - @doc """ Returns a pair of integers whose ratio is exactly equal to the original float and with a positive denominator. From 400eb9c37a9bfb2feeb1bbc77a3094ed19eeefc4 Mon Sep 17 00:00:00 2001 From: Peter Ullrich Date: Thu, 30 Apr 2026 22:12:09 +0200 Subject: [PATCH 2/7] Add additional edge cases to tests --- lib/elixir/test/elixir/float_test.exs | 84 +++++++++++++++++++++++++++ 1 file changed, 84 insertions(+) diff --git a/lib/elixir/test/elixir/float_test.exs b/lib/elixir/test/elixir/float_test.exs index 751a482d1cb..df411741949 100644 --- a/lib/elixir/test/elixir/float_test.exs +++ b/lib/elixir/test/elixir/float_test.exs @@ -103,6 +103,36 @@ defmodule FloatTest do assert Float.floor(5.0e-324, precision) === 0.0 end end + + test "with already-exact floats does not bump" do + # The integer rounding step must not add 1 when the truncated remainder + # is exactly zero (e.g. -1.5 has no content below the 1st decimal place). + assert Float.floor(-1.5, 1) === -1.5 + assert Float.floor(-1.875, 3) === -1.875 + assert Float.floor(-3.0, 5) === -3.0 + end + + test "with extremely small floats" do + # `tiny_round`: ceil(+) and floor(-) must bump to ±10^-precision. + assert Float.floor(-1.0e-200, 5) === -1.0e-5 + assert Float.floor(-1.0e-200, 15) === -1.0e-15 + assert Float.floor(1.0e-200, 5) === 0.0 + end + + test "with already-integer floats" do + # |f| >= 2^52 — fast path returns the float unchanged. + assert Float.floor(:math.pow(2, 52), 5) === 4_503_599_627_370_496.0 + assert Float.floor(1.0e20, 3) === 1.0e20 + assert Float.floor(-1.0e20, 3) === -1.0e20 + end + + test "with very large floats hits the bignum slow path" do + # n = round(|f| * 10^p) >= 2^53 — exercises bignum_to_float / align. + # The slow path must round-trip representable floats back to themselves. + assert Float.floor(1.234e15, 3) === 1.234e15 + assert Float.floor(1.234567e11, 5) === 1.234567e11 + assert Float.floor(-1.234567e11, 5) === -1.234567e11 + end end test "ceil/1" do @@ -163,6 +193,30 @@ defmodule FloatTest do assert Float.ceil(-5.0e-324, precision) === -0.0 end end + + test "with already-exact floats does not bump" do + assert Float.ceil(1.5, 1) === 1.5 + assert Float.ceil(1.875, 3) === 1.875 + assert Float.ceil(3.0, 5) === 3.0 + end + + test "with extremely small floats" do + assert Float.ceil(1.0e-200, 5) === 1.0e-5 + assert Float.ceil(1.0e-200, 15) === 1.0e-15 + assert Float.ceil(-1.0e-200, 5) === -0.0 + end + + test "with already-integer floats" do + assert Float.ceil(:math.pow(2, 52), 5) === 4_503_599_627_370_496.0 + assert Float.ceil(1.0e20, 3) === 1.0e20 + assert Float.ceil(-1.0e20, 3) === -1.0e20 + end + + test "with very large floats hits the bignum slow path" do + assert Float.ceil(1.234e15, 3) === 1.234e15 + assert Float.ceil(1.234567e11, 5) === 1.234567e11 + assert Float.ceil(-1.234567e11, 5) === -1.234567e11 + end end describe "round/2" do @@ -203,6 +257,36 @@ defmodule FloatTest do assert Float.round(-5.0e-324, precision) === -0.0 end end + + test "rounds up across a digit boundary" do + # Rounding pushes the integer answer to 10^precision; the float-emission + # step must produce the next-magnitude value cleanly. + assert Float.round(0.9995, 3) === 1.0 + assert Float.round(0.9999, 3) === 1.0 + assert Float.round(99.9999, 3) === 100.0 + assert Float.round(-99.9999, 3) === -100.0 + end + + test "with already-integer floats" do + assert Float.round(:math.pow(2, 52), 5) === 4_503_599_627_370_496.0 + assert Float.round(1.0e20, 3) === 1.0e20 + assert Float.round(-1.0e20, 3) === -1.0e20 + end + + test "with very large floats hits the bignum slow path" do + assert Float.round(1.234e15, 3) === 1.234e15 + assert Float.round(1.234567e11, 5) === 1.234567e11 + assert Float.round(-1.234567e11, 5) === -1.234567e11 + end + + test "preserves documented tie behavior" do + # The actual binary representation of 5.5675 is 5.567499999..., so + # half-up rounding correctly gives 5.567 (not 5.568). + assert Float.round(5.5675, 3) === 5.567 + assert Float.round(-5.5675, 3) === -5.567 + assert Float.round(12.5, 0) === 13.0 + assert Float.round(-12.5, 0) === -13.0 + end end describe "ratio/1" do From f8c5ca79fdef76e4a1f551e703947696bfbe92ad Mon Sep 17 00:00:00 2001 From: Peter Ullrich Date: Thu, 30 Apr 2026 22:17:39 +0200 Subject: [PATCH 3/7] Make code more readable --- lib/elixir/lib/float.ex | 147 +++++++++++++++++++++------------------- 1 file changed, 76 insertions(+), 71 deletions(-) diff --git a/lib/elixir/lib/float.ex b/lib/elixir/lib/float.ex index aa6da2688f6..77b7d15a108 100644 --- a/lib/elixir/lib/float.ex +++ b/lib/elixir/lib/float.ex @@ -366,131 +366,136 @@ defmodule Float do # Float rounding via Russ Cox's "unrounded scaling" algorithm specialised for # fixed-precision rounding. Reference: https://research.swtch.com/fp # - # 1. Decompose `f = (-1)^s * m * 2^e` (m has 53 bits incl. implicit leading 1). - # 2. Bounded multiply `prod = m * 10^p` (≤103 bits, exact for p in 0..15). - # 3. Round `prod / 2^-e` to integer `n` per the requested rounding mode. - # 4. Emit float closest to `n / 10^p`: - # - fast path: when n < 2^53, IEEE float division is correctly rounded. + # 1. Decompose float = (-1)^sign * mantissa * 2^binary_exp + # (mantissa has 53 bits incl. implicit leading 1). + # 2. Bounded multiply `product = mantissa * 10^precision` (≤103 bits, exact). + # 3. Round `product / 2^-binary_exp` to an integer per the requested mode. + # 4. Emit float closest to `rounded_int / 10^precision`: + # - fast path: when rounded_int < 2^53, IEEE float division is correctly rounded. # - slow path: bignum scaling + manual mantissa extraction. defp round(num, _precision, _rounding) when is_float(num) and num == 0.0, do: num defp round(float, precision, mode) do - <> = <> + <> = <> cond do # Subnormal — tiny but non-zero; treat per-mode (ceil(+) and floor(-) bump # to 10^-precision; everything else rounds to signed zero). - exp_field == 0 -> + exp == 0 -> tiny_round(sign, precision, mode) - # |f| >= 2^52 — already integer-valued at any precision >= 1. - exp_field - 1075 >= 0 -> + # |float| >= 2^52 — already integer-valued at any precision >= 1. + exp - 1075 >= 0 -> float true -> - do_round(sign, @power_of_2_to_52 ||| mant, 1075 - exp_field, precision, mode) + mantissa = @power_of_2_to_52 ||| mantissa + shift = 1075 - exp + do_round(sign, mantissa, shift, precision, mode) end end - # |f * 10^p| < 0.5 — integer round is 0; ceil/floor still bump per sign. - defp do_round(sign, _m, shift, precision, mode) when shift >= 104 do + # |float * 10^precision| < 0.5 — integer round is 0; ceil/floor still bump per sign. + defp do_round(sign, _mantissa, shift, precision, mode) when shift >= 104 do tiny_round(sign, precision, mode) end - defp do_round(sign, m, shift, precision, mode) do - pow = power_of_10(precision) - prod = m * pow + defp do_round(sign, mantissa, shift, precision, mode) do + power = power_of_10(precision) + product = mantissa * power half = 1 <<< (shift - 1) - q = prod >>> shift - rem = prod - (q <<< shift) - n = round_step(mode, sign, q, rem, half) + quotient = product >>> shift + remainder = product - (quotient <<< shift) + rounded_int = round_step(mode, sign, quotient, remainder, half) cond do - n == 0 -> + rounded_int == 0 -> signed_zero(sign) - n < @power_of_2_to_52 <<< 1 -> - # Both n and pow fit in 53 bits, so IEEE float division is correctly rounded. - r = :erlang.float(n) / :erlang.float(pow) - if sign == 1, do: -r, else: r + rounded_int < @power_of_2_to_52 <<< 1 -> + # Both rounded_int and power fit in 53 bits, so IEEE float division + # is correctly rounded. + result = :erlang.float(rounded_int) / :erlang.float(power) + if sign == 1, do: -result, else: result true -> - bignum_to_float(sign, n, pow) + bignum_to_float(sign, rounded_int, power) end end - defp round_step(:half_up, _sign, q, rem, half) do - if rem >= half, do: q + 1, else: q + defp round_step(:half_up, _sign, quotient, remainder, half) do + if remainder >= half, do: quotient + 1, else: quotient end - defp round_step(:floor, 0, q, _rem, _half), do: q - defp round_step(:floor, 1, q, rem, _half) when rem > 0, do: q + 1 - defp round_step(:floor, 1, q, _rem, _half), do: q - defp round_step(:ceil, 0, q, rem, _half) when rem > 0, do: q + 1 - defp round_step(:ceil, 0, q, _rem, _half), do: q - defp round_step(:ceil, 1, q, _rem, _half), do: q + defp round_step(:floor, 0, quotient, _remainder, _half), do: quotient + defp round_step(:floor, 1, quotient, remainder, _half) when remainder > 0, do: quotient + 1 + defp round_step(:floor, 1, quotient, _remainder, _half), do: quotient + + defp round_step(:ceil, 0, quotient, remainder, _half) when remainder > 0, do: quotient + 1 + defp round_step(:ceil, 0, quotient, _remainder, _half), do: quotient + defp round_step(:ceil, 1, quotient, _remainder, _half), do: quotient defp signed_zero(0), do: 0.0 defp signed_zero(1), do: -0.0 - # Result of rounding a non-zero float whose |f * 10^precision| < 0.5. + # Result of rounding a non-zero float whose |float * 10^precision| < 0.5. # ceil(+) → +10^-precision, floor(-) → -10^-precision, others → signed 0. defp tiny_round(0, precision, :ceil), do: 1.0 / :erlang.float(power_of_10(precision)) defp tiny_round(1, precision, :floor), do: -1.0 / :erlang.float(power_of_10(precision)) defp tiny_round(sign, _precision, _mode), do: signed_zero(sign) - # Slow path: emit float closest to (sign * n / pow) when n >= 2^53. - # The binary emission step is always IEEE round-to-nearest-even, regardless - # of the integer-rounding mode used for `n`. - defp bignum_to_float(sign, n, pow) do - s = bit_length(n) - bit_length(pow) - 53 - {num, den, exp} = align(n, pow, s) + # Slow path: emit float closest to `sign * rounded_int / power` when + # rounded_int >= 2^53. The binary emission step is always IEEE + # round-to-nearest-even, regardless of the integer-rounding mode. + defp bignum_to_float(sign, rounded_int, power) do + shift_adjust = bit_length(rounded_int) - bit_length(power) - 53 + {numerator, denominator, exp} = align(rounded_int, power, shift_adjust) - quo = div(num, den) - rem = num - quo * den - half = den >>> 1 + quotient = div(numerator, denominator) + remainder = numerator - quotient * denominator + half = denominator >>> 1 - mant = + mantissa = cond do - rem > half -> quo + 1 - rem < half -> quo - (quo &&& 1) === 1 -> quo + 1 - true -> quo + remainder > half -> quotient + 1 + remainder < half -> quotient + (quotient &&& 1) === 1 -> quotient + 1 + true -> quotient end - <> = <> - f + <> = <> + result end - # Pick (num, den, exp) so that num/den ∈ [2^52, 2^53) and the resulting - # float = num/den * 2^(exp-52). - defp align(n, pow, s) when s >= 0 do - new_pow = pow <<< s + # Pick (numerator, denominator, exp) so that numerator/denominator ∈ [2^52, 2^53) + # and the resulting float = numerator/denominator * 2^(exp-52). + defp align(rounded_int, power, shift_adjust) when shift_adjust >= 0 do + new_power = power <<< shift_adjust - if n < new_pow <<< 53, - do: {n, new_pow, 52 + s}, - else: {n, new_pow <<< 1, 53 + s} + if rounded_int < new_power <<< 53, + do: {rounded_int, new_power, 52 + shift_adjust}, + else: {rounded_int, new_power <<< 1, 53 + shift_adjust} end - defp align(n, pow, s) do - new_n = n <<< -s - boundary = pow <<< 52 + defp align(rounded_int, power, shift_adjust) do + shifted = rounded_int <<< -shift_adjust + boundary = power <<< 52 - if new_n >= boundary, - do: {new_n, pow, 52 + s}, - else: {new_n <<< 1, pow, 51 + s} + if shifted >= boundary, + do: {shifted, power, 52 + shift_adjust}, + else: {shifted <<< 1, power, 51 + shift_adjust} end defp bit_length(0), do: 0 - defp bit_length(n) when n > 0, do: bit_length(n, 0) - defp bit_length(n, acc) when n >= 1 <<< 64, do: bit_length(n >>> 64, acc + 64) - defp bit_length(n, acc) when n >= 1 <<< 16, do: bit_length(n >>> 16, acc + 16) - defp bit_length(n, acc) when n >= 1 <<< 4, do: bit_length(n >>> 4, acc + 4) - defp bit_length(n, acc) when n >= 1, do: bit_length(n >>> 1, acc + 1) - defp bit_length(_, acc), do: acc - - Enum.reduce(0..15, 1, fn x, acc -> - defp power_of_10(unquote(x)), do: unquote(acc) + defp bit_length(integer) when integer > 0, do: bit_length(integer, 0) + defp bit_length(integer, acc) when integer >= 1 <<< 64, do: bit_length(integer >>> 64, acc + 64) + defp bit_length(integer, acc) when integer >= 1 <<< 16, do: bit_length(integer >>> 16, acc + 16) + defp bit_length(integer, acc) when integer >= 1 <<< 4, do: bit_length(integer >>> 4, acc + 4) + defp bit_length(integer, acc) when integer >= 1, do: bit_length(integer >>> 1, acc + 1) + defp bit_length(_integer, acc), do: acc + + Enum.reduce(0..15, 1, fn exponent, acc -> + defp power_of_10(unquote(exponent)), do: unquote(acc) acc * 10 end) From 12a691398986edfba708b60fa7bba967501783a3 Mon Sep 17 00:00:00 2001 From: Peter Ullrich Date: Thu, 30 Apr 2026 22:35:22 +0200 Subject: [PATCH 4/7] Remove deprecated `:erlang.float/1` --- lib/elixir/lib/float.ex | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/lib/elixir/lib/float.ex b/lib/elixir/lib/float.ex index 77b7d15a108..892fa2325d9 100644 --- a/lib/elixir/lib/float.ex +++ b/lib/elixir/lib/float.ex @@ -349,7 +349,9 @@ defmodule Float do def round(float, 0) when float == 0.0, do: float def round(float, 0) when is_float(float) do - case float |> :erlang.round() |> :erlang.float() do + rounded = :erlang.round(float) * 1.0 + + case rounded do zero when zero == 0.0 and float < 0.0 -> -0.0 rounded -> rounded end @@ -415,7 +417,7 @@ defmodule Float do rounded_int < @power_of_2_to_52 <<< 1 -> # Both rounded_int and power fit in 53 bits, so IEEE float division # is correctly rounded. - result = :erlang.float(rounded_int) / :erlang.float(power) + result = rounded_int / power if sign == 1, do: -result, else: result true -> @@ -440,8 +442,8 @@ defmodule Float do # Result of rounding a non-zero float whose |float * 10^precision| < 0.5. # ceil(+) → +10^-precision, floor(-) → -10^-precision, others → signed 0. - defp tiny_round(0, precision, :ceil), do: 1.0 / :erlang.float(power_of_10(precision)) - defp tiny_round(1, precision, :floor), do: -1.0 / :erlang.float(power_of_10(precision)) + defp tiny_round(0, precision, :ceil), do: 1.0 / power_of_10(precision) + defp tiny_round(1, precision, :floor), do: -1.0 / power_of_10(precision) defp tiny_round(sign, _precision, _mode), do: signed_zero(sign) # Slow path: emit float closest to `sign * rounded_int / power` when From e2b4b468006cd199ddd186d6e8cbcc499aca5ef4 Mon Sep 17 00:00:00 2001 From: Peter Ullrich Date: Fri, 1 May 2026 12:27:07 +0200 Subject: [PATCH 5/7] Update function comment to reflect that this is not an implementation of Cox, but only a refactor of the existing exact rational scaling approach. Fix bug when mantissa is exactly `2^53`. Fix Float.round/2 slow-path mantissa overflow when shift_adjust < 0 and the integer quotient lands in [2^53, 2^54) --- lib/elixir/lib/float.ex | 56 ++++++++++++++++++++++++++++++----------- 1 file changed, 42 insertions(+), 14 deletions(-) diff --git a/lib/elixir/lib/float.ex b/lib/elixir/lib/float.ex index 892fa2325d9..28f17c9166a 100644 --- a/lib/elixir/lib/float.ex +++ b/lib/elixir/lib/float.ex @@ -365,16 +365,34 @@ defmodule Float do raise ArgumentError, invalid_precision_message(precision) end - # Float rounding via Russ Cox's "unrounded scaling" algorithm specialised for - # fixed-precision rounding. Reference: https://research.swtch.com/fp + # Decimal-place rounding via exact rational scaling. This is the bignum + # core used by reference implementations like David M. Gay's "Correctly + # Rounded Binary-Decimal and Decimal-Binary Conversions" (cited in the + # @doc above), Python's round(), and Java's BigDecimal.setScale. # - # 1. Decompose float = (-1)^sign * mantissa * 2^binary_exp - # (mantissa has 53 bits incl. implicit leading 1). - # 2. Bounded multiply `product = mantissa * 10^precision` (≤103 bits, exact). - # 3. Round `product / 2^-binary_exp` to an integer per the requested mode. - # 4. Emit float closest to `rounded_int / 10^precision`: - # - fast path: when rounded_int < 2^53, IEEE float division is correctly rounded. - # - slow path: bignum scaling + manual mantissa extraction. + # 1. Decompose float exactly: |float| = mantissa / 2^shift. + # 2. Scale exactly: |float| * 10^precision = mantissa * 10^precision / 2^shift. + # Because precision is bounded to 0..15, the product fits in ~103 bits + # (53-bit mantissa + ~50-bit power of ten) and BEAM bignums handle it + # directly without approximation. + # 3. Round the exact rational to an integer per the requested mode + # (half_up / floor / ceil) using quotient and remainder. + # 4. Emit the float closest to rounded_int / 10^precision: + # - fast path: when rounded_int < 2^53, both operands are exactly + # representable as floats and IEEE division is correctly rounded. + # - slow path: bignum alignment + manual mantissa extraction with + # round-to-nearest-even for the trailing bit. + # + # The integer-rounding decision (step 3) and the binary-emission decision + # (step 4) are deliberately independent: step 3 picks the exact rational + # the user asked for, step 4 picks the closest float to that rational. + # Conflating them is the classic source of double-rounding bugs. + # + # Faster algorithms exist (Cox 2026's table-based uscale; Ryū / Schubfach + # for round-trip printing) but target different problems or assume + # fixed-width machine arithmetic that BEAM doesn't expose efficiently. + # At precision <= 15, the exact path is small, easy to audit, and fast + # enough that a more complex algorithm has not been justified by benchmarks. defp round(num, _precision, _rounding) when is_float(num) and num == 0.0, do: num defp round(float, precision, mode) do @@ -386,7 +404,7 @@ defmodule Float do exp == 0 -> tiny_round(sign, precision, mode) - # |float| >= 2^52 — already integer-valued at any precision >= 1. + # |float| >= 2^52 — has no fractional bits, return unchanged. exp - 1075 >= 0 -> float @@ -465,6 +483,15 @@ defmodule Float do true -> quotient end + # Carry-bit normalization: `mantissa` lives in [2^52, 2^53]. The upper + # bound `2^53` is reachable when `align/3` returns an upper-bound quotient + # or when rounding carries. Rebalance into the canonical [2^52, 2^53) + # range so the 52-bit packing below doesn't silently truncate. + {mantissa, exp} = + if mantissa == @power_of_2_to_52 <<< 1, + do: {@power_of_2_to_52, exp + 1}, + else: {mantissa, exp} + <> = <> result end @@ -481,11 +508,12 @@ defmodule Float do defp align(rounded_int, power, shift_adjust) do shifted = rounded_int <<< -shift_adjust - boundary = power <<< 52 - if shifted >= boundary, - do: {shifted, power, 52 + shift_adjust}, - else: {shifted <<< 1, power, 51 + shift_adjust} + cond do + shifted >= power <<< 53 -> {shifted, power <<< 1, 53 + shift_adjust} + shifted >= power <<< 52 -> {shifted, power, 52 + shift_adjust} + true -> {shifted <<< 1, power, 51 + shift_adjust} + end end defp bit_length(0), do: 0 From ec2a0ce2a5ecfb28ee867e7e24e71f2ef6ec69e4 Mon Sep 17 00:00:00 2001 From: Peter Ullrich Date: Fri, 1 May 2026 12:27:19 +0200 Subject: [PATCH 6/7] Update tests to cover edge cases better --- lib/elixir/test/elixir/float_test.exs | 29 ++++++++++++++++++++++++--- 1 file changed, 26 insertions(+), 3 deletions(-) diff --git a/lib/elixir/test/elixir/float_test.exs b/lib/elixir/test/elixir/float_test.exs index df411741949..607accbd4b4 100644 --- a/lib/elixir/test/elixir/float_test.exs +++ b/lib/elixir/test/elixir/float_test.exs @@ -121,7 +121,6 @@ defmodule FloatTest do test "with already-integer floats" do # |f| >= 2^52 — fast path returns the float unchanged. - assert Float.floor(:math.pow(2, 52), 5) === 4_503_599_627_370_496.0 assert Float.floor(1.0e20, 3) === 1.0e20 assert Float.floor(-1.0e20, 3) === -1.0e20 end @@ -129,6 +128,9 @@ defmodule FloatTest do test "with very large floats hits the bignum slow path" do # n = round(|f| * 10^p) >= 2^53 — exercises bignum_to_float / align. # The slow path must round-trip representable floats back to themselves. + assert Float.floor(2_661_101_816_343_531.5, 1) === 2_661_101_816_343_531.5 + assert Float.floor(3.0e15, 1) === 3.0e15 + assert Float.floor(-3.0e15, 1) === -3.0e15 assert Float.floor(1.234e15, 3) === 1.234e15 assert Float.floor(1.234567e11, 5) === 1.234567e11 assert Float.floor(-1.234567e11, 5) === -1.234567e11 @@ -207,12 +209,14 @@ defmodule FloatTest do end test "with already-integer floats" do - assert Float.ceil(:math.pow(2, 52), 5) === 4_503_599_627_370_496.0 assert Float.ceil(1.0e20, 3) === 1.0e20 assert Float.ceil(-1.0e20, 3) === -1.0e20 end test "with very large floats hits the bignum slow path" do + assert Float.ceil(2_661_101_816_343_531.5, 1) === 2_661_101_816_343_531.5 + assert Float.ceil(3.0e15, 1) === 3.0e15 + assert Float.ceil(-3.0e15, 1) === -3.0e15 assert Float.ceil(1.234e15, 3) === 1.234e15 assert Float.ceil(1.234567e11, 5) === 1.234567e11 assert Float.ceil(-1.234567e11, 5) === -1.234567e11 @@ -268,15 +272,23 @@ defmodule FloatTest do end test "with already-integer floats" do - assert Float.round(:math.pow(2, 52), 5) === 4_503_599_627_370_496.0 assert Float.round(1.0e20, 3) === 1.0e20 assert Float.round(-1.0e20, 3) === -1.0e20 + assert Float.round(3.0e15, 3) === 3.0e15 + assert Float.round(-3.0e15, 3) === -3.0e15 end test "with very large floats hits the bignum slow path" do + assert Float.round(3.0e15, 1) === 3.0e15 + assert Float.round(-3.0e15, 1) === -3.0e15 assert Float.round(1.234e15, 3) === 1.234e15 assert Float.round(1.234567e11, 5) === 1.234567e11 assert Float.round(-1.234567e11, 5) === -1.234567e11 + + assert Float.round(2_251_799_813_685_248.5, 1) === 2_251_799_813_685_248.5 + assert Float.round(2_251_799_813_685_249.5, 1) === 2_251_799_813_685_249.5 + assert Float.round(-2_251_799_813_685_248.5, 1) === -2_251_799_813_685_248.5 + assert Float.round(-2_251_799_813_685_249.5, 1) === -2_251_799_813_685_249.5 end test "preserves documented tie behavior" do @@ -289,6 +301,17 @@ defmodule FloatTest do end end + test "round/2, floor/2, and ceil/2 preserve powers of two" do + # Powers of two are exactly representable and have no fractional content, + # so rounding at any precision must return the input unchanged. + for k <- 0..60 do + f = :math.pow(2, k) + assert Float.round(f, 1) === f + assert Float.floor(f, 1) === f + assert Float.ceil(f, 1) === f + end + end + describe "ratio/1" do test "with 0.0" do assert Float.ratio(0.0) == {0, 1} From e9b32690713e90bfa424ba1e17292d113ae907a9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jos=C3=A9=20Valim?= Date: Fri, 1 May 2026 12:31:37 +0200 Subject: [PATCH 7/7] Apply suggestion from @josevalim --- lib/elixir/lib/float.ex | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/lib/elixir/lib/float.ex b/lib/elixir/lib/float.ex index 28f17c9166a..8f778fed54d 100644 --- a/lib/elixir/lib/float.ex +++ b/lib/elixir/lib/float.ex @@ -349,9 +349,7 @@ defmodule Float do def round(float, 0) when float == 0.0, do: float def round(float, 0) when is_float(float) do - rounded = :erlang.round(float) * 1.0 - - case rounded do + case :erlang.round(float) * 1.0 do zero when zero == 0.0 and float < 0.0 -> -0.0 rounded -> rounded end