Implementing the Extended Euclidean Algorithm with Unsigned Inputs

As the previous post showed, it’s possible to correctly implement the Extended Euclidean Algorithm using one signed integral type for all input parameters, intermediate variables, and output variables. None of the calculations will overflow. The implementation was given as follows:

template <class T>
void extended_euclidean(const T a, const T b, T* pGcd, T* pX, T* pY)
{
   static_assert(std::numeric_limits<T>::is_integer, "");
   static_assert(std::numeric_limits<T>::is_signed, "");
   assert(a >= 0 && b >= 0);   // precondition
   T x0=1, y0=0, a0=a;
   T x1=0, y1=1, a1=b;

   while (a1 != 0) {
      T q = a0/a1;
      T a2 = a0 - q*a1;
      T x2 = x0 - q*x1;
      T y2 = y0 - q*y1;
      x0=x1; y0=y1; a0=a1;
      x1=x2; y1=y2; a1=a2;
   }
   *pX = x0;
   *pY = y0;
   *pGcd = a0;
}

Figure 1: C++ Extended Euclidean Function, Signed Integer Inputs.

The precondition a >= 0 && b >= 0 is required because the proof from the last post didn’t consider any cases where a or b was negative. I haven’t investigated whether the Extended Euclidean algorithm might be valid for any and all cases of negative inputs, so I won’t speculate about what might be possible for those cases. Instead let’s focus on the proven case a >= 0 && b >= 0.

It seems a bit wasteful to use signed integers for inputs a and b when the function precondition disallows a and b from ever being negative. Why not just use unsigned integer inputs that have the same precision as the signed integer inputs would have had? The one potential problem with doing this is that when we look at the bounds, it’s not quite as obvious as the all-signed case that overflow can’t occur. Displaying the bounds proven in the last post, throughout every loop iteration in Figure 1 we know the following:
0 <= q <= max(a,b)
0 <= (q*a1) <= max(a,b)
0 <= a2 < max(a,b)
abs(q*x1) <= b
abs(x2) <= b
abs(q*y1) <= max(1,a)
abs(y2) <= max(1,a)

Assuming we go ahead and change the type for inputs a and b to be unsigned, it should be immediately clear from the bounds and Figure 1 that in order to avoid overflow, q, a2, a1, a0, and *pGcd will all need to be changed to use the same unsigned integer type as the input parameters.

More problematic are x2 and y2 (and likewise x1, y1, x0, y0, *pX, and *pY). On the one hand, these variables need to remain signed types since the bounds above don’t restrict them from taking on negative values. But on the other hand, the bounds above allow them to be as large as a or b, and since a and b are now unsigned types, the largest possible value for a or b would overflow the corresponding signed type. Fortunately we can resolve this conflict, by using less generally applicable but tighter bounds for these variables. Let’s review the second set of proven bounds from the last post. For every loop iteration other than the final iteration (when a2 == 0), we know the following tighter bounds hold:
0 <= q <= max(a,b)/2
abs(q*x1) <= max(1,b/2)
abs(x2) <= max(1,b/2)
abs(q*y1) <= max(1,a/2)
abs(y2) <= max(1,a/2)

Superficially, it seems like these bounds should be enough to guarantee we can use signed types for x2 and y2 without any problems. The only loop iteration where they could overflow is the final iteration, and their values don’t affect anything on the last iteration, so seemingly overflow (if it happened) would be fine at that point. But in truth it’s unacceptable, because overflow of a signed integral type is undefined behavior in C and C++. Upon overflow on the last iteration, a reasonable possibly is everything would be okay, but another reasonable possibility is a trap would occur that crashes the program. Anything is allowed after invoking undefined behavior, so in the end it’s hopeless to try to reason about whether the program would be okay. The only safe thing to do is to view any undefined behavior as a fatal bug that must be fixed.

Since we can not allow undefined behavior, we must ensure signed integer overflow is completely impossible in the function. Fortunately this is easy. We can simply test if we are in the final loop iteration (by checking if a1 == 0), and break out if so. Most of the work in the final iteration is irrelevant since it’s never used later. The only exception is the reassignments (x0=x1; y0=y1; a0=a1), which still need to be done before causing a break.

