import java.io.File;
import java.io.FileNotFoundException;
import java.io.PrintWriter;

import org.apache.commons.math.linear.Array2DRowRealMatrix;
import org.apache.commons.math.linear.ArrayRealVector;
import org.apache.commons.math.linear.DecompositionSolver;
import org.apache.commons.math.linear.LUDecompositionImpl;
import org.apache.commons.math.linear.RealMatrix;
import org.apache.commons.math.linear.RealVector;


public class Solution {
	public final static int a=-1;
	public final static int b=1;
	public final static int INTEGRAL_STEP=100000;
	public final static int lambda=1;
	private double h;
	private int n=1;
	private int sum;
	private double B[][];
	private double d[];
	private double c[];
	private double eps_prev=0;
	public double f(double x){
		return Math.cos(4*Math.PI*x)-
		(lambda*2*x*Math.sin(Math.PI*x))/(Math.PI*(x*x-16));
	}
	public double K(double x, double s){
		return Math.cos(Math.PI*x*s);
		//return Math.abs(x+s);
	}
	public double u(double x){
		return Math.cos(4*Math.PI*x);
		//return Math.abs(x);
	}
	public double ksi_0(double x){
		if (x>=a && x<=(a+h)){
			return (a+h-x)/h;
		}
		else{
			return 0;
		}
	}

	public double ksi_i(double x, int i){
		if (x>=a+(i-1)*h && x<=(a+h*i)){
			return (x-(a+(i-1)*h))/h;
		}
		else
			if (x>=a+i*h && x<=(a+h*(i+1))){
				return (x+(i+1)*h-x)/h;
			}
		else{
			return 0;
		}
	}

	public double ksi_n(double x){
		if (x>=a+h*(n-1) && x<=(a+h*n)){
			return (a+h-x)/h;
		}
		else{
			return 0;
		}
	}
	public double intergrateBij(int i, int j){
		double sum=0;
		double h=this.h/INTEGRAL_STEP;;
		double x;

		if (i==0){
			for (int k=0; k<INTEGRAL_STEP;k++){
				x=a+k*h;
				sum+=ksi_0(x)*K(x,a+this.h*j);
			}
		}
		else{
			if (i==n){
				for (int k=0; k<INTEGRAL_STEP;k++){
					x=b-k*h;
					sum+=ksi_n(x)*K(x,a+this.h*j);
				}
			}
			else{
				double x0=a+(i-1)*this.h;
				for (int k=0; k<INTEGRAL_STEP;k++){
					x=x0+h*k;
					sum+=ksi_i(x,i)*K(x,a+this.h*j);
				}
				x0=a+i*this.h;
				for (int k=0; k<INTEGRAL_STEP;k++){
					x=x0+h*k;
					sum+=ksi_i(x,i)*K(x,a+this.h*j);
				}
			}
		}
		return sum*h;
	}
	public double integrateDi(int i){
		double sum=0;
		double h=this.h/INTEGRAL_STEP;;
		double x;
		if (i==0){
			for (int k=0; k<INTEGRAL_STEP;k++){
				x=a+k*h;
				sum+=ksi_0(x)*f(x);
			}
		}
		else{
			if (i==n){
				for (int k=0; k<INTEGRAL_STEP;k++){
					x=b-k*h;
					sum+=ksi_n(x)*f(x);
				}
			}
			else{
				double x0=a+(i-1)*this.h;
				for (int k=0; k<INTEGRAL_STEP;k++){
					x=x0+h*k;
					sum+=ksi_i(x,i)*f(x);
				}
				x0=a+i*this.h;
				for (int k=0; k<INTEGRAL_STEP;k++){
					x=x0+h*k;
					sum+=ksi_i(x,i)*f(x);
				}
			}
		}
		return sum*h;
	}
	public void initMartix(){
		d=new double[n+1];
		B=new double[n+1][];
		for (int i=0; i<n+1; i++){
			B[i]=new double[n+1];
		}
		for (int i=0; i<n+1; i++){
			for (int j=0; j<n+1; j++){
				B[i][j]=lambda*intergrateBij(i, j);
				if (i!=j){
					B[i][j]=-B[i][j];
				}
				else{
					B[i][j]=1-B[i][j];
				}
			}
			d[i]=integrateDi(i);
		}
		RealMatrix coefficients=new Array2DRowRealMatrix(B);
		DecompositionSolver solver=new LUDecompositionImpl(coefficients).getSolver();
		RealVector constants=new ArrayRealVector(d);
		RealVector solution=solver.solve(constants);
		c=solution.getData();
	}
	public double u_wave(double x){
		double sum=0;
		for (int i=0; i<n+1;i++){
			sum+=c[i]*K(x,a+i*h);
		}
		sum*=lambda;
		sum+=f(x);
		return sum;
	}

	public void solve() throws FileNotFoundException{
		PrintWriter writer;
		double x_i;
		double eps_max;
		double eps;
		double u_1,u_2;
		long start, finish;
		for (int i=1; i<=10; i++){
			System.out.println("I="+i);
			start=System.currentTimeMillis();
			sum=0;
			n*=2;
			h=(b-a)/(double)n;
			initMartix();
			eps_max=-1;
			for (int j=1; j<=n; j++){
				x_i=a+(j-1)*h;
				u_1=u(x_i);
				u_2=u_wave(x_i);
				System.out.println("with u="+u_1+" and u_wave="+u_2);
				eps=Math.abs(u_1-u_2);
				if (eps>eps_max){
					eps_max=eps;
				}
			}
			finish=System.currentTimeMillis();
			writer = new PrintWriter(new File("out" + i + ".txt"));
			writer.println("h_i = "+h);
			writer.println("eps = "+eps_max);
			writer.println("r_i = " + Math.log((eps_prev)/eps_max)/Math.log(2));
			eps_prev=eps_max;
			writer.println("t_i = " + (finish-start));
			writer.close();
		}


	}

}