30

The code below performs a fast inverse square root operation by some bit hacks. The algorithm was probably developed by Silicon Graphics in early 1990's and it's appeared in Quake 3 too. more info

However I get the following warning from GCC C++ compiler: dereferencing type-punned pointer will break strict-aliasing rules

Should I use static_cast, reinterpret_cast or dynamic_cast instead in such situations?

float InverseSquareRoot(float x)
{
    float xhalf = 0.5f*x;
    int32_t i = *(int32_t*)&x;
    i = 0x5f3759df - (i>>1);
    x = *(float*)&i;
    x = x*(1.5f - xhalf*x*x);
    return x;
}
5
  • The paper gives the value 0x5fe6ec85e7de30da as the optimal constant for doubles in the conclusion. Jul 22, 2013 at 14:49
  • "Lomont pointed out that the "magic number" for 64 bit IEEE754 size type double is 0x5fe6ec85e7de30da, but in fact it was shown to be exactly 0x5fe6eb50c7b537a9" (from Wikipedia) Maybe I'll do some tests.
    – plasmacel
    Jul 22, 2013 at 14:51
  • Oh interesting. I didn't know about McEniry's paper. Thanks :) Jul 22, 2013 at 14:53
  • 1
    This trick is completely free of undefined behavior, assuming the float and int bits are compatible to each other: union A { float x; int32_t y; }; int32_t value = A{3.14f}.y; the reason is, that the aliasing rule is never applied here. However, this trick only works up to C++11 (I think C++14 will change the rules so that this becomes UB). Jul 23, 2013 at 19:59
  • 1
    Your pseudo_cast in C++ should also static_assert that both types are trivially copiable (see std::is_trivially_copyable).
    – Ruslan
    Mar 12, 2018 at 16:53

8 Answers 8

40

Forget casts. Use memcpy.

float xhalf = 0.5f*x;
uint32_t i;
assert(sizeof(x) == sizeof(i));
std::memcpy(&i, &x, sizeof(i));
i = 0x5f375a86 - (i>>1);
std::memcpy(&x, &i, sizeof(i));
x = x*(1.5f - xhalf*x*x);
return x;

The original code tries to initialize the int32_t by first accessing the float object through an int32_t pointer, which is where the rules are broken. The C-style cast is equivalent to a reinterpret_cast, so changing it to reinterpret_cast would not make much difference.

The important difference when using memcpy is that the bytes are copied from the float into the int32_t, but the float object is never accessed through an int32_t lvalue, because memcpy takes pointers to void and its insides are "magical" and don't break the aliasing rules.

17
  • 6
    @Jan Why is it “obviously” slower? I don’t find that obvious at all. Jul 22, 2013 at 14:27
  • 3
    @plasmacel Mostly because it works. I wouldn't expect any discernible difference in performance from a reasonable compiler. (Yes, I am sure MSVC will prove me wrong...) Jul 22, 2013 at 14:31
  • 9
    It's only slower if the compiler does a horrible job of optimizing it. Here's proof that gcc 4.8 compiles the cast and memcpy version to the same assembly. And unknown version of clang.
    – Casey
    Jul 22, 2013 at 14:31
  • 12
    I fear that my "proof that it compiles to the same assembly" comment above will be construed by someone to indicate that it's reasonable to write code that breaks the aliasing rules. Realize that this is relying on undefined behavior, and just as in any other case of undefined behavior it may produce different results with different compilers, or different compiler versions, or on Tuesday instead of Monday, .... etc. The most dangerous possible outcome of undefined behavior is to do what you expect during testing.
    – Casey
    Jul 22, 2013 at 14:49
  • 2
    @MikeSeymour, C89 said, "An object shall have its stored value accessed only by an lvalue that has one of the following types: the declared type of the object, a qualified version of the declared type of the object, [signed or unsigned variants], an aggregate or union type that includes one of the aforementioned types among its members [], or a character type" in 6.3 Expressions. That's basically the same as the C99 wording and bans the original code. Jul 25, 2013 at 15:10
13

There are a few good answers here that address the type-punning issue.

I want to address the "fast inverse square-root" part. Don't use this "trick" on modern processors. Every mainstream vector ISA has a dedicated hardware instruction to give you a fast inverse square-root. Every one of them is both faster and more accurate than this oft-copied little hack.

These instructions are all available via intrinsics, so they are relatively easy to use. In SSE, you want to use rsqrtss (intrinsic: _mm_rsqrt_ss( )); in NEON you want to use vrsqrte (intrinsic: vrsqrte_f32( )); and in AltiVec you want to use frsqrte. Most GPU ISAs have similar instructions. These estimates can be refined using the same Newton iteration, and NEON even has the vrsqrts instruction to do part of the refinement in a single instruction without needing to load constants.

0
8

If you have access to C++20 or later then you can use std::bit_cast