Incorporating all the discussed changes, we get the following potential implementation:

template <class S, class U>
void extended_euclidean(const U a, const U b, U* pGcd, S* pX, S* pY)
{
   static_assert(std::numeric_limits<S>::is_integer, "");
   static_assert(std::numeric_limits<S>::is_signed, "");
   static_assert(std::numeric_limits<U>::is_integer, "");
   static_assert(!(std::numeric_limits<U>::is_signed), "");
   static_assert(std::is_same<std::make_signed<U>::type, S>::value, "");
   S x0=1, y0=0;
   U a0=a;
   S x1=0, y1=1;
   U a1=b;

   while (a1 != 0) {
      U q = a0/a1;
      U a2 = a0 - q*a1;
      if (a2 == 0) {
         x0=x1; y0=y1; a0=a1;
         break;
      }
      S x2 = x0 - q*x1;
      S y2 = y0 - q*y1;
      x0=x1; y0=y1; a0=a1;
      x1=x2; y1=y2; a1=a2;
   }
   *pX = x0;
   *pY = y0;
   *pGcd = a0;
}

Figure 2: Potential (But Broken) Implementation, Unsigned Integer Inputs.

But there’s still a problem. We determined that q must be unsigned and x1 must be signed, and Figure 2 reflects this. Therefore during the calculation in Figure 2 of (q*x1), the compiler will implicitly convert x1 to the unsigned type of q in order to carry out the calculation. Quoting from the C++17 standard section 8 [expr] (11.5.3):
“Many binary operators that expect operands of arithmetic type cause conversions… The following rules shall be applied… if the operand that has unsigned integer type has rank equal to the rank of the type of the other operand, the operand with signed integer type shall be converted to the type of the operand with unsigned integer type.”
[The C standard specifies very similar behavior.] The problem is we didn’t want to convert x1 to unsigned. And even putting that aside, the conversion causes further conversions which can eventually result in implementation defined behavior. What we really wanted and intended is to convert q from unsigned to signed, which we can do with an explicit cast.

Explicitly casting q to signed integral type works well in the code of Figure 2. We know from the tight bounds given above that on every loop iteration except the final iteration,
0 <= q <= max(a,b)/2
Figure 2 has a break on the last loop iteration before it gets to the calculations (q*x1 and q*y1) which would require type conversion of q to signed integral. So the conversions will never be performed on the final loop iteration, and therefore no conversion of q will ever overflow the signed result type.

Incorporating the explicit casts of q, we get the following implementation:

template <class S, class U>
void extended_euclidean(const U a, const U b, U* pGcd, S* pX, S* pY)
{
   static_assert(std::numeric_limits<S>::is_integer, "");
   static_assert(std::numeric_limits<S>::is_signed, "");
   static_assert(std::numeric_limits<U>::is_integer, "");
   static_assert(!(std::numeric_limits<U>::is_signed), "");
   static_assert(std::is_same<std::make_signed<U>::type, S>::value, "");
   S x0=1, y0=0;
   U a0=a;
   S x1=0, y1=1;
   U a1=b;

   while (a1 != 0) {
      U q = a0/a1;
      U a2 = a0 - q*a1;
      if (a2 == 0) {
         x0=x1; y0=y1; a0=a1;
         break;
      }
      S x2 = x0 - static_cast<S>(q)*x1;
      S y2 = y0 - static_cast<S>(q)*y1;
      x0=x1; y0=y1; a0=a1;
      x1=x2; y1=y2; a1=a2;
   }
   *pX = x0;
   *pY = y0;
   *pGcd = a0;
}

Figure 3: Correct Implementation, Unsigned Integer Inputs.

This implementation in Figure 3 should be correct, but the conditional branch inside the loop isn’t ideal from a performance standpoint. It would be nice if we could remove the conditional branch. Indeed we can, by re-arranging the code.

Optimized Implementation

