1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
|
/*
* ns.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;
}
|
