Download All

Code covered by the BSD License  

Highlights from
Fractions Toolbox

from Fractions Toolbox by Ben Petschel
create and manipulate fractions (K+N/D) using exact arithmetic

Fractions Toolbox

Fractions Toolbox

This demo file shows how one might use the Fractions Toolbox

The code was inspired by, and to a large extent based upon, John D'Errico's Variable Precision Integer Toolbox.

Fractions are stored in the form K + N / D where N>=0 and D>=0.

The toolbox was designed to be very flexible in the types of objects that can be used, for example doubles, integers and vpi numbers. The fraction components can in principle be any objects provided certain minimal requirements are satisfied.

Author: Ben Petschel 30/7/2009

Contents

The creator for fraction objects is fr

help fr

  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

Fractions are created by providing the numerator and denominator

FRAC = fr

FRAC =
   0

FRAC = fr(1)

FRAC =
   1

FRAC = fr(1,3)

FRAC =
   1 / 3

FRAC = fr(4,3)

FRAC =
   1 + 1 / 3

FRAC = fr(1,1,3)

FRAC =
   1 + 1 / 3

FRAC = fr(-4,3)

FRAC =
   -2 + 2 / 3

Non-integers are automatically converted (including NaN and Inf)

keep in mind the default 1e-6 relative accuracy of rat

FRAC = fr(1/3)

FRAC =
   1 / 3

FRAC = fr(pi)

FRAC =
   3 + 4703 / 33215

FRAC = fr(eps)

FRAC =
   1 / 4503599627370496

Fractions can be infinite

FRAC = fr(NaN)

FRAC =
   NaN

FRAC = fr(0,0)

FRAC =
   NaN

FRAC = fr(1,0)

FRAC =
   Inf

FRAC = fr(Inf)

FRAC =
   Inf

FRAC = fr(-1,0)

FRAC =
  -Inf

FRAC = fr(-Inf)

FRAC =
  -Inf

For very large numerators or denominators, use vpi numbers

doubles can only represent integers exactly up to 2^53-1. If you have John D'Errico's Variable Precision Toolbox, use vpi numbers. Although fr takes care to avoid integer overflow whenever possible, it is recommended that you use vpi numbers. 

FRAC = fr(1,2^100)

Warning: Any N, D or K that are doubles must be no larger than 2^53 - 1, otherwise roundoff error will result. 
Warning: one of the fraction components is a double exceeding bitmax; results may be inaccurate - use vpi-based fractions
instead 
FRAC =
   1 / 1267650600228229400000000000000

try
  FRAC = fr(1,vpi(2)^100)
catch
  % need VPI toolbox
end;

FRAC =
   1 / 1267650600228229401496703205376

Arithmetic manipulation of fractions

FRAC = fr(1,3)

FRAC =
   1 / 3

FRAC + 1

ans =
   1 + 1 / 3

FRAC + fr(2,3)

ans =
   1

FRAC*3

ans =
   1

FRAC/4

ans =
   1 / 12

FRAC^2

ans =
   1 / 9

Work with arrays of fractions

A = fr(1,1:10)

A =
 1x10 fraction array with elements (reading down columns)
   1
   1 / 2
   1 / 3
   1 / 4
   1 / 5
   1 / 6
   1 / 7
   1 / 8
   1 / 9
   1 / 10

B = A.^2

B =
 1x10 fraction array with elements (reading down columns)
   1
   1 / 4
   1 / 9
   1 / 16
   1 / 25
   1 / 36
   1 / 49
   1 / 64
   1 / 81
   1 / 100

Compute the sum

sum(A)

ans =
   2 + 2341 / 2520

Or the product

prod(B)

ans =
   1 / 13168189440000

Convolution

conv(A,B)

ans =
 1x19 fraction array with elements (reading down columns)
   1
   3 / 4
   41 / 72
   65 / 144
   8009 / 21600
   1127 / 3600
   190513 / 705600
   167101 / 705600
   13371157 / 63504000
   240427 / 1270080
   1106573 / 15240960
   10153 / 254016
   438401 / 17781120
   1517161 / 95256000
   529817 / 50803200
   42857 / 6350400
   9761 / 2332800
   19 / 8100
   1 / 1000

Using vpi fractions for larger arrays is recommended

%sum(fr(1,1:100))
try
  sum(fr(vpi(1),1:100))
catch
  % need VPI toolbox
end;

ans =
   5 + 522561233577855727314756256041670736351 / 2788815009188499086581352357412492142272

Relational operators

All of the standard operators are provided, <, >, <=, >=, ==, ~=

a=fr(1,2)

a =
   1 / 2

b=fr(2,3)

b =
   2 / 3

1/(a*b) == 3

ans =

     1


a > b

ans =

     0


a <= 1

ans =

     1

Concatenation using []

