MATLAB Central Newsreader - robust truncation at n decimal places

robust truncation at n decimal places

Felipe G. Nievinski — Sun, 07 Dec 2008 03:42:44 -0500

Hi. Given
    x = 0.99899998
    y = 0.99900000
    za = sscanf(sprintf('%.7f\n',x), '%f')
    zb = fix(x./10^-7).*10^-7
    za == y
    zb == y
then za == y is true, but zb == y is false.
As you see, it's possible to perform a robust truncation at n decimal
places through conversion to/from character string. I was trying to
achive the same using computationally cheaper double precision
arithmetic, without success. Is that it or am I missing something?

Thanks,
Felipe.

Re: robust truncation at n decimal places

Roger Stafford — Sun, 07 Dec 2008 05:53:02 -0500

"Felipe G. Nievinski" <fgnievinski@gmail.com> wrote in message <545015a5-3f51-41e3-a40b-1cf7db531b11@35g2000pry.googlegroups.com>...
> Hi. Given
> x = 0.99899998
> y = 0.99900000
> za = sscanf(sprintf('%.7f\n',x), '%f')
> zb = fix(x./10^-7).*10^-7
> za == y
> zb == y
> then za == y is true, but zb == y is false.
> As you see, it's possible to perform a robust truncation at n decimal
> places through conversion to/from character string. I was trying to
> achive the same using computationally cheaper double precision
> arithmetic, without success. Is that it or am I missing something?
>
> Thanks,
> Felipe.

Your 'fix' is doing truncation - that is to say, rounding toward zero - whereas sprintf('%.7f\n',x) is doing a round to nearest. It's not surprising you get different answers. Also for better accuracy you should be multiplying by 10^7 instead of dividing by 10^-7 and then dividing by 10^7 instead of multiplying by 10^-7, because the number 10^-7 cannot be represented precisely with binary floating point numbers.

Nevertheless even with these corrections, you probably cannot count on always getting precisely the same answers clear out to the 53rd binary place unless you know exactly how sprintf's and sscanf's algorithms work. You should be making comparison like abs(za-zb) < tol where tol is something like 10^-15*abs(x)

Roger Stafford

Re: robust truncation at n decimal places

Felipe G. Nievinski — Sun, 07 Dec 2008 09:26:53 -0500

On Dec 6, 10:53 pm, "Roger Stafford"
<ellieandrogerxy...@mindspring.com.invalid> wrote:
> "Felipe G. Nievinski" <fgnievin...@gmail.com> wrote in message <545015a5-=
3f51-41e3-a40b-1cf7db531...@35g2000pry.googlegroups.com>...
>
> > Hi. Given
> > x =3D 0.99899998
> > y =3D 0.99900000
> > za =3D sscanf(sprintf('%.7f\n',x), '%f')
> > zb =3D fix(x./10^-7).*10^-7
> > za =3D=3D y
> > zb =3D=3D y
> > then za =3D=3D y is true, but zb =3D=3D y is false.
> > As you see, it's possible to perform a robust truncation at n decimal
> > places through conversion to/from character string. I was trying to
> > achive the same using computationally cheaper double precision
> > arithmetic, without success. Is that it or am I missing something?
>
> > Thanks,
> > Felipe.
>
> Your 'fix' is doing truncation - that is to say, rounding toward zero -=
whereas sprintf('%.7f\n',x) is doing a round to nearest. It's not surpris=
ing you get different answers. Also for better accuracy you should be mult=
iplying by 10^7 instead of dividing by 10^-7 and then dividing by 10^7 inst=
ead of multiplying by 10^-7, because the number 10^-7 cannot be represented=
precisely with binary floating point numbers.
>
Thanks for indicating the round-like behavior of sprintf -- I was
overlooking that. Now I get the expected result for the particular x
value I gave:
    x =3D 0.99899998
    y =3D 0.99900000
    za =3D sscanf(sprintf('%.7f\n',x), '%f')
    zb =3D round(x./10^-7).*10^-7
    za =3D=3D y
    zb =3D=3D y

Then I ran the code for another test cases.
    x =3D 257104079.99899998
    y =3D 257104079.9990000
    za =3D sscanf(sprintf('%.7f\n',x), '%f')
    zb =3D round(x.*10^7)./10^7
    za =3D=3D y
    zb =3D=3D y
The test above would only pass after I've replaced x*10^-7 by x/(10^7)
-- it's interesting that the two operations are not equivalent.

But for the test below the algorithm fails:
    x =3D 257098540.00000003
    y =3D 257098540.0000000
    za =3D sscanf(sprintf('%.7f\n',x), '%f')
    zb =3D round(x.*10^7)./10^7
    za =3D=3D y
    zb =3D=3D y
:-(

> Nevertheless even with these corrections, you probably cannot count on =
always getting precisely the same answers clear out to the 53rd binary plac=
e unless you know exactly how sprintf's and sscanf's algorithms work.
Actually, sprintf can be dismissed in the comparisons above, because I
know the answer I expect, e.g.,
    in =3D 257098540.00000003
    out =3D 257098540.0000000
I've only included sprintf to demonstrate that what I need is not
impossible, even in a finite precision computer. You're right that if
sprintf can do it, there must be a way to do it -- I hope without
going to/from character strings.

Maybe dissecting the floating-point number into exponent and mantissa?
I tried this, using [f,e]=3Dlog2(x), where x =3D f * 2^e,:
    x =3D 257098540.00000003
    y =3D 257098540.0000000
    za =3D sscanf(sprintf('%.7f\n',x), '%f')
    zb =3D round(x.*10^7)./10^7
    [f,e] =3D log2(x);
    f =3D f*10^7;
    temp =3D pow2(f, e);
    temp =3D round(temp);
    [f,e] =3D log2(temp);
    f =3D f/10^7;
    zc =3D pow2(f, e)

    za =3D=3D y
    zb =3D=3D y
    zc =3D=3D y
But it failed the test, too. For the time being I have to keep using
the costly sprintf/scanf solution, simply because it's the only one
that I know of passing all test cases.

> You should be making comparison like abs(za-zb) < tol where tol is someth=
ing like 10^-15*abs(x)
>
Usually I compare two floating point numbers as you prescribe, but for
this particular application I need exact equality: I input two vectors
containing time tags into intersect() to find the common epochs, and
intersect internally relies on ismember(), which checks for exact
equality between vector elements. I could try to generalize ismember()
-- in the past, I've generalized unique(x, tol) -- but since my input
data is given at finite precision, I find it reasonable to expect
exact equality, after rounding/truncation. (My time tags are given in
the format year month day hour minute seconds -- all fields are
integers except seconds, which has seven decimal places. I end up with
decimals beyond the seventh place because I convert them through
datenum()).

Thanks a lot for your ideas.

Felipe.

Re: robust truncation at n decimal places

Martin Eisenberg — Sun, 07 Dec 2008 12:24:23 -0500

Felipe G. Nievinski wrote:

> I could try to generalize ismember() -- in the past, I've
> generalized unique(x, tol) -- but since my input data is given at
> finite precision, I find it reasonable to expect exact equality,
> after rounding/truncation. (My time tags are given in the format
> year month day hour minute seconds -- all fields are integers
> except seconds, which has seven decimal places. I end up with
> decimals beyond the seventh place because I convert them
> through datenum()).

Narrowly construed, what you want is impossible because for any given
number of decimal digits, many numbers thus representable don't have
terminating (but only periodic) binary expansions. What is your
seconds field's original representation? Does the path it takes
warrant the idea of containing exactly seven decimal places by the
time you get to the current point?

Martin

--
Quidquid latine scriptum est, altum videtur.

Re: robust truncation at n decimal places

Felipe G. Nievinski — Sun, 07 Dec 2008 16:51:54 -0500

On Dec 7, 5:24=A0am, Martin Eisenberg <martin.eisenb...@udo.edu> wrote:
> Felipe G. Nievinski wrote:
> > I could try to generalize ismember() -- in the past, I've
> > generalized unique(x, tol) -- but since my input data is given at
> > finite precision, I find it reasonable to expect exact equality,
> > after rounding/truncation. (My time tags are given in the format
> > year month day hour minute seconds -- all fields are integers
> > except seconds, which has seven decimal places. I end up with
> > decimals beyond the seventh place because I convert them
> > through datenum()).
>
> Narrowly construed, what you want is impossible
Well, what I want (as defined in the test cases above) is possible
with sprintf/sscanf -- I'm just looking for a computationally cheaper
algorithm, avoiding conversion to/from char strings.

> because for any given
> number of decimal digits, many numbers thus representable don't have
> terminating (but only periodic) binary expansions. What is your
> seconds field's original representation?
Seconds are given as a base 10 decimal number with seven significant
figures, e.g., 14.9999999 s.
Perhaps it's important to mention that my events are widely spaced in
time (e.g., usually 1 s apart) compared to the time tags precision
(1e-7 s).

> Does the path it takes
> warrant the idea of containing exactly seven decimal places by the
> time you get to the current point?
>
I don't understand what you mean by "path".

Thanks for your comments.

I'm still on the lookout for a robust rounding without sprintf/sscanf,
though.

Felipe.

Re: robust truncation at n decimal places

Martin Eisenberg — Sun, 07 Dec 2008 23:43:46 -0500

Felipe G. Nievinski wrote:
> On Dec 7, 5:24

Re: robust truncation at n decimal places

Keegan — Fri, 23 Jul 2010 22:43:05 -0400

Felipe, this is from:
http://home.online.no/~pjacklam/matlab/software/util/fullindex.html

function y = fixdec(x, n)
%FIXDEC Round towards zero with a specified number of decimals.
%
% Y = FIXDEC(X, N) rounds the elements of X to N decimals.
%
% For instance, fixdec(10*sqrt(2) + i*pi/10, 4) returns 14.1421 + 0.3141i
%
% See also: FIX, FLOOR, CEIL, ROUND, FIXDIG, ROUNDDEC, ROUNDDIG.

% Author: Peter J. Acklam
% Time-stamp: 2004-09-22 20:08:10 +0200
% E-mail: pjacklam@online.no
% URL: http://home.online.no/~pjacklam

   % Check number of input arguments.
   error(nargchk(2, 2, nargin));

   % Quick exit if either argument is empty.
   if isempty(x) || isempty(n)
      y = [];
      return
   end

   % Get size of input arguments.
   size_x = size(x);
   size_n = size(n);
   scalar_x = all(size_x == 1); % True if x is a scalar.
   scalar_n = all(size_n == 1); % True if n is a scalar.

   % Check size of input arguments.
   if ~scalar_x & ~scalar_n & ~isequal(size_x, size_n)
      error(['When both arguments are non-scalars they must have' ...
             ' the same size']);
   end

   f = 10.^n;
   y = fix(x .* f) ./ f;