problem with double's precision

Hello,

I've got a bit of experience in C++, but I'm writing my first app that
is dependent on relatively precise math functions. The app requires
that I get a time stamp based on s sample number, from a time series.

This seems liek an easy thing to do:

long lSample = 500; (for example)
double dSampleRate = 1000.0;
double dTimestamp = (double)lSample / dSampleRate;

The problem I am having is that the division by 1000 is not precise,
presumably because 1/(10^n) cannot be represented as a standard
floating point number. So, with each subsequent sample, the calculated
time stamp is further from the actual time of the sample.

The same problem occurs if I simply add the sample period to the
previous timestamp, and start at dTimestamp = 0.0. It also occurs when
trying to get the sample period, e.g.

double dSamplePeriod = 1.0 / dSampleRate

then dSamplePeriod = 0.00100000004749745

Currently, I am using a class off the web, called CLargeDouble, to get
around this, but this class uses string conversion, and as such, slows
the program immensely.

This app is being compiled using MS Visual C++ 7, and run on XP
machines. However, when I ran it on our terminal server (Win 2000
server), the 1/1000 operation did _not_ include the garbage starting at
the 6th or 7th siginificant digit.

My question is, is there a standard, fast workaround to this problem?
Or does it seem that this involves some problem in my code, which only
occurs when running on an XP platform? My hope is to avoid truncating
at some arbitrary number of significant digits, since I will be using
different sample rates (e.g. 1024hz, the inverse of which can be
respesented in floating point).

Thanks in advance for any help.
-D.T.

* towers@email.arizona.edu:[color=blue]
>
> I've got a bit of experience in C++, but I'm writing my first app that
> is dependent on relatively precise math functions. The app requires
> that I get a time stamp based on s sample number, from a time series.
>
> This seems liek an easy thing to do:
>
> long lSample = 500; (for example)
> double dSampleRate = 1000.0;
> double dTimestamp = (double)lSample / dSampleRate;
>
> The problem I am having is that the division by 1000 is not precise,
> presumably because 1/(10^n) cannot be represented as a standard
> floating point number. So, with each subsequent sample, the calculated
> time stamp is further from the actual time of the sample.[/color]

Nope. The calculated time stamp is approximately the hypothetical value
associated with the sample number, and that approximation is not significant,
and in particular it doesn't get worse for larger values. Rather, the
assumption that sample number can be used to compute a time stamp is probably
invalid, i.e., your samples do not actually come at a fixed rate.

[color=blue]
> The same problem occurs if I simply add the sample period to the
> previous timestamp, and start at dTimestamp = 0.0.[/color]

That is a different (and additional) problem, of accumulated error.

--
A: Because it messes up the order in which people normally read text.
Q: Why is it such a bad thing?
A: Top-posting.
Q: What is the most annoying thing on usenet and in e-mail?

towers@email.arizona.edu wrote:[color=blue]
> I've got a bit of experience in C++, but I'm writing my first app that
> is dependent on relatively precise math functions. The app requires
> that I get a time stamp based on s sample number, from a time series.
>
> This seems liek an easy thing to do:
>
> long lSample = 500; (for example)
> double dSampleRate = 1000.0;
> double dTimestamp = (double)lSample / dSampleRate;[/color]

So, supposedly your 'dTimestamp' should just be 0.5, right? Is it? The
value 0.5 is exactly representable in binary, BTW.
[color=blue]
> The problem I am having is that the division by 1000 is not precise,
> presumably because 1/(10^n) cannot be represented as a standard
> floating point number. So, with each subsequent sample, the calculated
> time stamp is further from the actual time of the sample.[/color]

That's why you should probably simply define 'dTimestamp' as 0.5. Also,
instead of accumulating your time value (and the error it carries),
calculate it: (step_index * 500.) / 1000
[color=blue]
> The same problem occurs if I simply add the sample period to the
> previous timestamp, and start at dTimestamp = 0.0. It also occurs when
> trying to get the sample period, e.g.
>
> double dSamplePeriod = 1.0 / dSampleRate
>
> then dSamplePeriod = 0.00100000004749745[/color]

That looks awfully imprecise. Are you sure you're using 'double'?
[color=blue]
> Currently, I am using a class off the web, called CLargeDouble, to get
> around this, but this class uses string conversion, and as such, slows
> the program immensely.[/color]

As soon as I see a class whose name begins with 'C', I shudder. It's
apparently written by a Microsoft-centric programmer, and that's
frightening. You should be able to find a better representation of
a larger precision FP number on the web than one that uses string as
its internal storage...
[color=blue]
> This app is being compiled using MS Visual C++ 7, and run on XP
> machines. However, when I ran it on our terminal server (Win 2000
> server), the 1/1000 operation did _not_ include the garbage starting at
> the 6th or 7th siginificant digit.[/color]

<shrug> Unless you post your code (without the CLargeDouble class), it
is impossible to comment on all of that.
[color=blue]
> My question is, is there a standard, fast workaround to this problem?[/color]

Yes, more than likely. Without seeing your code, however, we cannot
advise on using one technique or another.
[color=blue]
> Or does it seem that this involves some problem in my code, which only
> occurs when running on an XP platform?[/color]

Nothing can be said about that. 'comp.lang.c++' is platform-indifferent.
[color=blue]
> My hope is to avoid truncating
> at some arbitrary number of significant digits, since I will be using
> different sample rates (e.g. 1024hz, the inverse of which can be
> respesented in floating point).[/color]

Well, you're correct about this, whenever division is used, you should be
able to figure out a possible rounding error by replacing it with the
inverse multiplication. Dividing by 1000 carries potential rounding
errors, while dividing by 1024 shouldn't, at least on a system where
floating point numbers are IEEE-standard.

V

Thanks for the input. I've fiddles with the code a bit more, and now
have something reasonable - no iterative addition of the sample period,
just direct assignment (timestamp = period * sample number). For some
reason, this didn't work in an earlier version of the program - gave
the same results as the iterative solution (i.e. multiplying the period
by N gave the same value as adding the period N time). Now, seems to
work fine, putting aside the variability how precise each of products
turn out to be , e.g.

dSamplePeriod * (double)lIndex = dTimestamp;

(sample output)
0.001000000000 * 4 = 0.004000000190
0.001000000000 * 5 = 0.004999999888

If this does seem imprecise for a double, I'm guessing that is a
platform/compiler issue;

Anyhow, though the more subtle points of floating point math still
escape me, your suggestions got my code running a couple orders of
magnitude faster - so, many, many thanks!

-D.T.

Victor Bazarov wrote:[color=blue]
> towers@email.arizona.edu wrote:[color=green]
> > I've got a bit of experience in C++, but I'm writing my first app that
> > is dependent on relatively precise math functions. The app requires
> > that I get a time stamp based on s sample number, from a time series.
> >
> > This seems liek an easy thing to do:
> >
> > long lSample = 500; (for example)
> > double dSampleRate = 1000.0;
> > double dTimestamp = (double)lSample / dSampleRate;[/color]
>
> So, supposedly your 'dTimestamp' should just be 0.5, right? Is it? The
> value 0.5 is exactly representable in binary, BTW.
>[color=green]
> > The problem I am having is that the division by 1000 is not precise,
> > presumably because 1/(10^n) cannot be represented as a standard
> > floating point number. So, with each subsequent sample, the calculated
> > time stamp is further from the actual time of the sample.[/color]
>
> That's why you should probably simply define 'dTimestamp' as 0.5. Also,
> instead of accumulating your time value (and the error it carries),
> calculate it: (step_index * 500.) / 1000
>[color=green]
> > The same problem occurs if I simply add the sample period to the
> > previous timestamp, and start at dTimestamp = 0.0. It also occurs when
> > trying to get the sample period, e.g.
> >
> > double dSamplePeriod = 1.0 / dSampleRate
> >
> > then dSamplePeriod = 0.00100000004749745[/color]
>
> That looks awfully imprecise. Are you sure you're using 'double'?
>[color=green]
> > Currently, I am using a class off the web, called CLargeDouble, to get
> > around this, but this class uses string conversion, and as such, slows
> > the program immensely.[/color]
>
> As soon as I see a class whose name begins with 'C', I shudder. It's
> apparently written by a Microsoft-centric programmer, and that's
> frightening. You should be able to find a better representation of
> a larger precision FP number on the web than one that uses string as
> its internal storage...
>[color=green]
> > This app is being compiled using MS Visual C++ 7, and run on XP
> > machines. However, when I ran it on our terminal server (Win 2000
> > server), the 1/1000 operation did _not_ include the garbage starting at
> > the 6th or 7th siginificant digit.[/color]
>
> <shrug> Unless you post your code (without the CLargeDouble class), it
> is impossible to comment on all of that.
>[color=green]
> > My question is, is there a standard, fast workaround to this problem?[/color]
>
> Yes, more than likely. Without seeing your code, however, we cannot
> advise on using one technique or another.
>[color=green]
> > Or does it seem that this involves some problem in my code, which only
> > occurs when running on an XP platform?[/color]
>
> Nothing can be said about that. 'comp.lang.c++' is platform-indifferent.
>[color=green]
> > My hope is to avoid truncating
> > at some arbitrary number of significant digits, since I will be using
> > different sample rates (e.g. 1024hz, the inverse of which can be
> > respesented in floating point).[/color]
>
> Well, you're correct about this, whenever division is used, you should be
> able to figure out a possible rounding error by replacing it with the
> inverse multiplication. Dividing by 1000 carries potential rounding
> errors, while dividing by 1024 shouldn't, at least on a system where
> floating point numbers are IEEE-standard.
>
> V[/color]

towers@email.arizona.edu wrote:[color=blue]
> Thanks for the input. I've fiddles with the code a bit more, and now
> have something reasonable - no iterative addition of the sample period,
> just direct assignment (timestamp = period * sample number). For some
> reason, this didn't work in an earlier version of the program - gave
> the same results as the iterative solution (i.e. multiplying the period
> by N gave the same value as adding the period N time). Now, seems to
> work fine, putting aside the variability how precise each of products
> turn out to be , e.g.
>
> dSamplePeriod * (double)lIndex = dTimestamp;
>
> (sample output)
> 0.001000000000 * 4 = 0.004000000190
> 0.001000000000 * 5 = 0.004999999888
>
> If this does seem imprecise for a double, I'm guessing that is a
> platform/compiler issue;
>
> Anyhow, though the more subtle points of floating point math still
> escape me, your suggestions got my code running a couple orders of
> magnitude faster - so, many, many thanks![/color]

Why not use fixed point instead of floating point values? An easy way
to represent fixed point is simply to use integer values scaled
proportionately. If .001 is the precision needed than the integer value
1 represents 1/1000. Values will remain accurate within that margin:

const long kScaling = 1000;

long lSample = 500 * kScaling;
long dSampleRate = 1;
long dTimestamp = lSample / dSampleRate;

long dSamplePeriod = 1000 / dSampleRate

then dSamplePeriod will equal 1 (i.e. 1/1000th)

1 * 4 = 4 (i.e .004)
1 * 5 = 5 (i.e .005)

The values need be divided by 1000 only when they are to be displayed.
Internally, the program simply performs integer arithmetic.

Greg

doubles are not that accurate (normally a lot better than float)
If you need accuracy (and not speed) you might want to look at MAPM

http://www.tc.umn.edu/~ringx004/mapm-main.html