Computing high 64 bits of a 64x64 int product in C

Asked 2009-10-09T01:40:47.383

I would like my C function to efficiently compute the high 64 bits of the product of two 64 bit signed ints. I know how to do this in x86-64 assembly, with imulq and pulling the result out of %rdx. But I'm at a loss for how to write this in C at all, let alone coax the compiler to do it efficiently.

Does anyone have any suggestions for writing this in C? This is performance sensitive, so "manual methods" (like Russian Peasant, or bignum libraries) are out.

This dorky inline assembly function I wrote works and is roughly the codegen I'm after:

static long mull_hi(long inp1, long inp2) {
    long output = -1;
    __asm__("movq %[inp1], %%rax;"
            "imulq %[inp2];"
            "movq %%rdx, %[output];"
            : [output] "=r" (output)
            : [inp1] "r" (inp1), [inp2] "r" (inp2)
            :"%rax", "%rdx");
    return output;
}

c 64-bit math

Edit Report

2 Answers

Since you did a pretty good job solving your own problem with the machine code, I figured you deserved some help with the portable version. I would leave an ifdef in where you do just use the assembly if in gnu on x86.

Anyway, here is an implementation based on my general answer. I'm pretty sure this is correct, but no guarantees, I just banged this out last night. You probably should get rid of the statics positive_result[] and result_negative - those are just artefacts of my unit test.

#include <stdlib.h>
#include <stdio.h>

// stdarg.h doesn't help much here because we need to call llabs()

typedef unsigned long long uint64_t;
typedef   signed long long  int64_t;

#define B32 0xffffffffUL

static uint64_t positive_result[2]; // used for testing
static int result_negative;         // used for testing

static void mixed(uint64_t *result, uint64_t innerTerm)
{
  // the high part of innerTerm is actually the easy part

    result[1] += innerTerm >> 32;

  // the low order a*d might carry out of the low order result

    uint64_t was = result[0];

    result[0] += (innerTerm & B32) << 32;

    if (result[0] < was) // carry!
      ++result[1];
}


static uint64_t negate(uint64_t *result)
{
  uint64_t t = result[0] = ~result[0];
  result[1] = ~result[1];
  if (++result[0] < t)
    ++result[1];
  return result[1];
}

uint64_t higherMul(int64_t sx, int64_t sy)
{
    uint64_t x, y, result[2] = { 0 }, a, b, c, d;

    x = (uint64_t)llabs(sx);
    y = (uint64_t)llabs(sy);

    a = x >> 32;
    b = x & B32;
    c = y >> 32;
    d = y & B32;

  // the highest and lowest order terms are easy

    result[1] = a * c;
    result[0] = b * d;

  // now have the mixed terms ad + bc to worry about

    mixed(result, a * d);
    mixed(result, b * c);

  // now deal with the sign

    positive_result[0] = result[0];
    positive_result[1] = r


                        
                        
                            answered 2009-10-09T21:34:11.513


            
            
                
                    
                        
                            
                            
                        
                        1
                        
                            
                            
                        
                        
                    
                    
                        
                        Wait, you have a perfectly good, optimized assembly solution already
working for this, and you want to back it out and try to write it in
an environment that doesn't support 128 bit math?  I'm not following.

As you're obviously aware, this operation is a single instruction on
x86-64.  Obviously nothing you do is going to make it work any better.
If you really want portable C, you'll need to do something like
DigitalRoss's code above and hope that your optimizer figures out what
you're doing.

If you need architecture portability but are willing to limit yourself
to gcc platforms, there are __int128_t (and __uint128_t) types in the
compiler intrinsics that will do what you want.

                        
                        
                            answered 2009-10-09T02:32:43.050



        
        
            Your Answer