Gray codes addition

Question

As some of you know, I have been working with Gray codes and am trying to design a good gray_code class in modern C++. I will end up posting several parts of the implementation, but what I want to be reviewed is mostly the addition function. Anyway, here is the code of the class essentials:

template<typename Unsigned>
struct gray_code
{
    static_assert(std::is_unsigned<Unsigned>::value,
                  "gray code only supports built-in unsigned integers");

    // Underlying unsigned integer
    using value_type = Unsigned;

    // Variable containing the gray code
    value_type value;

    ////////////////////////////////////////////////////////////
    // Constructors operations

    // Default constructor
    constexpr gray_code() noexcept:
        value(0u)
    {}

    /**
     * @brief Construction from an unsigned integer.
     *
     * The integer is converted to gray code. The value
     * is preserved. The representation is not.
     *
     * @param value Unsigned integer to convert
     */
    constexpr explicit gray_code(value_type value) noexcept:
        value( (value >> 1) ^ value )
    {}

    ////////////////////////////////////////////////////////////
    // Assignment operations

    /**
     * @brief Assigns an unsigned integer.
     *
     * It works the same way as the equivalent
     * constructor does.
     */
    auto operator=(value_type other) & noexcept
        -> gray_code&
    {
        value = (other >> 1) ^ other;
        return *this;
    }

    ////////////////////////////////////////////////////////////
    // Conversion operations

    /**
     * @brief Conversion to the underlying type.
     *
     * See http://www.dspguru.com/dsp/tricks/gray-code-conversion
     */
    operator value_type() const noexcept
    {
        value_type res = value;
        for (value_type mask = std::numeric_limits<value_type>::digits / 2
             ; mask ; mask >>= 1)
        {
            res ^= res >> mask;
        }
        return res;
    }
};

As you can see, a gray_code<T> is just a wrapper around an unsigned built-in type T. The main operation are the construction and the conversion to and from the type T; the underlying value can be accessed thanks to the public variable value. The purpose of this class is to be easy to use. Thanks to the conversion operator, a gray_code can be used wherever the corresponding integer can be used. One of my goals is to provide operations specialized for gray_code only if the operation can be faster than converting the Gray code to the corresponding integer, performing the operation on the integer and converting the integer back to a Gray code. Conversions are really fast, so operations that are faster than the double conversion are hard to write. I did write some of them, not included here (comparison, bitwise operations, conversion to bool, increment and decrement; you can assume that these functions work between gray_codes for the rest of the question) and I am currently working on the addition.

My addition algorithm is based on the pseudo-code algorithm provided in this paper by R. W. Doran (am I the only one there thinking of Dorian Gray?), and attributed to Harold Lucal. Here is the original pseudo-code addition algorithm (⊕ means xor):

procedure add (n: integer; A,B:word; PA,PB:bit;
               var S:word; var PS:bit; var CE, CF:bit);
var i: integer; E, F, T: bit;
begin
   E := PA; F := PB;
   for i:= 0 to n-1 do begin {in parallel, using previous inputs}
       S[i] := (E and F) ⊕ A[i] ⊕ B[i];
       E := (E and (not F)) ⊕ A[i];
       F := ((not E) and F) ⊕ B[i];
   end;
   CE := E; CF := F;
end;

I translated the algorithm in C++. Here is the naive implementation (your can find the corresponding C version here):

template<typename Unsigned>
auto operator+(gray_code<Unsigned> lhs, gray_code<Unsigned> rhs) noexcept
    -> gray_code<Unsigned>
{
    // parity of lhs and rhs
    bool lhs_p = is_odd(lhs);
    bool rhs_p = is_odd(rhs);

    gray_code<Unsigned> res{};
    for (Unsigned i{} ; i < std::numeric_limits<Unsigned>::digits ; ++i)
    {
        // Get the ith bit of rhs and lhs
        bool lhs_i = (lhs.value >> i) & 1u;
        bool rhs_i = (rhs.value >> i) & 1u;

        // Copy lhs_p and rhs_p (see {in parallel} in the original algorithm)
        bool lhs_p_cpy = lhs_p;
        bool rhs_p_cpy = rhs_p;

        // Set the ith bit of res
        Unsigned res_i = (lhs_p_cpy & rhs_p_cpy) ^ lhs_i ^ rhs_i;
        res |= (res_i << i);

        // Update e and f
        lhs_p = (lhs_p_cpy & (!rhs_p_cpy)) ^ lhs_i;
        rhs_p = ((!lhs_p_cpy) & rhs_p_cpy) ^ rhs_i;
    }
    return res;
}

