#include <stdio.h>
#include <math.h>
// 計算する値の個数
#define N 20
// 円周率
#define M_PI acos(-1.0)
// 微分係数を計算する点
#define x0 (M_PI / 8)
// 真の値 = sqrt(2)
#define d0 sqrt(2)
// 関数 f
double f
(double x
) { return sin(2 * x
); }
// 相対誤差を計算する関数
double err
(double x
) { return fabs(x
- d0
) / d0
; }
// D_f
double dif_f(double x, double h) { return (f(x + h) - f(x)) / h; }
// D_m
double dif_m(double x, double h) { return (f(x + h) - f(x - h)) / (2 * h); }
// 結果は引数として受け取った配列 b に格納する
// 配列 a の長さは L で与えられ、結果は b[0] から b[L-2] までに入る
void accel(double a[], int L, int p, double b[]) {
for (int k = 0; k < L - 1; k++)
b
[k
] = (pow(2, p
) * a
[k
+ 1] - a
[k
]) / (pow(2, p
) - 1);}
// オーダーを計算
// 危険:算術エラーに対する処理を行わないので、サンドボックスでの実験以外で使ってはならない
void order(double a[], double L) {
for (int k = 0; k < L - 2; k++)
printf("%2d: %lf\n", k
, -log2
((a
[k
+ 2] - a
[k
+ 1]) / (a
[k
+ 1] - a
[k
]))); }
// 計算結果を表示
void print_result(double a[], double h[], int L) {
for(int k = 0; k < L; k++)
printf("%2d: %e %.16lf %e\n", k
, h
[k
], a
[k
], err
(a
[k
])); }
int main(void) {
double h[N], a[4][N];
for(int k = 0; k < N; k++) {
// 刻み幅
// D_m の計算値
a[0][k] = dif_m(x0, h[k]);
}
//order(a[0], N);
accel(a[0], N, 2, a[1]);
order(a[1], N - 1);
//print_result(a[1], h, N - 1);
return 0;
}
I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxtYXRoLmg+CgovLyDoqIjnrpfjgZnjgovlgKTjga7lgIvmlbAKI2RlZmluZSBOIDIwCgovLyDlhoblkajnjocKI2RlZmluZSBNX1BJIGFjb3MoLTEuMCkKCi8vIOW+ruWIhuS/guaVsOOCkuioiOeul+OBmeOCi+eCuQojZGVmaW5lIHgwIChNX1BJIC8gOCkKCi8vIOecn+OBruWApCA9IHNxcnQoMikKI2RlZmluZSBkMCBzcXJ0KDIpCgovLyDplqLmlbAgZgpkb3VibGUgZihkb3VibGUgeCkgeyByZXR1cm4gc2luKDIgKiB4KTsgfQoKLy8g55u45a++6Kqk5beu44KS6KiI566X44GZ44KL6Zai5pWwCmRvdWJsZSBlcnIoZG91YmxlIHgpIHsgcmV0dXJuIGZhYnMoeCAtIGQwKSAvIGQwOyB9CgovLyBEX2YKZG91YmxlIGRpZl9mKGRvdWJsZSB4LCBkb3VibGUgaCkgeyByZXR1cm4gKGYoeCArIGgpIC0gZih4KSkgLyBoOyB9CgovLyBEX20KZG91YmxlIGRpZl9tKGRvdWJsZSB4LCBkb3VibGUgaCkgeyByZXR1cm4gKGYoeCArIGgpIC0gZih4IC0gaCkpIC8gKDIgKiBoKTsgfQoKLy8g57WQ5p6c44Gv5byV5pWw44Go44GX44Gm5Y+X44GR5Y+W44Gj44Gf6YWN5YiXIGIg44Gr5qC857SN44GZ44KLCi8vIOmFjeWIlyBhIOOBrumVt+OBleOBryBMIOOBp+S4juOBiOOCieOCjOOAgee1kOaenOOBryBiWzBdIOOBi+OCiSBiW0wtMl0g44G+44Gn44Gr5YWl44KLCnZvaWQgYWNjZWwoZG91YmxlIGFbXSwgaW50IEwsIGludCBwLCBkb3VibGUgYltdKSB7CiAgICBmb3IgKGludCBrID0gMDsgayA8IEwgLSAxOyBrKyspCiAgICAgICAgYltrXSA9IChwb3coMiwgcCkgKiBhW2sgKyAxXSAtIGFba10pIC8gKHBvdygyLCBwKSAtIDEpOwp9CgovLyDjgqrjg7zjg4Djg7zjgpLoqIjnrpcKLy8g5Y2x6Zm677ya566X6KGT44Ko44Op44O844Gr5a++44GZ44KL5Yem55CG44KS6KGM44KP44Gq44GE44Gu44Gn44CB44K144Oz44OJ44Oc44OD44Kv44K544Gn44Gu5a6f6aiT5Lul5aSW44Gn5L2/44Gj44Gm44Gv44Gq44KJ44Gq44GECnZvaWQgb3JkZXIoZG91YmxlIGFbXSwgZG91YmxlIEwpIHsKICAgIGZvciAoaW50IGsgPSAwOyBrIDwgTCAtIDI7IGsrKykKICAgICAgICBwcmludGYoIiUyZDogJWxmXG4iLCBrLCAtbG9nMigoYVtrICsgMl0gLSBhW2sgKyAxXSkgLyAoYVtrICsgMV0gLSBhW2tdKSkpOwp9CgovLyDoqIjnrpfntZDmnpzjgpLooajnpLoKdm9pZCBwcmludF9yZXN1bHQoZG91YmxlIGFbXSwgZG91YmxlIGhbXSwgaW50IEwpIHsKICAgIGZvcihpbnQgayA9IDA7IGsgPCBMOyBrKyspCiAgICAgICAgcHJpbnRmKCIlMmQ6ICVlICAlLjE2bGYgICVlXG4iLCBrLCBoW2tdLCBhW2tdLCBlcnIoYVtrXSkpOwp9CgppbnQgbWFpbih2b2lkKSB7CiAgICBkb3VibGUgaFtOXSwgYVs0XVtOXTsKCiAgICBmb3IoaW50IGsgPSAwOyBrIDwgTjsgaysrKSB7CiAgICAgICAgLy8g5Yi744G/5bmFCiAgICAgICAgaFtrXSA9IHBvdygyLCAtayAtIDMpOwoKICAgICAgICAvLyBEX20g44Gu6KiI566X5YCkCiAgICAgICAgYVswXVtrXSA9IGRpZl9tKHgwLCBoW2tdKTsKICAgIH0KCiAgICAvL29yZGVyKGFbMF0sIE4pOwogICAgYWNjZWwoYVswXSwgTiwgMiwgYVsxXSk7CgogICAgb3JkZXIoYVsxXSwgTiAtIDEpOwogICAgLy9wcmludF9yZXN1bHQoYVsxXSwgaCwgTiAtIDEpOwoKICAgIHJldHVybiAwOwp9