diff options
Diffstat (limited to '')
| -rw-r--r-- | buch/chapters/110-elliptisch/agm/sn.cpp | 52 |
1 files changed, 52 insertions, 0 deletions
diff --git a/buch/chapters/110-elliptisch/agm/sn.cpp b/buch/chapters/110-elliptisch/agm/sn.cpp new file mode 100644 index 0000000..9e1b047 --- /dev/null +++ b/buch/chapters/110-elliptisch/agm/sn.cpp @@ -0,0 +1,52 @@ +/* + * sn.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> +#include <gsl/gsl_sf_elljac.h> + +static const int N = 10; + +inline long double sqrl(long double x) { + return x * x; +} + +int main(int argc, char *argv[]) { + long double u = 0.6; + long double k = 0.9; + long double kprime = sqrt(1 - sqrl(k)); + + long double a[N], b[N], x[N+1]; + a[0] = 1; + b[0] = kprime; + + for (int n = 0; n < N-1; n++) { + printf("a[%d] = %22.18Lf b[%d] = %22.18Lf\n", n, a[n], n, b[n]); + a[n+1] = (a[n] + b[n]) / 2; + b[n+1] = sqrtl(a[n] * b[n]); + } + + x[N] = sinl(u * a[N-1]); + printf("x[%d] = %22.18Lf\n", N, x[N]); + + for (int n = N - 1; n >= 0; n--) { + x[n] = 2 * a[n] * x[n+1] / (a[n] + b[n] + (a[n] - b[n]) * sqrl(x[n+1])); + printf("x[%2d] = %22.18Lf\n", n, x[n]); + } + + printf("sn(%7.4Lf, %7.4Lf) = %20.24Lf\n", u, k, x[0]); + + double sn, cn, dn; + double m = sqrl(k); + gsl_sf_elljac_e((double)u, m, &sn, &cn, &dn); + printf("sn(%7.4Lf, %7.4Lf) = %20.24f\n", u, k, sn); + printf("cn(%7.4Lf, %7.4Lf) = %20.24f\n", u, k, cn); + printf("dn(%7.4Lf, %7.4Lf) = %20.24f\n", u, k, dn); + + return EXIT_SUCCESS; +} |