We’ll remove the conditional branch inside the loop of Figure 3 by rotating statements from the bottom of the loop to the top of the loop. Each rotation will require changing the code that precedes and follows the loop (and sometimes the loop condition too), so that the meaning of the function stays the same. Aside from the rotation, nothing inside the loop should change. The end goal of the rotations is to get the position of the nested conditional branch to be at the end of the loop, so that we can merge it into the loop condition.

Let’s first rotate the reassignments from the end of the loop in Figure 3 to the top of the loop. [If this description is hard to follow, it may help to compare Figure 3 with Figure 4.] The rotation affects the code after the loop, since x0, y0, and a0 won’t be reassigned upon the end of the final loop iteration. So to compensate, *pX, *pY, and *pGcd will now need to be set to the values of x1, y1, and a1. This change to *pX, *pY, and *pGcd will need to also be valid if the loop is never taken at all, so we’ll need to change the initial values of x1, y1, and a1 to the initial values in Figure 3 given to x0, y0, and a0. Fortunately, since we will have already done the rotation that places the loop reassignments at the top of the loop, this last change also ensures the correct values will get reassigned upon the first entry to the loop. At least in part. We haven’t discussed a second part of the effect of having the reassignments at the top of the loop: the reassignments use x2, y2, and a2, which are variables that haven’t been declared or initialized yet. We’ll fix this by adding declarations and initializations to the start of the function for x2, y2, and a2, setting them to the initial values that had formerly been given to x1, y1, and a1 in Figure 3. As a result, upon first entry of the loop, the reassignments at the top of the loop will re-set the correct initial values, leaving the meaning of the function unchanged.

There’s one last detail. We’ll need to change the loop condition to refer to a2 instead of a1, again due to the rotation of the reassignment statements. Since we’ll have already declared and initialized a2 at the top of the function (to the initial value formerly used in Figure 3 by a1), this preserve the meaning of the function perfectly.

It’s easiest to see all this in actual code. By rotating the loop’s reassignments in Figure 3, we’ll have changed the implementation to the following code:

template <class S, class U>
void extended_euclidean(const U a, const U b, U* pGcd, S* pX, S* pY)
{
   static_assert(std::numeric_limits<S>::is_integer, "");
   static_assert(std::numeric_limits<S>::is_signed, "");
   static_assert(std::numeric_limits<U>::is_integer, "");
   static_assert(!(std::numeric_limits<U>::is_signed), "");
   static_assert(std::is_same<std::make_signed<U>::type, S>::value, "");
   S x1=1, y1=0;
   U a1=a;
   S x2=0, y2=1;
   U a2=b;

   while (a2 != 0) {
      S x0=x1;
      S y0=y1;
      U a0=a1;
      x1=x2; y1=y2; a1=a2;

      U q = a0/a1;
      a2 = a0 - q*a1;
      if (a2 == 0)
         break;
      x2 = x0 - static_cast<S>(q)*x1;
      y2 = y0 - static_cast<S>(q)*y1;
   }
   *pX = x1;
   *pY = y1;
   *pGcd = a1;
}

Figure 4: Optimization Step 1, Unsigned Integer Inputs.

We haven’t yet removed the conditional branch inside the loop, but we’re closer to the goal. We still need to rotate the calculations of x2 and y2 from the bottom of the loop to the top, and again fix the code before and after the loop to preserve the meaning of the function. As it turns out, there’s no need to change any code after the function. But now that the calculations are at the top of the loop, the variable q gets used before it has been declared or initialized. In order to preserve the meaning of the function, we’ll need to add a declaration and initialization of q at the start of the function. Initializing q to zero works quite well. On the first entry of the loop, the calculations of x2 and y2 will thereby contain a multiplication by zero (since q == 0), effectively reducing those calculations on the first iteration to reassignments from x0 and y0. This means we’ll need to add to the start of the function declarations and initializations for x0 and y0, setting them to the initial values that Figure 4 had used for x2 and y2.

After making these changes to Figure 4, we have the following code:

template <class S, class U>
void extended_euclidean(const U a, const U b, U* pGcd, S* pX, S* pY)
{
   static_assert(std::numeric_limits<S>::is_integer, "");
   static_assert(std::numeric_limits<S>::is_signed, "");
   static_assert(std::numeric_limits<U>::is_integer, "");
   static_assert(!(std::numeric_limits<U>::is_signed), "");
   static_assert(std::is_same<std::make_signed<U>::type, S>::value, "");
   S x1=1, y1=0;
   U a1=a;
   S x0=0, y0=1;
   U a2=b;
   U q=0;

   while (a2 != 0) {
      S x2 = x0 - static_cast<S>(q)*x1;
      S y2 = y0 - static_cast<S>(q)*y1;
      x0=x1;
      y0=y1;
      U a0=a1;
      x1=x2; y1=y2; a1=a2;

      q = a0/a1;
      a2 = a0 - q*a1;
      if (a2 == 0)
         break;
   }
   *pX = x1;
   *pY = y1;
   *pGcd = a1;
}

Figure 5: Optimization Step 2, Unsigned Integer Inputs.

We can now get rid of the conditional branch at the end of the function, since it’s redundant. Making this change and cleaning up the code a little, we get the following optimized function for the Extended Euclidean Algorithm with unsigned integer inputs:

template <class S, class U>
void extended_euclidean(const U a, const U b, U* pGcd, S* pX, S* pY)
{
   static_assert(std::numeric_limits<S>::is_integer, "");
   static_assert(std::numeric_limits<S>::is_signed, "");
   static_assert(std::numeric_limits<U>::is_integer, "");
   static_assert(!(std::numeric_limits<U>::is_signed), "");
   static_assert(std::is_same<std::make_signed<U>::type, S>::value, "");
   S x1=1, y1=0;
   U a1=a;
   S x0=0, y0=1;
   U a2=b, q=0;

   while (a2 != 0) {
      S x2 = x0 - static_cast<S>(q)*x1;
      S y2 = y0 - static_cast<S>(q)*y1;
      x0=x1; y0=y1;
      U a0=a1;
      x1=x2; y1=y2; a1=a2;

      q = a0/a1;
      a2 = a0 - q*a1;
   }
   *pX = x1;
   *pY = y1;
   *pGcd = a1;
}

Figure 6: Optimized Implementation, Unsigned Integer Inputs.

Conclusion

The Extended Euclidean algorithm can be correctly and very efficiently implemented using unsigned integer inputs and signed integer outputs, so long as all the function parameters have the same precision. Figure 6 provides the optimized implementation. There is a header file with the complete implementation, and a CMake file that will create a test project for the implementation if you wish to run it.

This entry was posted in Uncategorized. Bookmark the permalink.

