add numerical computation of lambert W function

2 files changed, 110 insertions, 0 deletions

diff --git a/buch/chapters/020-exponential/code/Makefile b/buch/chapters/020-exponential/code/Makefile
index b189dae..f7c1e9c 100644
--- a/buch/chapters/020-exponential/code/Makefile
+++ b/buch/chapters/020-exponential/code/Makefile
@@ -4,6 +4,14 @@
 # 
 # (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
 #
+all:	lambertw-test
 
 xxl:	xxl.c
 	gcc -o xxl -W -O2 -lgsl xxl.c
+
+lambertw:	lambertw.c
+	gcc -o lambertw -W -O2 lambertw.c
+
+lambertw-test:	lambertw
+	./lambertw -l l.txt -- -0.35 -0.367 -0.3678 -0.36787 -0.367879 -0.2 1 2 3 10 100 1000 -0.1 -0.01 -0.001 -0.0001
+
diff --git a/buch/chapters/020-exponential/code/lambertw.c b/buch/chapters/020-exponential/code/lambertw.c
new file mode 100644
index 0000000..14ee190
--- /dev/null
+++ b/buch/chapters/020-exponential/code/lambertw.c
@@ -0,0 +1,102 @@
+/*
+ * lambertw.c
+ *
+ * (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
+ */
+#include <stdlib.h>
+#include <stdio.h>
+#include <getopt.h>
+#include <math.h>
+
+static double	precision = 1e-16;
+static FILE	*logfile = NULL;
+
+double	f(double x) {
+	return x * exp(x);
+}
+
+double	lambertW0(double y) {
+	double	x = log(1 + y);
+	double	delta = 1;
+	int	count = 0;
+	while ((count < 20) && (delta > precision)) {
+		double	error = f(x) - y;
+		if (logfile) {
+			fprintf(logfile, "%20.16f    %20.16f\n", x, error);
+		}
+		double	xnew = x - error / ((1 + x) * exp(x));
+		delta = fabs(xnew - x);
+		x = xnew;
+		count++;
+	}
+	if (count == 0) {
+		fprintf(stderr, "count exceeded\n");
+		return -1;
+	}
+	if (logfile) {
+		fprintf(logfile, "%d\n", count);
+	}
+	return x;
+}
+
+double	lambertW1(double y) {
+	double	x = log(-y);
+	double	delta = 1;
+	int	count = 0;
+	while ((count < 20) && (delta > precision)) {
+		double	error = f(x) - y;
+		if (logfile) {
+			fprintf(logfile, "%20.16f    %20.16f\n", x, error);
+		}
+		double	xnew = x - error / ((1 + x) * exp(x));
+		delta = fabs(xnew - x);
+		x = xnew;
+		count++;
+	}
+	if (count == 0) {
+		fprintf(stderr, "count exceeded\n");
+		return -1;
+	}
+	if (logfile) {
+		fprintf(logfile, "%d\n", count);
+	}
+	return x;
+}
+
+struct option	options[] = {
+{	"precision",	required_argument,	NULL,	'p'	},
+{	"logfile",	required_argument,	NULL,	'l'	},
+{	"minus",	no_argument,		NULL,	'm'	},
+{	NULL,		0,			NULL,    0	}
+};
+
+int	main(int argc, char *argv[]) {
+	int	c;
+	int	l;
+	int	minusone = 0;
+	while (EOF != (c = getopt_long(argc, argv, "p:l:m", options, &l)))
+		switch (c) {
+		case 'l':
+			logfile = fopen(optarg, "w");
+			break;
+		case 'p':
+			precision = atof(optarg);
+			break;
+		case 'm':
+			minusone = 1;
+			break;
+		}
+
+	for (; optind < argc; optind++) {
+		double	y = atof(argv[optind]);
+		double	x = (y >= 0) ? lambertW0(y) :
+				((minusone) ? lambertW1(y) : lambertW0(y));
+		printf("%20.16f -> %20.16f  (%20.16f)\n", x, y, y - f(x));
+	}
+
+	if (logfile) {
+		fclose(logfile);
+	}
+
+	return EXIT_SUCCESS;
+}