var fPrec = 4; // floating output precision
var gIterMax = 100; // solving
var EPS = 0.0001 ;
var gIter=0; // dbg

/*----------------------------
mylog
-----------------------------------*/
function mylog(x)
{
if(x > 0) return Math.log(x);
return 0;
}

/*---------------------------
MATH LIB
----------------------------------*/
sin = Math.sin ;
log = mylog;
abs=Math.abs;
acos=Math.acos;
asin=Math.asin;
atan=Math.atan;
atan2=Math.atan2;
ceil=Math.ceil;
cos=Math.cos;
exp=Math.exp;
floor=Math.floor;
max=Math.max;
min=Math.min;
pow=Math.pow;
random=Math.random;
round=Math.round;
sqrt=Math.sqrt;
tan=Math.tan;
PI = Math.PI ;
E = Math.E ;
GAMMA = 0.57721566490153286060651209008240243104215933593992 ; 

/*---------------------------------
nprompt
-----------------------------*/
function nprompt(q,ndef)
{
var nret = prompt(q,ndef);
if(nret == null) return ndef;
nret = parseInt(nret);
if(isNaN(nret)) return ndef;
return nret;
}

/*----------------------------
format
-------------------------------*/
function format(x)
{
if(x >=0 || x < 0)
	{
	if(isint(x)) return x;
	if(fPrec < 0) fPrec = 4;
	return x.toFixed(fPrec);
	}
return x;
}

var _U = 	new Array(); // internal
/*---------------------------------
use(v)
-----------------------------------*/
function use (v)
{
 _U[v]=1;
}

/*-------------------------
least not used > s
usage:  for(n=notused(0);;n=notused(n)) {...}
-----------------------------*/
function notused(s)
{
for(var v=s+1;;v++) if(_U[v] == undefined) return v;
}

/*--------------------------------
digitsum
------------------------------------*/
function digitsum (n)
{
var s= 0; var d;
	while(n > 0) {d= n%10; s+=d ; n-=d; n/= 10;}
	return s;
}

/*--------------------------
indigits (a,b) b visible in a = a shows b
------------------------------*/
function indigits(a,b)
{
return (''+a).indexOf(''+b) >=0  ;
}

/*---------------------------------
isint
--------------------------------*/
function isint(n)
{
// if(n < 0) n=-n;
return (Math.floor(n) == n );
}
var isInt = isint ;
/*-------------------------------
powmod
-----------------------------------*/
function powmod(b,e,m){
var r = 1;
    while (e > 0) {
        if (e & 1) {
            r = (r * b) % m;
        }
        e >>= 1;
        b = (b * b) % m;
    }
    return r;
}

/*-----------------------------------------
numdiv
---------------------------------------------*/
function numdiv(n)
{
var nd= 2; // 1 & n
if(n <= 1) return 1;
for(var d=2;(n/d) >= d;d++) if(n % d == 0) nd++;
return nd;
}

/*-----------------------------------------
sumdiv
---------------------------------------------*/
function sumdiv(n)
{
var sd= 1+n; // 1 & n
if(n <= 1) return 1;
for(var d=2;(n/d) >= d;d++) if(n % d == 0) sd += d ;
return sd;
}

var _nth = 2;
var _fnth = 3;

/*-------------------------------------------
nthprime(n)
-----------------------------------------------*/
function nthprime(n)
{
if(n < 1) return undefined;
if(n == 1) return 2;
if(_nth > n) {_nth = 2; _fnth=3;}
if(_nth == n ) return _fnth;
for(var p =_fnth+2;;p+=2)
	{
	if(isprime(p)) 
		{ 
		_nth++;
		_fnth = p; 
		if(_nth == n) return p;
		}
	}
} // nthprime

var _np = 1;
var _fnp = 0;
/*----------------------------------------------------
numprimes (<= n)
use a cache array _ffnp[x] = numprimes(10**x) NYI
------------------------------------------------------*/
function numprimes(n)
{
	if(_np > n) {_np = 1; _fnp=0;} // reset cache
	if(_np == n) return _fnp;
	_np++;
	for(; _np <=n; _np++) if(isprime(_np)) _fnp++;
	_np--;
	return _fnp;
}
/*------------------------------------------
 isprime
-----------------------------------------------*/
function isprime(n) // brute force for small numbers
	{
	if(2 > n) return 0;
	if(2 == n) return 1;
	if (! (n & 1)) return 0;
	for (var div = 3; (n / div) >= div; div += 2) if (0 == n % div) return 0;
	return 1;
	} // isprime


