diff --git a/benches/bigdecimal/bench.mojo b/benches/bigdecimal/bench.mojo index 1daab2d..697e964 100644 --- a/benches/bigdecimal/bench.mojo +++ b/benches/bigdecimal/bench.mojo @@ -4,6 +4,7 @@ from bench_bigdecimal_multiply import main as bench_multiply from bench_bigdecimal_divide import main as bench_divide from bench_bigdecimal_sqrt import main as bench_sqrt from bench_bigdecimal_exp import main as bench_exp +from bench_bigdecimal_ln import main as bench_ln from bench_bigdecimal_scale_up_by_power_of_10 import main as bench_scale_up @@ -19,6 +20,7 @@ mul: Multiply div: Divide (true divide) sqrt: Square root exp: Exponential +ln: Natural logarithm all: Run all benchmarks q: Exit ========================================= @@ -39,6 +41,8 @@ scaleup: Scale up by power of 10 bench_sqrt() elif command == "exp": bench_exp() + elif command == "ln": + bench_ln() elif command == "all": bench_add() bench_sub() @@ -46,6 +50,7 @@ scaleup: Scale up by power of 10 bench_divide() bench_sqrt() bench_exp() + bench_ln() elif command == "q": return elif command == "scaleup": diff --git a/benches/bigdecimal/bench_bigdecimal_ln.mojo b/benches/bigdecimal/bench_bigdecimal_ln.mojo new file mode 100644 index 0000000..79d738c --- /dev/null +++ b/benches/bigdecimal/bench_bigdecimal_ln.mojo @@ -0,0 +1,673 @@ +""" +Comprehensive benchmarks for BigDecimal logarithm function (ln). +Compares performance against Python's decimal module with diverse test cases. +""" + +from decimojo import BigDecimal, RoundingMode +from python import Python, PythonObject +from time import perf_counter_ns +import time +import os +from collections import List + + +fn open_log_file() raises -> PythonObject: + """ + Creates and opens a log file with a timestamp in the filename. + + Returns: + A file object opened for writing. + """ + var python = Python.import_module("builtins") + var datetime = Python.import_module("datetime") + + # Create logs directory if it doesn't exist + var log_dir = "./logs" + if not os.path.exists(log_dir): + os.makedirs(log_dir) + + # Generate a timestamp for the filename + var timestamp = String(datetime.datetime.now().isoformat()) + var log_filename = log_dir + "/benchmark_bigdecimal_ln_" + timestamp + ".log" + + print("Saving benchmark results to:", log_filename) + return python.open(log_filename, "w") + + +fn log_print(msg: String, log_file: PythonObject) raises: + """ + Prints a message to both the console and the log file. + + Args: + msg: The message to print. + log_file: The file object to write to. + """ + print(msg) + log_file.write(msg + "\n") + log_file.flush() # Ensure the message is written immediately + + +fn run_benchmark_ln( + name: String, + value: String, + iterations: Int, + log_file: PythonObject, + mut speedup_factors: List[Float64], +) raises: + """ + Run a benchmark comparing Mojo BigDecimal ln with Python Decimal ln. + + Args: + name: Name of the benchmark case. + value: String representation of the number to calculate the ln of. + iterations: Number of iterations to run. + log_file: File object for logging results. + speedup_factors: Mojo List to store speedup factors for averaging. + """ + log_print("\nBenchmark: " + name, log_file) + log_print("Value: " + value, log_file) + + # Set up Mojo and Python values + var mojo_value = BigDecimal(value) + var pydecimal = Python.import_module("decimal") + var py_value = pydecimal.Decimal(value) + + # Execute the operations once to verify correctness + try: + var mojo_result = mojo_value.ln() + var py_result = py_value.ln() + + # Display results for verification + log_print("Mojo result: " + String(mojo_result), log_file) + log_print("Python result: " + String(py_result), log_file) + + # Benchmark Mojo implementation + var t0 = perf_counter_ns() + for _ in range(iterations): + _ = mojo_value.ln() + var mojo_time = (perf_counter_ns() - t0) / iterations + if mojo_time == 0: + mojo_time = 1 # Prevent division by zero + + # Benchmark Python implementation + t0 = perf_counter_ns() + for _ in range(iterations): + _ = py_value.ln() + var python_time = (perf_counter_ns() - t0) / iterations + + # Calculate speedup factor + var speedup = python_time / mojo_time + speedup_factors.append(Float64(speedup)) + + # Print results with speedup comparison + log_print( + "Mojo ln: " + String(mojo_time) + " ns per iteration", + log_file, + ) + log_print( + "Python ln: " + String(python_time) + " ns per iteration", + log_file, + ) + log_print("Speedup factor: " + String(speedup), log_file) + except e: + log_print("Error occurred during benchmark: " + String(e), log_file) + log_print("Skipping this benchmark case", log_file) + + +fn main() raises: + # Open log file + var log_file = open_log_file() + var datetime = Python.import_module("datetime") + + # Create a Mojo List to store speedup factors for averaging later + var speedup_factors = List[Float64]() + + # Display benchmark header with system information + log_print("=== DeciMojo BigDecimal Logarithm (ln) Benchmark ===", log_file) + log_print("Time: " + String(datetime.datetime.now().isoformat()), log_file) + + # Try to get system info + try: + var platform = Python.import_module("platform") + log_print( + "System: " + + String(platform.system()) + + " " + + String(platform.release()), + log_file, + ) + log_print("Processor: " + String(platform.processor()), log_file) + log_print( + "Python version: " + String(platform.python_version()), log_file + ) + except: + log_print("Could not retrieve system information", log_file) + + var iterations = 100 + var pydecimal = Python().import_module("decimal") + + # Set Python decimal precision to match Mojo's + pydecimal.getcontext().prec = 28 + log_print( + "Python decimal precision: " + String(pydecimal.getcontext().prec), + log_file, + ) + log_print("Mojo decimal precision: 28", log_file) + + # Define benchmark cases + log_print( + "\nRunning decimal ln benchmarks with " + + String(iterations) + + " iterations each", + log_file, + ) + + # === BASIC LOGARITHM TESTS === + + # Case 1: ln(1) = 0 + run_benchmark_ln( + "ln(1)", + "1", + iterations, + log_file, + speedup_factors, + ) + + # Case 2: ln(e) = 1 + run_benchmark_ln( + "ln(e)", + "2.71828182845904523536028747135", + iterations, + log_file, + speedup_factors, + ) + + # Case 3: ln(2) + run_benchmark_ln( + "ln(2)", + "2", + iterations, + log_file, + speedup_factors, + ) + + # Case 4: ln(10) + run_benchmark_ln( + "ln(10)", + "10", + iterations, + log_file, + speedup_factors, + ) + + # Case 5: ln(0.5) + run_benchmark_ln( + "ln(0.5)", + "0.5", + iterations, + log_file, + speedup_factors, + ) + + # === VALUES CLOSE TO 1 (CHALLENGING FOR PRECISION) === + + # Case 6: ln(0.9) + run_benchmark_ln( + "ln(0.9)", + "0.9", + iterations, + log_file, + speedup_factors, + ) + + # Case 7: ln(0.99) + run_benchmark_ln( + "ln(0.99)", + "0.99", + iterations, + log_file, + speedup_factors, + ) + + # Case 8: ln(0.999) + run_benchmark_ln( + "ln(0.999)", + "0.999", + iterations, + log_file, + speedup_factors, + ) + + # Case 9: ln(1.001) + run_benchmark_ln( + "ln(1.001)", + "1.001", + iterations, + log_file, + speedup_factors, + ) + + # Case 10: ln(1.01) + run_benchmark_ln( + "ln(1.01)", + "1.01", + iterations, + log_file, + speedup_factors, + ) + + # Case 11: ln(1.1) + run_benchmark_ln( + "ln(1.1)", + "1.1", + iterations, + log_file, + speedup_factors, + ) + + # === SMALL VALUE TESTS === + + # Case 12: ln(0.1) + run_benchmark_ln( + "ln(0.1)", + "0.1", + iterations, + log_file, + speedup_factors, + ) + + # Case 13: ln(0.01) + run_benchmark_ln( + "ln(0.01)", + "0.01", + iterations, + log_file, + speedup_factors, + ) + + # Case 14: ln(0.001) + run_benchmark_ln( + "ln(0.001)", + "0.001", + iterations, + log_file, + speedup_factors, + ) + + # Case 15: ln(0.0001) + run_benchmark_ln( + "ln(0.0001)", + "0.0001", + iterations, + log_file, + speedup_factors, + ) + + # Case 16: ln(1e-10) + run_benchmark_ln( + "ln(1e-10)", + "1e-10", + iterations, + log_file, + speedup_factors, + ) + + # Case 17: ln(1e-20) + run_benchmark_ln( + "ln(1e-20)", + "1e-20", + iterations, + log_file, + speedup_factors, + ) + + # === LARGE VALUE TESTS === + + # Case 18: ln(100) + run_benchmark_ln( + "ln(100)", + "100", + iterations, + log_file, + speedup_factors, + ) + + # Case 19: ln(1000) + run_benchmark_ln( + "ln(1000)", + "1000", + iterations, + log_file, + speedup_factors, + ) + + # Case 20: ln(10000) + run_benchmark_ln( + "ln(10000)", + "10000", + iterations, + log_file, + speedup_factors, + ) + + # Case 21: ln(1e10) + run_benchmark_ln( + "ln(1e10)", + "1e10", + iterations, + log_file, + speedup_factors, + ) + + # Case 22: ln(1e20) + run_benchmark_ln( + "ln(1e20)", + "1e20", + iterations, + log_file, + speedup_factors, + ) + + # Case 23: ln(1e50) + run_benchmark_ln( + "ln(1e50)", + "1e50", + iterations, + log_file, + speedup_factors, + ) + + # Case 24: ln(1e100) + run_benchmark_ln( + "ln(1e100)", + "1e100", + iterations, + log_file, + speedup_factors, + ) + + # === SCIENTIFIC NOTATION TESTS === + + # Case 25: ln(1.234e-5) + run_benchmark_ln( + "ln(1.234e-5)", + "1.234e-5", + iterations, + log_file, + speedup_factors, + ) + + # Case 26: ln(5.678e12) + run_benchmark_ln( + "ln(5.678e12)", + "5.678e12", + iterations, + log_file, + speedup_factors, + ) + + # Case 27: ln(9.876e-10) + run_benchmark_ln( + "ln(9.876e-10)", + "9.876e-10", + iterations, + log_file, + speedup_factors, + ) + + # Case 28: ln(3.14159e20) + run_benchmark_ln( + "ln(3.14159e20)", + "3.14159e20", + iterations, + log_file, + speedup_factors, + ) + + # === SPECIAL VALUE TESTS === + + # Case 29: ln(PI) + run_benchmark_ln( + "ln(PI)", + "3.14159265358979323846264338328", + iterations, + log_file, + speedup_factors, + ) + + # Case 30: ln(e^2) = 2 + run_benchmark_ln( + "ln(e^2)", + "7.3890560989306502272304274606", + iterations, + log_file, + speedup_factors, + ) + + # Case 31: ln(sqrt(2)) + run_benchmark_ln( + "ln(sqrt(2))", + "1.4142135623730950488016887242", + iterations, + log_file, + speedup_factors, + ) + + # === PRECISE DECIMAL TESTS === + + # Case 32: ln(1.5) + run_benchmark_ln( + "ln(1.5)", + "1.5", + iterations, + log_file, + speedup_factors, + ) + + # Case 33: ln(2.5) + run_benchmark_ln( + "ln(2.5)", + "2.5", + iterations, + log_file, + speedup_factors, + ) + + # Case 34: ln(3.75) + run_benchmark_ln( + "ln(3.75)", + "3.75", + iterations, + log_file, + speedup_factors, + ) + + # Case 35: ln(12.34567890) + run_benchmark_ln( + "ln(12.34567890)", + "12.34567890", + iterations, + log_file, + speedup_factors, + ) + + # Case 36: ln(0.1234567890) + run_benchmark_ln( + "ln(0.1234567890)", + "0.1234567890", + iterations, + log_file, + speedup_factors, + ) + + # === POWERS OF 2 === + + # Case 37: ln(4) = ln(2^2) + run_benchmark_ln( + "ln(4)", + "4", + iterations, + log_file, + speedup_factors, + ) + + # Case 38: ln(8) = ln(2^3) + run_benchmark_ln( + "ln(8)", + "8", + iterations, + log_file, + speedup_factors, + ) + + # Case 39: ln(16) = ln(2^4) + run_benchmark_ln( + "ln(16)", + "16", + iterations, + log_file, + speedup_factors, + ) + + # Case 40: ln(32) = ln(2^5) + run_benchmark_ln( + "ln(32)", + "32", + iterations, + log_file, + speedup_factors, + ) + + # Case 41: ln(64) = ln(2^6) + run_benchmark_ln( + "ln(64)", + "64", + iterations, + log_file, + speedup_factors, + ) + + # === EXTREME VALUES === + + # Case 42: ln(1e200) + run_benchmark_ln( + "ln(1e200)", + "1e200", + iterations, + log_file, + speedup_factors, + ) + + # Case 43: ln(1e-200) + run_benchmark_ln( + "ln(1e-200)", + "1e-200", + iterations, + log_file, + speedup_factors, + ) + + # Case 44: ln(1.797693e308) - near max double + run_benchmark_ln( + "ln(1.797693e308)", + "1.797693e308", + iterations, + log_file, + speedup_factors, + ) + + # Case 45: ln(4.940656e-324) - near min double + run_benchmark_ln( + "ln(4.940656e-324)", + "4.940656e-324", + iterations, + log_file, + speedup_factors, + ) + + # === PRECISE SPECIAL CASES === + + # Case 46: ln(2.718281828459045) - ln(e) to high precision + run_benchmark_ln( + "ln(e precise)", + "2.718281828459045", + iterations, + log_file, + speedup_factors, + ) + + # Case 47: ln(1.414213562373095) - ln(√2) to high precision + run_benchmark_ln( + "ln(sqrt(2) precise)", + "1.414213562373095", + iterations, + log_file, + speedup_factors, + ) + + # Case 48: ln(1.1892071150027210667174999705...) - ln(2^(1/4)) + run_benchmark_ln( + "ln(2^(1/4))", + "1.189207115002721", + iterations, + log_file, + speedup_factors, + ) + + # Case 49: ln(1.0000000000000000000000001) - very close to 1 + run_benchmark_ln( + "ln(1.0000000000000000000000001)", + "1.0000000000000000000000001", + iterations, + log_file, + speedup_factors, + ) + + # Case 50: ln(123456789.123456789) + run_benchmark_ln( + "ln(123456789.123456789)", + "123456789.123456789", + iterations, + log_file, + speedup_factors, + ) + + # Calculate average speedup factor (ignoring any cases that might have failed) + if len(speedup_factors) > 0: + var sum_speedup: Float64 = 0.0 + for i in range(len(speedup_factors)): + sum_speedup += speedup_factors[i] + var average_speedup = sum_speedup / Float64(len(speedup_factors)) + + # Display summary + log_print( + "\n=== BigDecimal Logarithm (ln) Benchmark Summary ===", log_file + ) + log_print( + "Benchmarked: " + + String(len(speedup_factors)) + + " different ln cases", + log_file, + ) + log_print( + "Each case ran: " + String(iterations) + " iterations", log_file + ) + log_print( + "Average speedup: " + String(average_speedup) + "×", log_file + ) + + # List all speedup factors + log_print("\nIndividual speedup factors:", log_file) + for i in range(len(speedup_factors)): + log_print( + String("Case {}: {}×").format( + i + 1, round(speedup_factors[i], 2) + ), + log_file, + ) + else: + log_print("\nNo valid benchmark cases were completed", log_file) + + # Close the log file + log_file.close() + print("Benchmark completed. Log file closed.") diff --git a/mojoproject.toml b/mojoproject.toml index 389685b..5c1ebef 100644 --- a/mojoproject.toml +++ b/mojoproject.toml @@ -11,6 +11,7 @@ version = "0.2.0" [dependencies] max = "==25.2" +mpmath = ">=1.3.0,<2" [tasks] # format the code diff --git a/src/decimojo/bigdecimal/arithmetics.mojo b/src/decimojo/bigdecimal/arithmetics.mojo index b463941..0485740 100644 --- a/src/decimojo/bigdecimal/arithmetics.mojo +++ b/src/decimojo/bigdecimal/arithmetics.mojo @@ -331,32 +331,25 @@ fn true_divide( ) -fn true_divide_fast( - x1: BigDecimal, - x2: BigDecimal, - minimum_precision: Int = 28, +fn true_divide_inexact( + x1: BigDecimal, x2: BigDecimal, number_of_significant_digits: Int ) raises -> BigDecimal: - """Returns the quotient of two numbers with precision. + """Returns the quotient of two numbers with number of significant digits. This function is a faster version of true_divide, but it does not - return the exact number of digits equal to `minimum_precision`. - It returns 9 digits more than `minimum_precision` in the result. - It is recommended to use this function when the exact number of - digits is not critical. + return the exact result of the division since no extra buffer digits are + added to the dividend during calculation. + It is recommended to use this function when you already know the dividend + has enough digits to produce a result with the desired precision. Then + use rounding to get the result with the desired precision. Args: x1: The first operand (dividend). x2: The second operand (divisor). - minimum_precision: The minumum number of significant digits in the result. + number_of_significant_digits: The number of significant digits in the + result. Returns: The quotient of x1 and x2. - - Notes: - - To improve the speed: - - The dividend and divisor are scaled up always by whole words. - - The result is truncated to above the minimum precision. - - Trailing zeros are not checked and removed for exact division. """ # Check for division by zero @@ -367,7 +360,7 @@ fn true_divide_fast( if x1.coefficient.is_zero(): return BigDecimal( coefficient=BigUInt.ZERO, - scale=x1.scale - x2.scale, + scale=number_of_significant_digits, sign=x1.sign != x2.sign, ) @@ -376,12 +369,10 @@ fn true_divide_fast( var x1_digits = x1.coefficient.number_of_digits() var x2_digits = x2.coefficient.number_of_digits() - # Calculate how many buffer digits we need in the dividend - # We want: (x1_digits + buffer) - x2_digits > mininum_precision - # For quickly removing the buffer digits, it is selected to be multiple of 9. - var buffer_digits = minimum_precision - (x1_digits - x2_digits) + # Calculate how many digits we need in the dividend + # We want: x1_digits - x2_digits >= mininum_precision + var buffer_digits = number_of_significant_digits - (x1_digits - x2_digits) buffer_digits = max(0, buffer_digits) - buffer_digits = (buffer_digits // 9 + 1) * 9 # Scale up the dividend to ensure sufficient precision var scaled_x1 = x1.coefficient @@ -393,12 +384,11 @@ fn true_divide_fast( var result_scale = buffer_digits + x1.scale - x2.scale var result_digits = quotient.number_of_digits() - if result_digits > minimum_precision: - var digits_to_remove = result_digits - minimum_precision - digits_to_remove = digits_to_remove // 9 * 9 + if result_digits > number_of_significant_digits: + var digits_to_remove = result_digits - number_of_significant_digits quotient = quotient.remove_trailing_digits_with_rounding( digits_to_remove, - RoundingMode.ROUND_HALF_UP, + RoundingMode.ROUND_DOWN, remove_extra_digit_due_to_rounding=False, ) # Adjust the scale accordingly diff --git a/src/decimojo/bigdecimal/bigdecimal.mojo b/src/decimojo/bigdecimal/bigdecimal.mojo index 6b238dd..024b4a4 100644 --- a/src/decimojo/bigdecimal/bigdecimal.mojo +++ b/src/decimojo/bigdecimal/bigdecimal.mojo @@ -542,6 +542,11 @@ struct BigDecimal: """Returns the exponential of the BigDecimal number.""" return decimojo.bigdecimal.exponential.exp(self, precision) + @always_inline + fn ln(self, precision: Int = 28) raises -> Self: + """Returns the natural logarithm of the BigDecimal number.""" + return decimojo.bigdecimal.exponential.ln(self, precision) + @always_inline fn max(self, other: Self) raises -> Self: """Returns the maximum of two BigDecimal numbers.""" @@ -567,14 +572,14 @@ struct BigDecimal: ) @always_inline - fn true_divide_fast( - self, other: Self, mininum_precision: Int + fn true_divide_inexact( + self, other: Self, number_of_significant_digits: Int ) raises -> Self: - """Returns the result of true division with minimum precision. - See `arithmetics.true_divide_fast()` for more information. + """Returns the result of true division with inexact precision. + See `arithmetics.true_divide_inexact()` for more information. """ - return decimojo.bigdecimal.arithmetics.true_divide_fast( - self, other, mininum_precision + return decimojo.bigdecimal.arithmetics.true_divide_inexact( + self, other, number_of_significant_digits ) @always_inline diff --git a/src/decimojo/bigdecimal/exponential.mojo b/src/decimojo/bigdecimal/exponential.mojo index 1a1f957..826131f 100644 --- a/src/decimojo/bigdecimal/exponential.mojo +++ b/src/decimojo/bigdecimal/exponential.mojo @@ -121,7 +121,7 @@ fn sqrt(x: BigDecimal, precision: Int = 28) raises -> BigDecimal: while guess != prev_guess and iteration_count < 100: prev_guess = guess - var quotient = x.true_divide_fast(guess, working_precision) + var quotient = x.true_divide_inexact(guess, working_precision) var sum = guess + quotient guess = sum.true_divide(BigDecimal(BigUInt(2), 0, 0), working_precision) iteration_count += 1 @@ -191,7 +191,7 @@ fn exp(x: BigDecimal, precision: Int = 28) raises -> BigDecimal: - Precision tracking. """ # Extra working precision to ensure final result accuracy - alias BUFFER_DIGITS = 5 + alias BUFFER_DIGITS = 9 var working_precision = precision + BUFFER_DIGITS # Handle special cases @@ -230,7 +230,7 @@ fn exp(x: BigDecimal, precision: Int = 28) raises -> BigDecimal: k += 1 # Calculate exp(x/2^k) - var reduced_x = x.true_divide_fast(threshold, working_precision) + var reduced_x = x.true_divide_inexact(threshold, working_precision) # var t_after_range_reduction = time.perf_counter_ns() @@ -318,7 +318,7 @@ fn exp_taylor_series( for _ in range(1, max_number_of_terms): # Calculate next term: x^i/i! = x^{i-1} * x/i # We can use the previous term to calculate the next one - var add_on = x.true_divide_fast( + var add_on = x.true_divide_inexact( BigDecimal(n, 0, False), minimum_precision ) term = term * add_on @@ -346,3 +346,235 @@ fn exp_taylor_series( # print("DEBUG: final result", result) return result^ + + +# ===----------------------------------------------------------------------=== # +# Logarithmic functions +# ===----------------------------------------------------------------------=== # +fn ln(x: BigDecimal, precision: Int = 28) raises -> BigDecimal: + """Calculate the natural logarithm of x to the specified precision. + + Args: + x: The input value. + precision: Desired precision in significant digits. + + Returns: + The natural logarithm of x to the specified precision. + + Raises: + Error: If x is negative or zero. + """ + alias BUFFER_DIGITS = 9 # word-length, easy to append and trim + var working_precision = precision + BUFFER_DIGITS + + # Handle special cases + if x.sign: + raise Error( + "Error in `ln`: Cannot compute logarithm of negative number" + ) + if x.coefficient.is_zero(): + raise Error("Error in `ln`: Cannot compute logarithm of zero") + if x == BigDecimal(BigUInt.ONE, 0, False): + return BigDecimal(BigUInt.ZERO, 0, False) # ln(1) = 0 + + # Range reduction to improve convergence + # ln(x) = ln(m * 2^a * 5^b) = + # = ln(m) + a*ln(2) + b*ln(5) + # = ln(m) + a*ln(2) + b*(ln(5/4) + ln(4)) + # = ln(m) + a*ln(2) + b*(ln(1.25) + 2*ln(2)) + # = ln(m) + (a+b*2)*ln(2) + b*ln(1.25) + # where 0.5 <= m < 1.5 + # Use Taylor series for ln(m) = ln(1+z) + var m = x + var power_of_2: Int = 0 + var power_of_5: Int = 0 + # First, scale down to [0.1, 1) + var power_of_10 = m.exponent() + 1 + m.scale += power_of_10 + # Second, scale to [0.5, 1.5) + if m < BigDecimal(BigUInt(List[UInt32](135)), 3, False): + # [0.1, 0.135) * 10 -> [1, 1.35) + power_of_10 -= 1 + m.scale -= 1 + elif m < BigDecimal(BigUInt(List[UInt32](275)), 3, False): + # [0.135, 0.275) * 5 -> [0.675, 1.375)] + power_of_5 -= 1 + m = m * BigDecimal(BigUInt(List[UInt32](5)), 0, False) + elif m < BigDecimal(BigUInt(List[UInt32](65)), 2, False): + # [0.275, 0.65) * 2 -> [0.55, 1.3)] + power_of_2 -= 1 + m = m * BigDecimal(BigUInt(List[UInt32](2)), 0, False) + else: # [0.65, 1) -> no change + pass + # Replace 10 with 5 and 2 + power_of_5 += power_of_10 + power_of_2 += power_of_10 + + # print("Input: {} = {} * 2^{} * 5^{}".format(x, m, power_of_2, power_of_5)) + + # Use series expansion for ln(m) = ln(1+z) = z - z²/2 + z³/3 - ... + var result = ln_series_expansion( + m - BigDecimal(BigUInt.ONE, 0, False), working_precision + ) + + # print("Result after series expansion:", result) + + # Apply range reduction adjustments + # ln(m) + (a+b*2)*ln(2) + b*ln(1.25) + # TODO: Use precomputed ln(2) for better performance + # It is only calculated once and saved in the cache (global variable) + # Need Mojo to support global variables + if power_of_2 + power_of_5 * 2 != 0: + # if len(decimojo.cache_ln2.coefficient.words) != 0: + # var ln2 = decimojo.cache_ln2 + var ln2 = compute_ln2(working_precision) + result += ln2 * BigDecimal.from_int(power_of_2 + power_of_5 * 2) + if power_of_5 != 0: + # if len(decimojo.cache_ln5.coefficient.words) != 0: + # var ln1d25 = decimojo.cache_ln1d25 + var ln1d25 = compute_ln1d25(working_precision) + result += ln1d25 * BigDecimal.from_int(power_of_5) + + # Round to final precision + result.round_to_precision( + precision=precision, + rounding_mode=RoundingMode.ROUND_HALF_EVEN, + remove_extra_digit_due_to_rounding=True, + ) + + return result^ + + +fn ln_series_expansion( + z: BigDecimal, working_precision: Int +) raises -> BigDecimal: + """Calculate ln(1+z) using optimized series expansion. + + Args: + z: The input value, should be |z| < 1 for fast convergence. + working_precision: Desired working precision in significant digits. + + Returns: + The ln(1+z) computed to the specified working precision. + + Notes: + + The last few digits of result are not accurate as there is no buffer for + precision. You need to use a larger precision to get the last few digits + accurate. The precision is only used to determine the number of terms in + the series expansion, not for the final result. + """ + + # print("DEBUG: ln_series_expansion for z =", z) + + if z.is_zero(): + return BigDecimal(BigUInt.ZERO, 0, False) + + var max_terms = Int(working_precision * 2.5) + 1 + var result = BigDecimal(BigUInt.ZERO, working_precision, False) + var term = z + var k = BigUInt.ONE + + # Use the series ln(1+z) = z - z²/2 + z³/3 - z⁴/4 + ... + result += term # First term is just x + # print("DEBUG: term =", term, "result =", result) + # print("DEBUG: k =", k, "max_terms =", max_terms) + + for _ in range(2, max_terms): + # Update for next iteration - multiply by z and divide by k + term = term * z # z^k + k += BigUInt.ONE + + # Alternate sign: -1^(k+1) = -1 when k is even, 1 when k is odd + var sign = (k % BigUInt(2) == BigUInt.ZERO) + var next_term = term.true_divide_inexact( + BigDecimal(k, 0, False), working_precision + ) + + if sign: + result -= next_term + else: + result += next_term + + # print("DEBUG: k =", k, "max_terms =", max_terms) + # print("DEBUG: term =", term, "next_term =", next_term) + # print("DEBUG: result =", result) + + # Check for convergence + if next_term.exponent() < -working_precision: + break + + # print("DEBUG: ln_series_expansion result:", result) + result.round_to_precision( + precision=working_precision, + rounding_mode=RoundingMode.ROUND_DOWN, + remove_extra_digit_due_to_rounding=False, + ) + return result^ + + +fn compute_ln2(working_precision: Int) raises -> BigDecimal: + """Compute ln(2) to the specified working precision. + + Args: + working_precision: Desired precision in significant digits. + + Returns: + The ln(2) computed to the specified precision. + + Notes: + + The last few digits of result are not accurate as there is no buffer for + precision. You need to use a larger precision to get the last few digits + accurate. The precision is only used to determine the number of terms in + the series expansion, not for the final result. + """ + # Directly using Taylor series expansion for ln(2) is not efficient + # Instead, we can use the identity: + # ln((1+x)/(1-x)) = 2*arcth(x) = 2*(x + x³/3 + x⁵/5 + ...) + # For x = 1/3: + # ln(2) = 2*(1/3 + (1/3)³/3 + (1/3)⁵/5 + ...) + + var max_terms = Int(working_precision * 2.5) + 1 + + var number_of_words = working_precision // 9 + 1 + var words = List[UInt32](capacity=number_of_words) + for _ in range(number_of_words): + words.append(UInt32(333_333_333)) + var x = BigDecimal(BigUInt(words), number_of_words * 9, False) # x = 1/3 + + var result = BigDecimal(BigUInt.ZERO, 0, False) + var term = x * BigDecimal(BigUInt(2), 0, False) # First term: 2*(1/3) + var k = BigDecimal(BigUInt.ONE, 0, False) + + # Add terms: 2*(x + x³/3 + x⁵/5 + ...) + for _ in range(1, max_terms): + result += term + var new_k = k + BigDecimal(BigUInt(List[UInt32](2)), 0, False) + term = (term * x * x * k).true_divide_inexact(new_k, working_precision) + k = new_k^ + + # print("DEBUG: k =", k, "max_terms =", max_terms) + + if term.exponent() < -working_precision: + break + result.round_to_precision( + precision=working_precision, + rounding_mode=RoundingMode.ROUND_DOWN, + remove_extra_digit_due_to_rounding=False, + ) + return result^ + + +fn compute_ln1d25(precision: Int) raises -> BigDecimal: + """Compute ln(1.25) to the specified precision. + + Args: + precision: Desired precision in significant digits. + + Returns: + The ln(1.25) computed to the specified precision. + """ + var z = BigDecimal(BigUInt(List[UInt32](25)), 2, False) + var result = ln_series_expansion(z^, precision) + return result^ diff --git a/src/decimojo/bigdecimal/rounding.mojo b/src/decimojo/bigdecimal/rounding.mojo index be7caad..d3dd053 100644 --- a/src/decimojo/bigdecimal/rounding.mojo +++ b/src/decimojo/bigdecimal/rounding.mojo @@ -48,9 +48,31 @@ fn round_to_precision( var ndigits_coefficient = number.coefficient.number_of_digits() var ndigits_to_remove = ndigits_coefficient - precision - number.coefficient = number.coefficient.remove_trailing_digits_with_rounding( - ndigits=ndigits_to_remove, - rounding_mode=rounding_mode, - remove_extra_digit_due_to_rounding=remove_extra_digit_due_to_rounding, + + if ndigits_to_remove == 0: + return + + if ndigits_to_remove < 0: + number = number.extend_precision(precision_diff=-ndigits_to_remove) + return + + number.coefficient = ( + number.coefficient.remove_trailing_digits_with_rounding( + ndigits=ndigits_to_remove, + rounding_mode=rounding_mode, + remove_extra_digit_due_to_rounding=False, + ) ) number.scale -= ndigits_to_remove + + if remove_extra_digit_due_to_rounding and ( + number.coefficient.number_of_digits() > precision + ): + number.coefficient = ( + number.coefficient.remove_trailing_digits_with_rounding( + ndigits=1, + rounding_mode=RoundingMode.ROUND_DOWN, + remove_extra_digit_due_to_rounding=False, + ) + ) + number.scale -= 1 diff --git a/src/decimojo/biguint/biguint.mojo b/src/decimojo/biguint/biguint.mojo index cae3e86..7755980 100644 --- a/src/decimojo/biguint/biguint.mojo +++ b/src/decimojo/biguint/biguint.mojo @@ -227,7 +227,7 @@ struct BigUInt(Absable, IntableRaising, Writable): return Self() if value < 0: - raise Error("Error in `from_int()`: The value is negative") + raise Error("Error in `BigUInt.from_int()`: The value is negative") var list_of_words = List[UInt32]() var remainder: Int = value @@ -277,13 +277,14 @@ struct BigUInt(Absable, IntableRaising, Writable): else: if value != value: # Check for NaN raise Error( - "Error in `from_scalar()`: Cannot convert NaN to BigUInt" + "Error in `BigUInt.from_scalar()`: Cannot convert NaN to" + " BigUInt" ) # Convert to string with full precision try: return Self.from_string(String(value)) except e: - raise Error("Error in `from_scalar()`: ", e) + raise Error("Error in `BigUInt.from_scalar()`: ", e) return Self() diff --git a/tests/bigdecimal/test_bigdecimal_exponential.mojo b/tests/bigdecimal/test_bigdecimal_exponential.mojo index d392b2c..bd45f72 100644 --- a/tests/bigdecimal/test_bigdecimal_exponential.mojo +++ b/tests/bigdecimal/test_bigdecimal_exponential.mojo @@ -1,5 +1,5 @@ """ -Test BigDecimal exponential operations including square root. +Test BigDecimal exponential operations including square root and natural logarithm. """ from python import Python @@ -108,6 +108,89 @@ fn test_negative_sqrt() raises: print("✓ Square root of negative number correctly raises an error") +fn test_ln() raises: + """Test BigDecimal natural logarithm with various test cases.""" + print("------------------------------------------------------") + print("Testing BigDecimal natural logarithm (ln)...") + + var pydecimal = Python.import_module("decimal") + + # Load test cases from TOML file + var test_cases = load_test_cases(exponential_file_path, "ln_tests") + print("Loaded", len(test_cases), "test cases for natural logarithm") + + # Track test results + var passed = 0 + var failed = 0 + + # Run all test cases in a loop + for i in range(len(test_cases)): + var test_case = test_cases[i] + var input_value = BigDecimal(test_case.a) + var expected = BigDecimal(test_case.expected) + + # Calculate natural logarithm + var result = input_value.ln() + + try: + # Using String comparison for easier debugging + testing.assert_equal( + String(result), String(expected), test_case.description + ) + passed += 1 + except e: + print( + "=" * 50, + "\n", + i + 1, + "failed:", + test_case.description, + "\n Input:", + test_case.a, + "\n Expected:", + test_case.expected, + "\n Got:", + String(result), + "\n Python decimal result (for reference):", + String(pydecimal.Decimal(test_case.a).ln()), + ) + failed += 1 + + print("BigDecimal ln tests:", passed, "passed,", failed, "failed") + testing.assert_equal(failed, 0, "All natural logarithm tests should pass") + + +fn test_ln_invalid_inputs() raises: + """Test that natural logarithm with invalid inputs raises appropriate errors. + """ + print("------------------------------------------------------") + print("Testing BigDecimal natural logarithm with invalid inputs...") + + # Test 1: ln of zero should raise an error + var zero = BigDecimal("0") + var exception_caught = False + try: + _ = zero.ln() + exception_caught = False + except: + exception_caught = True + testing.assert_true(exception_caught, "ln(0) should raise an error") + print("✓ ln(0) correctly raises an error") + + # Test 2: ln of negative number should raise an error + var negative = BigDecimal("-1") + exception_caught = False + try: + _ = negative.ln() + exception_caught = False + except: + exception_caught = True + testing.assert_true( + exception_caught, "ln of negative number should raise an error" + ) + print("✓ ln of negative number correctly raises an error") + + fn main() raises: print("Running BigDecimal exponential tests") @@ -117,4 +200,10 @@ fn main() raises: # Test sqrt of negative number test_negative_sqrt() + # Run ln tests + test_ln() + + # Test ln with invalid inputs + test_ln_invalid_inputs() + print("All BigDecimal exponential tests passed!") diff --git a/tests/bigdecimal/test_data/bigdecimal_exponential.toml b/tests/bigdecimal/test_data/bigdecimal_exponential.toml index 6a9ac5c..130dc79 100644 --- a/tests/bigdecimal/test_data/bigdecimal_exponential.toml +++ b/tests/bigdecimal/test_data/bigdecimal_exponential.toml @@ -150,3 +150,140 @@ description = "Square root of π" input = "2.71828182845904523536" expected = "1.648721270700128146848563608" description = "Square root of e" + +# === BASIC NATURAL LOGARITHM TESTS === +[[ln_tests]] +input = "1" +expected = "0" +description = "Natural logarithm of 1 (exactly 0)" + +[[ln_tests]] +input = "2.718281828459045235360287471" +expected = "0.9999999999999999999999999999" +description = "Natural logarithm of e (approximately 1)" + +[[ln_tests]] +input = "10" +expected = "2.302585092994045684017991455" +description = "Natural logarithm of 10" + +[[ln_tests]] +input = "2" +expected = "0.6931471805599453094172321215" +description = "Natural logarithm of 2" + +[[ln_tests]] +input = "0.5" +expected = "-0.6931471805599453094172321215" +description = "Natural logarithm of 0.5 (negative result)" + +# === POWERS OF e === +[[ln_tests]] +input = "7.389056098930650227230427461" +expected = "2.000000000000000000000000000" +description = "Natural logarithm of e^2 (exactly 2)" + +[[ln_tests]] +input = "20.08553692318766774092852965" +expected = "3.000000000000000000000000000" +description = "Natural logarithm of e^3 (exactly 3)" + +[[ln_tests]] +input = "0.3678794411714423215955237702" +expected = "-0.9999999999999999999999999999" +description = "Natural logarithm of 1/e (exactly -1)" + +[[ln_tests]] +input = "0.1353352832366126918939994949" +expected = "-2.000000000000000000000000001" +description = "Natural logarithm of 1/e^2 (exactly -2)" + +# === POWERS OF 10 === +[[ln_tests]] +input = "100" +expected = "4.605170185988091368035982909" +description = "Natural logarithm of 10^2" + +[[ln_tests]] +input = "0.01" +expected = "-4.605170185988091368035982909" +description = "Natural logarithm of 10^-2" + +[[ln_tests]] +input = "1000" +expected = "6.907755278982137052053974364" +description = "Natural logarithm of 10^3" + +[[ln_tests]] +input = "0.001" +expected = "-6.907755278982137052053974364" +description = "Natural logarithm of 10^-3" + +# === SMALL NUMBERS === +[[ln_tests]] +input = "0.1" +expected = "-2.302585092994045684017991455" +description = "Natural logarithm of 0.1" + +[[ln_tests]] +input = "0.00001" +expected = "-11.51292546497022842008995727" +description = "Natural logarithm of 1e-5" + +[[ln_tests]] +input = "0.0000000001" +expected = "-23.02585092994045684017991455" +description = "Natural logarithm of 1e-10" + +# === LARGE NUMBERS === +[[ln_tests]] +input = "1000000" +expected = "13.81551055796427410410794873" +description = "Natural logarithm of 1e6" + +[[ln_tests]] +input = "1000000000000" +expected = "27.63102111592854820821589746" +description = "Natural logarithm of 1e12" + +[[ln_tests]] +input = "1e50" +expected = "115.1292546497022842008995727" +description = "Natural logarithm of 1e50" + +# === VALUES NEAR 1 === +[[ln_tests]] +input = "0.9" +expected = "-0.1053605156578263012275009808" +description = "Natural logarithm of 0.9 (near 1)" + +[[ln_tests]] +input = "1.1" +expected = "0.09531017980432486004395212328" +description = "Natural logarithm of 1.1 (near 1)" + +[[ln_tests]] +input = "0.999999" +expected = "-0.000001000000500000333333583333533" +description = "Natural logarithm very close to 0 (0.999999)" + +[[ln_tests]] +input = "1.000001" +expected = "9.999995000003333330833335333E-7" +description = "Natural logarithm very close to 0 (1.000001)" + +# === MATHEMATICAL CONSTANTS === +[[ln_tests]] +input = "3.141592653589793238462643383" +expected = "1.144729885849400174143427351" +description = "Natural logarithm of π" + +[[ln_tests]] +input = "1.618033988749894848204586834" +expected = "0.4812118250596034474977589132" +description = "Natural logarithm of φ (golden ratio)" + +[[ln_tests]] +input = "1.414213562373095048801688724" +expected = "0.3465735902799726547086160606" +description = "Natural logarithm of √2"