J
jacob navia
James said:
It is completely meaningless.
It is completely meaningless.
R
rembremading
Thanks for your replies. It helps a lot.
In order to conclude:
double any_function()
{
}
(1) double temp1, temp2;
(2) temp1 = any_function();
(3) temp2 = temp1 - any_function();
can lead to a result temp2 != 0.0 because temp1 is stored in the memory and
the result of any_function() in (2) is therefore rounded to the closest 64
bit number, whereas the result of any_function() in (3) is stored the
register only, such that the difference in (3) can have a finite outcome.
Because of the precision discussion I would like to remind that in terms of
C numbers this result ~ 10e-17 *is* equal to 0.0 (in the context of the
subtraction of two order 1 numbers).
In this sense it is not an issue of missing precision of the '-' operation.
I think that actually the comparison
if(floating_point_number_1==floating_point_number_2)
is nonsense, if the floating point variables are the result of some
calculation. And I am doing this comparison only because of performance
issues.
On the other hand in
double temp1, temp2, temp3;
temp1 = any_function();
temp2 = any_function();
temp3 = temp1 - temp2;
temp3 == 0 because temp1 and temp2 are stored in the memory and therefore
the result of any_function is rounded in the same way.
Is this correct?
Unfortunately, when I define temp1 and temp2 as register variables
register double temp1, temp2;
Then, surprisingly, I still get a finite temp1. (of course register
variables are not guaranteed to end up in the register)
The solution for my concrete problem will probably be to choose the
numerical constants such, that the numbers occurring are exactly
representable numbers, which should be possible.
-Andreas
In order to conclude:
double any_function()
{
}
(1) double temp1, temp2;
(2) temp1 = any_function();
(3) temp2 = temp1 - any_function();
can lead to a result temp2 != 0.0 because temp1 is stored in the memory and
the result of any_function() in (2) is therefore rounded to the closest 64
bit number, whereas the result of any_function() in (3) is stored the
register only, such that the difference in (3) can have a finite outcome.
Because of the precision discussion I would like to remind that in terms of
C numbers this result ~ 10e-17 *is* equal to 0.0 (in the context of the
subtraction of two order 1 numbers).
In this sense it is not an issue of missing precision of the '-' operation.
I think that actually the comparison
if(floating_point_number_1==floating_point_number_2)
is nonsense, if the floating point variables are the result of some
calculation. And I am doing this comparison only because of performance
issues.
On the other hand in
double temp1, temp2, temp3;
temp1 = any_function();
temp2 = any_function();
temp3 = temp1 - temp2;
temp3 == 0 because temp1 and temp2 are stored in the memory and therefore
the result of any_function is rounded in the same way.
Is this correct?
Unfortunately, when I define temp1 and temp2 as register variables
register double temp1, temp2;
Then, surprisingly, I still get a finite temp1. (of course register
variables are not guaranteed to end up in the register)
The solution for my concrete problem will probably be to choose the
numerical constants such, that the numbers occurring are exactly
representable numbers, which should be possible.
-Andreas
J
jacob navia
rembremading wrote:
[snip]
Can you tell me if you use the following program the result is better?
This program sets the precision to 64 bits instead of
80 bits. (If I have no bugs... )
void SetPrecision64(void)
{
void (*fn)(void);
unsigned char code[] = {
0x50, // push %eax
0xd9,0x3c,0x24, // fnstcw (%esp,1)
0x8b,0x04,0x24, // mov (%esp,1),%eax
0x0f,0xba,0x34,0x24,0x08, // btrl $0x8,(%esp,1)
0x66,0x81,0x0c,0x24,0x00,0x02,// orw $0x200,(%esp,1)
0xd9,0x2c,0x24, // fldcw (%esp,1)
0x59, // pop %ecx
0xc3 // ret
};
fn = (void (*)(void))code;
fn();
}
[snip]
Can you tell me if you use the following program the result is better?
This program sets the precision to 64 bits instead of
80 bits. (If I have no bugs... )
void SetPrecision64(void)
{
void (*fn)(void);
unsigned char code[] = {
0x50, // push %eax
0xd9,0x3c,0x24, // fnstcw (%esp,1)
0x8b,0x04,0x24, // mov (%esp,1),%eax
0x0f,0xba,0x34,0x24,0x08, // btrl $0x8,(%esp,1)
0x66,0x81,0x0c,0x24,0x00,0x02,// orw $0x200,(%esp,1)
0xd9,0x2c,0x24, // fldcw (%esp,1)
0x59, // pop %ecx
0xc3 // ret
};
fn = (void (*)(void))code;
fn();
}
J
James Kuyper
jacob said:
3.9975222353 is a very precise but inaccurate answer to the question
"what is 2.0 + 2.0?".
The numerical limits section describes and puts minimum requirements on
the precision; while section 5.2.4.2.2 makes explicit the fact that the
standard imposes no requirements on the accuracy.
B
Ben Bacarisse
rembremading said:
<obsession type="favourite">
You mean non-zero rather than finite both here and later on.
</obsession>
D
Dik T. Winter
Have a look at the Ada standard, where that is precisely the case. But
the C standard only goes halfway there. They properly define a model
for floating-point representation, but give no guarantee about the
precision within that model, and that is something the Ada standard
does do. What does happen in Ada systems if the maximal model for
the floating-point hardware does not meet the requirements of precision,
the model is reduced.
T
Thad Smith
behaviour.rembremading said:
This is cross-posted to comp.std.c. See the final paragraph for a
standards question. This is the second posting attempt.
The original code temp3 may be non-zero because of extra precision:
5.2.4.2.2:
The value of operations with floating operands and values subject to the
usual arithmetic conversions and of floating constant are evaluated to a
format whose range and precision may be greater than required by the
type. The use of evaluation formats is characterized by the
implementation-defined value of FLT_EVAL_METHOD.
The standard requires that values be rounded to their specified
precision by storing into a variable of the type or cast to the type
(See 5.1.2.3 footnote 4). Neither of these requirements are in effect
for the use of GRID_POINT(2) within the statement computing temp3.
The compiler is probably passing the evaluated value back in a floating
point register with more precision than a double variable.
Regarding the effect of replacing (double)i by the constant 2, you
permit the compiler to evaluate the expression in GRID_POINT() at
compile time. The compiler probably is then propagating that constant
to the assignments in main. What is apparently different in the
compile-time evaluation is that the value computed in GRID_POINT() gets
rounded to double precision before it is used in main.
Is the compiler rounding of the constant value, compared with the
register passing with extra precision, conforming? 5.1.2.3 footnote 4
says "an implicit spilling of a register is not permitted to alter the
value." Is the compile time evaluation, vs. runtime evaluation, an
implicit spilling of a register? If not, what is an "implicit spilling
of a register"?
J
James Kuyper
jacob said:
While the original poster's problem was with the accuracy of the
calculations, your very first response was about the precision, not the
accuracy. An 80-bit number can hold an inaccurate number with much
greater precision than an 64-bit number; in itself that doesn't make the
number any more accurate.
The standard does say things about the precision. The values for
LDBL_EPSILON, DBL_EPSILON, and FLT_EPSILON depend very directly on the
precision of floating point operations, and the standard sets a maximum
value for the magnitude of that constant, which implies a minimum number
of bits of precision.
The standard also says things about when it's permissible to use a
higher precision type for intermediate steps in calculations, and it
seems to be precisely that issue which is the actual reason for the
discrepancies that the original poster reported. See 5.2.4.2.2p8. If
FLT_EVAL_METHOD were 0, then he shouldn't have seen those discrepancies.
For any other value of FLT_EVAL_METHOD, such discrepancies are always a
possibility.
M
Mark Bluemel
jacob said:
http://www.network-theory.co.uk/docs/gccintro/gccintro_70.html
(previously cited) gives a briefer version, which may be specific to
gcc, but that's what the OP has the issue with ....
J
jacob navia
Mark said:
Yes, the one there is better adapted for gcc.
That document CONFIRMS what I am saying since the start.
All this problem can probably be solved by
changing the precision of the machine!
R
rembremading
Yes, your assembler code works fine, although I don't understand what it
does
For my purposes, however, it will be better to choose the constants such
that they are exactly representable, as outlined below.
Every additional function call, at this place in the code, is unwanted.
(although an inlined version of this function would be quite fast, I guess)
I would not want to activate the 64 bit computations for the entire program.
Actually for the most part of it calculation with higher precision are very
wellcome, as long as it does not affect the speed)
Nevertheless it is an interesting piece of code, which might be usefull in
the future.
does
For my purposes, however, it will be better to choose the constants such
that they are exactly representable, as outlined below.
Every additional function call, at this place in the code, is unwanted.
(although an inlined version of this function would be quite fast, I guess)
I would not want to activate the 64 bit computations for the entire program.
Actually for the most part of it calculation with higher precision are very
wellcome, as long as it does not affect the speed)
Nevertheless it is an interesting piece of code, which might be usefull in
the future.
C
CBFalconer
James said:
Please describe how you implement infinite precision in a finite
number of bits.
J
jacob navia
rembremading said:Yes, your assembler code works fine, although I don't understand what it
does
For my purposes, however, it will be better to choose the constants such
that they are exactly representable, as outlined below.
Every additional function call, at this place in the code, is unwanted.
(although an inlined version of this function would be quite fast, I guess)
I would not want to activate the 64 bit computations for the entire program.
Actually for the most part of it calculation with higher precision are very
wellcome, as long as it does not affect the speed)
Nevertheless it is an interesting piece of code, which might be usefull in
the future.
You need it call it ONCE when your program STARTS.
Then, all calculations are done in 64 bits only
and not in 80.
That code sets the precision flag of the floating point unit
to 64 bits, i.e. double precision.
M
Mark Bluemel
rembremading said:Yes, your assembler code works fine, although I don't understand what it
does
Then read http://www.network-theory.co.uk/docs/gccintro/gccintro_70.html
It answers your questions fairly comprehensively and gives a simple
way of disabling extended precision for the duration of a program's
execution.
The mechanism outlined is the simplest way of "fixing" your problem.
R
rembremading
Yes, I understand that I would have to call it only once, but in the main
threefold for loop of my program I have to do something like
res = ( (1.0+f1)*(1.0+f2)*f3*f4 - f1*f2*(1.0+f3)*(1.0+f4) );
Where f1..f4 are function values (not of the GRID_POINT function), which can
take values from 1.0e-100 to 1.0e2.
In this case, presumably, the intermediate values are stored in the register
and it is obvious that res is highly sensitive to the numerical precision
of f1..f4 and the intermediate results of the computation.
Therefore, if I can get the computation with 80 bit precision instead of 64
bit (at same costs), I definitely prefere the former one.
Although the final result is rounded to 64 bit precision this extra digits
can change a lot.
Because my system has dimension ~1e3 and the loop is threefold I cannot use
arbitrary precision math libraries. I cannot even use long double type
variables on my machine.
Now, the GRID_POINT function has to be evaluated when f1..f4 are computed.
Therefore I would have to switch 64 bit precision on (and off again) in
every single step. That is what I meant.
I understood that, but I do not understand how the function does it
(possibly lack of assembler knowledge)
There might be also a lack of English knowledge, which could be the reason
for missunderstandings. Sorry for that.
-Andreas
threefold for loop of my program I have to do something like
res = ( (1.0+f1)*(1.0+f2)*f3*f4 - f1*f2*(1.0+f3)*(1.0+f4) );
Where f1..f4 are function values (not of the GRID_POINT function), which can
take values from 1.0e-100 to 1.0e2.
In this case, presumably, the intermediate values are stored in the register
and it is obvious that res is highly sensitive to the numerical precision
of f1..f4 and the intermediate results of the computation.
Therefore, if I can get the computation with 80 bit precision instead of 64
bit (at same costs), I definitely prefere the former one.
Although the final result is rounded to 64 bit precision this extra digits
can change a lot.
Because my system has dimension ~1e3 and the loop is threefold I cannot use
arbitrary precision math libraries. I cannot even use long double type
variables on my machine.
Now, the GRID_POINT function has to be evaluated when f1..f4 are computed.
Therefore I would have to switch 64 bit precision on (and off again) in
every single step. That is what I meant.
I understood that, but I do not understand how the function does it
(possibly lack of assembler knowledge)
There might be also a lack of English knowledge, which could be the reason
for missunderstandings. Sorry for that.
-Andreas
J
jacob navia
Your problem is a very complicated one Andreas. Maybe you will getrembremading said:
better advice
in sci.math.num-analysis, a group that is dedicated to this kind
of problems. I have seen discussions about problems of this type there,
and they have a math level that goes WAY beyond what we have here in
comp.lang.c.
I think you have to rearrange the terms in your equation to make them
less sensible to roundoff errors. I think in comp.math.num-analysis
you will find people that can help you to do that.
K
Keith Thompson
I was referring to jacob's earlier statement, which did not mention
precision:
| Since the standard doesn't force even IEEE754, there is *nothing*
| the C language by itself can guarantee the user.
Look again, I *agreed* with some of what jacob was saying; the
C standard doesn't make any guarantees about the precsion of
floating-point operations.
K
Keith Thompson
Keith Thompson said:
And, to be clear, we should have been referring to accuracy rather
than precision.
K
Keith Thompson
Keith Thompson said:
And, to be clear, we should have been referring to accuracy rather
than precision.
C
Chris Torek
It also notes, in passing, that setting the precision field in
the FPU control word does *not* affect the exponent range. This
means that the in-FPU computations are *still* not quite what
one might expect from an IEEE implementation.
I have the same note (with an example) at the
end of <http://web.torek.net/torek/c/numbers.html>.
As someone else noted else-thread, one can get "proper" IEEE floating
point on the x86 using gcc's "-ffloat-store" flag. This does,
however, slow down computation.