feat(stdlib): Reimplement Number.gamma and Number.factorial (#2182)

spotandjake · ospencer · web-flow · commit 1e5f921a436d · 2024-11-07T03:53:55.000Z
Co-authored-by: Oscar Spencer &lt;oscar.spen@gmail.com&gt;
diff --git a/compiler/test/stdlib/number.test.gr b/compiler/test/stdlib/number.test.gr
@@ -1057,3 +1057,36 @@ assert Number.isClose(
 assert Number.isNaN(Number.tan(Infinity))
 assert Number.isNaN(Number.tan(-Infinity))
 assert Number.isNaN(Number.tan(NaN))
+
+// Number.gamma
+// Note: Currently gamma will overflow the memory when using a bigint as such there are no tests for this
+assert Number.gamma(1) == 1
+assert Number.gamma(3) == 2
+assert Number.gamma(10) == 362880
+assert Number.isClose(Number.gamma(0.5), Number.sqrt(Number.pi))
+assert Number.isClose(Number.gamma(1/2), Number.sqrt(Number.pi))
+assert Number.isClose(Number.gamma(5/2), (3/4) * Number.sqrt(Number.pi))
+assert Number.isClose(Number.gamma(7/2), (15/8) * Number.sqrt(Number.pi))
+assert Number.isClose(Number.gamma(-0.5), -2 * Number.sqrt(Number.pi))
+assert Number.isClose(Number.gamma(-1/2), -2 * Number.sqrt(Number.pi))
+assert Number.isClose(Number.gamma(-3/2), (4/3) * Number.sqrt(Number.pi))
+assert Number.isClose(Number.gamma(-5/2), (-8/15) * Number.sqrt(Number.pi))
+assert Number.isNaN(Number.gamma(-1))
+assert Number.isNaN(Number.gamma(-10))
+assert Number.isNaN(Number.gamma(-100))
+assert Number.isNaN(Number.gamma(Infinity))
+assert Number.isNaN(Number.gamma(NaN))
+
+// Number.factorial
+// Note: Currently factorial will overflow the memory when using a bigint as such there are no tests for this
+assert Number.factorial(0) == 1
+assert Number.factorial(10) == 3628800
+assert Number.factorial(5) == 120
+assert Number.factorial(0) == 1
+assert Number.factorial(-5) == -120
+assert Number.factorial(-10) == -3628800
+assert Number.isClose(Number.factorial(0.5), (1/2) * Number.sqrt(Number.pi))
+assert Number.isClose(Number.factorial(1/2), (1/2) * Number.sqrt(Number.pi))
+assert Number.isClose(Number.factorial(5/2), (15/8) * Number.sqrt(Number.pi))
+assert Number.isNaN(Number.factorial(Infinity))
+assert Number.isNaN(Number.factorial(NaN))
diff --git a/stdlib/number.gr b/stdlib/number.gr
@@ -1118,3 +1118,86 @@ provide let tan = (radians: Number) => {
     WasmI32.toGrain(newFloat64(value)): Number
   }
 }
+
+// Math.gamma implemented using the Lanczos approximation
+// https://en.wikipedia.org/wiki/Lanczos_approximation
+/**
+ * Computes the gamma function of a value using the Lanczos approximation.
+ *
+ * @param z: The value to interpolate
+ * @returns The gamma of the given value
+ *
+ * @example Number.gamma(1) == 1
+ * @example Number.gamma(3) == 2
+ * @example Number.isClose(Number.gamma(0.5), Number.sqrt(Number.pi))
+ *
+ * @since v0.7.0
+ */
+provide let rec gamma = z => {
+  if (z == 0 || isInteger(z) && z < 0) {
+    NaN
+  } else if (isInteger(z) && z > 0) {
+    let mut output = 1
+    for (let mut i = 1; i < z; i += 1) {
+      output *= i
+    }
+    output
+  } else {
+    let mut z = z
+    let g = 7
+    let c = [>
+      0.99999999999980993,
+      676.5203681218851,
+      -1259.1392167224028,
+      771.32342877765313,
+      -176.61502916214059,
+      12.507343278686905,
+      -0.13857109526572012,
+      9.9843695780195716e-6,
+      1.5056327351493116e-7,
+    ]
+    let mut output = 0
+    if (z < 0.5) {
+      output = pi / sin(pi * z) / gamma(1 - z)
+    } else if (z == 0.5) {
+      // Handle this case separately because it is out of the domain of Number.pow when calculating
+      output = 1.7724538509055159
+    } else {
+      z -= 1
+      let mut x = c[0]
+      for (let mut i = 1; i < g + 2; i += 1) {
+        x += c[i] / (z + i)
+      }
+
+      let t = z + g + 0.5
+      output = sqrt(2 * pi) * (t ** (z + 0.5)) * exp(t * -1) * x
+    }
+    if (abs(output) == Infinity) Infinity else output
+  }
+}
+
+/**
+ * Computes the factorial of an integer input or the gamma function of a non-integer input.
+ *
+ * @param n: The value to factorialize
+ * @returns The factorial of the given value
+ *
+ * @throws InvalidArgument(String): When `n` is a negative integer
+ *
+ * @example Number.factorial(0) == 1
+ * @example Number.factorial(3) == 6
+ * @example Number.isClose(Number.factorial(0.5), (1/2) * Number.sqrt(Number.pi))
+ *
+ * @since v0.7.0
+ */
+provide let rec factorial = n => {
+  if (isInteger(n) && n < 0) {
+    gamma(abs(n) + 1) * -1
+  } else if (!isInteger(n) && n < 0) {
+    throw Exception.InvalidArgument(
+      "Cannot compute the factorial of a negative non-integer",
+    )
+  } else {
+    gamma(n + 1)
+  }
+}
diff --git a/stdlib/number.md b/stdlib/number.md
@@ -1659,3 +1659,87 @@ Examples:
 Number.tan(0) == 0
 ```
 
+### Number.**gamma**
+
+<details disabled>
+<summary tabindex="-1">Added in <code>next</code></summary>
+No other changes yet.
+</details>
+
+```grain
+gamma : (z: Number) => Number
+```
+
+Computes the gamma function of a value using the Lanczos approximation.
+
+Parameters:
+
+|param|type|description|
+|-----|----|-----------|
+|`z`|`Number`|The value to interpolate|
+
+Returns:
+
+|type|description|
+|----|-----------|
+|`Number`|The gamma of the given value|
+
+Examples:
+
+```grain
+Number.gamma(1) == 1
+```
+
+```grain
+Number.gamma(3) == 2
+```
+
+```grain
+Number.isClose(Number.gamma(0.5), Number.sqrt(Number.pi))
+```
+
+### Number.**factorial**
+
+<details disabled>
+<summary tabindex="-1">Added in <code>next</code></summary>
+No other changes yet.
+</details>
+
+```grain
+factorial : (n: Number) => Number
+```
+
+Computes the factorial of an integer input or the gamma function of a non-integer input.
+
+Parameters:
+
+|param|type|description|
+|-----|----|-----------|
+|`n`|`Number`|The value to factorialize|
+
+Returns:
+
+|type|description|
+|----|-----------|
+|`Number`|The factorial of the given value|
+
+Throws:
+
+`InvalidArgument(String)`
+
+* When `n` is a negative integer
+
+Examples:
+
+```grain
+Number.factorial(0) == 1
+```
+
+```grain
+Number.factorial(3) == 6
+```
+
+```grain
+Number.isClose(Number.factorial(0.5), (1/2) * Number.sqrt(Number.pi))
+```
+