This algorithm works. However, it is not as efficient as converting the Gray codes to integers, performing a regular addition, and converting the result back to a Gray code.

---

In case you want it, here is the implementation of is_odd, even though it does not add that much overhead to the addition since it is only performed once. The parity of a Gray code corresponds to the parity of the number of bits set. The function uses compiler intrinsics when possible (which may be as simple as taking a parity bit on some architectures) and a \$ O(\log n) \$ algorithm to count bits otherwise.

template<typename Unsigned>
auto is_odd(gray_code<Unsigned> code)
    -> bool
{
    #if defined(__GNUC__) || defined(__clang__)

        return bool(__builtin_parity(code.value));

    #else

        unsigned nb_bits{};
        for (auto val = code.value ; val ; ++nb_bits)
        {
            // clear the least significant bit set
            val &= val - 1;
        }
        return bool(nb_bits % 2);

    #endif
}

Note: if you found this Q&A interesting an are interesting in another algorithm to add Gray codes, I made another Q&A about a custom algorithm to add Gray codes based on a powers of \$2\$ decomposition.

Answer 1

Loop invariants code motion

This is the name of an optimization performed by most compilers: when they detect code that actually does not depend on the state of the loop, they move it out of the loop. It's not as obvious as the usual invariants, but let's have a look at this line:

Unsigned res_i = (lhs_p_cpy & rhs_p_cpy) ^ lhs_i ^ rhs_i;

lhs_i and rhs_i correspond to the \$ i_{th} \$ bit of lhs and rhs. These bits do not depend on the loop (the loop does not change them) and can all be computed at once by performing lhs ^ rhs before running the loop. Therefore, we can simplify the function by moving them out of the loop. If we store them directly in res, we can even simplify the assignment to res to a mere res ^= res_i << i:

template<typename Unsigned>
auto operator+(gray_code<Unsigned> lhs, gray_code<Unsigned> rhs) noexcept
    -> gray_code<Unsigned>
{
    // parity of lhs and rhs
    bool lhs_p = is_odd(lhs);
    bool rhs_p = is_odd(rhs);

    gray_code<Unsigned> res = lhs ^ rhs;
    for (Unsigned i{} ; i < std::numeric_limits<Unsigned>::digits ; ++i)
    {
        // Get the ith bit of rhs and lhs
        bool lhs_i = (lhs.value >> i) & 1u;
        bool rhs_i = (rhs.value >> i) & 1u;

        // Copy lhs_p and rhs_p (see {in parallel} in the original algorithm)
        bool lhs_p_cpy = lhs_p;
        bool rhs_p_cpy = rhs_p;

        // Set the ith bit of res
        Unsigned res_i = lhs_p_cpy & rhs_p_cpy;
        res ^= res_i << i;

        // Update e and f
        lhs_p = (lhs_p_cpy & (!rhs_p_cpy)) ^ lhs_i;
        rhs_p = ((!lhs_p_cpy) & rhs_p_cpy) ^ rhs_i;
    }
    return res;
}

Unfortunately, lhs_i and rhs_i are also used to compute the new values of lhs_p and rhs_p, so we can't remove them from the body of the loop anyway. If we are to optimize the algorithm even more, we will need something else, for example...

More efficient parity algorithm

It hardly changes anything since is_odd is only used at the very beginning of the algorithm, but its fallback implementation when __builtin_parity cannot be used is suboptimal. Actually, the parity of an unsigned integer can be obtained by xoring together every bit, without having to compute the number of set bits. The overall xoring of bits can be done "in parallel" instead of bit by bit with the following algorithm:

for (std::size_t i = std::numeric_limits<Unsigned>::digits / 2u ; i ; i >>= 1u)
{
    code.value ^= (code.value >> i);
}
return bool(code.value & 1u);

