/*
	jfdiff.java -- Compute difference scheme to solve p(E)x=0
	Written September 12, 2012 for Math 701 by Eric Olson

	Adapted from fdiff.c Version 2.
*/

import static java.lang.System.*;
import static java.lang.Math.*;

class jfdiff {

static double VERSION=2.0;
static int N=200;
static int L=2;
static double[] lambda={ 1.0/3.0, 3.0/2.0 };

static public class polynomial {
	public int d;
	public double[] c = new double[L+1];
}

static void mkeigvec(double[] x,double lambda,int n){
	double xn;
	int i;
	x[0]=xn=1.0;
	for(i=1;i<n;i++){
		xn*=lambda;
		x[i]=xn;
	}
}

static void pmulfact(polynomial p,double lambda){
	int i;
	p.c[p.d+1]=p.c[p.d];
	for(i=p.d++;i>0;i--) p.c[i]=p.c[i-1]-p.c[i]*lambda;
	p.c[0]*=-lambda;
}

static void printpoly(polynomial p,String s){
	int i;
	out.printf("%g",p.c[0]);
	for(i=1;i<=p.d;i++) out.printf("%+g%s^%d",p.c[i],s,i);
}

static void printvec(double x[],int n){
	int i;
	out.printf("%g",x[0]);
	for(i=1;i<n;i++) out.printf(",%g",x[i]);
}

static double[] x = new double[N];
static double[] x_lambda1 = new double[N];
static double[] c = new double[L+1];

static public void main(String[] args) {
	int i,k;
	out.printf(
		"jfdiff -- Compute difference scheme to solve p(E)x=0 "+
		"Version %3.1f\n"+
		"Written September 12 for Math 701 by Eric Olson\n",VERSION);
	mkeigvec(x_lambda1,lambda[0],N);
	polynomial p = new polynomial();
	p.d=0; p.c[0]=4.0;
	for(k=0;k<L;k++) pmulfact(p,lambda[k]);
	mkeigvec(x,lambda[0],L);
	
	out.printf("\nx=(");
	printvec(x,L);
	out.printf(",...)\n{lambda_n}={");
	printvec(lambda,L);
	out.printf("}\np(x)=");
	printpoly(p,"x");
	out.printf("\n\n");

	for(i=L;i<N;i++){
		int j=0; double s=0.0;
		for(;j<L;j++) s+=p.c[j]*x[j+i-L];
		x[i]=-s/p.c[L];
	}

	out.printf("%16s %16s %16s\n","x","x_lambda1","|x-x_lambda1|");
	for(i=0;i<N;i++)
		out.printf("%16.8e %16.8e %16.8e\n",x[i],x_lambda1[i],
			abs(x[i]-x_lambda1[i]));
}

} // class jfdiff