crepl

commit 5ec6c6b393a59ae9c4391c40da716c1451bf63f8
parent 859a5f67c167d89987dcf83f628748a8fc111880
Author: Demonstrandum <moi@knutsen.co>
Date:   Mon, 29 Jun 2020 16:52:17 +0100

Fix factorial fiasco.

Diffstat:
M
README.md
 | 
1
-
M
src/builtin.c
 | 
67
+++++++++++++++++++++++++++++--------------------------------------
M
src/builtin.h
 | 
2
++
M
src/defaults.h
 | 
7
+++++++

4 files changed, 38 insertions(+), 39 deletions(-)
diff --git a/README.md b/README.md
@@ -30,7 +30,6 @@ sudo make install  # Installs the program system wide.
  - [ ] Throw errors on overflows until we implement bignums.
  - [ ] Imaginary numbers (using `complex.h`).
  - [ ] User defined functions.
- - [ ] Extend factorial to positive reals and complex values using Gamma function.
  - [ ] Garbage collection.
  - [ ] Numerical equation solver (polynomial, simultaneous, &c.).
  - [ ] Add more functionality, notably for calculus.
diff --git a/src/builtin.c b/src/builtin.c
@@ -105,22 +105,34 @@ fsize nice_sin(fsize alpha)
 	return sin(alpha);
 }
 
-#define gamma_upper_bound 1000
-#define gamma_resolution 1000000
-#define gamma_h gamma_upper_bound/gamma_resolution  //upper bound / resolution
+fsize gamma_complete(fsize num) {
+	if (num == 0)
+		return INF;
 
-fsize gamma_integrad(float x, fsize num) {
-	return pow(x, num) * exp(-x);
-}
+	fsize integral = 0;
+	fsize fractional = modfl(num, &integral);
 
-fsize gammae(fsize num) {
-	fsize sum = 0;
+	if (fractional == 0) {
+		if (num < 0)
+			return INF;
+		fsize res = 1;
+		for (fsize i = num - 1; i > 1; --i)
+			res *= i;
+		return res;
+	}
 
-	for (int i = 1; i < gamma_resolution; i++) {
-    	sum += gamma_integrad(i*gamma_h, num);
-    }
+	if (num < 0) {
+		// Gamma(x) = (1/x) * Gamma(x + 1)
+		fsize acc = 1 / num;
+		fsize i = num + 1;
+		do {
+			acc *= 1 / i;
+			i += 1;
+		} while (i < 0);
 
-    return 0.5 * gamma_h * (gamma_integrad(0, num) + gamma_integrad(gamma_upper_bound, num) + 2 * sum);
+		return acc * tgammal(i);
+	}
+	return tgammal(num);
 }
 
 MATH_WRAPPER(sin, nice_sin)
@@ -145,6 +157,7 @@ MATH_WRAPPER(atanh, atanhl)
 // TODO: atan2, hypot
 MATH_WRAPPER(ceil, ceil)
 MATH_WRAPPER(floor, floor)
+MATH_WRAPPER(Gamma, gamma_complete)
 
 DataValue *builtin_neg(DataValue input)
 {
@@ -176,37 +189,15 @@ DataValue *builtin_factorial(DataValue input)
 {
 	NumberNode *num = type_check("!", LHS, T_NUMBER, &input);
 
+	if (num == NULL)
+		return NULL;
+
 	NumberNode tmp = num_to_float(*num);
 	NumberNode *new_num = malloc(sizeof(NumberNode));
 	memcpy(new_num, &tmp, sizeof(NumberNode));
+	new_num->value.f = gamma_complete(new_num->value.f + 1);
 
 	DataValue *result = wrap_data(T_NUMBER, new_num);
-	if (num == NULL)
-		return NULL;
-
-	fsize integral = 0;
-	fsize fractional = modfl(num->value.f, &integral);
-	if ((num->type == FLOAT && (num->value.f < 0))
-	||  (num->type == INT && num->value.i < 0)) {
-		ERROR_TYPE = EXECUTION_ERROR;
-		strcpy(ERROR_MSG,
-			"factorial (`!') is only defined for positve real numbers.");
-	} else if ((num->type == FLOAT && (fractional != 0))) {
-		new_num->value.f = gammae(num->value.f);
-		result->value = new_num;
-		return result;
-		}
-
-	if (new_num->value.f == 0) {
-		new_num->value.f = 1;
-		result->value = new_num;
-		return result;
-	}
-	ssize i = new_num->value.f - 1;
-	while (i > 1) {
-		new_num->value.f *= i;
-		--i;
-	}
 	result->value = new_num;
 	return result;
 }
diff --git a/src/builtin.h b/src/builtin.h
@@ -36,6 +36,7 @@ DataValue *builtin_ceil(DataValue);
 DataValue *builtin_floor(DataValue);
 DataValue *builtin_factorial(DataValue);
 DataValue *builtin_neg(DataValue);
+DataValue *builtin_Gamma(DataValue);
 
 NumberNode *num_add(NumberNode, NumberNode);
 NumberNode *num_sub(NumberNode, NumberNode);
@@ -83,5 +84,6 @@ static const struct _func_name_pair builtin_fns[] = {
 	{ "!", { builtin_factorial } },
 	FUNC_PAIR(neg),
 	{ "-", { builtin_neg } },
+	FUNC_PAIR(Gamma),
 };
 
diff --git a/src/defaults.h b/src/defaults.h
@@ -7,11 +7,18 @@
 #include <stdbool.h>
 #include <string.h>
 #include <ctype.h>
+#include <math.h>
 
 #include "error.h"
 
 #define len(array) (sizeof(array) / sizeof((array)[0]))
 
+#ifdef INFINITY
+	#define INF INFINITY
+#else
+	#define INF (1.0 / 0.0)
+#endif
+
 typedef uint8_t u8 ;
 typedef uint16_t u16;
 typedef uint32_t u32;