Re: Fractions

[Paul Kienzle]

You can get a rational approximation to a number with the function
rat.m that I have written.  And you can display a matrix of fractions
as a string with the function rats.m.  See below.  

E.g.,

octave:1> rats(hilb(9))
ans = 
  1   1/2  1/3  1/4  1/5  1/6  1/7  1/8  1/9
 1/2  1/3  1/4  1/5  1/6  1/7  1/8  1/9 1/10
 1/3  1/4  1/5  1/6  1/7  1/8  1/9 1/10 1/11
 1/4  1/5  1/6  1/7  1/8  1/9 1/10 1/11 1/12
 1/5  1/6  1/7  1/8  1/9 1/10 1/11 1/12 1/13
 1/6  1/7  1/8  1/9 1/10 1/11 1/12 1/13 1/14
 1/7  1/8  1/9 1/10 1/11 1/12 1/13 1/14 1/15
 1/8  1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16
 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 

Any non-standard octave functions can be found in matcompat
        http://users.powernet.co.uk/kienzle/octave/matcompat

Hope this helps.

Paul Kienzle
address@hidden

On Thu, Mar 29, 2001 at 11:20:02PM +1000, Menaka Lashitha Bandara wrote:
> Hi people,
> 
> can Octave work purely with fractions???
> 
> I'm getting too many decimals in my work.
> 
> Thanks.
> 
> Lashi.
> 
> 
> 
> -------------------------------------------------------------
> Octave is freely available under the terms of the GNU GPL.
> 
> Octave's home on the web:  http://www.octave.org
> How to fund new projects:  http://www.octave.org/funding.html
> Subscription information:  http://www.octave.org/archive.html
> -------------------------------------------------------------
> 
> 

## Copyright (C) 2001 Paul Kienzle
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software
## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

## [n,d] = rat(x,tol)
## Find a rational approximation to x within tolerance using a continued
## fraction expansion. E.g,
##
##    rat(pi) = 3 + 1/(7 + 1/16) = 355/113
##    rat(e) = 3 + 1/(-4 + 1/(2 + 1/(5 + 1/(-2 + 1/(-7))))) = 1457/536

function [n,d] = rat(x,tol)

  if (nargin != [1,2] || nargout != 2)
    usage("[n,d] = rat(x,tol)");
  endif

  y = x(:);
  if (nargin < 2)
    tol = 1e-6 * norm(y,1);
  endif

  ## First step in the approximation is the integer portion
  n = round(y);  # first element in the continued fraction
  d = ones(size(y));
  frac = y-n;
  lastn = ones(size(y));
  lastd = zeros(size(y));

  ## grab new factors until all continued fractions converge
  while (1)
    ## determine which fractions have not yet converged
    idx = find (abs(y-n./d) >= tol);
    if (isempty(idx)) break; endif

    ## grab the next step in the continued fraction
    flip = 1./frac(idx);
    step = round(flip); # next element in the continued fraction
    frac(idx) = flip-step;

    ## update the numerator/denominator
    nextn = n;
    nextd = d;
    n(idx) = n(idx).*step + lastn(idx);
    d(idx) = d(idx).*step + lastd(idx);
    lastn = nextn;
    lastd = nextd;
  endwhile

  ## move the minus sign to the top
  n = n.*sign(d);
  d = abs(d);

  ## return the same shape as you receive
  n = reshape(n, size(x));
  d = reshape(d, size(x));

endfunction


opyright (C) 2001 Paul Kienzle
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software
## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

## S = rat(x,tol)
## Convert x into a rational approximation represented as a string.  You
## can convert the string back into a matrix as follows:
##
##    eval(["[",rats(hilb(4)),"];"])

function T = rats(x)
  if nargin != 1
    usage("S = rats(x)");
  endif

  [n, d] = rat(x');
  [nr, nc] = size(x);

  len = nr*nc;
  S = sprintf("%d/%d ", [n(:), d(:)]');
  if (len == 1)
    T = S;
  else
    index = findstr(S, "/1 ");
    if (index)
      S = toascii(S);
      S([index, index+1]) = [];
      S = setstr(S);
    endif
    index = find(S == " ");
    shift = [index(1), diff(index)];
    cellsize = max(shift);
    shift = cellsize - shift;
    assign = ones(size(S));
    index = index(1:len-1);
    assign([1,index]) = assign([1,index]) + ceil(shift/2);
    assign(index) = assign(index) + floor(shift(1:len-1)/2);
    assign = cumsum(assign);
    T = setstr(toascii(" ")*ones(1, nr*nc*cellsize+1));
    T(assign+1) = S;
    T(nc*cellsize+1:nc*cellsize:nr*nc*cellsize) = "\n";
    T(1)='\n';
  endif

endfunction



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

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