' ******************************************************************
' This program solves a series of Hilbert systems of increasing
' order by the method of Gauss-Jordan.
'
' Hilbert systems have ill-conditioned matrices (i.e. with
' determinants close to zero), so that the matrix is considered
' singular for an order which depends on the numerical precision
' of the software.
' ******************************************************************

#INCLUDE "fbmath.bi"

SUB Hilbert(N AS INTEGER, A() AS DOUBLE, B() AS DOUBLE)
' ------------------------------------------------------------------
' Setup of the Hilbert system of order N
' ------------------------------------------------------------------
' Matrix A is the system matrix:
'       (                                             )
'       ( 1      1/2     1/3     1/4     ... 1/N      )
'       ( 1/2    1/3     1/4     1/5     ... 1/(N+1)  )
'   A = ( 1/3    1/4     1/5     1/6     ... 1/(N+2)  )
'       ( ........................................... )
'       ( 1/N    1/(N+1) 1/(N+2) 1/(N+3) ... 1/(2N-1) )
'       (                                             )
'
' Vector B is the vector of constant terms:
'
'           N
'   B(i) = Sum  A(i,j)
'          j=1
'
' The solution vector is [1 1 1 ... 1]
' ------------------------------------------------------------------

  DIM AS INTEGER I, J

  REDIM AS DOUBLE A(1 TO N, 1 TO N), B(1 TO N)

  A(1, 1) = 1
  FOR J = 2 TO N                 ' First row of matrix
    A(1, J) = 1 / J
  NEXT J

  FOR I = 2 TO N
    A(I, N) = 1 / (N + I - 1)    ' Last column of matrix
    FOR J = 1 TO N - 1
      A(I, J) = A(I - 1, J + 1)  ' Fill matrix
    NEXT J
  NEXT I

  FOR I = 1 TO N                 ' Constant terms
    B(I) = 0
    FOR J = 1 TO N
      B(I) = B(I) + A(I, J)
    NEXT J
  NEXT I
END SUB

SUB WriteSolution(ErrCode AS INTEGER, B() AS DOUBLE, Det AS DOUBLE)
' ------------------------------------------------------------------
' Result output
' ------------------------------------------------------------------

  DIM AS INTEGER I, N

  N = UBOUND(B)

  CLS
  PRINT "Hilbert system of order"; N; " - ";

  IF ErrCode = MatSing THEN
    PRINT "Singular matrix !"
  ELSE
    PRINT "Solution vector :"
    PRINT
    FOR I = 1 TO N
      PRINT USING "#####.######"; B(I)
    NEXT I
    PRINT
    PRINT "Determinant = "; Det
  END IF

  LOCATE 25, 1
  PRINT "Press a key ...";
  SLEEP
END SUB

DIM AS INTEGER N        ' Order of the system
DIM AS DOUBLE  Det      ' Determinant
DIM AS INTEGER ErrCode  ' Error code

N = 1
ErrCode = MatOk

' Initialize matrix and vector
DIM AS DOUBLE A(1 TO N, 1 TO N), B(1 TO N)

' Main loop
DO WHILE ErrCode = MatOk
  N = N + 1
  Hilbert N, A(), B()              ' Generates Hilbert system
  LinEq A(), B(), Det              ' Solves Hilbert system
  ErrCode = MathErr()              ' Get error code
  WriteSolution ErrCode, B(), Det  ' Writes solution
LOOP                               ' Until singular matrix found
