[decimojo] Improve the division operation for BigDecimal and others

README.md

@@ -2,7 +2,7 @@
 
 A comprehensive decimal and integer mathematics library for [Mojo](https://www.modular.com/mojo).
 
-**[中文·漢字»](https://zhuyuhao.com/decimojo/docs/readme_zht.html)**　|　**[Repository on GitHub»](https://github.com/forfudan/decimojo)**
+**[中文·漢字»](https://zhuyuhao.com/decimojo/docs/readme_zht.html)**　|　**[Repository on GitHub»](https://github.com/forfudan/decimojo)**  | **[Changelog](https://zhuyuhao.com/decimojo/docs/changelog.html)**
 
 ## Overview
 

docs/changelog.md

@@ -1,6 +1,6 @@
-# DeciMojo released changelog
+# DeciMojo changelog
 
-This is a list of RELEASED changes for the DeciMojo Package.
+This is a list of RELEASED changes for the DeciMojo Package. For the unreleased changes, please refer to **[changelog_unreleased](https://zhuyuhao.com/decimojo/docs/changelog_unreleased.html)**.
 
 ## 01/04/2025 (v0.2.0)
 

docs/changelog_unreleased.md

@@ -0,0 +1,11 @@
+# DeciMojo unreleased changelog
+
+This is a list of UNRELEASED changes for the DeciMojo Package. For the released changes, please refer to **[changelog](https://zhuyuhao.com/decimojo/docs/changelog.html)**.
+
+### ⭐️ New
+
+- Add comprehensive `BigDecimal` implementation with unlimited precision arithmetic.
+
+### 🛠️ Fixed
+
+- Fix a bug in `BigUInt` multiplication where the calcualtion of carry is mistakenly skipped if a word of x2 is zero (PR #70).

mojoproject.toml

@@ -10,7 +10,7 @@ readme = "README.md"
 version = "0.2.0"
 
 [dependencies]
-max = ">=25.2"
+max = "==25.2"
 
 [tasks]
 # format the code
@@ -31,7 +31,11 @@ test = "magic run package && magic run mojo test tests --filter"
 t = "clear && magic run package && magic run mojo test tests --filter"
 
 # benches
-bench_dec = "clear && magic run package && cd benches/decimal && magic run mojo -I ../ bench.mojo && cd ../.. && magic run clean"
-bench_bint = "clear && magic run package && cd benches/bigint && magic run mojo -I ../ bench.mojo && cd ../.. && magic run clean"
-bench_buint = "clear && magic run package && cd benches/biguint && magic run mojo -I ../ bench.mojo && cd ../.. && magic run clean"
-bench_bdec = "clear && magic run package && cd benches/bigdecimal && magic run mojo -I ../ bench.mojo && cd ../.. && magic run clean"
+bench_decimal = "clear && magic run package && cd benches/decimal && magic run mojo -I ../ bench.mojo && cd ../.. && magic run clean"
+bench_bigint = "clear && magic run package && cd benches/bigint && magic run mojo -I ../ bench.mojo && cd ../.. && magic run clean"
+bench_biguint = "clear && magic run package && cd benches/biguint && magic run mojo -I ../ bench.mojo && cd ../.. && magic run clean"
+bench_bigdecimal = "clear && magic run package && cd benches/bigdecimal && magic run mojo -I ../ bench.mojo && cd ../.. && magic run clean"
+bench_dec = "magic run bench_decimal"
+bench_bint = "magic run bench_bigint"
+bench_buint = "magic run bench_biguint"
+bench_bdec = "magic run bench_bigdecimal"

src/decimojo/bigdecimal/arithmetics.mojo

@@ -188,31 +188,32 @@ fn multiply(x1: BigDecimal, x2: BigDecimal) raises -> BigDecimal:
     )
 
 
+# TODO: Optimize when divided by power of 10
 fn true_divide(
-    x1: BigDecimal, x2: BigDecimal, max_precision: Int = 28
+    x1: BigDecimal, x2: BigDecimal, precision: Int = 28
 ) raises -> BigDecimal:
     """Returns the quotient of two numbers.
 
     Args:
         x1: The first operand (dividend).
         x2: The second operand (divisor).
-        max_precision: The maximum precision for the result. It should be
-            non-negative.
+        precision: The number of significant digits in the result.
 
     Returns:
-        The quotient of x1 and x2, with precision up to max_precision.
+        The quotient of x1 and x2, with precision up to `precision`
+        significant digits.
 
     Notes:
 
     - If the coefficients can be divided exactly, the number of digits after
         the decimal point is the difference of the scales of the two operands.
     - If the coefficients cannot be divided exactly, the number of digits after
-        the decimal point is max_precision.
+        the decimal point is precision.
     - If the division is not exact, the number of digits after the decimal
-        point is calcuated to max_precision + BUFFER_DIGITS, and the result is
-        rounded to max_precision according to the specified rules.
+        point is calcuated to precision + BUFFER_DIGITS, and the result is
+        rounded to precision according to the specified rules.
     """
-    alias BUFFER_DIGITS = 2  # Buffer digits for rounding
+    alias BUFFER_DIGITS = 9  # Buffer digits for rounding
 
     # Check for division by zero
     if x2.coefficient.is_zero():
@@ -221,87 +222,89 @@ fn true_divide(
     # Handle dividend of zero
     if x1.coefficient.is_zero():
         return BigDecimal(
-            coefficient=BigUInt(UInt32(0)),
-            scale=max(0, x1.scale - x2.scale),
+            coefficient=BigUInt.ZERO,
+            scale=x1.scale - x2.scale,
             sign=x1.sign != x2.sign,
         )
 
-    # TODO: Divided by power of 10
+    # First estimate the number of significant digits needed in the dividend
+    # to produce a result with precision significant digits
+    var x1_digits = x1.coefficient.number_of_digits()
+    var x2_digits = x2.coefficient.number_of_digits()
 
-    # Check whether the coefficients can be divided exactly
+    # Check whether the coefficients can already be divided exactly
     # If division is exact, return the result immediately
-    if len(x1.coefficient.words) >= len(x2.coefficient.words):
-        # Check if x1 is divisible by x2
+    if x1_digits >= x2_digits:
         var quotient: BigUInt
         var remainder: BigUInt
         quotient, remainder = x1.coefficient.divmod(x2.coefficient)
+        # Calculate the expected result scale
+        var result_scale = x1.scale - x2.scale
         if remainder.is_zero():
-            return BigDecimal(
-                coefficient=quotient,
-                scale=x1.scale - x2.scale,
-                sign=x1.sign != x2.sign,
-            )
+            # For exact division, calculate significant digits in result
+            var num_sig_digits = quotient.number_of_digits()
+            # If the significant digits are within precision, return as is
+            if num_sig_digits <= precision:
+                return BigDecimal(
+                    coefficient=quotient^,
+                    scale=result_scale,
+                    sign=x1.sign != x2.sign,
+                )
+            else:  # num_sig_digits > precision
+                # Otherwise, need to truncate to max precision
+                var digits_to_remove = num_sig_digits - precision
+                var quotient = quotient.remove_trailing_digits_with_rounding(
+                    digits_to_remove,
+                    RoundingMode.ROUND_HALF_EVEN,
+                    remove_extra_digit_due_to_rounding=True,
+                )
+                result_scale -= digits_to_remove
+                return BigDecimal(
+                    coefficient=quotient^,
+                    scale=result_scale,
+                    sign=x1.sign != x2.sign,
+                )
+
+    # Calculate how many additional digits we need in the dividend
+    # We want: (x1_digits + additional) - x2_digits ≈ precision
+    # Scale factor needs to account for the scales of the operands?
+    var additional_digits = precision + BUFFER_DIGITS - (x1_digits - x2_digits)
+    additional_digits = max(0, additional_digits)
+
+    # Scale up the dividend to ensure sufficient precision
+    var scaled_x1 = x1.coefficient
+    if additional_digits > 0:
+        scaled_x1 = scaled_x1.scale_up_by_power_of_10(additional_digits)
+
+    # Perform division
+    var quotient: BigUInt
+    var remainder: BigUInt
+    quotient, remainder = scaled_x1.divmod(x2.coefficient)
+    var result_scale = additional_digits + x1.scale - x2.scale
 
-    # Calculate how many extra digits we need to scale x1 by
-    # We want (max_precision + BUFFER_DIGITS) decimal places in the result
-    var desired_result_scale = max_precision + BUFFER_DIGITS
-    var current_result_scale = x1.scale - x2.scale
-    var scale_factor = max(0, desired_result_scale - current_result_scale)
-
-    # Scale the dividend coefficient
-    var scaled_x1_coefficient = x1.coefficient
-    if scale_factor > 0:
-        scaled_x1_coefficient = x1.coefficient.scale_up_by_power_of_10(
-            scale_factor
-        )
+    # Check if division is exact
+    var is_exact = remainder.is_zero()
 
-    # Perform the division and get remainder
-    var result_coefficient: BigUInt
-    var remainder: BigUInt
-    result_coefficient, remainder = scaled_x1_coefficient.divmod(x2.coefficient)
-    var result_scale = x1.scale + scale_factor - x2.scale
+    # Check total digits in the result
+    var result_digits = quotient.number_of_digits()
 
     # If the division is exact
     # we may need to remove the extra trailing zeros.
     # TODO: Think about the behavior, whether division should always return the
-    # maximum precision even if the result scale is less than max_precision.
-    # Example: 1 / 1 = 1.0000000000000000000000000000
-    if remainder.is_zero():
-        # result_scale == scale_factor + (x1.scale - x2.scale)
-        var number_of_trailing_zeros = result_coefficient.number_of_trailing_zeros()
-
-        # If number_of_trailing_zeros <= scale_factor:
-        #   Just remove the trailing zeros, the scale is larger than expected
-        #   scale (x1.scale - x2.scale) because the division is exact but with
-        #   fractional part.
-        # If number_of_trailing_zeros > scale_factor:
-        #   We can remove at most scale_factor digits because the result scale
-        #   should be no less than expected scale
-        var number_of_zeros_to_remove = min(
-            number_of_trailing_zeros, scale_factor
-        )
-        result_coefficient = result_coefficient.scale_down_by_power_of_10(
-            number_of_zeros_to_remove
-        )
-
-        return BigDecimal(
-            coefficient=result_coefficient^,
-            scale=result_scale - number_of_zeros_to_remove,
-            sign=x1.sign != x2.sign,
-        )
+    # `precision` even if the result scale is less than precision.
+    # Example: 10 / 4 = 2.50000000000000000000000000000
+    # If exact division, remove trailing zeros
+    if is_exact:
+        var num_trailing_zeros = quotient.number_of_trailing_zeros()
+        if num_trailing_zeros > 0:
+            quotient = quotient.scale_down_by_power_of_10(num_trailing_zeros)
+            result_scale -= num_trailing_zeros
+            # Recalculate digits after removing trailing zeros
+            result_digits = quotient.number_of_digits()
 
     # Otherwise, the division is not exact or have too many digits
-    # round to max_precision
-    # TODO: Use round() function when it is available
-    var digits_to_remove = result_scale - max_precision
-    if digits_to_remove > BUFFER_DIGITS:
-        print(
-            "Warning: Remove (={}) more than BUFFER_DIGITS digits (={}), the"
-            " algorithm may not be optimal.".format(
-                digits_to_remove, BUFFER_DIGITS
-            )
-        )
-
+    # round to precision
+    # If we have too many significant digits, reduce to precision
     # Extract the digits to be rounded
     # Example: 2 digits to remove
     # divisor = 100
@@ -312,26 +315,18 @@ fn true_divide(
     # If rounding_digits == half_divisor, round up if the last digit of
     # result_coefficient is odd
     # If rounding_digits < half_divisor, round down
-    var divisor = BigUInt.ONE.scale_up_by_power_of_10(digits_to_remove)
-    var half_divisor = divisor // BigUInt(2)
-    var rounding_digits: BigUInt
-    result_coefficient, rounding_digits = result_coefficient.divmod(divisor)
-
-    # Apply rounding rules
-    var round_up = False
-    if rounding_digits > half_divisor:
-        round_up = True
-    elif rounding_digits == half_divisor:
-        round_up = result_coefficient.words[0] % 2 == 1
-
-    if round_up:
-        result_coefficient += BigUInt(1)
-
-    # Update scale
-    result_scale = max_precision
+    if result_digits > precision:
+        var digits_to_remove = result_digits - precision
+        quotient = quotient.remove_trailing_digits_with_rounding(
+            digits_to_remove,
+            RoundingMode.ROUND_HALF_EVEN,
+            remove_extra_digit_due_to_rounding=True,
+        )
+        # Adjust the scale accordingly
+        result_scale -= digits_to_remove
 
     return BigDecimal(
-        coefficient=result_coefficient^,
+        coefficient=quotient^,
         scale=result_scale,
         sign=x1.sign != x2.sign,
     )