/*------------------------------------------
factorial - 
----------------------------------------*/
function fact(n)
{
var f= 1;
	if(n > 25) return Infinity; 
	if(n <= 1) return 1;
	for(var i=2;i<=n;i++) f*=i;
	return f;
}

/*--------------------------------------------------------------
binomial  bin(n,k) = 0 if k > n -  bin(n,k) = n*n-1*..n-k+1 / k!
-------------------------------------------------------------------*/
function cbin(n,k)
{
var num=1; var den = 1; var g;
if(k > n) return 0;
for(var i=1;i <=k ; i++)
	{
	den *= i ;
	num *= n+1-i;
	g = gcd(num,den); if(g != 1) { num/=g ; den /=g;}
	}
	return num/den ;
} // cbin

var _ndiff = -1; // cache
var _fdiff = 0;
/*-------------------------------------
diff(n,f) : f(n+1)- f(n) : using cache
-------------------------------*/
function diff(n,f,k)
{

if(k <=1 || arguments.length == 2) k= 1;
if(k == 1)
	{
	var fnext,delta;

	fnext  = f(n+1) ; ;
	if(n != _ndiff) _fdiff = f(n);
	delta = fnext - _fdiff ;
	_ndiff = n+1; _fdiff = fnext;
	return delta;
	}
	return diff(n+1,f,k-1) - diff(n,f,k-1);
}

	
function sigma(nmin,nmax,f){var s= 0;for(var i=nmin;i<=nmax;i++) s += f(i);return s;}
function product(nmin,nmax,f){var s= 1;for(var i=nmin;i<=nmax;i++) s *= f(i);return s;}
function collatz(n) {if(n&1) return n*3+1; return n/2;}
function gcd(a, b) {var t;if(b <=0 || a <= 0) return 1; while (b != 0) {t=b; b= a % b; a=t;} return a;}
function nextprime(n) {if(n <=1) return 2; if(n&1)n+=2; else n++; while(!isprime(n)) n+=2; return n;}

/*---------------------------------------------
resetcaches()
--------------------------------------*/
function resetcaches()
{
	_U = new Array();
	_ndiff = -1;
	_delta=0; // ??? check me
}

/*---------------------------------------
opposite
-----------------------------------*/
function opposite(a,b)
{
if(isNaN(a) || isNaN(b) || a== undefined || b==undefined) return false;
return ((a<0 && b>0) || (a>0 && b<0)) ;
}

/*---------------------------------------------
_inv(fun, x)
returns y such as f(y) = x
------------------------------------------*/
function _inv(f_,x,ymin,ymax)
{
var y,fy,ffrom ; 
var yfrom,yto;
var _delta = (ymax-ymin)/2 ;
if(gInt) _delta= Math.floor(_delta);

	if(f_(ymin) == x) return ymin;
	
// try small ints
	if(gInt)
	{
	yto = Math.ceil(min(ymin + 1000,ymax));
	for(yfrom = Math.floor(ymin); yfrom <= yto ;yfrom++)
		if(f_(yfrom) == x) return yfrom;
	if(yfrom >= ymax) return undefined ; // no hope if gInt
	}
	
	
// general
	yfrom = ymin; yto = ymax;

while (!opposite(f_(yfrom)-x,f_(yto)-x))
	{
	if(f_(yfrom) == x) return yfrom;
	yfrom += _delta;
	if(yfrom > yto) {_delta /= 2; yfrom = ymin; if(gInt) _delta= Math.floor(_delta);}
	if(_delta < EPS) return undefined;
	}
	
	ffrom = f_(yfrom);
while(1)
	{
	y = yfrom + (yto-yfrom)/2 ;
	if(gInt) y = Math.floor(y+0.5);
	fy = f_(y);
	if(abs(fy - x) < EPS) break;
	if(opposite(fy - x,ffrom -x )) yto = y; else {yfrom = y; ffrom = f_(yfrom);}
	if(yto - yfrom < EPS) break;
	}
				
	if(y - Math.floor(y) < EPS) y = Math.floor(y);
	else
	if(Math.ceil(y) - y < EPS) y = Math.ceil(y);
	return y;		
} // _inv


