fix agm

3 files changed, 72 insertions, 0 deletions

diff --git a/buch/chapters/110-elliptisch/agm/Makefile b/buch/chapters/110-elliptisch/agm/Makefile
new file mode 100644
index 0000000..e7975e1
--- /dev/null
+++ b/buch/chapters/110-elliptisch/agm/Makefile
@@ -0,0 +1,10 @@
+#
+# Makefile
+#
+# (c) 2022 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
+#
+
+agm:	agm.cpp
+	g++ -O -Wall -g -std=c++11 agm.cpp -o agm `pkg-config --cflags gsl` `pkg-config --libs gsl`
+	./agm
+
diff --git a/buch/chapters/110-elliptisch/agm/agm.cpp b/buch/chapters/110-elliptisch/agm/agm.cpp
new file mode 100644
index 0000000..fdb0441
--- /dev/null
+++ b/buch/chapters/110-elliptisch/agm/agm.cpp
@@ -0,0 +1,42 @@
+/*
+ * agm.cpp
+ *
+ * (c) 2022 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
+ */
+#include <cstdlib>
+#include <cstdio>
+#include <cmath>
+#include <iostream>
+#include <gsl/gsl_sf_ellint.h>
+
+
+
+int	main(int argc, char *argv[]) {
+	long double	a = 1;
+	long double	b = sqrtl(2.)/2;
+	if (argc >= 3) {
+		a = std::stod(argv[1]);
+		b = std::stod(argv[2]);
+	}
+
+	{
+		long double	an = a;
+		long double	bn = b;
+		for (int i = 0; i < 10; i++) {
+			printf("%d  %24.18Lf  %24.18Lf  %24.18Lf\n",
+				i, an, bn, a * M_PI / (2 * an));
+			long double A = (an + bn) / 2;
+			bn = sqrtl(an * bn);
+			an = A;
+		}
+	}
+
+	{
+		double	k = b/a;
+		k = sqrt(1 - k*k);
+		double	K = gsl_sf_ellint_Kcomp(k, GSL_PREC_DOUBLE);
+		printf("                             %24.18f  %24.18f\n", k, K);
+	}
+
+	return EXIT_SUCCESS;
+}
diff --git a/buch/chapters/110-elliptisch/agm/agm.m b/buch/chapters/110-elliptisch/agm/agm.m
new file mode 100644
index 0000000..dcb3ad8
--- /dev/null
+++ b/buch/chapters/110-elliptisch/agm/agm.m
@@ -0,0 +1,20 @@
+#
+# agm.m
+#
+# (c) 2022 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
+#
+format long
+
+n = 10;
+a = 1;
+b = sqrt(0.5);
+
+for i = (1:n)
+	printf("%20.16f %20.16f\n", a, b);
+	A = (a+b)/2;
+	b = sqrt(a*b);
+	a = A;
+end	
+
+K = pi / (2 * a)
+