[3 2]*[a,b;b,a]*[a^3 b^4]'

ans =
   409 / 432

Conversion from fraction back to double

double(fr(1,3))

ans =

    0.3333

Extracting the numerator and denominator

[n,d]=rat(fr(4,3))

n =

     4


d =

     3


[k,n,d]=rat(fr(4,3))

k =

     1


n =

     1


d =

     3

Extracting digits after the decimal point (including repeating sequence)

base 10: 1/7=0.142857(142857)

[dig,rep]=digits(fr(1,7))

dig =

     1     4     2     8     5     7


rep =

     1     4     2     8     5     7

base 7: 1/7=0.1(0)

[dig,rep]=digits(fr(1,7),[],7)

dig =

     1


rep =

     0

Extracting a fixed number of digits after the decimal point

base 10: 1/7=0.142+(3/3500)

[dig,rep]=digits(fr(1,7),3)

dig =

     1     4     2

rep =
   3 / 3500

base 7: 1/7=0.100+(0)

[dig,rep]=digits(fr(1,7),3,7)

dig =

     1     0     0

rep =
   0

fr(142,1000)+fr(3,3500)

ans =
   1 / 7

Partial fractions

partial expansion of 1/12 as sum(a/p^k) where 0<a<p

FPART=partial(fr(1,12))

FPART =
 1x4 fraction array with elements (reading down columns)
   -1
   1 / 2
   1 / 4
   1 / 3

sum(FPART)

ans =
   1 / 12

partial expansion of 1/12 as sum(a/p^k) where 0<a<p^k

FPART=partial(fr(1,12),false)

FPART =
 1x3 fraction array with elements (reading down columns)
   -1
   3 / 4
   1 / 3

sum(FPART)

ans =
   1 / 12

Continued fraction expansions and best rational approximations

continued fraction representation of a fraction:

cf=cfrac(fr(3,8))

cf =

     0     2     1     2

continued fraction expansion of a square root (including repeating term)

[cf,rep]=cfracsqrt(fr(13))

cf =

     3


rep =

     1     1     1     1     6

best rational approximations with a given denominator limit

[r1,r2] = bestrat(fr(cf),rep,100)

r1 =
   3 + 23 / 38
r2 =
   3 + 43 / 71

Factoring numerator and denominator

F=factor(fr(35,24))

F =
 1x6 fraction array with elements (reading down columns)
   5
   7
   1 / 2
   1 / 2
   1 / 2
   1 / 3

prod(F)

ans =
   1 + 11 / 24

Solving linear equations

use of fractions eliminates many issues due to machine precision

A=fr(ones(2),[2,3;5,7]);
B=fr(ones(2,1),[11;13]);

rref([A,B])

ans =
 2x3 fraction array with elements (reading down columns)
   1
   0
   0
   1
   -3 + 49 / 143
   4 + 37 / 143

A\B

ans =
 2x1 fraction array with elements (reading down columns)
   -3 + 49 / 143
   4 + 37 / 143

double(A)\double(B)-double(A\B)

ans =

  1.0e-014 *

   -0.1332
    0.2665

Solving rank-deficient linear equations

fr/mldivide ("\") can handle singular and non-square A

A=[1,1;0,0];
B=[1;0];

the built-in version does not always find a particular solution

A\B

Warning: Matrix is singular to working precision. 

ans =

   NaN
   NaN

the 1-output version of "\" gives NaN values to indicate when the system has multiple solutions however substituting 0 for NaN still gives a solution 

X=fr(A)\fr(B)
X(isnan(X))=0;

Warning: Coefficient matrix is singular and the system has indeterminate components.  Recommend using [X,N]=A\B 
X =
 2x1 fraction array with elements (reading down columns)
   1
   NaN

A*X-B

ans =
 2x1 fraction array with elements (reading down columns)
   0
   0

the 2-output version of "\" provides the general solution X+N*V where V is an arbitrary array or fraction array of compatible size 

[X,N]=fr(A)\fr(B)
V=fr(5,13);

X =
 2x1 fraction array with elements (reading down columns)
   1
   0
N =
 2x1 fraction array with elements (reading down columns)
   -1
   1

A*(X+N*V)-B

ans =
 2x1 fraction array with elements (reading down columns)
   0
   0

when the equations have no solution, "\" returns an empty matrix and a warning (this differs from the built-in version) e.g. {x=0,x=1} 

fr([1;1])\[0;1]

Warning: Coefficient matrix is singular and the system has inconsistent components.  Recommend using least squares
X=lsq(A,B) 

ans =

   Empty matrix: 1-by-0

to find the least squares solution, use lsq:

lsq(fr([1;1]),[0;1])

ans =
   1 / 2

See the help files for more details


Published with MATLAB® 7.8

Contact us at files@mathworks.com