7 Responses to Implementing the Extended Euclidean Algorithm with Unsigned Inputs

  1. ecirhana says:

    Hello!
    I asked a question about modular inverses and unsigned integers on StackOverflow [1] and was pointed to your article.
    Please, if I wanted to calculate the modular inverse, i.e. unsigned pX, would it be possible to get rid of `static_cast<S>` in `S x2 = x0 – static_cast<S>(q)*x1;` ?
    In other words, for modular inverses, would it be possible to use only `U`?
    And if that was possible, how would the final `if (x1 < 0) x1 += b0;` for correcting the result look like?
    Thank you!
    [1] https://stackoverflow.com/questions/53560302/modular-inverses-and-unsigned-integers

    Like

    • hurchalla says:

      Hi Ecir,
      The simplest way to implement the modular multiplicative inverse using only unsigned inputs and outputs would be to call the unsigned extended euclidean function from this post – you can use either the code of Figure 6 (the optimized version) or Figure 3 (the unoptimized but more straightforward version). You shouldn’t need to change that function code at all. More specifically, answering your first question, I don’t think you need to make any changes to the unsigned extended euclidean function that would get rid of the internal signed variables or internal signed casts. (Getting rid of the signed variables and casts probably isn’t even possible).
      After calling the unsigned input extended euclidean function, in order to get a result that is >= 0, you’ll need to add the modulus to the result if the result would otherwise be less than 0. The addition won’t have any problems with overflow. Without going into detail on the C/C++ standards, the addition will work like you want even though it has signed and unsigned addends.
      I haven’t compiled or tested the following code yet, but I expect it is fine and it should provide you an easy implementation of the modular multiplicative inverse, using all unsigned inputs and outputs.

      I hope this helps. Thanks for your comment, If I have some time in the next few days, I’ll write a blog post here with details on calculating the modular multiplicative inverse and why it works.

      // Returns false if the inverse doesn't exist.  Returns true if the inverse exists, and places
      // the result in *pInverseValue.
      template <class U>
      bool modular_multiplicative_inverse(const U modulus, const U value, U* pInverseValue)
      {
          if (modulus == 0)  // ordinarily operations modulo 0 are undefined
              return false;
          if (modulus == 1) {
              // without this "if" clause, when modulus == 1 and value == 1,
              // this function would calculate the result to be 1.  That result
              // wouldn't be completely wrong, but it isn't reduced.  We always
              // want a fully reduced result.  When modulus == 1, the fully
              // reduced result will always be 0.
              *pInverseValue = 0;
              return true;
          }
          typedef std::make_signed<U>::type S;
          U gcd;
          S x, y;
          // the code from Figure 6 or 3 could implement the next function
          extended_euclidean_unsigned_inputs<S,U>(modulus, value, &gcd, &x, &y);
          if (gcd != 1)
              return false;
          else {
              if (y >= 0)
                  *pInverseValue == static_cast<U>(y);
              else
                  *pInverseValue == static_cast<U>(y) + modulus;
              return true;
          }
      }
      

      Like

      • ecirhana says:

        Thank you very much for you answer! Please, one more question: wouldn’t it work if I was using only unsigned integers in the main loop and only at the end I would do something like `if (static_cast(y1) < 0) y1 += modulus`? And thank you for the article, it was very informative!

        Like

      • ecirhana says:

        Thank you very much for you answer! Please, one more question: wouldn’t it work if I was using only unsigned integers in the main loop and only at the end I would do something like `if (static_cast(y1) < 0) y1 += modulus`? And thank you for the article, it was very informative!

        (Sorry, the previous reply probably got formatted badly.)

        Like

      • ecirhana says:

        Thank you very much for you answer! Please, one more question: wouldn’t it work if I was using only unsigned integers in the main loop and only at the end I would do something like `if (static_cast<S>(y1) < 0) y1 += modulus`? And thank you for the article, it was very informative!

        (Sorry!! The comments seems to be eating the < sign!)

        Like

  2. hurchalla says:

    That’s an interesting question. For the sake of simplicity and code that is easy to understand, I still think best solution is to use the unsigned extended euclidean function, with code similar to my last reply.

    But thinking about your question with the curiosity of what is possible, I believe you can implement the modular multiplicative inverse using all unsigned types. The extended euclidean algorithm unfortunately requires some variables with signed values. But we can work-around this requirement by representing signed values within an unsigned type: we would just need to interpret those (unsigned type) variables’ bit pattern as a two’s compliment signed value. Throughout all the arithmetic operations, there would be no changes needed despite using these unsigned types. But when we initialize those variables we would need to do a manual conversion step into this representation, and we would need to do a manual conversion step out of this representation when we want to convert one of those variables back to a normally-interpreted type like one of the standard unsigned or signed integral types. Converting into this representation is easy – we just do a normal cast to the unsigned type – for example, unsigned int y = static_cast<unsigned int>(x). The C/C++ standards essentially require that this conversion results in a two’s complement bit representation. Converting out of this representation is slightly harder. In our particular application of the modular multiplicative inverse, we will want to convert out of this representation and into just a normally interpreted unsigned type. So long as at the point of the conversion the variable represents a positive or zero value, we can simply change how we interpret the variable (we no longer view it as two’s complement and instead view it as a normal unsigned value) – we don’t even need to do a cast, since the source and destination types are both the same unsigned type. At the point of conversion, the variable can not be allowed to be represent a negative value though, since in that case there would be no way to convert out of our representation and into a normally-interpreted unsigned type.

    As you noted, and as I said in my earlier reply too, code for the modular multiplicative inverse should add the modulus to the result when the result is negative, and this will achieve a positive or zero value for the final result – which we would return from the function. Since we have a positive or zero value final result, if we were using the two’s complement interpretation I’ve discussed, at this point we would be safe to convert into normally-interpreted unsigned type! The only question we need to answer is how do we know if the (two’s complement representation) result variable is negative? The literal C or C++ type we would be using for that variable is unsigned, so obviously comparing that variable for less than zero wouldn’t do what we want. Unsigned variables are never less than zero, and the compiler doesn’t know anything about our unusual two’s complement interpretation of the variable. In order to know if our variable represents a negative value, we need to check if the the first bit of the variable is a 1. We can do this by shifting right by one less than the bit depth of our type, and then checking if the shifted value == 1. For example, if our type is unsigned int and the bit depth of that type is 32, then we would shift right by 31. And if the shifted value == 1, then we need to add the modulus to our (unshifted) result to get the final non-negative value.

    Here’s what this idea looks like when implemented. I tested this following code only very lightly, but nevertheless it passed the basic tests, and the idea seems to work as desired.

    
    #include <limits.h>
    
    // Returns false if the inverse doesn't exist.  Returns true if the inverse exists, and places
    // the result in *pInverseValue.
    template <class U>
    bool modular_multiplicative_inverse(const U modulus, const U value, U* pInverseValue)
    {
        static_assert(std::numeric_limits<U>::is_integer, "");
        static_assert(!(std::numeric_limits<U>::is_signed), "");
    
        if (modulus == 0)  // ordinarily operations modulo 0 are undefined
            return false;
        if (modulus == 1) {
            // without this "if" clause, when modulus == 1 and value == 1,
            // this function would calculate the result to be 1.  That result
            // wouldn't be completely wrong, but it isn't reduced.  We always
            // want a fully reduced result.  When modulus == 1, the fully
            // reduced result will always be 0.
            *pInverseValue = 0;
            return true;
        }
        U gcd;
        U inv;    // inv should be interpreted as a signed value in two's complement representation.
        {       // This entire clause is adapted from Figure 3.
            U a0=modulus, a1=value;
            // y0 and y1 should be interpreted as signed values in two's complement representation.
            U y0=0, y1=1;
            while (a1 != 0) {
                U q = a0/a1;
                U a2 = a0 - q*a1;
                // I've removed the check of (a2 == 0) followed by a break, since there is no longer
                // any possibility of undefined behavior from signed integral overflow.
                // y2 should be interpreted as a signed value in two's complement representation
                U y2 = y0 - q*y1;
                y0=y1; a0=a1;
                y1=y2; a1=a2;
            }
            inv = y0;
            gcd = a0;
        }
    
        if (gcd != 1)
            return false;
        else {
            constexpr std::size_t numbits = sizeof(U) * CHAR_BIT;
            if ((inv >> (numbits  - 1)) == 0)    // Check if the highest bit is 0 or 1
                *pInverseValue = inv;            // Highest bit == 0 indicates a non-negative value
            else                                            // Highest bit == 1 indicates a negative value
                *pInverseValue = inv + modulus;    // Adding modulus will result in a value >= 0
            return true;
        }
    }

    Like

    • ecirhana says:

      Thank you very much yet again, very interesting!

      One very last question (I hope!), instead of checking if the highest bit is set with the shift, shouldn’t temporarily casting `inv ` to signed type work as well? I mean, if unsigned int has highest bit set and I cast it into signed int it becomes negative so theoretically I could use `if (static_cast(inv) < 0)`, instead of `if ((inv >> (numbits – 1)) == 1)` or am I mistaken?
      The advantage would be saving the shift and just interpreting `pInverseValue` as signed. The only question then is whether the magnitude of `inv` is less than `numbits – 1`, and I think you proved before that indeed it is.
      But I understand that now it uses just unsigned types. Thanks once more!

      Like

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s