This algorithm will xor the first half of an unsigned integer with the second half, then repeat the operation with the result... until there are only two bits left to xor, leading to a more efficient operation than a population count. This page even provides a more subtle implementation of the algorithm to "skip" the last iterations. I didn't really dig the proposed implementations with magic numbers since I want my code to still work with integers of any size.

Getting rid of the temporaries

The part of the Gray addition which is likely to be be the slowest is the body of the for loop. To find a possible optimization, I decided to see wether I could get rid of lhs_p_cpy and rhs_p_cpy in it, which are copies of lhs_p and rhs_p required to compute the new value of rhs_p. While it is easy to get rid of rhs_p_cpy by simply replacing it by rhs_p, getting rid of lhs_p_cpy is harder because lhs_p has already been updated meanwhile. Without any trickery, the following loop body is the best we can get:

Unsigned res_i = lhs_p & rhs_p;
res ^= res_i << i;

bool lhs_p_cpy = lhs_p;
bool lhs_i = (lhs.value >> i) & 1u;
bool rhs_i = (rhs.value >> i) & 1u;
lhs_p = (lhs_p_cpy & not rhs_p) ^ lhs_i;
rhs_p = (rhs_p & not lhs_p_cpy) ^ rhs_i;

Therefore, I decided to build some truth tables. In the following table, \$ lhs_{old} \$ and \$ rhs_{old} \$ are possible values of lhs_p and rhs_p before the update while \$ lhs_{new} \$ and \$ rhs_{new} \$ are the values of the same variables after the update. I did not take any interest in lhs_i and rhs_i since none of them is used to compute both \$ lhs_{new} \$ and \$ rhs_{new} \$; therefore, I simply ignore them in the rest of this reflection.

\begin{array} {|cc|cc|} \hline lhs_{old} & rhs_{old} & lhs_{new} & rhs_{new} \\ \hline 0 & 0 & 0 & 0\\ 0 & 1 & 0 & 1\\ 1 & 0 & 1 & 0\\ 1 & 1 & 0 & 0\\ \hline \end{array}

We can easily replace lhs_p_cpy by \$ lhs_{old} \$ in the computation of \$ lhs_{new} \$ but we can't use it to compute \$ rhs_{new} \$ since the update already occured. My first thought to get rid of lhs_p_cpy was "can we compute \$ rhs_{new} \$ with only \$ rhs_{old} \$ and \$ lhs_{new} \$?". Looking at the truth table above, it appears that it is not possible. Trying to compute \$ rhs_{new} \$ before \$ lhs_{new} \$ doesn't solve the problem either.

However, we already have computed another value in the loop: res_i. Therefore, I injected \$ res_i \$ in the table and looked at what could be done with it:

\begin{array} {|cc|c|cc|} \hline lhs_{old} & rhs_{old} & res_i = lhs_{old} \land rhs_{old} & lhs_{new} & rhs_{new} \\ \hline 0 & 0 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 & 1\\ 1 & 0 & 0 & 1 & 0\\ 1 & 1 & 1 & 0 & 0\\ \hline \end{array}

We can infer from this "truth table" that \$ rhs_{old} \land \lnot lhs_{old} = rhs_{old} \land \lnot res_i \$. If we inject this new discovery back in the code, we can use it to totally get rid of lhs_p_cpy:

Unsigned res_i = lhs_p & rhs_p;
res ^= res_i << i;

bool lhs_i = (lhs.value >> i) & 1u;
bool rhs_i = (rhs.value >> i) & 1u;
lhs_p = (lhs_p & not rhs_p) ^ lhs_i;
rhs_p = (rhs_p & not res_i) ^ rhs_i;

By the way res_i can also be used to compute \$ lhs_{new} \$ so that the code looks more symmetrical (even though it does not improve the performance):

Unsigned res_i = lhs_p & rhs_p;
res ^= res_i << i;

bool lhs_i = (lhs.value >> i) & 1u;
bool rhs_i = (rhs.value >> i) & 1u;
lhs_p = (lhs_p & not res_i) ^ lhs_i;
rhs_p = (rhs_p & not res_i) ^ rhs_i;

With all optimizations turned on, this code is a little bit faster than the previous version, but it's still far from being as efficient as the version which converts the Gray code back to a regular integer, performs an integer addition and converts the result back to a Gray code. Hey, we can't optimize such an algorithm without twisting it in every way possible! This is why, we will also try to...

