/*  -*- MaTX -*-
 *
 *  NAME
 *      ihilbert() - Inverse Hilbert matrix
 *
 *  SYNOPSIS
 *      H = ihilbert(n)
 *         Matrix H;
 *         Integer n;
 *
 *  DESCRIPTION
 *      ihilbert(n) returns the n-by-n inverse hilbert matrix.
 *
 *  SEE ALSO
 *      hilbert
 */

Func Matrix ihilbert(n)
	Integer n;
{
	Integer i, j;
	Real a, b;
	Matrix H;

	a = n;
	H = Z(n);

	for (i = 1; i <= n; i++) {
		if (i > 1) {
			a = ((n-i+1.0)*a*(n+i-1.0))/(i-1.0)^2;
		}
		b = a*a;
		H(i,i) = b/(2*i-1.0);
		for (j = i+1; j <= n; j++) {
			b = -((n-j+1.0)*b*(n+j-1.0))/(j-1.0)^2;
			H(i,j) = b/(i+j-1.0);
			H(j,i) = b/(i+j-1.0);
		}
	}
	return H;
}