Code covered by the BSD License  

Highlights from
The Computation of Pi by Archimedes

image thumbnail

The Computation of Pi by Archimedes

by

 

23 Nov 2010 (Updated )

Archimedes wrote 3 1/7 > pi > 3 10/71. This is how he did it.

Editor's Notes:

This file was selected as MATLAB Central Pick of the Week

fr(varargin)
function FRAC = fr(varargin)
% fr: Creates a fraction object
% usage: F = fr
%        F = fr(K)     % K is a numeric variable
%        F = fr(N,D)   % N and D are numeric variables
%        F = fr(K,N,D) % K, N and D are numeric variables
%        F = fr(...)   % K, N and D are compatible objects (see below)
%
%
% Arguments:
%  K - the whole number part of the fraction
%  N - the numerator of the fraction
%  D - the denominator of the fraction
%
%  fr(K,N,D) returns a fraction object representing K + N/D.
%  K, N, D can be integers (double, single, int*) or any "compatible"
%  objects for which certain arithmetic and comparison operations are
%  defined (see below).  K, N, D are assumed to be 0, 0, 1 if not provided.
%
%
%  The fraction will be automatically reduced to lowest common form, N>0.
%  - non-integer floats N (single/double) are converted to fractions using
%    RAT however this has a default tolerance of 1e-6*abs(N).
%    For better accuracy it is recommended that N and D are precomputed
%    with rats using a smaller tolerance.
%  - Integers N are represented by N + 0 / 1
%  - Negative fractions -N/D are represented by -1 + (D-N) / D (if N<D)
%  - [-inf,nan,inf] are represented by 0 + [-1,0,1] / 0
%
%  The following object types are known to be compatible:
%    vpi (John D'Errico's Variable Precision Integer toolbox,
%         available on MATLAB File Exchange)
%
%  Non-standard objects must include 0, 1, -1 and require the following
%  operations to be defined in order to create a fraction object:
%    gcd
%    rem
%    sign
%    abs
%    +, - , .*, ./
%    ==, <, <=, >, >=, ~=
%  
%  The following additional operation definitions are recommended:
%    *, .^
%    sort
%    floor
%    factor
%    gcd (3-output form)
%    rat (if floor(x) or mod(x,1) is not always equal to x)
%
%  Examples:
%  f = fr      % creates a zero fraction (0+0/1)
%  f = fr([])  % creates an empty fraction object
%  f = fr(1,3) % creates a fraction representing 1/3 (0+1/3)
%  f = fr(1,vpi(3)^100) % creates a fraction representing 1/3^100
%                       % (requires VPI toolbox)
%
%  See also: double, single

% Author: Ben Petschel 25/7/2009
%
% Version history:
%   25/7/2009 - first release (using vpi/vpi as a template)
%   28/7/2009 - added checks to ensure fraction parts are all the same class


% process the input arguments
if nargin==0
  % creates a zero variable
  K = 0;
  N = 0;
  D = 1;
elseif nargin==1
  K=varargin{1};
  if isa(K,'fr'),
    % already a fraction object
    FRAC=K;
    return;
  else
    % convert to a fraction later
    N=0;
    D=1;
  end;
elseif nargin==2
  K=0;
  N=varargin{1};
  D=varargin{2};
elseif nargin==3
  K=varargin{1};
  N=varargin{2};
  D=varargin{3};
else
  error('fr:nargin','fr cannot have more than 3 input arguments');
end;


% create the fraction object
isfrac=false;
isreduced=true;
if isempty(N) || isempty(D) || isempty(K)
  % create an empty fraction
  FRAC.numer = 0;
  FRAC.denom = 1;
  FRAC.whole = 0;
  FRAC(1) = [];