float InverseSquareRoot(float x)
{
    float xhalf = 0.5f*x;
    int32_t i = std::bit_cast<int32_t>(x);
    i = 0x5f3759df - (i>>1);
    x = std::bit_cast<float>(i);
    x = x*(1.5f - xhalf*x*x);
    return x;
}

At the moment std::bit_cast is only supported by MSVC. See demo on Godbolt

While waiting for the implementation, if you're using Clang you can try __builtin_bit_cast. Just change the casts like this

int32_t i = __builtin_bit_cast(std::int32_t, x);
x = __builtin_bit_cast(float, i);

Demo

7

Update

I no longer believe this answer is correct, due to feedback I've gotten from the committee. But I want to leave it up for informational purposes. And I am purposefully hopeful that this answer can be made correct by the committee (if it chooses to do so). I.e. there's nothing about the underlying hardware that makes this answer incorrect, it is just the judgement of a committee that makes it so, or not so.


I'm adding an answer not to refute the accepted answer, but to augment it. I believe the accepted answer is both correct and efficient (and I've just upvoted it). However I wanted to demonstrate another technique that is just as correct and efficient:

float InverseSquareRoot(float x)
{
    union
    {
        float as_float;
        int32_t as_int;
    };
    float xhalf = 0.5f*x;
    as_float = x;
    as_int = 0x5f3759df - (as_int>>1);
    as_float = as_float*(1.5f - xhalf*as_float*as_float);
    return as_float;
}

Using clang++ with optimization at -O3, I compiled plasmacel's code, R. Martinho Fernandes code, and this code, and compared the assembly line by line. All three were identical. This is due to the compiler's choice to compile it like this. It had been equally valid for the compiler to produce different, broken code.