Find a better end condition

Currently, the algorithm loops over every bit in res before ending. It may be possible to end the loop sooner. Let's have a look at a big old truth table that we will use for the rest of the answer. For the sake of brevity, I used the names from the original pseudocode to name the columns. I hope that the parallel with the code I wrote isn't too hard to make (\$ A = lhs_i \$, \$ B = rhs_i \$, \$ E_{old} = lhs_p \$, \$ F_{old} = rhs_p \$, \$ S = A \oplus B \oplus res_i \$, \$ E_{new} = lhs_p' \$, \$ F_{new} = rhs_p' \$)

\begin{array} {|cccc|cc|cc|} \hline A & B & E_{old} & F_{old} & A \oplus B & S & E_{new} & F_{new}\\ \hline 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1\\ 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0\\ 0 & 0 & 1 & 1 & 0 & 1 & 0 & 0\\ 0 & 1 & 0 & 0 & 1 & 1 & 0 & 1\\ 0 & 1 & 0 & 1 & 1 & 1 & 0 & 0\\ 0 & 1 & 1 & 0 & 1 & 1 & 1 & 1\\ 0 & 1 & 1 & 1 & 1 & 0 & 0 & 1\\ \hline 1 & 0 & 0 & 0 & 1 & 1 & 1 & 0\\ 1 & 0 & 0 & 1 & 1 & 1 & 1 & 1\\ 1 & 0 & 1 & 0 & 1 & 1 & 0 & 0\\ 1 & 0 & 1 & 1 & 1 & 0 & 1 & 0\\ 1 & 1 & 0 & 0 & 0 & 0 & 1 & 1\\ 1 & 1 & 0 & 1 & 0 & 0 & 1 & 0\\ 1 & 1 & 1 & 0 & 0 & 0 & 0 & 1\\ 1 & 1 & 1 & 1 & 0 & 1 & 1 & 1\\ \hline \end{array}

Looking at the very first line of that table, we can see that when \$ A = B = E_{old} = F_{old} = 0 \$, then as a result, we have \$ S = E_{new} = F_{new} = 0 \$, which means that since our result was filled with \$ A \oplus B \$ from the start, then we can stop the loop when the aforementioned condition is met (everything equals \$0\$) since the following bits S will always be \$0\$. An easy way to do so is to "consume" lhs and rhs in the code with right shifts so that they eventually reach zero and to change the loop condition with bitwise ORs (here, only the for loop for the sake of simplicity):

for (Unsigned i{} ;
     lhs | rhs | lhs_p | rhs_p ;
     ++i, lhs >>= 1u, rhs >>= 1u)
{
    Unsigned res_i = lhs_p & rhs_p;
    res ^= res_i << i;

    bool lhs_i = lhs.value & 1u;
    bool rhs_i = rhs.value & 1u;
    lhs_p = (lhs_p & not res_i) ^ lhs_i;
    rhs_p = (rhs_p & not res_i) ^ rhs_i;
}

While this can reduce the complexity of the algorithm and greatly improve things for small values, lhs | rhs | lhs_p | rhs_p is still a bunch of operations by itself for a loop condition. Our goal is to stop the loop when all of these variables are \$0\$. One thing we can do is find before the loop which one between lhs and rhs has the highest set bit and loop until that one is consumed. Since they are both "consumed" at the same speed, when the one with the highest set bit is consumed, then it means that the other one has already been consumed, so we only have to check whether the biggest of the two is \$0\$. We can easily find which one has the highest bit by simply checking which one has the biggest value. Here is the algorithm once modified:

template<typename Unsigned>
auto operator+(gray_code<Unsigned> lhs, gray_code<Unsigned> rhs) noexcept
    -> gray_code<Unsigned>
{
    // Make sure lhs.value has the highest set bit
    if (lhs.value < rhs.value)
    {
        std::swap(lhs.value, rhs.value);
    }

    // parity of lhs and rhs
    bool lhs_p = is_odd(lhs);
    bool rhs_p = is_odd(rhs);

    gray_code<Unsigned> res = lhs ^ rhs;
    for (Unsigned i{} ;
         lhs | lhs_p | rhs_p ;
         ++i, lhs >>= 1u, rhs >>= 1u)
    {
        Unsigned res_i = lhs_p & rhs_p;
        res ^= res_i << i;

        bool lhs_i = lhs.value & 1u;
        bool rhs_i = rhs.value & 1u;
        lhs_p = (lhs_p & not res_i) ^ lhs_i;
        rhs_p = (rhs_p & not res_i) ^ rhs_i;
    }
    return res;
}

Have a better look at the truth table

From the big truth table we built, we can see that when \$ A = B = 0 \$, the only time when \$ S = 1 \$ is when \$ E_{old} = F_{old} = 1 \$. However, when \$ A = B = 0 \$, we can also see that \$ E_{new} \$ and \$ F_{new} \$ can never be \$1\$ at the same time. That means that once lhs and rhs are consumed, res can change only one more time and only if (lhs_p & rhs_p) is true. Therefore, we can extract the lhs_p | rhs_p part from the loop condition and apply exactly one res assignment after the loop (I used a while loop since the for loop was becoming a bit messy for the job):

template<typename Unsigned>
auto operator+(gray_code<Unsigned> lhs, gray_code<Unsigned> rhs) noexcept
    -> gray_code<Unsigned>
{
    // Make sure lhs.value has the highest set bit
    if (lhs.value < rhs.value)
    {
        std::swap(lhs, rhs);
    }

    // parity of lhs and rhs
    bool lhs_p = is_odd(lhs);
    bool rhs_p = is_odd(rhs);

    gray_code<Unsigned> res = lhs ^ rhs;
    Unsigned i{};
    while (lhs)
    {
        Unsigned res_i = lhs_p & rhs_p;
        res ^= res_i << i;

        bool lhs_i = lhs.value & 1u;
        bool rhs_i = rhs.value & 1u;
        lhs_p = (lhs_p & not res_i) ^ lhs_i;
        rhs_p = (rhs_p & not res_i) ^ rhs_i;

        ++i;
        lhs >>= 1u;
        rhs >>= 1u;
    }
    // Last value in case lhs_p and rhs_p are 1
    res ^= (lhs_p & rhs_p) << i;
    return res;
}

Sometimes, more is less

Currently, we use the fact that at some point lhs and rhs will eventually become \$0\$ and that only one modification can happen right after this point. However, we currently don't use the fact that at some point rhs or lhs is \$0\$. We can actually write two loops instead of one: one loop that runs until the smallest of lhs.value and rhs.value is \$0\$ and another that runs until the other one reaches \$0\$. The second loop can be simplified since we know that one rhs (the smallest value from the std::swap) is \$0\$ during the iteration:

template<typename Unsigned>
auto operator+(gray_code<Unsigned> lhs, gray_code<Unsigned> rhs) noexcept
    -> gray_code<Unsigned>
{
    // Make sure lhs.value has the highest set bit
    if (lhs.value < rhs.value)
    {
        std::swap(lhs, rhs);
    }

    // parity of lhs and rhs
    bool lhs_p = is_odd(lhs);
    bool rhs_p = is_odd(rhs);

    gray_code<Unsigned> res = lhs ^ rhs;
    Unsigned i{};
    while (rhs)
    {
        Unsigned res_i = lhs_p & rhs_p;
        res ^= res_i << i;

        bool lhs_i = lhs.value & 1u;
        bool rhs_i = rhs.value & 1u;
        lhs_p = (lhs_p & not res_i) ^ lhs_i;
        rhs_p = (rhs_p & not res_i) ^ rhs_i;

        ++i;
        lhs >>= 1u;
        rhs >>= 1u;
    }

    // We know that rhs is 0 now, let's get rid of it
    while (lhs)
    {
        Unsigned res_i = lhs_p & rhs_p;
        res ^= res_i << i;

        bool lhs_i = lhs.value & 1u;
        lhs_p = (lhs_p & not res_i) ^ lhs_i;
        rhs_p = rhs_p & not res_i;

        ++i;
        lhs >>= 1u;
    }

    // Last value in case lhs_p and rhs_p are 1
    res ^= (lhs_p & rhs_p) << i;
    return res;
}

Rethink the logic

