Rounding double

[Marco Manfredini]

double roundto(long double value, unsigned digits)
{
long double fv=fabsl(value);
if (fv > powl(10.0L,DBL_DIG) || digits > DBL_DIG)
return value; // Out of range
long long p = powl(10.0L, digits);
return roundl(value * p ) / (double)p;
}

Supposing long double is extended precision, why should this not
work?

Suppose long double is 128 bit and long long 64 bit. Why should this work?

 

[Richard Heathfield]

jacob navia said:

Richard Heathfield wrote:


This is perfectly OK since double has only 16 decimal digits,
and the first ones in your result appear at 1e-21 (if I counted
correctly). Those values can't be accurately represented with
double precision.

Nor can 0.3.

You can't go beyond 1e-16, since DBL_DIG is 15.

Writing it in in a more sensible way:

double roundto(long double value, unsigned digits)
{
long double fv=fabsl(value);
if (fv > powl(10.0L,DBL_DIG) || digits > DBL_DIG)
return value; // Out of range
long long p = powl(10.0L, digits);
return roundl(value * p ) / (double)p;
}

Could someone else please share with us the result they get when they drive
this function with inputs of 0.33 and 1 ? I'm getting some very, very
strange results here.

 

[Mark McIntyre]

OK That's much better tone Heathfield. Let's calm down.

Can I also suggest we start being polite too? In most of the Western
world, referring to someone by their surname alone is considered either
rude or authoritarian. Nobody has called me just plain "Mcintyre" since
I left school.

I use more precision than double,

I think this is the tricky part - the code will of necessity be
nonportable, since C doesn't constrain the maximum precision of doubles.

 

[Mark McIntyre]

Who cares?
That works in lcc-win.
Lcc-win gives you the best of C and C++

This is a very bad answer, and precisely what you should /not/ have said
if you wanted to be taken seriously.

 

[jacob navia]

jacob navia said:


Nor can 0.3.


Could someone else please share with us the result they get when they drive
this function with inputs of 0.33 and 1 ? I'm getting some very, very
strange results here.

d:\lcc\mc68\test>type tdouble3.c
#include <stdio.h>
#include <float.h>
#include <math.h>

double roundto(long double value, unsigned digits)

{
long double fv = fabs(value);
if (fv > powl(10.0L,DBL_DIG) || digits > DBL_DIG)
return value;
long long p = powl(10.0L, digits);
return roundl(p*value)/p;
}

int main(int argc,char *argv[])
{
long double v = atof(argv[1]);
short int ndp = 0;
while(ndp < DBL_DIG)
{
long double newv = v;
newv = roundto(newv, ndp);
printf("%.19Lf: %hd decimals %.16Lf\n",
v,
ndp++,
newv);
}

return 0;
}



d:\lcc\mc68\test>tdouble3 0.33
0.3300000000000000160: 0 decimals 0.0000000000000000
0.3300000000000000160: 1 decimals 0.3000000000000000
0.3300000000000000160: 2 decimals 0.3300000000000000
0.3300000000000000160: 3 decimals 0.3300000000000000
0.3300000000000000160: 4 decimals 0.3300000000000000
0.3300000000000000160: 5 decimals 0.3300000000000000
0.3300000000000000160: 6 decimals 0.3300000000000000
0.3300000000000000160: 7 decimals 0.3300000000000000
0.3300000000000000160: 8 decimals 0.3300000000000000
0.3300000000000000160: 9 decimals 0.3300000000000000
0.3300000000000000160: 10 decimals 0.3300000000000000
0.3300000000000000160: 11 decimals 0.3300000000000000
0.3300000000000000160: 12 decimals 0.3300000000000000
0.3300000000000000160: 13 decimals 0.3300000000000000
0.3300000000000000160: 14 decimals 0.3300000000000000

d:\lcc\mc68\test>tdouble3 1
1.0000000000000000000: 0 decimals 1.0000000000000000
1.0000000000000000000: 1 decimals 1.0000000000000000
1.0000000000000000000: 2 decimals 1.0000000000000000
1.0000000000000000000: 3 decimals 1.0000000000000000
1.0000000000000000000: 4 decimals 1.0000000000000000
1.0000000000000000000: 5 decimals 1.0000000000000000
1.0000000000000000000: 6 decimals 1.0000000000000000
1.0000000000000000000: 7 decimals 1.0000000000000000
1.0000000000000000000: 8 decimals 1.0000000000000000
1.0000000000000000000: 9 decimals 1.0000000000000000
1.0000000000000000000: 10 decimals 1.0000000000000000
1.0000000000000000000: 11 decimals 1.0000000000000000
1.0000000000000000000: 12 decimals 1.0000000000000000
1.0000000000000000000: 13 decimals 1.0000000000000000
1.0000000000000000000: 14 decimals 1.0000000000000000

d:\lcc\mc68\test>

 

[Richard Heathfield]

jacob navia said:

Obviously numbers are written in
binary, but we can do the dismissed quantity very
small. That means rounding to the nth digit.

Yes, but your code *doesn't* round to the nth digit. If it did, only n
non-zero digits would remain. Yes, you can get very, very close to the
right value, but you can't always hit it.

If we use higher precision in the calculations, and only round at the
end, we get the best possible result for that rounding.

Fine, but we can't store the *correct* result, the result that the OP asked
for, in a double - not in general, anyway, although obviously there are
various values that *do* work.

Unrepresentable numbers in floating point will not go away

Precisely. And therefore the OP's problem is insoluble in general terms, so
he must look for a workaround.

but we can get the BEST approximation.

Sure, but that is not what the OP asked for.

 

[Mark McIntyre]

This is perfectly OK since double has only 16 decimal digits,

This is incorrect. The ISO standard requires a double to have at least
10 digits, but specifies no upper bound.

5.2.4.2.2 Characteristics of floating types
8. The values given in the following list shall be replaced by constant
expressions with implementation-defined values that are greater or equal
in magnitude (absolute value) to those shown, with the same sign:
FLT_DIG 6
DBL_DIG 10
LDBL_DIG 10

There are a couple of worked examples showing 15 digits, and one
footnote that references the ISO/IEC floating point standard, but
nothing else normative.

 

[Richard Heathfield]

jacob navia said:

d:\lcc\mc68\test>type tdouble3.c

No, I said "someone else". I wouldn't trust you to make a cup of coffee,
let alone a computer program.

 

[Keith Thompson]

jacob navia wrote:
[...]

Excuse me nbut you REALLY believe that it is impossible to write in
standard C

double roundto(double n,int digits);

????

Of course it's possible to write a function with that declaration.

The OP asked "how to round a double value with a specific number
of digits after the decimal points". What the OP is asking for is in
fact not possible in most C implementations.

For example, presumably roundto(3.14159, 2) would yield 3.14. The
problem is that 3.14 is not exactly representable in binary floating-point.

(C does allow for decimal floating-point; see C99 5.2.4.2.2, which
requires the radix to be an integer > 1, not necessarily a power of 2.
But few if any real-world C implementations do this. Of course, an
implementation with decimal floating-point would not define
__STDC_IEC_559__, and any code that depends on decimal floating-point is
extremely non-portable.)

You can, of course, write a function that gives you a close
approximation of the mathematically correct rounded result; for example,
where roundto(3.14159, 2) yields the nearest representable approximation
to 3.14. And this might even be good enough for the OP's purposes. But
it doesn't *quite* satisfy what the OP asked for.

Consider this: how many "digits after the decimal point" does the
nearest double approximation of 3.14 have? On one system, the result I
get is 3.140000000000000124344978758017532527446746826171875, so the
answer is 51 digits rather than the desired 2 digits.