17
  • 1
    @jogojapan: Good comment. I thought it important to show code. English can be so vague and ambiguous. On your second point, I'm glad you brought it up. I think the C++ committee needs to address this issue and recommend safe techniques for type punning. There have been discussions in the Library Working Group specifically requesting direction from the Core or Evolution Working Group on how this should best be accomplished. These questions, to the best of my knowledge, have so far gone unanswered. However, my belief at this time is that compiler writers will not violate the C99 contract. Jul 23, 2013 at 3:11
  • 4
    As far as I'm aware memcpy is the recommended approach for type punning in C++. C++ compilers may maintain the C99 union guarantee, however it's still undefined behavior in C++.
    – bames53
    Jul 23, 2013 at 16:04
  • 2
    For up to C++11, you can use my above union trick to be free of aliasing issues: union A { float x; int32_t y; }; int32_t value = A{3.14f}.y; (I don't really think that this is anymore safer than doing it without the temporary :D). The reason this "works" is that the initializer is a prvalue and hence is free of the aliasing rule restrictions. This however will change in C++14 because the initializer will be an xvalue :) Jul 23, 2013 at 20:05
  • 2
    @JohannesSchaub-litb I don't believe the strict aliasing rules are the only restriction on type punning with unions in C++. For example I believe your example runs afoul of 8.5 [dcl.init]/16 where it says "the initial value of the object being initialized is the (possibly converted) value of the initializer expression." since the object designated by A{3.14f}.y does not have any value. C++ omits specification of behavior in this situation and so the behavior is undefined.
    – bames53
    Jul 23, 2013 at 21:29
  • 1
    @JohannesSchaub-litb And by 'C99 guarantee' I mean the guarantee that reading a non-active member will reinterpret the stored representation of the active member. In C11 they added an explicit note to that effect, but I haven't looked at the C spec enough to know if the normative language of the C99 spec truly does imply what the C11 note claims.
    – bames53
    Jul 24, 2013 at 1:04
1

Take a look at this for more information on type punning and strict aliasing.

The only safe cast of a type into an array is into a char array. If you want one data address to be switchable to different types you will need to use a union

14
  • 1
    Using a union for type punning doesn't conform with the Standard either (although it typically works).
    – jogojapan
    Jul 22, 2013 at 14:40
  • 1
    @jogojapan - strictly speaking type-punning unions do violate the strict aliasing rule and are not defined but they are a common enough idiom for GCC to support them (and no generated the warning)
    – doron
    Jul 22, 2013 at 14:50
  • 2
    The strict aliasing rule doesn't mention them (at least not explicitly). But the rules for unions state that you must only read from the member of the union that was last set. For type punning you have to write one member and then read another.
    – jogojapan
    Jul 22, 2013 at 14:53
  • 1
    @jogojapan 3.10/10 does explicitly mention unions in bullet 6. Accessing the value stored in an object through a glvalue whose type is a union containing the dynamic type of the object as a member is defined behavior. You may violate the "last member set" rule, but it does not violate aliasing. </nitpick>
    – Casey
    Jul 22, 2013 at 15:03
  • 4
    I've opened a query about this on the Core Working Group mailing list. Jul 23, 2013 at 19:59
1

The cast invokes undefined behaviour. No matter what form of cast you use, it will still be undefined behaviour. It is undefined no matter what type of cast you use.

Most compilers will do what you expect, but gcc likes being mean and is probably going to assume you didn't assign the pointers despite all indication you did and reorder the operation so they give some strange result.

Casting a pointer to incompatible type and dereferencing it is an undefined behaviour. The only exception is casting it to or from char, so the only workaround is using std::memcpy (as per R. Martinho Fernandes' answer). (I am not sure how much it is defined using unions; It does stand a better chance of working though).

That said, you should not use C-style cast in C++. In this case, static_cast would not compile, nor would dynamic_cast, forcing you to use reinterpret_cast and reinterpret_cast is a strong suggestion you might be violating strict aliasing rules.

16
  • 2
    There are different reasons for undefined behavior. In this case, of course, there's no way the standards committee could define it, particularly in the face of possible trapping representations for float (usual) or int (rare, but can exist). On the other hand, the intend is clearly that it should have the behavior which someone familiar with the architecture would expect. A compiler which breaks this, when the cast is immediately visible, is simply broken. Jul 22, 2013 at 14:30
  • 3
    @JamesKanze: Compiler that breaks this when the cast is immediately visible is conforming to the specification. Gcc does just that (in fact I think it won't break it as written, because the int pointer is cast back, but it would almost certainly read the original value back from the original variable, because it deleted the dependency and would consider the value left in floating pointer register valid.)
    – Jan Hudec
    Jul 22, 2013 at 14:35
  • 2
    @JamesKanze: Better to break right away in ALL situations rather than work sometimes and not others. If you move the cast out of the function and it silently breaks would you be happy?
    – Zan Lynx
    Jul 22, 2013 at 14:52
  • 1
    @JamesKanze: No, the memcpy technique is defined. The exact values are implementation defined, but they must be the value represented by the same sequence of bytes in memory. On the other hand type punning is undefined and the compiler may do anything it pleases. For union I don't remember the specification exactly, but the rules are basically similar.
    – Jan Hudec
    Jul 23, 2013 at 6:02
  • 1
    @JamesKanze: No. You are mixing two very different undefineds. The behaviour of reinterpret cast is "Undefined". That means the compiler is free to do anything ranging from what you expect to formatting your harddrive and making daemons flying from your nose. But the memcpy behaviour is "Implementation Defined". The specification says the implementation must copy the actual bytes and only does not say what the value represented by the same bits shall mean. Yes, it may cause fault if the value is a signalling NaN, but it may not involve any daemons. And must be deterministic.
    – Jan Hudec
    Jul 23, 2013 at 8:03
1

Based on the answers here I made a modern "pseudo-cast" function for ease of application.

C99 version (while most compilers support it, theoretically could be undefined behavior in some)

template <typename T, typename U>
inline T pseudo_cast(const U &x)
{
    static_assert(std::is_trivially_copyable<T>::value && std::is_trivially_copyable<U>::value, "pseudo_cast can't handle types which are not trivially copyable");

    union { U from; T to; } __x = {x};
    return __x.to;
}

Universal versions (based on the accepted answer)

Cast types with the same size:

#include <cstring>

template <typename T, typename U>
inline T pseudo_cast(const U &x)
{
    static_assert(std::is_trivially_copyable<T>::value && std::is_trivially_copyable<U>::value, "pseudo_cast can't handle types which are not trivially copyable");
    static_assert(sizeof(T) == sizeof(U), "pseudo_cast can't handle types with different size");
    
    T to;
    std::memcpy(&to, &x, sizeof(T));
    return to;
}

Cast types with any sizes:

#include <cstring>

template <typename T, typename U>
inline T pseudo_cast(const U &x)
{
    static_assert(std::is_trivially_copyable<T>::value && std::is_trivially_copyable<U>::value, "pseudo_cast can't handle types which are not trivially copyable");

    T to = T(0);
    std::memcpy(&to, &x, (sizeof(T) < sizeof(U)) ? sizeof(T) : sizeof(U));
    return to;
}

Use it like:

float f = 3.14f;
uint32_t u = pseudo_cast<uint32_t>(f);

Update for C++20

C++20 introduces constexpr std::bit_cast in header <bit> which is functionally equivalent for types with the same size. Nevertheless, the above versions are still useful if you want to implement this functionality yourself (supposed that constexpr is not required), or if you want to support types with different sizes.

2
  • C++20 has this as std::bit_cast<T> in #include <bit> Apr 9, 2021 at 22:24
  • @PeterCordes Yeah, that makes this somewhat obsolete. However it is still useful if you want to implement that yourself - supposed that constexpr is not required. I will update the answer.
    – plasmacel
    Apr 9, 2021 at 23:07
0

The only cast that will work here is reinterpret_cast. (And even then, at least one compiler will go out of its way to ensure that it won't work.)

But what are you actually trying to do? There's certainly a better solution, that doesn't involve type punning. There are very, very few cases where type punning is appropriate, and they all are in very, very low level code, things like serialization, or implementing the C standard library (e.g. functions like modf). Otherwise (and maybe even in serialization), functions like ldexp and modf will probably work better, and certainly be more readable.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.