aboutsummaryrefslogtreecommitdiffstats
path: root/buch/chapters/110-elliptisch/agm
diff options
context:
space:
mode:
Diffstat (limited to '')
-rw-r--r--buch/chapters/110-elliptisch/agm/Makefile10
-rw-r--r--buch/chapters/110-elliptisch/agm/agm.cpp42
-rw-r--r--buch/chapters/110-elliptisch/agm/agm.m (renamed from buch/chapters/110-elliptisch/agm.m)2
3 files changed, 53 insertions, 1 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.m b/buch/chapters/110-elliptisch/agm/agm.m
index 2f0a1ea..dcb3ad8 100644
--- a/buch/chapters/110-elliptisch/agm.m
+++ b/buch/chapters/110-elliptisch/agm/agm.m
@@ -16,5 +16,5 @@ for i = (1:n)
a = A;
end
-E = 2 / (pi * a)
+K = pi / (2 * a)