Path: news.mathworks.com!newsfeed-00.mathworks.com!newsfeed2.dallas1.level3.net!news.level3.com!postnews.google.com!v5g2000prm.googlegroups.com!not-for-mail
From: "Felipe G. Nievinski" <fgnievinski@gmail.com>
Newsgroups: comp.soft-sys.matlab,sci.math.num-analysis
Subject: Re: robust truncation at n decimal places
Date: Sun, 7 Dec 2008 01:26:53 -0800 (PST)
Organization: http://groups.google.com
Lines: 109
Message-ID: <fac1f966-62ff-427d-b67c-7cce5297f210@v5g2000prm.googlegroups.com>
References: <545015a5-3f51-41e3-a40b-1cf7db531b11@35g2000pry.googlegroups.com> 
	<ghfobu$blk$1@fred.mathworks.com>
NNTP-Posting-Host: 128.138.43.113
Mime-Version: 1.0
Content-Type: text/plain; charset=ISO-8859-1
Content-Transfer-Encoding: quoted-printable
X-Trace: posting.google.com 1228642013 13697 127.0.0.1 (7 Dec 2008 09:26:53 GMT)
X-Complaints-To: groups-abuse@google.com
NNTP-Posting-Date: Sun, 7 Dec 2008 09:26:53 +0000 (UTC)
Complaints-To: groups-abuse@google.com
Injection-Info: v5g2000prm.googlegroups.com; posting-host=128.138.43.113; 
	posting-account=iCBMwwoAAAARtxHOK9kIhRHsTae4FBZX
User-Agent: G2/1.0
X-HTTP-UserAgent: Mozilla/5.0 (Windows; U; Windows NT 6.0; en-US) 
	AppleWebKit/525.19 (KHTML, like Gecko) Chrome/0.4.154.29 Safari/525.19,gzip(gfe),gzip(gfe)
Xref: news.mathworks.com comp.soft-sys.matlab:505443 sci.math.num-analysis:105694

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.