elseif any([numel(N),numel(D),numel(K)]>1)
  % create an array: check the sizes agree and then replicate frac objects
  
  % check sizes agree and find common sizes
  sn=size(N);
  sd=size(D);
  sk=size(K);
  scell={sn,sd,sk};
  sf=scell{find(cellfun(@(x)any(x>1),scell),1,'first')}; % get first non-1x1 size
  sok=all(cellfun(@(x)(all(x==1)||isequal(x,sf)),scell)); % check all sizes are 1x1 or sf
  
  if sok,
    % sizes are ok, so proceed
    FRAC.numer = 0;
    FRAC.denom = 1;
    FRAC.whole = 0;
  
    FRAC = repmat(FRAC,sf);
    
    if all(sn==1)
      N=repmat(N,sf);
    end;
    if all(sd==1)
      D=repmat(D,sf);
    end;
    if all(sk==1)
      K=repmat(K,sf);
    end;
    for i=1:prod(sf)
      FRAC(i)=fr(K(i),N(i),D(i));
    end;
    
  else
    error('fr:inputsize','array sizes of N, D and K do not match');
  end;

else
  % All scalars
  
  % FK, FN, FD are only converted to fractions if they are are not integers
  [tfk,FK]=tofrac(K);
  [tfn,FN]=tofrac(N);
  [tfd,FD]=tofrac(D);
  if tfk || tfn || tfd,
    if ~tfn && ~tfd,
      % adding K part was fraction, N and D are integers
      if N==0 && D~=0,
        % N/D = 0
        FRAC = FK;
      else
        FRAC = FK + fr(N,D);
      end;
    else
      % N and/or D are fractions
      FRAC = FK + FN/FD;
    end;
    isfrac=true;
    
  else
    if (isa(N,'double') && (abs(N)>=(2^53))) || (isa(D,'double') && (abs(D)>=(2^53))) || (isa(K,'double') && (abs(K)>=(2^53)))
      warning('fr:fr:largedouble','Any N, D or K that are doubles must be no larger than 2^53 - 1, otherwise roundoff error will result.');
    end
    
    % ensure that N,D,K are the same type
    try
      zeroterm = (N-N)+(D-D)+(K-K);
    catch
      % types are incompatible
      error('fr:fr:inputclass','Classes of N,D,K are incompatible for addition');
    end;
    
    FRAC.numer = N+zeroterm;
    FRAC.denom = D+zeroterm;
    FRAC.whole = K+zeroterm;

    isreduced=false;
  end;
end

% set the class for this variable
if ~isfrac
  FRAC = class(FRAC,'fr');
end;

% reduce the fraction to lowest form
if ~isreduced
  FRAC = freduce(FRAC);
end;


% make sure that the appropriate fraction method is
% used whenever one of the arguments is a fraction
% and any numeric type as the other argument.
superiorto('double','single','int8','uint8','int16', ...
  'uint16','int32','uint32','int64','uint64','logical')
try
  superiorto('vpi')
catch
  % in case VPI toolbox has never been used
end;


end % main function fr(...)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [tf,F]=tofrac(X)
% determines whether the input is a whole number or not,
% and convert to a fraction.  tf returns true if it was converted to a
% fraction.

tf=false;

try
  if isnan(X)
    tf=true;
    F=fr(0,0,0);
    return
  elseif isinf(X)
    tf=true;
    F=fr(0,sign(X),0);
    return
  end;
catch
  % for objects where isnan or isinf is not defined assume value is finite
end;

% otherwise see if X has a remainder mod 1
try
  K=floor(X);
  rem=X-K;
catch
  % for object types that don't have a "floor" function
  rem=mod(X,1);
  K=X-rem;
end;

if rem==0,
  F=X;
else
  try
    % assume double, otherwise error will result
    [N,D]=rat(rem);
    tf=true;
    F=fr(K,N,D);
  catch
    warning('fr:convertfraction','input value has a remainder mod 1, however function "rat" is not defined for this object; rounding downwards instead');
    F=X;
  end;
end;

end % helper function iswhole(X)

Contact us