Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- //#include <bits/stdc++.h>
- #include <iostream>
- #include <iomanip>
- using namespace std;
- typedef long long ll;
- typedef unsigned long long ull;
- #define all(x) x.begin(), x.end()
- #define rall(x) x.rbegin(), x.rend()
- #define endl '\n'
- #define boostIO() ios_base::sync_with_stdio(0);cin.tie(0);cout.tie(0);
- ll gcd(ll a, ll b) { return (b == 0 ? a : gcd(b, a % b)); }
- double f(double x, double y) {
- return x * x * x - 3 * x * y * y - x * x + y * y + x - 1;
- }
- double g(double x, double y) {
- return 3 * x * x * y - y * y * y - 2 * x * y + y;
- }
- double eps = 1e-9;
- double step = 0.15;
- //double Xm = -10, Xp = 10;
- //double Ym = -10, Yp = 10;
- //double Xm = -10, Xp = 10;
- //double Ym = -10, Yp = 0.5;
- double Xm = -10, Xp = 0.5;
- double Ym = -10, Yp = 0.5;
- double x = 0;
- double y = 0;
- double x0 = 0;
- double Y0 = 0;
- double x1 = 0;
- double Y1 = 0;
- double x2 = 0;
- double y2 = 0;
- void initial_points0() {
- double dl = 0.1;
- double curmin = 1000000.;
- while (curmin > 999999) {
- for (x = Xm; x < Xp; x = x + step)
- for (y = Ym; y < Yp; y = y + step) {
- double fc = f(x, y);
- double gc = g(x, y);
- if ((fc > eps) && (gc > eps)) {
- double dv = fabs(fc) + fabs(gc);
- if (dv > curmin) continue;
- x0 = x;
- Y0 = y;
- curmin = dv;
- }
- }
- dl = dl / 2;
- }
- cout << "x0 = " << x0 << endl;
- cout << "Y0 = " << Y0 << endl;
- }
- void initial_points1() {
- double dl = 0.1;
- double curmin = 1000000.;
- while (curmin > 999999) {
- for (x = Xm; x < Xp; x = x + step)
- for (y = Ym; y < Yp; y = y + step) {
- double fc = f(x, y);
- double gc = g(x, y);
- if ((fc < eps) && (gc < eps)) {
- double dv = fabs(fc) + fabs(gc);
- if (dv > curmin) continue;
- x1 = x;
- Y1 = y;
- curmin = dv;
- }
- }
- dl = dl / 2;
- }
- cout << "x1 = " << x1 << endl;
- cout << "Y1 = " << Y1 << endl;
- }
- void initial_points2() {
- double dl = 0.1;
- double curmin = 1000000.;
- while (curmin > 999999) {
- for (x = Xm; x < Xp; x = x + step)
- for (y = Ym; y < Yp; y = y + step) {
- double fc = f(x, y);
- double gc = g(x, y);
- if ((fc > eps) && (gc < eps)) {
- double dv = fabs(fc) + fabs(gc);
- if (dv > curmin) continue;
- x2 = x;
- y2 = y;
- curmin = dv;
- }
- }
- dl = dl / 2;
- }
- //cout << "x0 = " << x0 << endl;
- //cout << "Y0 = " << Y0 << endl;
- //cout << "x1 = " << x1 << endl;
- //cout << "Y1 = " << Y1 << endl;
- cout << "x2 = " << x2 << endl;
- cout << "y2 = " << y2 << endl;
- }
- double stri(double x0_, double Y0_, double x1_, double Y1_, double x2_, double y2_) {
- return (0.5 * (x0_ * (Y1_ - y2_) + x1_ * (y2_ - Y0_) + x2_ * (Y0_ - Y1_)));
- }
- signed main() {
- cout << fixed << setprecision(6);
- initial_points0();
- cout << endl;
- initial_points1();
- cout << endl;
- initial_points2();
- cout << endl;
- //defining of the last point before the first execution of the main cycle; any values far from the initial three
- //points:
- double xls = 1000000.;
- double yls = 1000000.;
- double dv = 1;
- double x, y;
- while (dv > 1e-9) //main cycle corresponding to steps 1 β 3 of the brief description
- {
- // points on two surfaces calculation:
- double zf0 = f(x0, Y0);
- double zg0 = g(x0, Y0);
- double zf1 = f(x1, Y1);
- double zg1 = g(x1, Y1);
- double zf2 = f(x2, y2);
- double zg2 = g(x2, y2);
- // the first plane coefficients: AX + BY + CZ + D = 0
- double Af = (Y1 - Y0) * (zf2 - zf0) - (y2 - Y0) * (zf1 - zf0);
- double Bf = (x2 - x0) * (zf1 - zf0) - (x1 - x0) * (zf2 - zf0);
- double Cf = (x1 - x0) * (y2 - Y0) - (x2 - x0) * (Y1 - Y0);
- double Df = -x0 * Af - Y0 * Bf - zf0 * Cf;
- // the second plane coefficients:
- double Ag = (Y1 - Y0) * (zg2 - zg0) - (y2 - Y0) * (zg1 - zg0);
- double Bg = (x2 - x0) * (zg1 - zg0) - (x1 - x0) * (zg2 - zg0);
- double Cg = (x1 - x0) * (y2 - Y0) - (x2 - x0) * (Y1 - Y0);
- double Dg = -x0 * Ag - Y0 * Bg - zg0 * Cg;
- // the stability checking:
- double dt = Af * Bg - Ag * Bf;
- if (fabs(dt) < 0.000001) break;
- // three planes intersection (z=0 by the construction):
- x = (-Df * Bg + Dg * Bf) / dt;
- y = (-Af * Dg + Ag * Df) / dt;
- // value of both functions in the new point:
- double zfc = f(x, y);
- double zgc = g(x, y);
- // proximity calculation (the current new point (x,y) is compared with the last one (xls,yls):
- dv = fabs(x - xls) + fabs(y - yls);
- // βnewβ last point:
- xls = x;
- yls = y;
- // square of triangles corresponding to the discarding of one of the old points are calculated (function stri):
- double s0 = stri(x1, Y1, x2, y2, x, y);
- double s1 = stri(x0, Y0, x2, y2, x, y);
- double s2 = stri(x0, Y0, x1, Y1, x, y);
- // choice of the maximal triangle and corresponding denotation of considered points (x0,Y0), (x1,Y1), //(x2,y2) :
- if ((s0 >= s1) && (s0 >= s2)){
- x0 = x1;
- Y0 = Y1;
- x1 = x2;
- Y1 = y2;
- x2 = x;
- y2 = y;
- continue;
- }
- if ((s1 >= s0) && (s1 >= s2)){
- x1 = x2;
- Y1 = y2;
- x2 = x;
- y2 = y;
- continue;
- }
- if ((s2 >= s0) && (s2 >= s1)){
- x2 = x;
- y2 = y;
- continue;
- }
- }//end of the main cycle
- cout << "x = " << xls << endl;
- cout << "y = " << yls << endl;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement