[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
rat and rats ported from octave-forge
|
From: |
David Bateman |
|
Subject: |
rat and rats ported from octave-forge |
|
Date: |
Fri, 20 Jul 2007 12:34:46 +0200 |
|
User-agent: |
Thunderbird 1.5.0.7 (X11/20060921) |
I hesitated in the past to port the rat and rats functions from
octave-forge as we don't have the "format rat" command available. Also
the octave-forge rats command used a tolerance rather than a maximum
string length. Attached is a patch that ports the rat function from
octave-forge, implements the "format rat" command and reimplements rats
in terms of the "format rat" display code. Note that the results are
slightly different that matlab in a couple of cases..
1) Something like 1.233 will be represented as 1233/1000 in matlab,
whereas in octave it will be represented as 127/103
2) Matlab seems to stop expanding the representation earlier than it
needs to. For example in matlab 1.0014 is represented as 700/699, where
as in Octave 7003/6993 which equally well fits in the maximum string
length.
Regards
David
--
David Bateman address@hidden
Motorola Labs - Paris +33 1 69 35 48 04 (Ph)
Parc Les Algorithmes, Commune de St Aubin +33 6 72 01 06 33 (Mob)
91193 Gif-Sur-Yvette FRANCE +33 1 69 35 77 01 (Fax)
The information contained in this communication has been classified as:
[x] General Business Information
[ ] Motorola Internal Use Only
[ ] Motorola Confidential Proprietary
*** ./doc/interpreter/io.txi.orig49 2007-07-20 11:41:43.350214563 +0200
--- ./doc/interpreter/io.txi 2007-07-20 12:05:14.863613548 +0200
***************
*** 25,30 ****
--- 25,31 ----
* Terminal Output::
* Terminal Input::
* Simple File I/O::
+ * Rational Approximations::
@end menu
@node Terminal Output
***************
*** 213,218 ****
--- 214,225 ----
@DOCSTRING(octave_core_file_name)
+ @node Rational Approximations
+ @subsection Rational Approximations
+
+ @DOCSTRING(rat)
+
+ @DOCSTRING(rats)
@node C-Style I/O Functions
@section C-Style I/O Functions
*** ./scripts/general/Makefile.in.orig49 2007-07-20 11:41:52.399755004
+0200
--- ./scripts/general/Makefile.in 2007-07-20 11:42:05.352097121 +0200
***************
*** 28,34 ****
isdefinite.m isdir.m isequal.m isequalwithequalnans.m isscalar.m \
issquare.m issymmetric.m isvector.m logical.m logspace.m lookup.m mod.m \
nargchk.m nextpow2.m nthroot.m num2str.m perror.m pol2cart.m \
! polyarea.m postpad.m prepad.m quadl.m randperm.m rem.m \
repmat.m rot90.m rotdim.m shift.m shiftdim.m sortrows.m \
sph2cart.m strerror.m sub2ind.m trapz.m tril.m triu.m
--- 28,34 ----
isdefinite.m isdir.m isequal.m isequalwithequalnans.m isscalar.m \
issquare.m issymmetric.m isvector.m logical.m logspace.m lookup.m mod.m \
nargchk.m nextpow2.m nthroot.m num2str.m perror.m pol2cart.m \
! polyarea.m postpad.m prepad.m quadl.m randperm.m rat.m rem.m \
repmat.m rot90.m rotdim.m shift.m shiftdim.m sortrows.m \
sph2cart.m strerror.m sub2ind.m trapz.m tril.m triu.m
*** ./scripts/general/rat.m.orig49 2007-07-20 11:02:58.559702196 +0200
--- ./scripts/general/rat.m 2007-07-20 12:08:13.121477530 +0200
***************
*** 0 ****
--- 1,92 ----
+ ## Copyright (C) 2001 Paul Kienzle
+ ##
+ ## This file is part of Octave.
+ ##
+ ## Octave 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, or (at your option)
+ ## any later version.
+ ##
+ ## Octave 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 Octave; see the file COPYING. If not, write to the Free
+ ## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+ ## 02110-1301, USA.
+
+ ## -*- texinfo -*-
+ ## @deftypefn {Function File} address@hidden =} rat (@var{x}, @var{tol})
+ ## @deftypefnx {Function File} address@hidden, @var{d}] =} rat (@var{x},
@var{tol})
+ ##
+ ## Find a rational approximation to @var{x} within tolerance defined
+ ## by @var{tol} using a continued fraction expansion. E.g,
+ ##
+ ## @example
+ ## rat(pi) = 3 + 1/(7 + 1/16) = 355/113
+ ## rat(e) = 3 + 1/(-4 + 1/(2 + 1/(5 + 1/(-2 + 1/(-7))))) = 1457/536
+ ## @end example
+ ##
+ ## Called with two arguments returns the numerator and deniminator seperately
+ ## as two matrices.
+ ## @end deftypefn
+ ## @seealso{rats}
+
+ function [n,d] = rat(x,tol)
+
+ if (nargin != [1,2] || nargout != 2)
+ usage("[n,d] = rat(x,tol)");
+ endif
+
+ y = x(:);
+
+ ## replace inf with 0 while calculating ratios
+ y(isinf(y)) = 0;
+
+ ## default norm
+ 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));
+
+ ## use 1/0 for Inf
+ n(isinf(x)) = sign(x(isinf(x)));
+ d(isinf(x)) = 0;
+
+ endfunction
*** ./src/pr-output.cc.orig49 2007-07-19 17:26:53.282692954 +0200
--- ./src/pr-output.cc 2007-07-20 12:10:47.254644239 +0200
***************
*** 90,95 ****
--- 90,101 ----
// First char for > 0, second for < 0, third for == 0.
static std::string plus_format_chars = "+ ";
+ // TRUE means always print in a rational approximation
+ static bool rat_format = false;
+
+ // Used to force the length of the rational approximation string for Frats
+ static int rat_string_len = -1;
+
// TRUE means always print like dollars and cents.
static bool bank_format = false;
***************
*** 112,117 ****
--- 118,124 ----
static bool print_big_e = false;
class pr_formatted_float;
+ class pr_rational_float;
static int
current_output_max_field_width (void)
***************
*** 170,175 ****
--- 177,185 ----
friend std::ostream& operator << (std::ostream& os,
const pr_formatted_float& pff);
+ friend std::ostream& operator << (std::ostream& os,
+ const pr_rational_float& pff);
+
private:
// Field width. Zero means as wide as necessary.
***************
*** 221,226 ****
--- 231,351 ----
return os;
}
+ static inline std::string
+ rational_approx (double val, int len)
+ {
+ std::string s;
+
+ if (len <= 0)
+ len = 10;
+
+ if (xisinf (val))
+ s = "1/0";
+ else if (xisnan (val))
+ s = "0/0";
+ else if (val < INT_MIN || val > INT_MAX || D_NINT (val) == val)
+ {
+ std::ostringstream buf;
+ buf.flags (std::ios::fixed);
+ buf << std::setprecision (0) << xround(val);
+ s = buf.str ();
+ }
+ else
+ {
+ double lastn = 1.;
+ double lastd = 0.;
+ double n = xround (val);
+ double d = 1.;
+ double frac = val - n;
+ int m = 0;
+
+ std::ostringstream buf2;
+ buf2.flags (std::ios::fixed);
+ buf2 << std::setprecision (0) << static_cast<int>(n);
+ s = buf2.str();
+
+ while (1)
+ {
+ double flip = 1. / frac;
+ double step = xround (flip);
+ double nextn = n;
+ double nextd = d;
+ frac = flip - step;
+ n = n * step + lastn;
+ d = d * step + lastd;
+ lastn = nextn;
+ lastd = nextd;
+
+ std::ostringstream buf;
+ buf.flags (std::ios::fixed);
+ buf << std::setprecision (0) << static_cast<int>(n)
+ << "/" << static_cast<int>(d);
+ m++;
+
+ if (buf.str().length() >= static_cast<unsigned int>(len) && m > 1)
+ break;
+
+ s = buf.str();
+ if (m > 1000 || frac == 0.)
+ {
+ lastn = n;
+ lastd = d;
+ break;
+ }
+ }
+
+ if (lastd < 0.)
+ {
+ // Move sign to the top
+ lastd = - lastd;
+ lastn = - lastn;
+ std::ostringstream buf;
+ buf.flags (std::ios::fixed);
+ buf << std::setprecision (0) << static_cast<int>(lastn)
+ << "/" << static_cast<int>(lastd);
+ s = buf.str();
+ }
+ }
+
+ return s;
+ }
+
+ class
+ pr_rational_float
+ {
+ public:
+
+ const float_format& f;
+
+ double val;
+
+ pr_rational_float (const float_format& f_arg, double val_arg)
+ : f (f_arg), val (val_arg) { }
+ };
+
+ std::ostream&
+ operator << (std::ostream& os, const pr_rational_float& prf)
+ {
+ int fw = (rat_string_len > 0 ? rat_string_len : prf.f.fw);
+ std::string s = rational_approx (prf.val, fw);
+
+ if (fw >= 0)
+ os << std::setw (fw);
+
+ std::ios::fmtflags oflags =
+ os.flags (static_cast<std::ios::fmtflags>
+ (prf.f.fmt | prf.f.up | prf.f.sp));
+
+ if (fw > 0 && s.length() > static_cast<unsigned int>(fw))
+ os << "*";
+ else
+ os << s;
+
+ os.flags (oflags);
+
+ return os;
+ }
+
// Current format for real numbers and the real part of complex
// numbers.
static float_format *curr_real_fmt = 0;
***************
*** 299,305 ****
int ld, rd;
! if (bank_format)
{
fw = digits < 0 ? 4 : digits + 3;
if (inf_or_nan && fw < 4)
--- 424,435 ----
int ld, rd;
! if (rat_format)
! {
! fw = 0;
! rd = 0;
! }
! else if (bank_format)
{
fw = digits < 0 ? 4 : digits + 3;
if (inf_or_nan && fw < 4)
***************
*** 343,349 ****
fw = 4;
}
! if (! (bank_format || hex_format || bit_format)
&& (fw > Voutput_max_field_width || print_e || print_g))
{
if (print_g)
--- 473,479 ----
fw = 4;
}
! if (! (rat_format || bank_format || hex_format || bit_format)
&& (fw > Voutput_max_field_width || print_e || print_g))
{
if (print_g)
***************
*** 412,418 ****
int ld, rd;
! if (bank_format)
{
int digits = x_max > x_min ? x_max : x_min;
fw = digits <= 0 ? 4 : digits + 3;
--- 542,553 ----
int ld, rd;
! if (rat_format)
! {
! fw = 10;
! rd = 0;
! }
! else if (bank_format)
{
int digits = x_max > x_min ? x_max : x_min;
fw = digits <= 0 ? 4 : digits + 3;
***************
*** 483,489 ****
fw = 4;
}
! if (! (bank_format || hex_format || bit_format)
&& (print_e
|| print_g
|| (! Vfixed_point_format && fw > Voutput_max_field_width)))
--- 618,624 ----
fw = 4;
}
! if (! (rat_format || bank_format || hex_format || bit_format)
&& (print_e
|| print_g
|| (! Vfixed_point_format && fw > Voutput_max_field_width)))
***************
*** 564,570 ****
int ld, rd;
! if (bank_format)
{
int digits = r_x;
i_fw = 0;
--- 699,711 ----
int ld, rd;
! if (rat_format)
! {
! i_fw = 0;
! r_fw = 0;
! rd = 0;
! }
! else if (bank_format)
{
int digits = r_x;
i_fw = 0;
***************
*** 639,645 ****
}
}
! if (! (bank_format || hex_format || bit_format)
&& (r_fw > Voutput_max_field_width || print_e || print_g))
{
if (print_g)
--- 780,786 ----
}
}
! if (! (rat_format || bank_format || hex_format || bit_format)
&& (r_fw > Voutput_max_field_width || print_e || print_g))
{
if (print_g)
***************
*** 749,755 ****
int ld, rd;
! if (bank_format)
{
int digits = r_x_max > r_x_min ? r_x_max : r_x_min;
i_fw = 0;
--- 890,902 ----
int ld, rd;
! if (rat_format)
! {
! i_fw = 10;
! r_fw = 10;
! rd = 0;
! }
! else if (bank_format)
{
int digits = r_x_max > r_x_min ? r_x_max : r_x_min;
i_fw = 0;
***************
*** 835,841 ****
}
}
! if (! (bank_format || hex_format || bit_format)
&& (print_e
|| print_g
|| (! Vfixed_point_format && r_fw > Voutput_max_field_width)))
--- 982,988 ----
}
}
! if (! (rat_format || bank_format || hex_format || bit_format)
&& (print_e
|| print_g
|| (! Vfixed_point_format && r_fw > Voutput_max_field_width)))
***************
*** 949,955 ****
int ld, rd;
! if (bank_format)
{
int digits = x_max > x_min ? x_max : x_min;
fw = sign + digits < 0 ? 4 : digits + 3;
--- 1096,1107 ----
int ld, rd;
! if (rat_format)
! {
! fw = 10;
! rd = 0;
! }
! else if (bank_format)
{
int digits = x_max > x_min ? x_max : x_min;
fw = sign + digits < 0 ? 4 : digits + 3;
***************
*** 1012,1018 ****
fw = sign + 1 + ld + 1 + rd;
}
! if (! (bank_format || hex_format || bit_format)
&& (print_e
|| print_g
|| (! Vfixed_point_format && fw > Voutput_max_field_width)))
--- 1164,1170 ----
fw = sign + 1 + ld + 1 + rd;
}
! if (! (rat_format || bank_format || hex_format || bit_format)
&& (print_e
|| print_g
|| (! Vfixed_point_format && fw > Voutput_max_field_width)))
***************
*** 1211,1216 ****
--- 1363,1377 ----
}
}
}
+ else if (octave_is_NA (d))
+ {
+ if (fw > 0)
+ os << std::setw (fw) << "NA";
+ else
+ os << "NA";
+ }
+ else if (rat_format)
+ os << pr_rational_float (*fmt, d);
else if (xisinf (d))
{
const char *s;
***************
*** 1224,1236 ****
else
os << s;
}
- else if (octave_is_NA (d))
- {
- if (fw > 0)
- os << std::setw (fw) << "NA";
- else
- os << "NA";
- }
else if (xisnan (d))
{
if (fw > 0)
--- 1385,1390 ----
***************
*** 1695,1701 ****
double scale = 1.0;
set_format (cm, r_fw, i_fw, scale);
int column_width = i_fw + r_fw;
! column_width += (bank_format || hex_format|| bit_format) ? 2 : 7;
octave_idx_type total_width = nc * column_width;
octave_idx_type max_width = command_editor::terminal_cols ();
--- 1849,1856 ----
double scale = 1.0;
set_format (cm, r_fw, i_fw, scale);
int column_width = i_fw + r_fw;
! column_width += (rat_format || bank_format || hex_format
! || bit_format) ? 2 : 7;
octave_idx_type total_width = nc * column_width;
octave_idx_type max_width = command_editor::terminal_cols ();
***************
*** 2363,2369 ****
fw = digits + isneg;
}
! int column_width = fw + (bank_format ? 5 : 2);
octave_idx_type total_width = nc * column_width;
int max_width = command_editor::terminal_cols () - extra_indent;
octave_idx_type inc = nc;
--- 2518,2524 ----
fw = digits + isneg;
}
! int column_width = fw + (rat_format ? 0 : (bank_format ? 5 : 2));
octave_idx_type total_width = nc * column_width;
int max_width = command_editor::terminal_cols () - extra_indent;
octave_idx_type inc = nc;
***************
*** 2570,2575 ****
--- 2725,2771 ----
panic_impossible ();
}
+ DEFUN (rats, args, nargout,
+ "-*- texinfo -*-\n\
+ @deftypefn {Built-in Function} {} rats (@var{x}, @var{len})\n\
+ Convert @var{x} into a rational approximation represented as a string.\n\
+ You can convert the string back into a matrix as follows:\n\
+ \n\
+ @example\n\
+ eval(['[',rats(hilb(4)),'];'])\n\
+ @end example\n\
+ \n\
+ The optional second argument defines the maximum length of the string\n\
+ representing the elements of @var{x}. By default @var{len} is 13.\n\
+ @seealso{format, rat}\n\
+ @end deftypefn")
+ {
+ octave_value retval;
+ int nargin = args.length ();
+
+ rat_string_len = 10;
+ if (nargin == 2)
+ rat_string_len = args(1).nint_value () - 3;
+
+ if (!error_state)
+ {
+ if (nargin < 3 && nargout < 2)
+ {
+ bool save_rat_format = rat_format;
+ rat_format = true;
+ std::ostringstream buf;
+ args(0).print (buf);
+ retval = buf.str ();
+ rat_format = save_rat_format;
+ }
+ else
+ print_usage ();
+ }
+
+ rat_string_len = -1;
+ return retval;
+ }
+
DEFUN (disp, args, nargout,
"-*- texinfo -*-\n\
@deftypefn {Built-in Function} {} disp (@var{x})\n\
***************
*** 2659,2664 ****
--- 2855,2861 ----
{
free_format = false;
plus_format = false;
+ rat_format = false;
bank_format = false;
hex_format = 0;
bit_format = 0;
***************
*** 2804,2809 ****
--- 3001,3011 ----
init_format_state ();
plus_format = true;
}
+ else if (arg == "rat")
+ {
+ init_format_state ();
+ rat_format = true;
+ }
else if (arg == "bank")
{
init_format_state ();
***************
*** 2981,2986 ****
--- 3183,3191 ----
@item loose\n\
Insert blank lines above and below column number labels (this is the\n\
default).\n\
+ @item rat\n\
+ Print a rational approximation. That is the values are approximated\n\
+ by one small integer divided by another.\n\
@end table\n\
\n\
By default, Octave will try to print numbers with at least 5 significant\n\
2007-07-20 David Bateman <address@hidden>
* interpreter/io.txi: Document rat and rats in new sub-section.
2007-07-20 David Bateman <address@hidden>
* general/rat.m: New function for ration approximation imported
from octave-forge.
* general/Makefile.in: Include rat.m in SOURCES.
2007-07-20 David Bateman <address@hidden>
* pr-output.cc (rat_format, rat_string_len): Global variable
controlling behavior of rational approximation. Use throughout.
(class pr_rational_float): New class for rational approximation of
floats, specifically with the << operator defined.
(std::ostream& operator << (std::ostream&, const
pr_rational_float&)): Operator to print rational approximations of
double values.
(std::string rational_approx (double, int)): Function to convert a
double value to a string of maximum length giving the rational
approximation.
(pr_any_float): Include the output of rational approximations.
(Fformat): Add the "rat" format as an option.
(Frats): New function.
- rat and rats ported from octave-forge,
David Bateman <=