#include <math.h>
#include <stdlib.h>
#include <stdio.h>

/*
 * moore.c implementation of flow chart on p. 142 of Ramon E Moore
 * "Interval analysis", Prentice-Hall, 1966. The numbering of the 
 * labels B01, ..., B18 is for each box in the flow chart, starting
 * from top to bottom with the left column (B01 to B08) and 
 * continuing from top to bottom with the right column (B09 to B18).
 */


class iv {
	public:
		double l;
		double r;
		char t[100];
		iv () {
		}
		iv (double _l, double _r) {
			l = _l;
			r = _r;
		}
		iv clone () {
			return iv(l,r);
		}		
		void add (iv a) {
			l = l+a.l;
			r = r+a.r;
		}
		double abs() {
			double al,ar;
			al = fabs(l);
			ar = fabs(r);	
			return al > ar ? al : ar;
		}
		void mult (iv a) { //assumption: positive values
			l = l*a.l;
			r = r*a.r;
		}
		void pow (int a) { //assumption: positive values
			double accl = 1.;
			double accr = 1.;
			for (int i=0; i<a; i++) {
				accl = l*accl;
				accr = r*accr;
			}
			l = accl;
			r = accr;
		}
		void repr() {
			sprintf(t,"[%.16f,%.16f]", l, r);
		}
		void scale (double lambda) {
			l = l*lambda;
			r = r*lambda;
		}		
};

char s[256];
int xcounter;

void v() {
	printf("[%02d] %s\n",xcounter,s);
}

void vv() {
	//printf("[%02d] %s\n",xcounter,s);
}

double eps;
double range;
double xinterval;
double h,H;
int alternate;
iv x,A0,A1,Ar,Y0hAr,Yh,Y0,Yx;
int K;

int fac(int j) { //factorial
	int res = 1;
	for (int i=1;i<=j;i++)
		res = res * i;
	return res;
}

iv F(int j, iv in) {
	//p. 103 F^(j)(Y)=(j+1)!Y^(j+2)
	//factorial	
	iv store = in.clone();
	in = iv(1.,1.);	
	for (int i=1;i<=j+2;i++)
		in.mult(store);
	in.scale(fac(j+1));
	return in;
}

double f1038() { //determination of h
	double facc = pow(eps*Y0.abs()*fac(K)/F(K,Y0).abs(),(double)1/(K+1));
	if (facc < Y0.abs()/F(0,Y0).abs()) 
		h = facc;
	else
		h = pow(eps*F(0,Y0).abs()*fac(K)/F(K,Y0).abs(),(double)1/K);	
	//double h = pow(eps/(K+1.),1/(K+1.)); //p. 103 / 10-42
	return h;
}

iv f1039() { //first-order approx used for determination of A0/Ar
	iv n = F(0,Y0);
	n.mult(iv(0,h));
	n.add(Y0);
	return n;
}

iv f1040(iv x) { //k-th order approx
	iv res = Y0.clone();
	iv xpow;
	for (int j=1;j<=K-1;j++) {
		iv n = F(j-1,Y0);
		n.scale((double)1/(fac(j)));
		n.repr();
		sprintf(s,"f1040 coeff(%d): %s",j,n.t); vv();
		xpow = x.clone();
		xpow.pow(j);
		n.mult(xpow);
		res.add(n);
		res.repr();
		sprintf(s,"f1040 repr: %s",res.t); vv();
	}
	iv n = F(K-1,A1);
	A1.repr();
	sprintf(s,"f1040 A1: %s",A1.t); vv();
	n.scale((double)1/(fac(K)));
	n.repr();
	sprintf(s,"f1040 coeff(%d): %s",K,n.t); vv();
	xpow = x.clone();
	xpow.repr();
	sprintf(s,"f1040 x: %s",xpow.t); vv();
	xpow.pow(K);
	xpow.repr();
	sprintf(s,"f1040 xpow: %s",xpow.t); vv();
	n.mult(xpow);	
	n.repr();
	sprintf(s,"f1040 n: %s",n.t); vv();
	res.add(n);
	res.repr();
	sprintf(s,"f1040 res: %s",res.t); vv();
	return res;
}