The second loop can be rethought as "what happens when res_i is \$0\$ and what happens when it is \$1\$"? In terms of code, with some simplifications applied, it yields the following loop:

while (lhs)
{
    Unsigned res_i = lhs_p & rhs_p;
    bool lhs_i = lhs.value & 1u;
    if (res_i)
    {
        res ^= res_i << i;
        lhs_p = lhs_i;
        rhs_p = false;
    }
    else
    {
        lhs_p ^= lhs_i;
    }

    ++i;
    lhs >>= 1u;
}

While this algorithm is actually generally slower than the previous one because of branching, it permits to realize some interesting things: one of the branches always sets rhs_p to false and cannot be executed when rhs_p == false. The other branch does not update rhs_p which means that the loop can be uderstood as "do stuff while rhs_p == true", so we can rewrite the original second loop (not the new transformed second loop) as follows:

while (rhs_p)
{
    res ^= lhs_p << i;
    rhs_p = not lhs_p;
    lhs_p = lhs.value & 1u;

    ++i;
    lhs >>= 1u;
}

Many simplifications have been applied since we know during this loop that rhs_p == true. We end up with a second loop in the algorithm almost free compared to the first one, which means that when we add a number with a big highest set bit to one with a small highest set bit, it could perform significantly faster than the previous versions. Also note that since the last loop ends when rhs_p == false, then rhs_p & lhs_p is always false is the post-last-loop operation, so we can simply get rid of it altogether.

Now, we can observe that the loop updates res with lhs_p until rhs_p is false. However, rhs_p is set to not lhs_p which means that res is updated with lhs_p while lhs_p is \$0\$ except for the last iterations. In other words, res is meaningfully updated only once since when we reach lhs_p == 1, we leave the loop right after the next iteration. So we can tweak the loop to kick the res update out of it:

if (rhs_p)
{
    while (not lhs_p)
    {
        lhs_p = lhs.value & 1u;
        ++i;
        lhs >>= 1u;
    }
    res ^= lhs_p << i;
}

Combined one after the other, these optimizations leave us with the following algorithm:

template<typename Unsigned>
auto operator+(gray_code<Unsigned> lhs, gray_code<Unsigned> rhs) noexcept
    -> gray_code<Unsigned>
{
    if (lhs.value < rhs.value)
    {
        std::swap(lhs, rhs);
    }

    bool lhs_p = is_odd(lhs);
    bool rhs_p = is_odd(rhs);

    gray_code<Unsigned> res = lhs ^ rhs;
    Unsigned i{};
    // Algorithm until the smallest number is zero
    while (rhs)
    {
        Unsigned res_i = lhs_p & rhs_p;
        res ^= res_i << i;

        bool lhs_i = lhs.value & 1u;
        bool rhs_i = rhs.value & 1u;
        lhs_p = (lhs_p & not res_i) ^ lhs_i;
        rhs_p = (rhs_p & not res_i) ^ rhs_i;

        ++i;
        lhs >>= 1u;
        rhs >>= 1u;
    }
    // Algorithm until the largest number is zero
    if (rhs_p)
    {
        while (not lhs_p)
        {
            lhs_p = lhs.value & 1u;
            ++i;
            lhs >>= 1u;
        }
        res ^= lhs_p << i;
    }
    return res;
}

While this is the most versatile and optimized version of the algorithm since the beginning of this answer, it is still an order of magnitude slower than converting the Gray codes to regular integers, performing a regular addition and converting the result back to a Gray code. Nevertheless, it was fun to see how much this solution could be pushed and optimized. I hope that it may give ideas to people interested into adding two Gray codes in the future :)

Answer 2

All my comments are very minor:

Dry Code:

Your constructor/assignment operator are not DRY.

constexpr explicit gray_code(value_type value) noexcept:
    value( (value >> 1) ^ value )
{}

auto operator=(value_type other) & noexcept
    -> gray_code&
{
    value = (other >> 1) ^ other;
    return *this;
}

If things change you have to modify the code in multiple places. You should wrap the work in a function. That both these methods use:

constexpr explicit gray_code(value_type value) noexcept:
    value( convertValueToGray(value) )
{}

auto operator=(value_type other) & noexcept
    -> gray_code&
{
    value = convertValueToGray(other);
    return *this;
}

New Return Type:

Not sure I like the new style of return:
I personal don't think it is as clear (but that may be me and I need to get used to it). Personally I currently only use the new style when the return type is being derived.

auto operator=(value_type other) & noexcept
    -> gray_code&

When not deriving the return type I still prefer the old style.

gray_code& operator=(value_type other) & noexcept

But this is not a big deal. The project/company style guide will dictate what is prefered.

Public Member Variables:

Your value is public:
Not sure if this is on purpose? You have assignment that does conversion on the way in and a conversion operator to adjust the value on the way out. Do you really want to allow access to the raw value for manipulation without supervision?

struct gray_code
{
    value_type value;

Costly Conversion

The conversion back to value_type seems a bit heavy:

operator value_type() const noexcept
{
    value_type res = value;
    for (value_type mask = std::numeric_limits<value_type>::digits / 2
         ; mask ; mask >>= 1)
    {
        res ^= res >> mask;
    }
    return res;
}

If you are doing this more than a couple of times that may become expensive. Why not pre-prepare it and store it in a mutable member.

// Variable containing the gray code
value_type          value;
mutable value_type  convertedValue; // Created on construction/assignment.
                                    // Marked as mutable to indicate that
                                    // it is not part of the state of the
                                    // object but rather a cached value.

operator value_type() const noexcept  {return convertedValue;}

Personally I think your code is a little heavy on useless comments.

Answer 3

Some operations can be shaved from the main loop
(there are machines without "AndNot-instruction"):

// process bits from least significant to most in shorter operand
while (0 < rhs) {
    Unsigned res_i = lhs_p & rhs_p;
    res ^= res_i << i;

    lhs_p ^= res_i ^ lhs;
    rhs_p = (rhs_p ^ res_i ^ rhs) & 1;

    ++i;
    lhs >>= 1u;
    rhs >>= 1u;
}
// remaining bits of longer operand
if (0 == rhs_p)
    return res;
lhs_p &= 1;

(While expected run lengths for uniformly distributed bit patterns are (1 +) smaller than 1,) There is no need to process the remaining bits in a loop:

    if (0 != lhs_p)
        return res ^= lhs_p << i;
    Unsigned bit = lhs ^ (lhs & (lhs-1));
    return res ^ (bit << (i+1));

I tried to shift the parities left instead of shifting the operands right: unsurprisingly, the differences are negligible.

    Unsigned bit = 1;
    for ( ; bit <= rhs && bit != 0 ; bit <<= 1) {
        long both = lhs_p&rhs_p;
        res ^= both;
        lhs_p = (lhs_p ^ both ^ lhs) << 1;
        rhs_p = ((rhs_p ^ both ^ rhs) & bit) << 1;  
    }
    if (0 == rhs_p)
        return res;
    lhs_p &= bit;
    if (0 != lhs_p)
        return res ^= lhs_p;
    bit = lhs ^ (lhs & (lhs-bit));
    return res ^ (bit << 1);

With determination, the bit-hack variants of population count and parity can be coded without undue regard to type sizes:

/** bit masks */
const uintmax_t
    ALL_ONES = ~(uintmax_t)0,
    O_ONES = ALL_ONES/0xff, //0x...10101
    FIVES  = ALL_ONES/3,    //0x...33333
    THREES = ALL_ONES/5,    //0x...55555
    OFOFS  = ALL_ONES/17;   //0x...F0F0F

/** return number of bits set to one in <code>bits</code> */
// fast multiplication variant
int popCount(uintmax_t bits) {
    bits -= (bits >> 1) & FIVES;
    bits = (bits & THREES) + ((bits >> 2) & THREES);
// of course, O7O7s would do - but 255 doesn't divide by 7
    bits = bits + (bits >> 4) & OFOFS;
// needs one more "fold operation" and "ones spaced 16"
//  for 256 <= #bits < 65536
    return (bits * O_ONES) >> (sizeof bits * 8);
}

/** return (odd) parity of <code>bits</code> */
int parity(uintmax_t bits) {
    uintmax_t shifted;
    for (int shift = 4 ; 0 < ((shifted = bits) >>= shift) ; shift *= 2)
        bits ^= shifted;
    return (0x6996 >> ((int)bits & 0xf)) & 1;
}