/*---------------------------------------------
_zero (fun(n,y),n0 ,ymin,ymax)
find y0 : fun(n0,y0) == 0
uses global gInt
--------------------------------------------*/
function _zero(f_,n,ymin,ymax)
{
var y,fy,ffrom,m ; 
var yfrom,yto;
var _delta = (ymax-ymin)/2 ;
if(gInt) _delta= Math.floor(_delta);

	if(f_(n,ymin) == 0) return ymin;
	
// try small ints
	if(gInt)
	{
	yto = Math.ceil(min(ymin + 1000,ymax));
	for(yfrom = Math.floor(ymin); yfrom <= yto ;yfrom++)
		if(f_(n,yfrom) == 0) return yfrom;
	if(yfrom >= ymax) return undefined ; // no hope if gInt
	}
	
	
// general
	yfrom = ymin; yto = ymax;

while (!opposite(f_(n,yfrom),f_(n,yto)))
	{
	gIter++;// dbg
	if(f_(n,yfrom) == 0) return yfrom;
	yfrom += _delta;
	if(yfrom > yto) {_delta /= 2; yfrom = ymin; if(gInt) _delta= Math.floor(_delta);}
	if(_delta < EPS) return undefined;
	}
	
	m=gIterMax ;
	ffrom = f_(n,yfrom);
while(m--)
	{
	gIter++;
	y = yfrom + (yto-yfrom)/2 ;
	if(gInt) y = Math.floor(y+0.5);
	fy = f_(n,y);
	if(abs(fy) < EPS) break;
	if(opposite(fy,ffrom)) yto = y; else {yfrom = y; ffrom = f_(n,yfrom);}
	if(yto - yfrom < EPS) break;
	}
				
	if(y - Math.floor(y) < EPS) y = Math.floor(y);
	else
	if(Math.ceil(y) - y < EPS) y = Math.ceil(y);
	return y;		
} // _zero

var gm,gf1,gfn,gdgyn,gyn_1;
/*-----------------------------------
find a zero using secant method
works only if (yn_1,yn) are good approx
-------------------------------------*/
function _secant(f_,n,yn_1,yn)
{
var m = gIterMax ;
var d,fn,f1;
var yfrom = yn_1;
var yto = yn; // save for _zero

	if(gInt) return _zero(f_,n,yn_1,yn);
	
	f1 = f_(n,yn_1);
	if(abs(f1) < EPS) return yn_1 ;
	while(m--)
	{
gyn = yn;
gyn_1=yn_1;
gm= m;
gIter++;
	fn = f_(n,yn);
	if(abs(fn) < EPS) return yn ;
	d = fn-f1;
gfn=fn; gf1=f1;
	if(abs(d) < EPS) return _zero(f_,n,yfrom,yto);
	d =((yn - yn_1) * fn) / d;
gd=d;
	if(abs(d) < EPS) return _zero(f_,n,yfrom,yto); // original return yn;
	yn_1 = yn;
	f1 = fn;
	yn = yn - d;
	}
gm= m;
	return yn;
} // _secant

/*-----------------------------
derivative

        -3*y(x) + 4*y(x+h) - y(x+2h)
y'(x) = ------------------------------
                 2*h

----------------------------*/
function _deriv(f,x,h)
{
return (-3*f(x) + 4*f(x+h) - f(x+2*h)) /(2*h) ;
}

/*---------------------------------
deriv
------------------------------*/
function deriv(f,x)
{
var m = gIterMax;
var d1,d2;
var h= EPS;
var fx,fh,f2h,fh_,f2h_;

	fx=  f(x);
	fh=  f(x+h);
	f2h= f(x+2*h);
	fh_= f(x-h);
	f2h_ = f(x-2*h);
	
while(m--)
	{
	gIter++; // call f()
	d1= (-3*fx +4*fh - f2h) /(2*h);
	d2= (-3*fx +4*fh_ -f2h_)/(-2*h);
	if(abs(d1-d2) < EPS) return d1;
	h /=2 ;
	f2h =  fh;
	f2h_ = fh_;
	fh =  f(x+h);
	fh_ = f(x-h);
	}
	
return 0; // undef
} // deriv
	

/*-------------------------
adaptive integrate
---------------------------------*/

function adaptive(f,a,b,eps,S,fa,fb,fc,depth)
{
var c = (a+b)/2  ;
var h = b-a;
var d =(a+c)/2;
var e = (c+b)/2 ;
var fd = f(d);
var fe = f(e);
var Sleft =  (h/12)*(fa + 4*fd + fc); 
var Sright = (h/12)*(fc + 4*fe + fb);                                                          
var S2 = Sleft + Sright;    
if (depth <= 0 || abs(S2 - S) <= 15*eps )                                                   
    return S2 + (S2 - S)/15;                                                                        

return adaptive(f, a, c, eps/2, Sleft,  fa, fc, fd, depth-1) +                    
       adaptive(f, c, b, eps/2, Sright, fc, fb, fe, depth-1);    

}

/*-------------------------------
integrate
---------------------------------------*/
function integrate(f,a,b,depth)
{
var c = (a+b)/2 ;
var h = b-a ;
var fa = f(a);
var fb = f(b);
var fc = f(c);
var S  = (h/6)*(fa + 4*fc + fb);   
if(arguments.length < 4) depth = 10 ;
return adaptive(f,a,b,EPS,S,fa,fb,fc,depth);
}

var _gnorm = 1./sqrt(PI*2) ;
/*----------------------
Gauss
------------------------------*/
function gauss(x,m,s)
{
if(arguments.length < 3) s = 1;
if(arguments.length < 2) m = 0;
return _gnorm /s  * exp(- ((x-m)/s) * ((x-m)/s)  /2) ;
} // gauss

/*--------------------------------------
regression  
ldata --> odata,bdata

y = b1*x + b0
b1 = sum((x_i-x_avg)*(y_i-y_avg)) / sum((x_i-x_avg)^2)
b0 = y_avg - b1*x_avg
----------------------------*/
function _regression(ldata)
{
var b0=0
var b1=0;
var i;
var y_avg = 0;
var n = ldata.length ; // x in o..n-1
var x_avg = (n-1)/2;

for(i=0;i<n;i++) y_avg += ldata[i] ; y_avg /= n;
for(i=0;i<n;i++) {b1 += (i - x_avg)*(ldata[i]-y_avg) ; b0+= (i-x_avg)*(i-x_avg);}
b1/=b0;
b0=y_avg - b1*x_avg;

odata= new Array();
bdata =new Array();

for(i=0;i<n;i++)
	{
	odata.push (b1*i + b0);
	bdata.push (b1*i + b0);
	}
} // regression


/*----------------------------------------
http://en.wikipedia.org/wiki/Spline_(mathematics)#Algorithm_for_computing_natural_cubic_splines
clamped cubic splines
in array x [0...n] array y[0...n]
---------------------------------*/
var _sa,_sb,_sc,_sd,_sx ;

function ccsplines(x,y)
{
var i;
var l,mu; 
var h,alf,z ;

var n = x.length - 1;
_sa = new Array(n+1); for(i=0;i<=n;i++) _sa[i]= y[i];
_sb = new Array(n);
_sd = new Array(n);
h =  new Array(n); for(i=0;i<n;i++) h[i] =x[i+1] - x[i];
alf= new Array(n); for(i=1;i<n;i++) alf[i] =3/h[i]*(_sa[i+1] - _sa[i]) - 3/h[i-1]*(_sa[i] - _sa[i-1]) ;
_sc = new Array(n+1);
l = new Array(n+1);
mu = new Array(n+1);
z = new Array(n+1);
l[0]= 1.
mu[0]= z[0]=0;

for(i=1;i<=n-1;i++)
{
l[i] = 2*(x[i+1]-x[i-1]) - h[i-1]*mu[i-1];
mu[i] = h[i]/l[i];
z[i]= (alf[i]-h[i-1]*z[i-1])/l[i] ;
}

l[n]=1; z[n]=_sc[n]=0;

for(i=n-1;i>=0;i--)
{
_sc[i] = z[i] - mu[i]*_sc[i+1];
_sb[i] =(_sa[i+1]- _sa[i])/h[i] - h[i]*(_sc[i+1]+2*_sc[i]) / 3 ;
_sd[i] = (_sc[i+1]-_sc[i])/ (3 * h[i]);
}

_sx = new Array(n);
for(i=0;i<n;i++) _sx[i] = x[i];
} 

