52 lines
853 B
C++
52 lines
853 B
C++
#include "stdafx.h"
|
|
|
|
//-----------------------------------------------------------------------------
|
|
//
|
|
// 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))
|
|
//
|
|
//-----------------------------------------------------------------------------
|
|
|
|
#include "defs.h"
|
|
|
|
#define A p1
|
|
#define N p2
|
|
|
|
#define AELEM(i, j) A->u.tensor->elem[i * n + j]
|
|
|
|
void
|
|
hilbert(void)
|
|
{
|
|
int i, j, n;
|
|
save();
|
|
N = pop();
|
|
push(N);
|
|
n = pop_integer();
|
|
if (n < 2) {
|
|
push_symbol(HILBERT);
|
|
push(N);
|
|
list(2);
|
|
restore();
|
|
return;
|
|
}
|
|
push_zero_matrix(n, n);
|
|
A = pop();
|
|
for (i = 0; i < n; i++) {
|
|
for (j = 0; j < n; j++) {
|
|
push_integer(i + j + 1);
|
|
inverse();
|
|
AELEM(i, j) = pop();
|
|
}
|
|
}
|
|
push(A);
|
|
restore();
|
|
}
|