hilbert.cpp

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