/*--------------------------------------------
call a Spline i in 0....n-1 at  _sx[i] < x <  _sx[i+1]
-----------------------------------------*/
function spline (i,x)
{
return _sa[i] + _sb[i] * (x - _sx[i]) + _sc[i] * (x - _sx[i])* (x - _sx[i]) + _sd[i] *  (x - _sx[i])* (x - _sx[i])* (x - _sx[i]) ;
}

/*-------------------------------------------------
interface with GBgraph
input adata
what = 2 : into adata
what = 1 : into bdata
-----------------------------------------------------*/
var nSlices = 10; // number of slices for spline interpolation
function dospline (what)
{
var i,j;
var n = adata.length;
var xdata = new Array(n);
var ydata = new Array(n);
var x,dx ;

for(i=0;i<n;i++) {xdata[i] = i; ydata[i] = ldata[i];}
ccsplines(xdata,ydata);

if((n-1)*nSlices <= 1000) 
	slices= nSlices;
	else slices = Math.floor (1000 / (n-1));

adata = new Array((n-1)*slices+1);
ldata = new Array((n-1)*slices+1);
if(what == 1)
{
odata = new Array((n-1)*slices+1);
bdata = new Array((n-1)*slices+1);


// linear interpol : need same steps numbers on a and b
for(i=0;i<n-1;i++)
{
dy = (ydata[i+1] -ydata[i]) / slices ;
for(j=0;j<slices;j++)
	{
	adata[slices*i+j] = ydata[i] + j * dy;
	ldata[slices*i+j] = scale(adata[slices*i+j]);
	}
}
	adata[(n-1)*slices] = ydata[n-1] ;
	ldata[(n-1)*slices] = scale(adata[(n-1)*slices]);
	
// spline interpolate
for(i=0;i<n-1;i++) 
for(j=0;j<slices;j++)
	{
	bdata[slices*i +j ] = spline(i,i + j/slices);
	odata[slices*i +j ] = scale(bdata[slices*i +j ]);
	}
        bdata[(n-1)*slices] = ydata[n-1] ;
        odata[(n-1)*slices] = scale(bdata[(n-1)*slices]);
} // what = 1

if(what == 2)
{
// spline interpolate
for(i=0;i<n-1;i++) 
for(j=0;j<slices;j++)
	{
	adata[slices*i +j ] = spline(i,i + j/slices);
	ldata[slices*i +j ] = scale(adata[slices*i +j ]);
	}
        adata[(n-1)*slices] = ydata[n-1] ;
        ldata[(n-1)*slices] = scale(adata[(n-1)*slices]);
 } // what == 2


nMin = 0;
nSteps = (n-1)*slices  ;
nMult= 1 ;
gInt = 0;
gFrame=1; // no scrolling

} // dospline

/*----------------------------------------
tangent (f,x0) ---> sets bdata, odata
uses nSteps
equation is y = y0 +(x-x0) * f'(x0)
            x in [nMin, nMin+nSteps*nMult]
---------------------------------------------*/
function tangent(f,x0)
	{
	var fp = deriv(f,x0);
	var y0 = f(x0);
	var x; var y;
	odata = new Array();
	bdata = new Array();
	
	for(var i= 0;i<= nSteps;i++)
	{
	x = nMin + i*nMult;
	y = y0 + (x-x0) * fp ;
	odata.push(y);
	bdata.push(y);
	}
wantsBn = 0 ; // autocomputed
gFrame=1;
} // tangent