Don't like ads? PRO users don't see any ads ;-)
Public Pastes
Guest

Untitled

By: a guest on Feb 26th, 2012  |  syntax: C  |  size: 4.57 KB  |  hits: 23  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4.  
  5. const double EPS1 = 1e-5;
  6. const double EPS2 = 1e-5;
  7.  
  8. typedef struct Slice {
  9.     int n;
  10.     double *keys, *values;
  11. } Slice;
  12.  
  13. double f(int k, double x) {
  14.     switch(k) {
  15.         case 0: return exp(x) + 2;
  16.         case 1: return -1/x;
  17.         case 2: return -2*(x + 1)/3;
  18.         default: return 0;
  19.     }
  20. }
  21.  
  22. double fp(int k, double x) {
  23.      switch(k) {
  24.         case 0: return exp(x);
  25.         case 1: return 1/(x * x);
  26.         case 2: return -2*x/3;
  27.         default: return 0;
  28.     }
  29. }
  30.  
  31. double rootChord(int f1, int f2, double l, double r, double lv, double rv) {
  32.     double m = (l * rv - r * lv) / (rv - lv);
  33.     if (fabs(m - l) < EPS1 || fabs(r - m) < EPS1) {
  34.         return m;
  35.     }
  36.     double mv = f(f1, m) - f(f2, m);
  37.     if ((lv < 0) != (mv < 0)) {
  38.         return rootChord(f1, f2, l, m, lv, mv);
  39.     } else {
  40.         return rootChord(f1, f2, m, r, mv, rv);
  41.     }
  42. }
  43.  
  44. double rootTangent(int f1, int f2, double l, double r, double lv, double rv, double lp, double rp) {
  45.     //fprintf(stderr, "%.4lf, %.4lf\n", l, r);
  46.     //system("pause");
  47.     double m1 = l - lv / lp;
  48.     double m2 = r - rv / rp;
  49.     double m = (l < m1 && m1 < r) ? m1 : m2;
  50.     if (fabs(m - l) < EPS1 || fabs(r - m) < EPS1) {
  51.         return m;
  52.     }
  53.     double mv = f(f1, m) - f(f2, m);
  54.     double mp = fp(f1, m) - fp(f2, m);
  55.     if ((lv < 0) != (mv < 0)) {
  56.         return rootTangent(f1, f2, l, m, lv, mv, lp, mp);
  57.     } else {
  58.         return rootTangent(f1, f2, m, r, mv, rv, mp, rp);
  59.     }
  60. }
  61.  
  62. double root(int f1, int f2, double l, double r) {
  63.     //fprintf(stderr, "Root: %.4lf, %.4lf\n", l, r);
  64.     double lv = f(f1, l) - f(f2, l);
  65.     double rv = f(f1, r) - f(f2, r);
  66.     double lp = fp(f1, l) - fp(f2, l);
  67.     double rp = fp(f1, r) - fp(f2, r);
  68.     //fprintf(stderr, "C\n");
  69.     double root1 = rootChord(f1, f2, l, r, lv, rv);
  70.     //fprintf(stderr, "T\n");
  71.     double root2 = rootTangent(f1, f2, l, r, lv, rv, lp, rp);
  72.     //fprintf(stderr, "OK\n");
  73.     return root1;
  74. }
  75.  
  76. void expandArray(int fn, Slice *s) {
  77.     s->n *= 2;
  78.     s->keys = realloc(s->keys, (s->n + 1) * sizeof(double));
  79.     s->values = realloc(s->values, (s->n + 1) * sizeof(double));
  80.     for (int i = s->n; i >= 2; i -= 2) {
  81.         s->keys[i] = s->keys[i / 2];
  82.         s->values[i] = s->values[i / 2];
  83.     }
  84.     for (int i = 1; i < s->n; i += 2) {
  85.         s->keys[i] = (s->keys[i - 1] + s->keys[i + 1]) / 2.0;
  86.         s->values[i] = f(fn, s->keys[i]);
  87.     }
  88. }
  89.  
  90. Slice newSlice(int fn, double l, double r) {
  91.     Slice res;
  92.     res.n = 1;
  93.     res.keys = malloc(2 * sizeof(double));
  94.     res.values = malloc(2 * sizeof(double));
  95.     res.keys[0] = l;
  96.     res.values[0] = f(fn, l);
  97.     res.keys[1] = r;
  98.     res.values[1] = f(fn, r);
  99.     expandArray(fn, &res);
  100.     return res;
  101. }
  102.  
  103. double integrateArray(Slice s) {
  104.     double res = s.values[0] + s.values[s.n];
  105.     for (int i = 1; i < s.n; ++i) {
  106.         if (i % 2) {
  107.             res += 2 * s.values[i];
  108.         } else {
  109.             res += 4 * s.values[i];
  110.         }
  111.     }
  112.     res *= fabs(s.keys[0] - s.keys[s.n]) / (s.n * 3);
  113.     return res;
  114. }
  115.  
  116. double integrate(int fn, double l, double r) {
  117.     Slice sl = newSlice(fn, l, r);
  118.     double prev, cur = integrateArray(sl);
  119.     do {
  120.         prev = cur;
  121.         expandArray(fn, &sl);
  122.         cur = integrateArray(sl);
  123.     } while (fabs(cur - prev) > EPS2);
  124.     return prev;
  125. }
  126.  
  127. int main(void) {
  128.     double RANGES[3][2];
  129.     RANGES[0][0] = -3.0;
  130.     RANGES[0][1] = -1.0;
  131.     RANGES[1][0] = -5.0;
  132.     RANGES[1][1] = -3.0;
  133.    // fprintf(stderr, "To send: %.4lf, %.4lf\n", RANGES[1][0], RANGES[1][1]);
  134.     RANGES[2][0] = -1.0;
  135.     RANGES[2][1] = -0.1;
  136.     double cuts[3];
  137.     //fprintf(stderr, "---\n");
  138.     cuts[0] = root(1, 2, RANGES[0][0], RANGES[0][1]);
  139.     cuts[1] = root(0, 2, RANGES[1][0], RANGES[1][1]);
  140.     cuts[2] = root(0, 1, RANGES[2][0], RANGES[2][1]);
  141.     //fprintf(stderr, "---\n");
  142.     double ints[3];
  143.     ints[0] = integrate(0, cuts[1], cuts[2]);
  144.     ints[1] = integrate(1, cuts[0], cuts[2]);
  145.     ints[2] = integrate(2, cuts[0], cuts[1]);
  146.     //fprintf(stderr, "---\n");
  147.     double result;
  148.     if ((cuts[0] < cuts[1]) == (cuts[1] < cuts[2])) {
  149.         result = ints[1] - ints[0] - ints[2];
  150.     } else if ((cuts[0] < cuts[2]) == (cuts[2] < cuts[1])) {
  151.         result = ints[2] - ints[0] - ints[1];
  152.     } else {
  153.         result = ints[0] - ints[1] - ints[2];
  154.     }
  155.     //fprintf(stderr, "---\n");
  156.     printf("Result: %.4lf\n", fabs(result));
  157.     printf("Points: %.4lf %.4lf %.4lf\n", cuts[0], cuts[1], cuts[2]);
  158.     return 0;
  159. }
create a new version of this paste RAW Paste Data