pow and exp

Math/Core.scope

@@ -61,6 +61,25 @@ func dec powInt(dec base, int exponent) {
 	}
 }
 
+/%
+@approx: This function could have some inaccuracies at higher values.
+	Use @"powInt" if the @"exponent" is an integer for a more accurate result.
+@todo: Add support for negative bases.
+
+Computes @"base" ^ @"exponent".
+
+This function uses @"exp" and @"ln" to compute the power.
+%/
+func dec pow(dec base, dec exponent) {
+	// log10(base ^ exponent) = log10(base) * exponent 
+
+	if (base < 0.0) {
+		ret nan;
+	}
+
+	ret exp(ln(base) * exponent);
+}
+
 /%
 Computes the log (base 10) of @"x".
 
@@ -90,8 +109,8 @@ func dec log10(dec x) {
 		ret infinity;
 	}
 
-	// The Taylor series below only works between 0 and 1
-	// Transform x to fit in the range with a summand (b)
+	// The Taylor series below only works between 0 and 1.
+	// Transform x to fit in the range with a summand (b).
 	// log10(x : 0..inf) = log10(x : 0..=1) + b
 
 	// Get "b"
@@ -122,7 +141,7 @@ func dec log10(dec x) {
 	// Finally, return
 	// log10(x) = ln(x) * log10(e) + b
 
-	// log10(e) = 0.434294481903251
+	// log10(e)   = 0.434294481903251
 	ret seriesSum * 0.434294481903251 + b;
 }
 
@@ -149,7 +168,73 @@ func dec log(dec x, dec base) {
 		// TODO: prevent x <= 0.0
 		ret nan;
 	}
-
+	
 	// log_base(x) = log10(x) / log10(base)
 	ret log10(x) / log10(base);
+}
+
+/%
+Computes e ^ @"x".
+
+This function uses the Taylor series to compute the result.
+%/
+func dec exp(dec x) {
+	// Temp until global constants
+	dec E = 2.718281828459045;
+
+	// Special cases
+
+	if (x == -infinity) {
+		ret 0.0;
+	}
+
+	if (x == 0.0) {
+		ret 1.0;
+	}
+
+	if (x == infinity) {
+		ret infinity;
+	}
+
+	// Deal with negatives later (see below)
+
+	bool neg = x < 0.0;
+	if (neg) {
+		x = -x;
+	}
+
+	// The Taylor series below is only accurate between -2 and 2
+	// Transform x to fit in the range.
+	// e ^ 8.2 = e^3.2 * e^6.0
+	// Note that the exponent in the second power is an integer
+	// which we can compute already. Let b be the second power.
+	// This does not work if x is negative. We can deal with that
+	// by doing `1.0 / result` at the end.
+
+	// Get "b"
+
+	dec b = 1.0;
+	while (x > 2.0) {
+		x -= 1.0;
+		b *= E;
+	}
+
+	// Taylor series (for seriesSum)
+
+	// exp(x) = (x^0 / 0!) + (x^1 / 1!) + (x^2 / 2!) + ...
+	dec seriesSum = 1.0;
+	int factorial = 1;
+	for (int k : 1..23) {
+		factorial *= k;
+		seriesSum += powInt(x, k) / factorial -> dec;
+	}
+
+	// Finally, return
+
+	dec result = seriesSum * b;
+	if (neg) {
+		ret 1.0 / result;
+	} else {
+		ret result;
+	}
 }