# 
# Approximating the exponential function correctly (Python)
# 
# Copyright (c) 2020 Project Nayuki
# All rights reserved. Contact Nayuki for licensing.
# https://www.nayuki.io/page/approximating-eulers-number-correctly
# 


# Runs a demo that prints out some numbers
def main():
	for i in range(31):  # x = 0.0, 0.1, 0.2, ..., 3.0
		for j in range(1, 11):  # Accuracy 1 to 10
			x = i * 10 ** (j - 1)
			print(f"exp({format_decimal(x, j):12}) = {compute_exp(x, j)}")
		print()


# For example: compute_exp(20000, 4) = "7.3891"
def compute_exp(x, accuracy):
	if accuracy < 0:
		raise ValueError()
	if x < 0:
		raise ValueError("Negative numbers not supported")
	if x == 0:
		return format_decimal(10 ** accuracy, accuracy)
	
	extra_precision = x * 4343 // 10 ** (accuracy + 4) + 10  # Initial estimate based on x / log(10)
	while True:
		result = compute_exp_internal(x, accuracy, extra_precision)
		if result is not None:
			return result
		extra_precision += 2


def compute_exp_internal(x, accuracy, extra_precision):
	accuracy_scaler = 10 ** accuracy
	extra_scaler    = 10 ** extra_precision
	full_scaler = accuracy_scaler * extra_scaler
	
	sum_low  = 0
	sum_high = 0
	term_low  = full_scaler
	term_high = full_scaler
	floor_x = x // accuracy_scaler
	i = 0
	while term_low > 0:
		sum_low  += term_low
		sum_high += term_high
		term_low  = term_low  * x // accuracy_scaler
		term_high = term_high * x // accuracy_scaler + 1
		
		if i > floor_x and term_high < extra_scaler:
			sum_upper_bound = sum_high + term_high
			temp = divide_and_round(sum_low, extra_scaler)
			if divide_and_round(sum_upper_bound, extra_scaler) == temp:
				# Note: The number of terms used is i+1
				return format_decimal(temp, accuracy)
		
		i += 1
		term_low  = term_low  // i
		term_high = term_high // i + 1


# Any rounding mode works correctly with compute_eulers_number_internal().
# Round-half-to-even is implemented here, but truncation, flooring, etc. are acceptable too.
def divide_and_round(num, div):
	quot = num // div
	rem = num % div
	if rem * 2 > div or (rem * 2 == div and quot & 1 == 1):
		quot += 1
	return quot


def format_decimal(num, accuracy):
	s = str(num).rjust(accuracy + 1, '0')
	i = len(s) - accuracy  # Do not try to simplify this to '-accuracy' because it fails when accuracy == 0
	return s[ : i] + "." + s[i : ]


if __name__ == "__main__":
	main()