[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: rat and rats ported from octave-forge
|
From: |
David Bateman |
|
Subject: |
Re: rat and rats ported from octave-forge |
|
Date: |
Mon, 23 Jul 2007 14:08:54 +0200 |
|
User-agent: |
Thunderbird 1.5.0.7 (X11/20060921) |
David Bateman wrote:
> John W. Eaton wrote:
>
>> On 20-Jul-2007, David Bateman wrote:
>>
>> | 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.
>>
>> My guess is that someone will eventually complain that the results are
>> not precisely the same, but I'm not sure that matters at this point.
>>
>>
> For case 1) I think matlab is probably right to do what it does, I just
> have figured out a good way to duplicate it yet.. For point 2) matlab is
> not doing what it says it is doing in its documentation. It says it
> gives a string representation of the rational approximation upto a
> certain string length, where as it gives the first it finds of the
> maximum string length... I'd stick with the behavior I proposed..
>
Ok, I figured out what was wrong with case 1)... Basically the
fractional representation can have a negative numerator and denominator
and the signs are fixed at the end... For the representation of 1.233
this was the case and as "-1233/-1000" is longer than 9 characters it
was dropped as a valid representation. To get this right I had to check
the sign of the numerator and the denominator and check for string up to
two characters longer than actually desired.. The attached version of
the patch fixes this issue (and a small convergence test issue)..
Basically I think this version pretty much duplicates matlabs behavior
expect for point two above and that "rats" doesn't inlcude the leading
four whitespaces that matlab appends in the string length, as Octave
doesn't add them...
D.
--
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-22 20:59:52.613157119 +0200
***************
*** 1,4 ****
! @c Copyright (C) 1996, 1997 John W. Eaton
@c This is part of the Octave manual.
@c For copying conditions, see the file gpl.texi.
--- 1,4 ----
! @c Copyright (C) 1996, 1997, 2007 John W. Eaton
@c This is part of the Octave manual.
@c For copying conditions, see the file gpl.texi.
***************
*** 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 18:23:11.717055082 +0200
***************
*** 0 ****
--- 1,135 ----
+ ## 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)
+ print_usage ();
+ 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));
+
+ nd = ndims(y);
+ nsz = prod (size (y));
+ steps = zeros([nsz, 0]);
+
+ ## 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
+
+ if (nargout < 2)
+ tsteps = NaN .* ones (nsz, 1);
+ tsteps (idx) = step;
+ steps = [steps, tsteps];
+ endif
+
+ 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
+
+ if (nargout == 2)
+ ## 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;
+
+ ## reshape the output
+ n = reshape (n, size (x));
+ d = reshape (d, size (x));
+ else
+ n = "";
+ nsteps = size(steps, 2);
+ for i = 1: nsz
+ s = [int2str(y(i))," "];
+ j = 1;
+
+ while (true)
+ step = steps(i, j++);
+ if (isnan (step))
+ break;
+ endif
+ if (j > nsteps || isnan (steps(i, j)))
+ if (step < 0)
+ s = [s(1:end-1), " + 1/(", int2str(step), ")"];
+ else
+ s = [s(1:end-1), " + 1/", int2str(step)];
+ endif
+ break;
+ else
+ s = [s(1:end-1), " + 1/(", int2str(step), ")"];
+ endif
+ endwhile
+ s = [s, repmat(")", 1, j-2)];
+ n = cat (1, n, s);
+ endfor
+ endif
+
+ endfunction
*** ./src/pr-output.cc.orig49 2007-07-19 17:26:53.282692954 +0200
--- ./src/pr-output.cc 2007-07-23 13:59:04.389931274 +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,361 ----
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 (n < 0 && d < 0)
+ {
+ // Double negative, string can be two characters longer..
+ if (buf.str().length() > static_cast<unsigned int>(len + 2) &&
+ m > 1)
+ break;
+ }
+ else if (buf.str().length() > static_cast<unsigned int>(len) &&
+ m > 1)
+ break;
+
+ s = buf.str();
+
+ // Have we converged to 1/intmax ?
+ if (m > 100 || fabs (frac) < 1 / static_cast<double>(INT_MAX))
+ {
+ 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)
--- 434,445 ----
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)
--- 483,489 ----
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;
--- 552,563 ----
int ld, rd;
! if (rat_format)
! {
! fw = 9;
! 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)))
--- 628,634 ----
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;
--- 709,721 ----
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)
--- 790,796 ----
}
}
! 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;
--- 900,912 ----
int ld, rd;
! if (rat_format)
! {
! i_fw = 9;
! r_fw = 9;
! 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)))
--- 992,998 ----
}
}
! 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;
--- 1106,1117 ----
int ld, rd;
! if (rat_format)
! {
! fw = 9;
! 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)))
--- 1174,1180 ----
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 ****
--- 1373,1387 ----
}
}
}
+ 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)
--- 1395,1400 ----
***************
*** 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 ();
--- 1859,1866 ----
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;
--- 2528,2534 ----
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 ****
--- 2735,2781 ----
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 9.\n\
+ @seealso{format, rat}\n\
+ @end deftypefn")
+ {
+ octave_value retval;
+ int nargin = args.length ();
+
+ rat_string_len = 9;
+ if (nargin == 2)
+ rat_string_len = args(1).nint_value ();
+
+ 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 ****
--- 2865,2871 ----
{
free_format = false;
plus_format = false;
+ rat_format = false;
bank_format = false;
hex_format = 0;
bit_format = 0;
***************
*** 2804,2809 ****
--- 3011,3021 ----
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 ****
--- 3193,3201 ----
@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.
glpk.h, John W. Eaton, 2007/07/26