| 1 | #-----------------------------------------------------------------------------
|
| 2 | #
|
| 3 | # Create a Hilbert matrix
|
| 4 | #
|
| 5 | # Input: Dimension on stack
|
| 6 | #
|
| 7 | # Output: Hilbert matrix on stack
|
| 8 | #
|
| 9 | # Example:
|
| 10 | #
|
| 11 | # > hilbert(5)
|
| 12 | # ((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))
|
| 13 | #
|
| 14 | #-----------------------------------------------------------------------------
|
| 15 |
|
| 16 |
|
| 17 |
|
| 18 | #define A p1
|
| 19 | #define N p2
|
| 20 |
|
| 21 | #define AELEM(i, j) A->u.tensor->elem[i * n + j]
|
| 22 |
|
| 23 | hilbert = ->
|
| 24 | i = 0
|
| 25 | j = 0
|
| 26 | n = 0
|
| 27 | save()
|
| 28 | p2 = pop()
|
| 29 | push(p2)
|
| 30 | n = pop_integer()
|
| 31 | if (n < 2)
|
| 32 | push_symbol(HILBERT)
|
| 33 | push(p2)
|
| 34 | list(2)
|
| 35 | restore()
|
| 36 | return
|
| 37 | push_zero_matrix(n, n)
|
| 38 | p1 = pop()
|
| 39 | for i in [0...n]
|
| 40 | for j in [0...n]
|
| 41 | push_integer(i + j + 1)
|
| 42 | inverse()
|
| 43 | p1.tensor.elem[i * n + j] = pop()
|
| 44 | push(p1)
|
| 45 | restore()
|