import static java.lang.Math.*;
import static java.util.Arrays.fill;

public class GlodusExpo{

private double s, sih, q, a, b, a0, b0, a1, b1,a2, b2;
private int N;

public double[] re, re1, re2, uf0, ub0, uf1, ub1;
/** @shw - array for show */
double[] shw, shw1, shw2;

public double sigma(){ return s; }
// exponential filter x exp( -|x|/sigma );
// filtruojami duomenys tiesiskai interpoliuojami
// N - dimension of data to be filtered by the filter

public GlodusExpo( double sigma, int N )
{
 setSigma( sigma );
 init(N);
}

public void init( int N )
{
 this.N = N;
 uf0 = new double[N]; ub0 = new double[N];
 uf1 = new double[N]; ub1 = new double[N];
 re = new double[N]; re1 = new double[N]; re2 = new double[N];
}

public void setSigma( double sigma )
{
 this.s = sigma;
 q = exp(-1/s);
 a = 1/(1-q);
 b = q/((1-q)*(1-q));
 sih = sinh(1/s);
 double[][] A = new double[][]{{q0(),q1()},{q2(),q3()}};
 double[] a = linsolve2( A, new double[]{1,0});
 a0 = a[0]; b0 = a[1];  
 a1 = -1/q2();
 b1 = -a1;
 a = linsolve2( A, new double[]{0,2});
 a2 = a[0]; b2 = a[1];
// System.out.println( "\n==== sigma=" + sigma + " ========\na0=" + a0 + " b0=" + b0 );
 a0 = 1/(q0() + q1()/sigma);
 b0 = a0/sigma;
}

public void filter( int[] data, int N )
{
 if ( this.N < N ) init(N);
 for ( int n = 0; n != N; n++ ) re[n] = data[n];
 filter( re, N );
}

public void filter( double[] data, int N )
{
	if ( this.N < N ) init( N );
	final int N1 = N - 1;
	uf0[0] = a*data[0]; uf1[0] = b*data[0];
	for ( int n = 1, n0 = n-1; n < N; )
	{
		uf0[n] = q * uf0[n0] + data[n];
		uf1[n++] = q * ( uf1[n0] + uf0[n0++] );
	}
	ub0[N1] = a*data[N1]; ub1[N1] = b*data[N1];
	for ( int n = N1-1, n1 = n+1; n >= 0; )
	{
		ub0[n] = q * ub0[n1] + data[n];
		ub1[n--] = q * ( ub1[n1] + ub0[n1--] );
	}
	for ( int n = 0; n != N; n++ )
	{
		double d = ub0[n] + uf0[n] - data[n], e = ub1[n] + uf1[n];
		re[n] = a0*d + b0*e;
		re1[n] = a1*( uf1[n] - ub1[n] );
		//re1[n] = ( ( ub0[n] + ub1[n] ) - (uf0[n] + uf1[n]) );
		//if ( n==N/2 ) util.my.pause( q2() + " " + q3() + " " + d + " " + e );
		re2[n] = a2*d + b2*e;
	}

}

public void vidur( double[] data, int N )
{
	if ( this.N < N ) init( N );
	final int N1 = N - 1;
	uf0[0] = a*data[0]; uf1[0] = b*data[0];
	for ( int n = 1, n0 = n-1; n < N; )
	{
		uf0[n] = q * uf0[n0] + data[n];
		uf1[n++] = q * ( uf1[n0] + uf0[n0++] );
	}
	ub0[N1] = a*data[N1]; ub1[N1] = b*data[N1];
	for ( int n = N1-1, n1 = n+1; n >= 0; )
	{
		ub0[n] = q * ub0[n1] + data[n];
		ub1[n--] = q * ( ub1[n1] + ub0[n1--] );
	}
	for ( int n = 0; n != N; n++ )
		re[n] = a0*(ub0[n] + uf0[n] - data[n]) + b0*(ub1[n] + uf1[n]);
}

// int q^(|x|) dx = int exp(-|x|/s ) dx,  q = exp(-1/s) -> s = -1/ln(q)
public double integ0( double q )
{
 double s = -1/log(q);
 return 2*s;
}
// int |x|q^(|x|) dx = int |x|exp(-|x|/s ) dx = s * int q^(|x|) dx
public double integ1( double q )
{
 return s*integ0(q);
}
// int |x|^2q^(|x|) dx = int |x|^2exp(-|x|/s ) dx = 2 * s * int |x|q^(|x|) dx
public double integ2( double q )
{
 return 2*s*integ1(q);
}

// int_0^b e^(-x/sigma) dx
public double int0( double b )
{
 return -s * ( exp(-b/s) - 1 );
}

// int_0^b x e^(-x/sigma) dx
public double int1( double b )
{
 return -s * ( b*exp(-b/s) - int0(b) );
}

// int_0^b x^2 e^(-x/sigma) dx
public double int2( double b )
{
 return -s * ( b*b*exp(-b/s) - 2*int1(b) );
}

// sum_m=-infty^infty q^|m| = (1+q)/(1-q)
public static double q0( double q )
{
 return (1+q)/(1-q);
}

// sum_m=-infty^infty |m|q^|m| = 2*q/(1-q)^2
public static double q1( double q )
{
 return 2*q/((1-q)*(1-q));
}
// sum_m=-infty^infty m^2q^|m| = 2q(1+q)/(1-q)^3
public static double q2( double q )
{
 return 2*q*(1+q)/((1-q)*(1-q)*(1-q));
}
// sum_m=-infty^infty m^3q^|m| = 2q(1+4*q+q^2)/(1-q)^4
public static double q3( double q )
{
 return 2*q*(1+4*q+q*q)/((1-q)*(1-q)*(1-q)*(1-q));
}

// sum_m=-infty^infty q^|m| = (1+q)/(1-q)
public double q0()
{
 return (1+q)/(1-q);
}
// sum_m=-infty^infty |m|q^|m| = 2*q/(1-q)^2
public double q1()
{
 return 2*q/((1-q)*(1-q));
}
// sum_m=-infty^infty m^2q^|m| = 2q(1+q)/(1-q)^3
public double q2()
{
 return 2*q*(1+q)/((1-q)*(1-q)*(1-q));
}
// sum_m=-infty^infty m^3q^|m| = 2q(1+4*q+q^2)/(1-q)^4
public double q3()
{
 return 2*q*(1+4*q+q*q)/((1-q)*(1-q)*(1-q)*(1-q));
}
public static double[] linsolve2(double[][] a, double[] b){ // ax = b
// if ((a.length!=2)|(b.length!=2)) print_exit("linsolve2 parashyta tik 2x2 lygciu sistemai");
 double D = a[0][0]*a[1][1]-a[0][1]*a[1][0];
 return new double[]{(b[0]*a[1][1]-b[1]*a[0][1])/D,(a[0][0]*b[1]-a[1][0]*b[0])/D};
}

public static void main( String[] args )
{
	int N = 300;
	double sigma = 1., u[] = new double[N];
	fill( u, 2 );
	fill( u, N/2, N, -0.5 );
	GlodusExpo e = new GlodusExpo( sigma, N );
	e.vidur( u, N );
}

}