bits for floating value and converting

[myfem.analyzer]

Hi,

I saw a line of codes in the "Fast Inverse Square Root" like this:

Float InvSqrt(float x)
{
.......
int i = *(int *) &x; //get bits for floating value
.....
x = * (float *) &i; //convert bit back to float


If x =1, &x = 0x0013fe54 and *(int *) &x = 1065353216. Can anyone
help me why I got these numbers? If I want to do the same operation in
Fortran, how can I do it?


Thanks a lot.

Haven

 

[Walter Roberson]

I saw a line of codes in the "Fast Inverse Square Root" like this:

Float InvSqrt(float x)
{
......
int i = *(int *) &x; //get bits for floating value
....
x = * (float *) &i; //convert bit back to float

If x =1, &x = 0x0013fe54 and *(int *) &x = 1065353216. Can anyone
help me why I got these numbers?

&x is the memory address of the variable x. The exact address you
get may vary from invocation to invocation. That you got
0x0013fe54 was more or less chance. The exact value that you got
has no particular significance.

*(int *) &x means to reinterpret the memory location x
as an int. This operation was not certain to succeed: you -might-
have been running on a system in which int requires stricter
alignment than float. And in C in general, the length of a float
is not necessarily the same as the length of an int -- it is not
uncommon for the two to be the same length, but sometimes int
is shorter than float, and sometimes int is longer than float.
It looks like it happened that on your compiler, the two are
the same length.

The exact value you got from the reinterpretation, 1065353216,
is equivilent to hex 0x3f800000 . Through experience, I recognize
that as the IEEE 754 single precision floating point
representation of 1.0 . You can find out more about IEEE 754 at
http://en.wikipedia.org/wiki/IEEE_floating-point_standard

If I want to do the same operation in
Fortran, how can I do it?

That's a Fortran question rather than a C question. The answer
is going to depend upon which Fortran standard you are using,
and possibly upon system extensions. Historically, you could
often play tricks like that using the %VAL and %REF FORTRAN
extensions; I see a hint that this might have become VAL()
in later versions, but you'd need to check that in your documentation
and in comp.lang.fortran

You might find the most straight-forward way in Fortran to get at the
bit representations is to push the value into a common block declared
one way, and then call a function that declares the common block
a different way (e.g., with an INTEGER*4 variable in that location.)

 

[CBFalconer]

I saw a line of codes in the "Fast Inverse Square Root" like this:

Float InvSqrt(float x) {
......
int i = *(int *) &x; //get bits for floating value
....
x = * (float *) &i; //convert bit back to float

If x =1, &x = 0x0013fe54 and *(int *) &x = 1065353216. Can anyone
help me why I got these numbers? If I want to do the same operation
in Fortran, how can I do it?

The code is illegal, or at least system dependent. Conversion of
float* to int and of int to float* is undefined in C.

 

[Keith Thompson]

The code is illegal, or at least system dependent. Conversion of
float* to int and of int to float* is undefined in C.

I suggest being careful with the word "illegal". I don't think the C
standard uses that term.

The quoted code does not violate any syntax rule or constraint, so in
that sense it's legal; an implementation is not required to reject it,
or even to issue any kind of diagnostic.

Conversion between int* and float* isn't undefined, but it can result
in Bad Things Happening. Specifically (C99 6.3.2.3p7):

A pointer to an object or incomplete type may be converted to a
pointer to a different object or incomplete type. If the resulting
pointer is not correctly aligned for the pointed-to type, the
behavior is undefined. Otherwise, when converted back again,
the result shall compare equal to the original pointer. When a
pointer to an object is converted to a pointer to a character
type, the result points to the lowest addressed byte of the
object. Successive increments of the result, up to the size of
the object, yield pointers to the remaining bytes of the object.

So the above can certainly cause problems if int and float have
different alignment requirements, but it *can* work as expected on
some systems. But it's highly non-portable.

 

[Pierre Asselin]

Float InvSqrt(float x)
{
......
int i = *(int *) &x; //get bits for floating value
....
x = * (float *) &i; //convert bit back to float

[ ... ] If I want to do the same operation in
Fortran, how can I do it?

Fortran has a "transfer" statement for this kind
of type punning. It's going to be machine-dependent
in Fortran just as in C.

 

[myfem.analyzer]

Hi, Walter Roberson:

Thank you for your quick reply and details. The whole source codes of
the fast square root algorithm are:

float InvSqrt(float x)
{
float xhalf = 0.5f*x;
int i = *(int*)&x; // get bits for floating value
i = 0x5f3759df - (i>>1); // gives initial guess y0
x = *(float*)&i; // convert bits back to float
x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases
accuracy
return x;
}

Based on your explanation, *(int*)&x will get the "IEEE 754 single
precision floating point Representation" of x, and *(float*)&i will
convert that representation back to a floating number?

Also, what is the algorithm to covert between "IEEE 754 single
precision floating point Representation" and normal floating number?

Thanks a lot.

Sincerely yours,

Haven

 

[ymuntyan]

I suggest being careful with the word "illegal". I don't think the C
standard uses that term.

The quoted code does not violate any syntax rule or constraint, so in
that sense it's legal; an implementation is not required to reject it,
or even to issue any kind of diagnostic.

Conversion between int* and float* isn't undefined, but it can result
in Bad Things Happening. Specifically (C99 6.3.2.3p7):

A pointer to an object or incomplete type may be converted to a
pointer to a different object or incomplete type. If the resulting
pointer is not correctly aligned for the pointed-to type, the
behavior is undefined. Otherwise, when converted back again,
the result shall compare equal to the original pointer. When a
pointer to an object is converted to a pointer to a character
type, the result points to the lowest addressed byte of the
object. Successive increments of the result, up to the size of
the object, yield pointers to the remaining bytes of the object.

So the above can certainly cause problems if int and float have
different alignment requirements, but it *can* work as expected on
some systems. But it's highly non-portable.

It violates aliasing rules, and not just in theory:

muntyan@munt10:/tmp$ cat file.c
#include <stdio.h>

int main (void)
{
float a = 2.;
int b = *(int*)&a;
b = ~b;
a = *(float*)&b;
printf ("%f\n", a);
return 0;
}
muntyan@munt10:/tmp$ gcc file.c # no optimizations
muntyan@munt10:/tmp$ ./a.out
-2.000000
muntyan@munt10:/tmp$ gcc file.c -O2
muntyan@munt10:/tmp$ ./a.out
0.000000

Best regards,
Yevgen

 

[ymuntyan]

Hi, Walter Roberson:

Thank you for your quick reply and details. The whole source codes of
the fast square root algorithm are:

float InvSqrt(float x)
{
float xhalf = 0.5f*x;
int i = *(int*)&x; // get bits for floating value
i = 0x5f3759df - (i>>1); // gives initial guess y0
x = *(float*)&i; // convert bits back to float
x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases
accuracy
return x;

}

It's broken, you need to do something like

float InvSqrt(float x)
{
union {
float f;
int i;
} u = {x};
float xhalf = 0.5f*u.f;
u.i = 0x5f3759df - (u.i>>1);
u.f = u.f*(1.5f-xhalf*u.f*u.f);
return u.f;
}

It is broken too, but at least it's guaranteed to work with gcc, for
what it's worth.

Based on your explanation, *(int*)&x will get the "IEEE 754 single
precision floating point Representation" of x, and *(float*)&i will
convert that representation back to a floating number?

Also, what is the algorithm to covert between "IEEE 754 single
precision floating point Representation" and normal floating number?

No idea about IEEE 754 part.

Best regards,
Yevgen

 

[Walter Roberson]

Hi, Walter Roberson:

Based on your explanation, *(int*)&x will get the "IEEE 754 single
precision floating point Representation" of x, and *(float*)&i will
convert that representation back to a floating number?

Also, what is the algorithm to covert between "IEEE 754 single
precision floating point Representation" and normal floating number?

What I said was,

That is, on *your* system, *(int*)&x gave you a particular result, and I
recognize that result as being IEEE 754 single precision representation.
In other words, your "normal floating point numbers" -happen- to
be internally represented in IEEE 754 single precision format.
They already -are- in that format on your machine. On other machines
that did not happen to use IEEE 754 but which did happen to have
sizeof(int) == sizeof(float), then *(int*)&x would get at whatever
representation happened to be used on -that- machine. Perhaps VAX-F
format.

There is no question about algorithms to convert between IEEE 754
and "normal floating number" because "normal floating numbers" do not
exist as an independant entity. There are merely a whole bunch of
different standard (and not so standard) ways of representing floating
point numbers, and whatever system you are running on will use -one-
of them (actually, it might even use several different ones).
IEEE 754 is as "normal" as you get in the general purpose computing world.

*(int*)&x did not change anything about the bits themselves at the
location of the variable x: it just allowed you to examine them in
a different format. *(float*)&i would not change anything about
the bits themselves at the locaton of the variable i: it would just
allow you to exmaine them in a different format. Though if you
used *(float*)&i = SomeValue and i happened to be properly
aligned and the same size as a float, then the floating point
representation of SomeValue would be stored at the location of the
variable i .

 

[David Thompson]

Float InvSqrt(float x)
{
......
int i = *(int *) &x; //get bits for floating value
....
x = * (float *) &i; //convert bit back to float

[ ... ] If I want to do the same operation in
Fortran, how can I do it?

Fortran has a "transfer" statement for this kind
of type punning. It's going to be machine-dependent
in Fortran just as in C.

TRANSFER intrinsic _function_; since F90 (if that matters).
Or, since forever, the EQUIVALENCE statement (declaration).

Fortran (including FORTRAN) does require single-precision float,
called REAL, and 'default-kind' INTEGER, to be the same size,
the equivalent of which C does not, as already noted elsethread.
But like C it does not require any particular representation of
floats, and actually less than C for integers.

But as already noted, the place to ask about the Fortran part is
comp.lang.fortran . The what-does-the-C-do part is (was) ontopic here,
but in practice would also be accepted in c.l.f even though it's not
Fortran; first they have much lower volume and can thus afford to be
less strict about topicality; and second most Fortran implementations
have specific features for C interoperability (though not standardized
until recently) whereas most C implementations don't do anything for
interoperability with other languages.

- formerly david.thompson1 || achar(64) || worldnet.att.net

 

[Tim Prince]

Hi,

I saw a line of codes in the "Fast Inverse Square Root" like this:

Float InvSqrt(float x)
{
......
int i = *(int *) &x; //get bits for floating value
....
x = * (float *) &i; //convert bit back to float


If x =1, &x = 0x0013fe54 and *(int *) &x = 1065353216. Can anyone
help me why I got these numbers? If I want to do the same operation in
Fortran, how can I do it?

Presumably, this was meant to be K&R (pre-89 standard) C. Not that many
C programmers caught up for the first decade. If you used unions to
represent data which could be interpreted in either int or float
representation, the analogy to Fortran EQUIVALENCE might be clearer.
Since 1990, TRANSFER is the recommended (and less arcane) Fortran syntax
for copying a bit pattern between data types.
Most common CPU architectures include a fast approximate inverse square
root instruction. It is not accessed by type punning or similar tricks,
which defeat the description "fast" by requiring moves between data
types which aren't supported efficiently. The need for it is greatly
reduced in CPUs introduced this year (all brands of which I know) which
support more efficient sqrt() without approximations.
So, it appears the code of which you speak is of mainly historic
interest, and may not even work the same on any post-1990 CPU as it may
have done when it was written.

 

[Tim Prince]

It's broken, you need to do something like

float InvSqrt(float x)
{
union {
float f;
int i;
} u = {x};
float xhalf = 0.5f*u.f;
u.i = 0x5f3759df - (u.i>>1);
u.f = u.f*(1.5f-xhalf*u.f*u.f);
return u.f;
}

It is broken too, but at least it's guaranteed to work with gcc, for
what it's worth.

The only improvement here is the demonstration that the same thing can
be done under the 1989 standard as in the earlier form of C.
Even if your CPU happens to have the same correspondence between int and
float bit patterns as the CPU for which this was designed, there is
absolutely no point in this on a CPU which includes some form of
firmware sqrt(). I was going to say there might be academic interest,
but it is broken in that respect, in the absence of documentation of the
CPU for which it was designed.