int main () {

	eps = 0.000000001;
	//eps =   0.0000000074505806;
	K = 9;
	range = 0.9;
	xinterval = 0.01;
	xcounter = 0;
	goto B01;

B01:
	//yprime = f(y);
	Y0 = iv(1.,1.);
	Y0.repr();
	sprintf(s,"B01 Y0: %s",Y0.t); v();
	goto B09;

B02:
	alternate = 1;
	sprintf(s,"B02 alternate:=%d", alternate); v();
	goto B03;

B03:
	h = h/2;
	sprintf(s,"B03 h %f", h); v();
	goto B15;

B04:
	//substitute Y([0,h],A(r)) for A(r)
	Ar = Y0hAr.clone(); Ar.repr();
	sprintf(s,"B04 Ar: %s",Ar.t); v();
	goto B15;

B05:
	alternate = 2;
	sprintf(s,"B05 alternate:=%d", alternate); v();
	goto B04;

B06:
	sprintf(s,"B06 alternate==%d", alternate); v();
	if (alternate == 1) goto B05;
	else goto B02;

B07:
	//Is Y([0,h],A(r)) subseteq A(r)?
	Y0hAr.repr(); Ar.repr();	
	sprintf(s,"B07 Y0hAr: %s Ar: %s", Y0hAr.t, Ar.t); v();
	if (Y0hAr.l >= Ar.l && Y0hAr.r <= Ar.r) goto B08;
	else goto B06;

B08:
	//Determine Y(x) from (10-40) for x in [0,h] with Y([0,h],A(r))
	//substituted for A(1) and print Y(x) at desired values of
	//x in [0,h]. Then y(x+H) in Y(x).
	A1 = Y0hAr.clone();
	A1.repr();
	xcounter++;
	sprintf(s,"B08 xcounter: %d, h: %f, A1: %s",xcounter,h,A1.t); v();
	Yx = f1040(iv(h,h)); Yx.repr();
	sprintf(s,"B08 approx sol: %s, exact sol: %.16f", Yx.t, (double)1./(1.-(H+h))); v();
	for (double xf=0; xf<=h; xf+=xinterval) {
		Yx = f1040(iv(xf,xf));
		Yx.repr();
		printf("At Y(%f) is %s (exact:%.16f)\n", xf+H, Yx.t, (double)1./(1.-(H+xf)));
	}
	goto B17;

B09:
	//Generate subroutine for interval extensions of Taylor
	//coefficients (F(j)).
	//p. 103 F^(j)(Y)=(j+1)!Y^(j+2)
	sprintf(s,"B09"); v();
	goto B10;
	
B10:
	//Set H=0.	
	H=0;
	sprintf(s,"B10 H: %f",H); v();
	goto B11;

B11:
	//Set "alternate" at first.
	alternate = 1;	
	sprintf(s,"B11 alternate:=%d", alternate); v();
	goto B12;

B12:
	//Get h from (10-38).
	h = f1038();
	sprintf(s,"B12 h: %.16f",h); v();
	goto B13;

B13:
	//Get A(0) from (10-39).
	A0 = f1039(); A0.repr();
	sprintf(s,"B13 A0: %s", A0.t); v();
	goto B14;

B14:
	//Substitute A(0) for A(r).
	Ar = A0.clone(); Ar.repr();
	sprintf(s,"B14 Ar: %s", Ar.t); v();
	goto B15;

B15:
	//Substitute A(r) for A(1) and [0,h] for x and compute
	//Y([0,h],A(r)) from (10-40).
	A1 = Ar.clone();
	x = iv(0,h).clone();
	Y0hAr = f1040(iv(0,h)); Y0hAr.repr();
	sprintf(s,"B15 Y0hAr: %s", Y0hAr.t); v();
	goto B07;

B16:
	//Y(h)=Y0 and H=H+h
	Y0 = f1040(iv(h,h)); Y0.repr();
	H = H+h;
	sprintf(s,"B16 Y0: %s H: %f", Y0.t, H); v();
	goto B11;

B17:
	//Solution desired at x > h?
	sprintf(s,"B17 h: %f, range: %f", h, range); v();
	if (h<range) goto B16;
	else goto B18;	

B18: 
	//Stop
	sprintf(s,"B18 stop"); v();
	exit(0);	

}