package com.company;


import javax.script.ScriptEngine;
import javax.script.ScriptEngineManager;
import java.beans.Expression;
import java.io.Reader;

public class Main {

    /*  Aproksymacja
    *
    * - aproksymacja punktowa
    * - wielomiany ortogonalne
    * - aproksymacja średniokwadratowa ciągła
    *
    * */




    private static double[] gauss(double s[][]) {

        int n = s.length;
        for (int i = 0; i < n; i++) {
            double maxEl = Math.abs(s[i][i]);
            int maxRow = i;
            for (int k = i + 1; k < n; k++) {
                if (Math.abs(s[k][i]) > maxEl) {
                    maxEl = Math.abs(s[k][i]);
                    maxRow = k;
                }
            }

            for (int k = i; k < n + 1; k++) {
                double tmp = s[maxRow][k];
                s[maxRow][k] = s[i][k];
                s[i][k] = tmp;
            }

            for (int k = i + 1; k < n; k++) {
                double c = -s[k][i] / s[i][i];
                for (int j = i; j < n + 1; j++) {
                    if (i == j) {
                        s[k][j] = 0;
                    } else {
                        s[k][j] += c * s[i][j];
                    }
                }
            }
        }

        double[] var = new double[n];
        for (int i = n - 1; i > -1; i--) {
            var[i] = s[i][n] / s[i][i];
            for (int k = i - 1; k > -1; k--) {
                s[k][n] -= s[k][i] * var[i];
            }
        }

        return var;
    }



    private static void aproksymacjaSredniokwadratowa(int n, double a, double b, double px, double _x) {

        double [][]A = new double[n][n+1];

        int precyzja = 1000;
        double []fx = new double[precyzja+1];
        double []x = new double[precyzja+1];


        for (int i=0; i<n; i++) {
            for (int j=0; j<n+1; j++) {
                if(j<n) {
                    for(int k=0; k<precyzja; k++) {
                        x[k] = a + (k*(b-a))/precyzja;
                        fx[k] = Math.pow(x[k], i) * Math.pow(x[k], j) * px;
                    }
                } else {
                    for(int k=0; k<precyzja; k++) {
                        x[k] = a + (k*(b-a))/precyzja;
                        fx[k] = Math.pow(x[k], i) * Math.sqrt( x[k]*x[k] + 3*x[k] +7 ) * px;
                    }
                }
                A[i][j] = trapez(precyzja,a,b,fx);
            }
        }


        double []W = gauss(A);
        double result = 0;
        for (int i=0; i<W.length; i++) {
            result += W[i] * Math.pow(_x,i);
        }
        System.out.println(result);
    }

    private static void wielomianyOrtogonalny(int n, double a, double b, double px, double x){
        int precyzja = 1000;

        double []fx = new double[precyzja+1];
        double []xi = new double[precyzja+1];

        double []Fi = new double[n+1];
        double []_Fi = new double[n+1];
        double []Lambda = new double[n+1];
        double []C = new double[n+1];


        double gx = 0;

        for (int i=0; i<=n; i++) {

            for(int j=0; j<=precyzja; j++) {
                xi[j] = a + (j*(b-a))/precyzja;

                Fi[0] = 1;
                Fi[1] = xi[j];
                Fi[2] = 1.0/2 * (3*xi[j]*xi[j] - 1);

                if(i>2) Fi[i] = 1.0/i *((2*(i-1)+1)*xi[j] * Fi[i-1] - (i-1)*Fi[i-2] );

                fx[j] = px * Math.pow(Fi[i],2);
            }

            Lambda[i] =  trapez(precyzja,a,b,fx);

            for(int j=0; j<=precyzja; j++) {
                xi[j] = a + (j*(b-a))/precyzja;

                Fi[0] = 1;
                _Fi[0] = 1;
                Fi[1] = xi[j];
                _Fi[1] = x;
                Fi[2] = 1.0/2 * (3*xi[j]*xi[j] - 1);
                _Fi[2] = 1.0/2 * (3*x*x - 1);

                if(i>2) {
                    Fi[i] = 1.0/i *((2*(i-1)+1)*xi[j] * Fi[i-1] - (i-1)*Fi[i-2] );
                    _Fi[i] = 1.0/i *((2*(i-1)+1)*x * Fi[i-1] - (i-1)*Fi[i-2] );
                }

                fx[j] = px * Fi[i] * Math.sqrt( xi[j]*xi[j] + 3*xi[j] +7 );
            }

            C[i] = (1.0/Lambda[i]) * trapez(precyzja,a,b,fx);
            gx += C[i]*_Fi[i];

        }

        System.out.println(gx);
    }

    private static void aproksymacjaPunktowa(int n, double _x){
        double []x = {-1, -0.5, 0, 0.5, 1};
        double []fx = new double[x.length];

        // wyliczanie fx
        for(int i=0; i<x.length; i++) {
            fx[i] = Math.sqrt(x[i] * x[i] + 3 * x[i] + 7);
        }

        double s[][] = new double[x.length][n + 2*n + 2];
        double sum[] = new double[n + 2*n + 2];

        // tworzenie macierzy
        for(int i=0; i < x.length; i++) {
            for(int j=0, k=0; j <= n + 2*n +1 ; j++) {
                // x^n * y
                if(j>2*n)  {
                    s[i][j] = Math.pow(x[i], k) * fx[i];
                    k++;
                } else {
                    s[i][j] = Math.pow(x[i],j);
                }
            }
        }

        // wyliczanie sumy
        for(int i=0; i <= n + 2*n +1  ; i++) {
            sum[i] = 0;

            for(int j=0; j < x.length; j++) {
               sum[i] += s[j][i];
            }
        }

        double _s[][] = new double[n+1][n+2];
        //tworzenie macierzy z sum
        for(int i=0, k=0; i <= n; i++) {
            for(int j=0, z=i; j <= n+1; j++, z++) {
                if (j == n + 1) {
                    _s[i][j] = sum[2 * n + k + 1];
                    k++;
                } else {
                    _s[i][j] = sum[z];
                }
            }
        }

        double aGauss[] = gauss(_s);
        double result = 0;
        for(int i=0; i<n; i++) {
            result += aGauss[i] * Math.pow(_x,i);
        }

        System.out.println(result);
    }

    private static double trapez(int precyzja, double a, double b, double []fx) {

        double h = (b-a)/precyzja;
        double result = 0;

        for(int i=1; i<=precyzja; i++) {
            result+= fx[i-1] + fx[i];
        } result += h/2;

        return result/precyzja;
    }



    public static void main(String[] args) {

        int n = 50;
        double a = -1,
               b = 1;

        aproksymacjaSredniokwadratowa(n,a,b,1,-0.25);
        wielomianyOrtogonalny(n,a,b,1, -0.25);
        aproksymacjaPunktowa(n,-0.25);
    }
}