From octave-sources-request at bevo dot che dot wisc dot edu  Mon Aug 14 03:40:55 2000
Subject: invhilb.m
From: Dirk Laurie <dirk at calvyn dot puk dot ac dot za>
To: octave-sources at bevo dot che dot wisc dot edu
Date: Mon, 14 Aug 2000 10:40:42 +0200

Here is a replacement for invhilb.m with the following features:
  1. Fast (only one level of for-loop)
  2. Robust against overflow: only returns Inf when the entry
     concerned cannot be represented in IEEE double precision  
  3. Accurate (produces closest machine number to actual value)

Dirk Laurie

--------
% Inverse of Hilbert matrix
% Author: Dirk Laurie

function retval = invhilb (n)

  if (nargin != 1)
    usage ("invhilb (n)");
  endif

  nmax = length (n);
  if (nmax == 1)
    retval = zeros (n); 
    k=1:n; p=k.*bincoeff(k+n-1,k-1).*bincoeff(n,k);
    p(2:2:n)=-p(2:2:n);
    if n<203,
      for l = 1:n
        retval(l,:)=(p(l)*p)./(l:l+n-1); 
      endfor
    else
      for l = 1:n
        retval(l,:)=p(l)*(p./(l:l+n-1));
      endfor
    endif
  else
    error ("invhilb: expecting scalar argument, found something else");
  endif

endfunction
--------



-----------------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.che.wisc.edu/octave/octave.html
How to fund new projects:  http://www.che.wisc.edu/octave/funding.html
Subscription information:  http://www.che.wisc.edu/octave/archive.html
-----------------------------------------------------------------------