#-----------------------------------------------------------------------------
#
#	Create a Hilbert matrix
#
#	Input:		Dimension on stack
#
#	Output:		Hilbert matrix on stack
#
#	Example:
#
#	> hilbert(5)
#	((1,1/2,1/3,1/4),(1/2,1/3,1/4,1/5),(1/3,1/4,1/5,1/6),(1/4,1/5,1/6,1/7))
#
#-----------------------------------------------------------------------------



#define A p1
#define N p2

#define AELEM(i, j) A->u.tensor->elem[i * n + j]

hilbert = ->
	i = 0
	j = 0
	n = 0
	save()
	p2 = pop()
	push(p2)
	n = pop_integer()
	if (n < 2)
		push_symbol(HILBERT)
		push(p2)
		list(2)
		restore()
		return
	push_zero_matrix(n, n)
	p1 = pop()
	for i in [0...n]
		for j in [0...n]
			push_integer(i + j + 1)
			inverse()
			p1.tensor.elem[i * n + j] = pop()
	push(p1)
	restore()