Advertisement
Georgiy031

Untitled

Oct 21st, 2020
123
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 4.95 KB | None | 0 0
  1. //#include <bits/stdc++.h>
  2. #include <iostream>
  3. #include <iomanip>
  4. using namespace std;
  5. typedef long long ll;
  6. typedef unsigned long long ull;
  7. #define all(x) x.begin(), x.end()
  8. #define rall(x) x.rbegin(), x.rend()
  9. #define endl '\n'
  10. #define boostIO() ios_base::sync_with_stdio(0);cin.tie(0);cout.tie(0);
  11. ll gcd(ll a, ll b) { return (b == 0 ? a : gcd(b, a % b)); }
  12.  
  13. double f(double x, double y) {
  14.     return x * x * x - 3 * x * y * y - x * x + y * y + x - 1;
  15. }
  16. double g(double x, double y) {
  17.     return 3 * x * x * y - y * y * y - 2 * x * y + y;
  18. }
  19.  
  20. double eps = 1e-9;
  21. double step = 0.15;
  22.  
  23. //double Xm = -10, Xp = 10;
  24. //double Ym = -10, Yp = 10;
  25.  
  26. //double Xm = -10, Xp = 10;
  27. //double Ym = -10, Yp = 0.5;
  28.  
  29. double Xm = -10, Xp = 0.5;
  30. double Ym = -10, Yp = 0.5;
  31.  
  32. double x = 0;
  33. double y = 0;
  34. double x0 = 0;
  35. double Y0 = 0;
  36. double x1 = 0;
  37. double Y1 = 0;
  38. double x2 = 0;
  39. double y2 = 0;
  40.  
  41. void initial_points0() {
  42.     double dl = 0.1;
  43.     double curmin = 1000000.;
  44.  
  45.     while (curmin > 999999) {
  46.         for (x = Xm; x < Xp; x = x + step)
  47.             for (y = Ym; y < Yp; y = y + step) {
  48.  
  49.                 double fc = f(x, y);
  50.                 double gc = g(x, y);
  51.  
  52.                 if ((fc > eps) && (gc > eps)) {
  53.                     double dv = fabs(fc) + fabs(gc);
  54.  
  55.                     if (dv > curmin) continue;
  56.  
  57.                     x0 = x;
  58.                     Y0 = y;
  59.  
  60.                     curmin = dv;
  61.                 }
  62.             }
  63.         dl = dl / 2;
  64.     }
  65.     cout << "x0 = " << x0 << endl;
  66.     cout << "Y0 = " << Y0 << endl;
  67. }
  68.  
  69. void initial_points1() {
  70.     double dl = 0.1;
  71.  
  72.     double curmin = 1000000.;
  73.  
  74.     while (curmin > 999999) {
  75.         for (x = Xm; x < Xp; x = x + step)
  76.             for (y = Ym; y < Yp; y = y + step) {
  77.  
  78.                 double fc = f(x, y);
  79.                 double gc = g(x, y);
  80.  
  81.                 if ((fc < eps) && (gc < eps)) {
  82.  
  83.                     double dv = fabs(fc) + fabs(gc);
  84.  
  85.                     if (dv > curmin) continue;
  86.  
  87.                     x1 = x;
  88.                     Y1 = y;
  89.  
  90.                     curmin = dv;
  91.                 }
  92.             }
  93.         dl = dl / 2;
  94.     }
  95.     cout << "x1 = " << x1 << endl;
  96.     cout << "Y1 = " << Y1 << endl;
  97. }
  98. void initial_points2() {
  99.  
  100.     double dl = 0.1;
  101.  
  102.     double curmin = 1000000.;
  103.  
  104.     while (curmin > 999999) {
  105.         for (x = Xm; x < Xp; x = x + step)
  106.             for (y = Ym; y < Yp; y = y + step) {
  107.  
  108.                 double fc = f(x, y);
  109.                 double gc = g(x, y);
  110.  
  111.                 if ((fc > eps) && (gc < eps)) {
  112.  
  113.                     double dv = fabs(fc) + fabs(gc);
  114.  
  115.                     if (dv > curmin) continue;
  116.  
  117.                     x2 = x;
  118.                     y2 = y;
  119.  
  120.                     curmin = dv;
  121.                 }
  122.             }
  123.         dl = dl / 2;
  124.     }
  125.     //cout << "x0 = " << x0 << endl;
  126.     //cout << "Y0 = " << Y0 << endl;
  127.     //cout << "x1 = " << x1 << endl;
  128.     //cout << "Y1 = " << Y1 << endl;
  129.     cout << "x2 = " << x2 << endl;
  130.     cout << "y2 = " << y2 << endl;
  131. }
  132.  
  133. double stri(double x0_, double Y0_, double x1_, double Y1_, double x2_, double y2_) {
  134.     return (0.5 * (x0_ * (Y1_ - y2_) + x1_ * (y2_ - Y0_) + x2_ * (Y0_ - Y1_)));
  135. }
  136. signed main() {
  137.     cout << fixed << setprecision(6);
  138.     initial_points0();
  139.     cout << endl;
  140.     initial_points1();
  141.     cout << endl;
  142.     initial_points2();
  143.     cout << endl;
  144.  
  145.      //defining of the last point before the first execution of the main cycle; any values far from the initial three
  146.      //points:  
  147.  
  148.     double xls = 1000000.;
  149.     double yls = 1000000.;
  150.     double dv = 1;
  151.     double x, y;
  152.  
  153.     while (dv > 1e-9) //main cycle corresponding to steps 1 – 3 of the brief description
  154.     {
  155.         // points on two surfaces calculation:
  156.         double zf0 = f(x0, Y0);
  157.         double zg0 = g(x0, Y0);
  158.         double zf1 = f(x1, Y1);
  159.         double zg1 = g(x1, Y1);
  160.         double zf2 = f(x2, y2);
  161.         double zg2 = g(x2, y2);
  162.         // the first plane coefficients: AX + BY + CZ + D = 0
  163.         double Af = (Y1 - Y0) * (zf2 - zf0) - (y2 - Y0) * (zf1 - zf0);
  164.         double Bf = (x2 - x0) * (zf1 - zf0) - (x1 - x0) * (zf2 - zf0);
  165.         double Cf = (x1 - x0) * (y2 - Y0) - (x2 - x0) * (Y1 - Y0);
  166.         double Df = -x0 * Af - Y0 * Bf - zf0 * Cf;
  167.         // the second plane coefficients:
  168.         double Ag = (Y1 - Y0) * (zg2 - zg0) - (y2 - Y0) * (zg1 - zg0);
  169.         double Bg = (x2 - x0) * (zg1 - zg0) - (x1 - x0) * (zg2 - zg0);
  170.         double Cg = (x1 - x0) * (y2 - Y0) - (x2 - x0) * (Y1 - Y0);
  171.         double Dg = -x0 * Ag - Y0 * Bg - zg0 * Cg;
  172.             // the stability checking:
  173.         double dt = Af * Bg - Ag * Bf;
  174.  
  175.         if (fabs(dt) < 0.000001) break;
  176.         // three planes intersection (z=0 by the construction):
  177.         x = (-Df * Bg + Dg * Bf) / dt;
  178.         y = (-Af * Dg + Ag * Df) / dt;
  179.         // value of both functions in the new point:
  180.         double zfc = f(x, y);
  181.         double zgc = g(x, y);
  182.         // proximity calculation (the current new point (x,y) is compared with the last one (xls,yls):
  183.         dv = fabs(x - xls) + fabs(y - yls);
  184.         // β€œnew” last point:
  185.         xls = x;
  186.         yls = y;
  187.         // square of triangles corresponding to the discarding of one of the old points are calculated (function stri):
  188.         double s0 = stri(x1, Y1, x2, y2, x, y);
  189.         double s1 = stri(x0, Y0, x2, y2, x, y);
  190.         double s2 = stri(x0, Y0, x1, Y1, x, y);
  191.         // choice of the maximal triangle and corresponding denotation of considered points (x0,Y0), (x1,Y1), //(x2,y2) :
  192.         if ((s0 >= s1) && (s0 >= s2)){
  193.             x0 = x1;
  194.             Y0 = Y1;
  195.             x1 = x2;
  196.             Y1 = y2;
  197.             x2 = x;
  198.             y2 = y;
  199.             continue;
  200.         }
  201.         if ((s1 >= s0) && (s1 >= s2)){
  202.             x1 = x2;
  203.             Y1 = y2;
  204.             x2 = x;
  205.             y2 = y;
  206.             continue;
  207.         }
  208.         if ((s2 >= s0) && (s2 >= s1)){
  209.             x2 = x;
  210.             y2 = y;
  211.             continue;
  212.         }
  213.     }//end of the main cycle
  214.     cout << "x = " << xls << endl;
  215.     cout << "y = " << yls << endl;
  216. }
  217.  
  218.  
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement