begin	comment	program  apv 300969  5238,
			berekeningen voor 9-uurs-waarden;
	integer	i,j,n,m;
	real		s,p,a,b,c;

	comment	SP 143, see descr 27;
	real procedure Cgamma2(n);value n;integer n;
	begin	real g;

		g:=1.77245;if n>=7 then go to L;if n%2*2=n then g:=1;
	K:	n:=n-2;if n<1 then go to M;g:=g*n/2;go to K;
	L:	g:=g*exp(-n/2)*n**(n/2-.5)*(1+(.166667+.0138889/n)/n)/2**(n/2-1);
	M:	Cgamma2:=g
	end;

	comment	SP 144, see descr 27;
	real procedure Qgamma2(n,x);value n,x;integer n;real x;
	begin	integer k;
		real a,s,t,u;

		s:=0;if x<=0 then go to out;a:=n/2;if n>=100 then go to WH;
		if n<=6 then go to LE;if x<=a then go to I else go to II;
	LE:	if n%2*2!=n then go to L1;t:=exp(-x);s:=1-t;if n=2 then go to out;
		t:=t*x;s:=s-t;if n=4 then go to out;s:=s-t*x/2;go to out;
	L1:	t:=sqrt(x);s:=1-2*PHI(-t*1.4142135624);if n=1 then go to out;
		t:=exp(-x)*t/.8862269255;s:=s-t;if n=3 then go to out;
		s:=s-t*x*2/3;go to out;
	I :	s:=t:=exp(-x)*x**a/a/Cgamma2(n);u:=x/(a+1);k:=1;
	L2:	if u<=.001 then go to out;t:=t*x/(a+k);s:=s+t;
		k:=k+1;u:=(x/(a+k))**k;go to L2;
	II:	t:=exp(-x)*x**(a-1)/Cgamma2(n);s:=1-t;u:=(a-1)/x;k:=1;
	L3:	if u<=.001 then go to out;t:=t*(a-k)/x;s:=s-t;
		k:=k+1;u:=u*(a-k)/x;go to L3;
	WH:	s:=PHI(((x/a)**(1/3)+1/a/9-1)*sqrt(9*a));
	out:	Qgamma2:=s
	end;

	comment	SP 135, see descr 31;
	real procedure Chisq(n,x);value n,x;integer n;real x;
	Chisq:= Qgamma2(n,x/2);

	comment 	mca 2310;

	boolean procedure zeroin(x, y, fx, tolx); real x, y, fx, tolx; 
	begin real a, fa, b, fb, c, fc, tol, m, p, q; 
	    a:= x; fa:= fx; b:= x:= y; fb:= fx; 
	interpolate: c:= a; fc:= fa; 
	extrapolate: if abs(fc) < abs(fb) then 
	    begin a:= b; fa:= fb; x:= b:= c; fb:= fc; c:= a; fc:= fa 
	    end interchange;
	    tol:= tolx; m:= (c + b) * .5; if abs(m - b) > tol then 
	    begin p:= (b - a) * fb; if p >= 0 then q:= fa - fb else 
		begin q:= fb - fa; p:= - p end;
		a:= b; fa:= fb; 
		x:= b:= if p <= abs(q) * tol then sign(c - b) * tol + b else 
		if p < (m - b) * q then p / q + b else m; fb:= fx; 
		goto if sign(fb) = sign(fc) then interpolate else 
		extrapolate
	    end;
	    y:= c; zeroin:= sign(fb) * sign(fc) <= 0
	end zeroin;