1) I enclose a small but hopefully useful submission but for theoretical
more than data analytic calculations. It is loosely modelled on a
similar facility in MATLAB.
rational(x) is like round(x,k) except that the rounding is done to
a nearby rational approximation to x rather than to a fixed
number of decimal places. It can remove round-off error in some
cases and is useful, like round(), for enhancing displays, etc.
fractions(x) returns the same result as round(x) but as a character
string in ordinary fraction form. fractions(0.6666667) = "2/3"
Useful for levels in categories and tables and for revealing
patterns in highly structured arrays with rational entries.
rat(x,...) is a low level function used by rational() and fractions().
2) The main submission is given as a shar file below. Rather than a
formal Makefile I have included a short shell script, `rat.make' that
does a similar job. It should work in most UNIX environments without
change. The README also includes instructions for hand installation,
if preferred.
Author's Name: Bill Venables.
Dr. W. N. Venables, | Phone: +61 8 228 5418
Department of Statistics, | Fax: +61 8 224 0464
University of Adelaide, | Telex: UNIVAD AA 89141
GPO Box 498, Adelaide, S.A. 5000 | ACSnet: wvenable@spam.ua.oz.au
Email address: wvenable@spam.ua.oz.au [ACSnet, Australia - may need
some research for your site.]
I am willing to help if there are any difficulties, but email to and
from Australia can be expensive and the service is rather flaky.
3) This tiny piece of code is PUBLIC DOMAIN and redistributions, changes,
etc are unrestricted. Formal copyright is maintained by the author to
inhibit the unlikely event of someone trying to sell it.
--
#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
# README
# rat.make
# rat.d
# rational.d
# fractions.d
# rat.fun
# rational.c
#
if test -f 'README'
then
echo shar: will not over-write existing file "'README'"
else
echo x - 'README'
sed 's/^X//' >'README' << 'SHAR_EOF'
XThese are a few small functions to find numerical rational approximations
Xusing an obvious continued fraction method.
X
XThe function rational(..) does a similar job to round(..) except that the
Xadjustment made is determined by the closeness of the number to a rational
Xnumber with a "small" denominator rather than by a fixed number of decimal
Xplaces. If the rational number is the correct value, accuracy is
Xactually enhanced.
X
Xfractions() finds the same approximation as rational() but returns a
Xcharacter array with the quantity in traditional fractional form. e.g.
Xfractions(0.66666667) = "2/3".
X
XBoth functions are useful to improve the readability of displays of the
Xresults of theoretical calculations where simple patterns in the answers
Xare obscured by round-off error.
X
Xrat() is a low level function used by both rational() and fractions().
X
X--
XTo set up, place the contents of the archive in your intended working
Xdirectory. Then look at the shell script `rat.make' and check that it
Xwill work on your system. No major changes are likely and if none are
Xneeded execute it by
X
X$ sh rat.make
X
X NOTE: this sh file assumes that the command 'S' runs new S.
X ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
X
XAlternatively install by hand by
X
Xa) Creating subdirectories .Data and .Data/.Help if need by,
X
Xb) Copying `rat.d' to .Data/.Help/rat (Note: no `.d') and similarly
X `rational.d' and `fractions.d'
X
Xc) Create an object file in the working directory (i.e. this one) by
X
X $ cc -O -c rational.c
X
Xd) Install the functions using
X
X $ S < rat.fun
X
Xe) Test it.
X
XOccasionally rational() or fractions() is required for large arrays to
Xenhance the readability of displays. Since the continued fraction
Xalgorithm must be applied to an array essentially element by element it
Xis much more efficient to do the main calculation directly in C rather
Xthan S. For this reason the function uses the dynamic loading facility to
Xload the object file `rational.o', which, as it stands, is assumed to
Xreside in the working directory. There is no difficulty in changing this
Xto a more convenient arrangement.
X
X---
XStandard disclaimers - you get what you paid for.
X
XCopyright W. Venables, August 1989.
XThis software is placed in the public domain.
X
XENQUIRIES TO: wvenable@spam.ua.oz.au
SHAR_EOF
if test 2308 -ne "`wc -c < 'README'`"
then
echo shar: error transmitting "'README'" '(should have been 2308 characters)'
fi
fi
if test -f 'rat.make'
then
echo shar: will not over-write existing file "'rat.make'"
else
echo x - 'rat.make'
sed 's/^X//' >'rat.make' << 'SHAR_EOF'
X#!/bin/sh
X# This is in lieu of a proper makefile.
X# Modify if necessary,
X# Execute the file with /bin/sh (not csh) to set up rat(..),
X# rational() and fractions
X#
Xexport PATH; PATH=/bin:$PATH
Xhome=`pwd`
X#
X# check that the .Data and .Data/.Help directories are in place
X#
Xif test ! -d '.Data'
Xthen
X echo 'rat.make: creating directory' "'.Data'"
X mkdir '.Data'
Xfi
Xecho 'rat.make: entering directory' "'.Data'"
Xcd '.Data'
Xif test ! -d '.Help'
Xthen
X echo 'rat.make: creating directory' "'.Help'"
X mkdir '.Help'
Xfi
Xecho shar: entering directory "'.Help'"
Xcd '.Help'
Xecho 'rat.make: setting up .Help files'
Xif test -f 'rat'
Xthen
X echo 'rat.make: over-writing existing file' "'.Data/.Help/rat'"
Xfi
Xcp ${home}/rat.d ./rat
Xif test -f 'rational'
Xthen
X echo 'rat.make: over-writing existing file' "'.Data/.Help/rational'"
Xfi
Xcp ${home}/rational.d ./rational
Xif test -f 'fractions'
Xthen
X echo 'rat.make: over-writing existing file' "'.Data/.Help/fractions'"
Xfi
Xcp ${home}/fractions.d ./fractions
Xcd ${home}
Xecho 'rat.make: re-entering directory' "`pwd`"
X#
X# Compile "rational.o"
X#
Xcc -O -c rational.c
X#
X# Install S functions
X#
Xecho 'This could take some time...'
XS 'rat.d' << 'SHAR_EOF'
X.BG
X.FN rat
X.TL
Xrat: calculate a rational approximation to `x' by continued fractions.
X.CS
Xrat(x, leng=6, maxm=100, name="rat.o")
X.AG x
Xnumeric data object for which the rational approximation is needed.
X.AG leng
Xmaximum length of the continued fraction used.
X.AG maxm
Xmaximum partial denominator. If any partial denominator exceeds
X`maxm' the continued fraction terminates at that point.
X.AG name
Xname of object file containing the dynamically loaded C function.
X.RT
Xlist of two components, $a and $b, containing the numerators and the
Xdenominators of the rational approximations respectively.
X
X`rat()' can be viewed as providing a speculative way of removing
Xroundoff error if the correct answer should be rational numbers
Xwith "small" denominators.
X
X.EX
X# solve() accuracy check using a small Hilbert Matrix
X x <- matrix(0,5,5); x <- 1/(row(x) + col(x) -1)
X# confirm:
X rat(x)
X
X x1 <- solve(solve(x)); x2 <- rat(x1)
X# compare:
X x-x1
X# with:
X x-x2$num/x2$den
X
X.KW rational approximation, fractions
X.WR
X
XSee also `rational()' and `fractions()' which use `rat()'.
X
SHAR_EOF
if test 1080 -ne "`wc -c < 'rat.d'`"
then
echo shar: error transmitting "'rat.d'" '(should have been 1080 characters)'
fi
fi
if test -f 'rational.d'
then
echo shar: will not over-write existing file "'rational.d'"
else
echo x - 'rational.d'
sed 's/^X//' >'rational.d' << 'SHAR_EOF'
X.BG
X.FN rational
X.TL
Xrational: rational approximation to real numbers.
X.CS
Xrational(x, ...)
X.AG x
XNumeric array to be approximated.
X.AG ...
XOptional additional parameters to be passed to rat(x,...)
X.RT
XNumeric array of the same shape as x of rational approximations.
X
XIf x contains the results of a calculation for which the correct
Xanswers rational with "small" denominators, this may remove
Xroundoff error.
X
XSimilar to round(..), except that full accuracy is maintained for
Xmany exactly rational quantities
X
X.EX
X# solve() accuracy check using a small Hilbert Matrix
X x <- matrix(0,5,5); x <- 1/(row(x) + col(x) - 1); rat(x)
X
X x1 <- solve(solve(x))
X#compare:
X x-x1
X#with:
X x-rational(x1)
X
X.KW rational approximation, fractions
X.WR
X
XSee also `fractions()' and the primitive function `rat()'.
SHAR_EOF
if test 804 -ne "`wc -c < 'rational.d'`"
then
echo shar: error transmitting "'rational.d'" '(should have been 804 characters)'
fi
fi
if test -f 'fractions.d'
then
echo shar: will not over-write existing file "'fractions.d'"
else
echo x - 'fractions.d'
sed 's/^X//' >'fractions.d' << 'SHAR_EOF'
X.BG
X.FN fractions
X.TL
Xfractions: rational approximations to x in fraction form
X.CS
Xfractions(x, ...)
X.AG x
XArray to be approximated.
X.AG ...
XOptional additional parameters to be passed to rat(x,...)
X.RT
XA character array of the same shape as x of rational approximations in
Xfraction form. e.g 0.28571429 -> "2/7".
X
X.EX
X# somewhat pointless..
X
X x <- rep(1:7,1:7)/7
X
X#compare:
X table(x)
X
X#with:
X table(fractions(x))
X
X
X.KW rational approximation, fractions
X.WR
X
XSee also rational(..) and the primitive function rat(..).
SHAR_EOF
if test 527 -ne "`wc -c < 'fractions.d'`"
then
echo shar: error transmitting "'fractions.d'" '(should have been 527 characters)'
fi
fi
if test -f 'rat.fun'
then
echo shar: will not over-write existing file "'rat.fun'"
else
echo x - 'rat.fun'
sed 's/^X//' >'rat.fun' << 'SHAR_EOF'
X"rat"<-
Xfunction(x, leng = 6, maxm = 100, name = "rational.o")
X{
X if(!exists(".RAT.IS.LOADED", frame = 0)) {
X dyn.load(name)
X assign(".RAT.IS.LOADED", 1, frame = 0)
X }
X hold <- attributes(x)
X i <- is.na(x)
X if(any(i))
X x[i] <- 0
X x <- as.double(x)
X storage.mode(x) <- "double"
X num <- den <- x
X h <- .C("rational",
X x,
X as.integer(length(x)),
X as.integer(leng),
X as.integer(maxm),
X num = num,
X den = den)
X if(any(i)) {
X h$den[i] <- NA
X h$num[i] <- NA
X }
X attributes(h$num) <- hold
X attributes(h$den) <- hold
X list(num = h$num, den = h$den)
X}
X"rational"<-
Xfunction(x, ...)
X{
X y <- rat(x, ...)
X y$num/y$den
X}
X"fractions"<-
Xfunction(x, ...)
X{
X hold <- attributes(x)
X y <- rat(x, ...)
X denc <- ifelse(y$den == 1, "", paste("/", y$den, sep = ""))
X y <- paste(y$num, denc, sep = "")
X attributes(y) <- hold
X y
X}
SHAR_EOF
if test 816 -ne "`wc -c < 'rat.fun'`"
then
echo shar: error transmitting "'rat.fun'" '(should have been 816 characters)'
fi
fi
if test -f 'rational.c'
then
echo shar: will not over-write existing file "'rational.c'"
else
echo x - 'rational.c'
sed 's/^X//' >'rational.c' << 'SHAR_EOF'
X#include
Xextern char *S_alloc();
X
Xvoid rational(x, n, length, maxdenom, a, b)
Xlong *n, *length, *maxdenom;
Xdouble *x, *a, *b;
X/*
X * a[i]/b[i] gives a continued fraction approximation to x[i]
X * using up to `length' terms, each less than `maxdenom'
X */
X{
Xlong m = *n, k, len = *length;
Xdouble bot, top, tmp, *d = (double *) S_alloc((unsigned) len, sizeof(double)),
X maxd = (double) *maxdenom;
X while (m--)
X {
X tmp = *x++;
X for (k = 0; k < len; k++)
X {
X d[k] = floor(tmp + 0.5e0);
X tmp = tmp - d[k];
X if (fabs(tmp) * maxd <= 1.0e0)
X break;
X tmp = 1.0e0 / tmp;
X };
X top = d[k];
X bot = 1.0e0;
X while (k-- > 0)
X {
X tmp = top;
X top = d[k] * top + bot;
X bot = tmp;
X };
X if (bot < 0.0e0)
X {
X top = -top;
X bot = -bot;
X }
X *a++ = top;
X *b++ = bot;
X }
X}
SHAR_EOF
if test 991 -ne "`wc -c < 'rational.c'`"
then
echo shar: error transmitting "'rational.c'" '(should have been 991 characters)'
fi
fi
echo Done
exit 0