math_mod_lab1

#include "stdafx.h"//vs2010

#include <bits\stdc++.h>
#include <iostream>
#include <complex>
#include <cassert>
#include <vector>
#include <conio.h>
#include <windows.h>
//#include "decimal.h"
#include "gnuplot_i.hpp"

using namespace std;
//using namespace dec;

#define M_PI 3.14159265358979323846
#define FOUT
#define cyrillic setlocale(LC_CTYPE, "rus")
#define eps 1e-2

inline string dtos(double d){ostringstream strs; strs << d; return strs.str();}

void wait(){
    FlushConsoleInputBuffer(GetStdHandle(STD_INPUT_HANDLE));
    #ifndef FOUT
        cout << endl << "Press any key to continue..." << endl;
    #endif
    _getch();
}


string parse_path(char* path){
    string s = path;
    size_t pos = s.find_last_of('\\');
    if(pos != string::npos){
        return s.substr(0 , pos + 1);
    }
    if(s[s.size() - 1] != '\\')s += '\\';
    return s;
}

inline void concat(vector<char> & v , string s){
    v.pop_back();
    for(int i = 0 ; i < s.size() ; ++i){
        v.push_back(s[i]);
    }
    v.push_back('\0');
}

vector<double> gauss(vector<vector<double> > & A , vector<double> & B){//Gauss
    int n = 6;
    vector<double> X(n);
    for(int j = 0 ; j < n ; ++j){
        for(int i = j + 1 ; i < n ; ++i){
            double k = A[i][j] / A[j][j];
            for(int r = 0 ; r < n ; ++r){
                A[i][r] -= k * A[j][r];
            }
            B[i] -= k * B[j];
        }
    }

    for(int i = n - 1 ; i >= 0 ; --i){
        double res = B[i];
        for(int j = i + 1 ; j < n ; ++j){
            res -= A[i][j] * X[j];
        }
        res /= A[i][i];
        X[i] = res;
        if(X[i] == -0.)X[i] = 0;
    }
    return X;
//    for(int i = 0 ; i < n ; ++i)cout << X[i] << " ";
//    cout << endl;
}

int main(int argc, char* argv[]){
    cyrillic;
    cout.setf(ios::fixed);
    freopen("input.txt" , "r" , stdin);
 //   freopen("output.txt" , "w" , stdout);
    cout.precision(20);

    string parse_p = parse_path(argv[0]);
    vector<char> path(parse_p.begin() , parse_p.end());
    path.push_back('\0');
    concat(path , "temp.txt");

    vector<double> x;
    double tmp;

    while(cin >> tmp){
        x.push_back(tmp);
    }

    int n = x.size();

    vector<double> dft;

    FILE * file = fopen(&path[0] , "w");
    fprintf(file , "# X   Y\n");


    for(int i = 0 ; i < n ; ++i){
        double real = double(0);
        double imag = double(0);
        for(int j = 0 ; j < n ; ++j){
            double angle = double((2 * M_PI * i * j) / n);
            real += x[j] * double(cos(angle));
            imag += x[j] * double(sin(angle));
        }
        double val = sqrt((real * real) + (imag * imag));
        fprintf(file , "  %lf   %lf\n" , (i * 0.01) , val);
        dft.push_back(val);
    }
    fclose(file);
    vector<pair<double , double>> mx; //{coordinate , value}

    for(int i = 2 ; i < n / 2 ; ++i){
        if(dft[i] - dft[i - 1] > eps && dft[i] - dft[i + 1] > eps){
            mx.push_back(make_pair(i / 5.0, dft[i]));
        }
    }

//    for(int i = 0 ; i < mx.size() ; ++i){
//        cout << mx[i].first << ":" << mx[i].second << endl;
//    }

    vector<vector<double> > M(6 , vector<double>(6));
    vector<double> y_(6);

    for(int i = 0 ; i < n ; ++i){
        double ii = double(i * 0.01);

        for(int j = 0 ; j < 3 ; ++j){
            for(int k = 0 ; k < 3 ; ++k){
                M[j][k] += pow(ii , 6 - (j + k));
            }
        }

        for(int j = 0 ; j < 3 ; ++j){
            M[j][3] += pow(ii , 3 - j) * sin(2 * M_PI * ii * mx[0].first);
            M[3][j] += pow(ii , 3 - j) * sin(2 * M_PI * ii * mx[0].first);
        }

        M[3][3] += pow(sin(2 * M_PI * ii * mx[0].first) , 2);
        M[4][4] += pow(sin(2 * M_PI * ii * mx[1].first) , 2);
        M[4][3] += sin(2 * M_PI * ii * mx[0].first) * sin(2 * M_PI * ii * mx[1].first);
        M[3][4] += sin(2 * M_PI * ii * mx[0].first) * sin(2 * M_PI * ii * mx[1].first);

        for(int j = 0 ; j < 3 ; ++j){
            M[j][4] += pow(ii , 3 - j) * sin(2 * M_PI * ii * mx[1].first);
            M[4][j] += pow(ii , 3 - j) * sin(2 * M_PI * ii * mx[1].first);
        }

        for(int j = 0 ; j < 3 ; ++j){
            M[5][j] += pow(ii , 3 - j);
            M[j][5] += pow(ii , 3 - j);
        }

        M[3][5] += sin(2 * M_PI * ii * mx[0].first);
        M[5][3] += sin(2 * M_PI * ii * mx[0].first);

        M[4][5] += sin(2 * M_PI * ii * mx[1].first);
        M[5][4] += sin(2 * M_PI * ii * mx[1].first);

        M[5][5] = n;

        for(int j = 0 ; j < 3 ; ++j){
            y_[j] += pow(ii , 3 - j) * x[i];
        }

        y_[3] += sin(2 * M_PI * ii * mx[0].first) * x[i];
        y_[4] += sin(2 * M_PI * ii * mx[1].first) * x[i];
        y_[5] += x[i];
    }

    vector<double> a_i = gauss(M , y_);
/*  for(int i = 0 ; i < 6 ; ++i){
        cout << a_i[i] << " ";
    }*/

    string res = "";
    res += dtos(a_i[0]) + "*x*x*x+";
    res += dtos(a_i[1]) + "*x*x+";
    res += dtos(a_i[2]) + "*x+";
    res += dtos(a_i[3]) + "*sin(2*pi*x*" + dtos(mx[0].first) +")+";
    res += dtos(a_i[4]) + "*sin(2*pi*x*" + dtos(mx[1].first) +")+";
    res += dtos(a_i[5]);
    cout << endl << res << endl;

    Gnuplot gnu("lines");
    string aa = "set style line 1 lc rgb '#0060ad' lt 2 lw 2 pt 0 ps 1.5";
    string bb = &path[0]; bb = "plot '" + bb; bb += "' with linespoints ls 1";
    gnu.set_xrange(0. , 5.);
    gnu.set_title("DFT");
    gnu.plot_points(bb , aa);
    wait();
    gnu.reset_plot();
    gnu.set_title("f(x)");
    gnu.set_xrange(0. , 5.);
    gnu.plot_equation(res);
    wait();
    remove(&path[0]);

    return 0;
}