fork(1) download
copy 
  1. #include <stdio.h>
  2. #include <math.h>
  3.  
  4. // 計算する値の個数
  5. #define N 20
  6.  
  7. // 円周率
  8. #define M_PI acos(-1.0)
  9.  
  10. // 微分係数を計算する点
  11. #define x0 (M_PI / 8)
  12.  
  13. // 真の値 = sqrt(2)
  14. #define d0 sqrt(2)
  15.  
  16. // 関数 f
  17. double f(double x) { return sin(2 * x); }
  18.  
  19. // 相対誤差を計算する関数
  20. double err(double x) { return fabs(x - d0) / d0; }
  21.  
  22. // D_f
  23. double dif_f(double x, double h) { return (f(x + h) - f(x)) / h; }
  24.  
  25. // D_m
  26. double dif_m(double x, double h) { return (f(x + h) - f(x - h)) / (2 * h); }
  27.  
  28. // 結果は引数として受け取った配列 b に格納する
  29. // 配列 a の長さは L で与えられ、結果は b[0] から b[L-2] までに入る
  30. void accel(double a[], int L, int p, double b[]) {
  31.     for (int k = 0; k < L - 1; k++)
  32.         b[k] = (pow(2, p) * a[k + 1] - a[k]) / (pow(2, p) - 1);
  33. }
  34.  
  35. // オーダーを計算
  36. // 危険：算術エラーに対する処理を行わないので、サンドボックスでの実験以外で使ってはならない
  37. void order(double a[], double L) {
  38.     for (int k = 0; k < L - 2; k++)
  39.         printf("%2d: %lf\n", k, -log2((a[k + 2] - a[k + 1]) / (a[k + 1] - a[k])));
  40. }
  41.  
  42. // 計算結果を表示
  43. void print_result(double a[], double h[], int L) {
  44.     for(int k = 0; k < L; k++)
  45.         printf("%2d: %e  %.16lf  %e\n", k, h[k], a[k], err(a[k]));
  46. }
  47.  
  48. int main(void) {
  49.     double h[N], a[4][N];
  50.  
  51.     for(int k = 0; k < N; k++) {
  52.         // 刻み幅
  53.         h[k] = pow(2, -k - 3);
  54.  
  55.         // D_m の計算値
  56.         a[0][k] = dif_m(x0, h[k]);
  57.     }
  58.  
  59.     //order(a[0], N);
  60.     accel(a[0], N, 2, a[1]);
  61.  
  62.     order(a[1], N - 1);
  63.     //print_result(a[1], h, N - 1);
  64.  
  65.     return 0;
  66. }
Success #stdin #stdout 0s 4556KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
 0: 3.997887
 1: 3.999472
 2: 3.999877
 3: 3.999706
 4: 4.001734
 5: 4.060754
 6: inf
 7: -inf
 8: -nan
 9: 0.320660
10: -nan
11: -nan
12: -nan
13: -nan
14: -nan
15: -nan
16: inf
https://ideone.com/ukFX89
language:
C (gcc 8.3)
created:
6 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!