But my suspicion, based on the OP's original question and his followups
in this thread, is that he doesn't understand this, and is looking for a
non-existent *exact* solution. I might be misinterpreting what he's
said, but he hasn't stated his problem precisely enough for anyone to be
sure of what solution would satisfy it.

To the original poster: Can you clarify what you're looking for? Why do
you want to round double values to a specified number of decimal digits?
The most common reason to do that is for output; the best way to do
that is to do the rounding when you do the output, not by rounding the
numbers themselves.

 

[Dik T. Winter]

> And it IS possible to round correctly a double precision
> number. One of the possible solutions is:
>
> #include <stdio.h>
> #include <math.h>
>
> void RoundMyDouble (long double *value, short numberOfPrecisions)
>
> {
> long double fv=fabsl(*value);
> if (fv > 1e17)
> return ;
> long long p = powl(10.0L, numberOfPrecisions);
> *value = (long long)(*value * p + 0.5L) / (double)p;
> }

This does *not* round correctly. If you round 6.123465e+23 to four
places you should get: 6.123000e+23.

 

[jacob navia]

jacob navia said:



Yes, but your code *doesn't* round to the nth digit. If it did, only n
non-zero digits would remain. Yes, you can get very, very close to the
right value, but you can't always hit it.


Fine, but we can't store the *correct* result, the result that the OP asked
for, in a double - not in general, anyway, although obviously there are
various values that *do* work.

You are saying in essence that 1/10 is not representable in
binary. I know that. But we are speaking about
FLOATING POINT here, that as everyone knows is an approximation
to the reals, and never anything else!

Of course there are unrepresentable numbers. Does that mean
that we can't use floating point or what?

Obviously we can't get any 100% correct representation of
all numbers in floating point, like for instance in base
ten with 1/3!

So what?

If I go to primary school and say to some kid:

Look round me 0.345345345 to 7 places he will tell
me 0.3453453, even if it is obvious that there is a
repeating fraction.

You CUT at some point, as accurately as possible and that was it!

Precisely. And therefore the OP's problem is insoluble in general terms, so
he must look for a workaround.

NO!

The existence of unrepresentable numbers in floating point
does not mean we can't use floating point or round numbers
as we wish. It means that all our results will be ALWAYS
incorrect from the mathematically true result by some epsilon
that should be as small as possible.

Sure, but that is not what the OP asked for.

It is what he asked for. In the original message
there is no mention of 100% accuracy or of any accuracy requirements
in the first place.

I insist that
1) Given the nature of floating point numbers any solution will
be an approximation in most cases.
2) We can get the best approximation for double precision by
using more precision internally and then rounding at the
last moment.
3) The function proposed is maybe not the best, but it works
for all values in range.

 

[jacob navia]

This is incorrect. The ISO standard requires a double to have at least
10 digits, but specifies no upper bound.

5.2.4.2.2 Characteristics of floating types
8. The values given in the following list shall be replaced by constant
expressions with implementation-defined values that are greater or equal
in magnitude (absolute value) to those shown, with the same sign:
FLT_DIG 6
DBL_DIG 10
LDBL_DIG 10

There are a couple of worked examples showing 15 digits, and one
footnote that references the ISO/IEC floating point standard, but
nothing else normative.

That "footnote" is Annex F, and it IS normative. I cite (again)

F.2 Types
1 The C floating types match the IEC 60559 formats as follows:
— The float type matches the IEC 60559 single format.
— The double type matches the IEC 60559 double format

 

[Marco Manfredini]

jacob navia said:


No, I said "someone else". I wouldn't trust you to make a cup of coffee,
let alone a computer program.

long double roundto(long double value, unsigned digits)
{
long double fv=fabsl(value);
if (fv > powl(10.0L,DBL_DIG) || digits > DBL_DIG)
return value; // Out of range
long long p = powl(10.0L, digits);
return roundl(value * p ) / (double)p;
}

int main()
{
long double r=roundto(0.33,1);
printf("%Lf\n",r);
r=1/(r-0.3);

printf("%Lf\n",r);
}

$ ./a.out
0.300000
89984117432729520.078125

Works great.

 

[jacob navia]

jacob navia wrote, On 22/11/07 23:01:

It does not meet the OPs requirements as specified.


<snip>

Which is precisely the problem. If it is wanted for display then there
are better ways, if it is wanted for further calculations then it is
very important that the OP understand why it is not possible in general.

OK from a teaching viewpoint maybe. But there IS a general solution
within a double precision framework. Of course due to the
nature of FP it is an approximation but so what? That is
floating point and there is no way around it.

Instead of saying "Impossible" we should give the best
approximation that is all.

 

[jacob navia]

long double roundto(long double value, unsigned digits)
{

This is not my function excuse me.

My function returned a double, NOT a long double.

AND THAT IS IMPORTANT!

Works great.

Yes. Works great.

 

[Mark McIntyre]

d:\lcc\mc68\test>type tdouble3.c

I copy-pasted this exact code, and compiled under gcc 4.1.2:

thelinux clc_tests $ ./a.out 0.33
1374389535.0000000000000000000: 0 decimals 1374389535.0000000000000000
1374389535.0000000000000000000: 1 decimals 1374389535.0000000000000000

I'm giving this one a 'strange' score of 'high'

 

[Dik T. Winter]

> Now, are we really such a C "heads" that we can't answer the OP
> question IN HIS terms?

You discarded the solution of James Kuyper out of hand. Did you even
think about that solution?

> double roundto(double d,int digits);

double roundto(double d, int digits) {
char number[digits+10];
sprintf(number, "%.*e", (digits > 0 ? digits - 1 : 0), d);
sscanf(number, "%f", &d);
return d;
}

Or somesuch.

 

[James Kuyper]

Have you given any more thought to my hint? If not, why did you quote it?

Note: Richard Heathfield has been saying this is impossible, I've been
saying it is possible. That's because we're referring to two different
things. It is impossible in binary floating point formats to represent
exactly an arbitrary floating point value rounded to a specified number
of decimal digits. That's simply because 10 is not a power of 2. That's
what Richard was talking about.

However, it is possible to calculate a result that is as close to the
ideally rounded value as is permitted by the precision of the floating
point type, which is what I'm talking about.

double roundto(long double value, unsigned digits)
{
long double fv=fabsl(value);
if (fv > powl(10.0L,DBL_DIG) || digits > DBL_DIG)
return value; // Out of range
long long p = powl(10.0L, digits);
return roundl(value * p ) / (double)p;
}

Supposing long double is extended precision, why should this not
work?

Well, it unnecessarily fails to work when digits is negative. I suppose
you could impose the arbitrary requirement digits >= 0, but the function
has an obvious meaning for digits<0, and that meaning is just as useful
as the meaning for digits >= 0 (which is not meant to imply that it's
particularly useful in either case).

Specifying that long double has extended precision is not sufficient,
even when digits>=0. A conforming implementation of C can use extended
precision for both double and long double, so long as __STDC_IEC559__ is
not #defined. In that case there's no guarantee that 'p' is big enough
to store the largest integer that can be stored in a double. As a
result, the rounded result will have less accuracy than it should.

 

[Dik T. Winter]

> double roundto(long double value, unsigned digits)
> {
> long double fv=fabsl(value);
> if (fv > powl(10.0L,DBL_DIG) || digits > DBL_DIG)
> return value; // Out of range
> long long p = powl(10.0L, digits);
> return roundl(value * p ) / (double)p;
> }
>
> Supposing long double is extended precision, why should this not
> work?

Because the return of roundto(6.12345e+50, 4) should be 6.12300e+50.

 

[James Kuyper]

jacob navia wrote:
....

This is perfectly OK since double has only 16 decimal digits,
and the first ones in your result appear at 1e-21 (if I counted

Citation please? Where does the C standard impose such requirements for
all implementations?