Stern's diatomic series fusc: replacing #define functions with two return values

Question

I have a function to review which implements Stern's diatomic series fusc for a single word (ulong = unsigned long) in place, as a part of an arbitrary-precision implementation. In fact, there's just one part of it that's bothering me, though I'd be happy to get feedback on any of it, of course.

I have a (pseudofunction) define fusc8bits which 'returns' two values. It works, but I'd rather do this some other way. As this is a performance-critical function any replacements will need to keep this requirement in mind. Alternatively, and unexpectedly, a reviewer might say that my approach is actually a good one.

Unfortunately the arrays are used by other functions, so they can't be pulled inside the function. (I wish they weren't so visually imposing, though.)

This can be found on GitHub here.

static const ulong fuscAA[] = {
  1,  8,  7,  13, 6,  17, 11, 16, 5,  19, 14, 23, 9,  22, 13, 17, 4,  19, 15,
  26, 11, 29, 18, 25, 7,  24, 17, 27, 10, 23, 13, 16, 3,  17, 14, 25, 11, 30,
  19, 27, 8,  29, 21, 34, 13, 31, 18, 23, 5,  22, 17, 29, 12, 31, 19, 26, 7,
  23, 16, 25, 9,  20, 11, 13, 2,  13, 11, 20, 9,  25, 16, 23, 7,  26, 19, 31,
  12, 29, 17, 22, 5,  23, 18, 31, 13, 34, 21, 29, 8,  27, 19, 30, 11, 25, 14,
  17, 3,  16, 13, 23, 10, 27, 17, 24, 7,  25, 18, 29, 11, 26, 15, 19, 4,  17,
  13, 22, 9,  23, 14, 19, 5,  16, 11, 17, 6,  13, 7,  8,  1,  7,  6,  11, 5,
  14, 9,  13, 4,  15, 11, 18, 7,  17, 10, 13, 3,  14, 11, 19, 8,  21, 13, 18,
  5,  17, 12, 19, 7,  16, 9,  11, 2,  11, 9,  16, 7,  19, 12, 17, 5,  18, 13,
  21, 8,  19, 11, 14, 3,  13, 10, 17, 7,  18, 11, 15, 4,  13, 9,  14, 5,  11,
  6,  7,  1,  6,  5,  9,  4,  11, 7,  10, 3,  11, 8,  13, 5,  12, 7,  9,  2,
  9,  7,  12, 5,  13, 8,  11, 3,  10, 7,  11, 4,  9,  5,  6,  1,  5,  4,  7,
  3,  8,  5,  7,  2,  7,  5,  8,  3,  7,  4,  5,  1,  4,  3,  5,  2,  5,  3,
  4,  1,  3,  2,  3,  1,  2,  1,  1
};
static const ulong fuscAB[] = {
  8,  7,  13, 6,  17, 11, 16, 5,  19, 14, 23, 9,  22, 13, 17, 4,  19, 15, 26,
  11, 29, 18, 25, 7,  24, 17, 27, 10, 23, 13, 16, 3,  17, 14, 25, 11, 30, 19,
  27, 8,  29, 21, 34, 13, 31, 18, 23, 5,  22, 17, 29, 12, 31, 19, 26, 7,  23,
  16, 25, 9,  20, 11, 13, 2,  13, 11, 20, 9,  25, 16, 23, 7,  26, 19, 31, 12,
  29, 17, 22, 5,  23, 18, 31, 13, 34, 21, 29, 8,  27, 19, 30, 11, 25, 14, 17,
  3,  16, 13, 23, 10, 27, 17, 24, 7,  25, 18, 29, 11, 26, 15, 19, 4,  17, 13,
  22, 9,  23, 14, 19, 5,  16, 11, 17, 6,  13, 7,  8,  1,  7,  6,  11, 5,  14,
  9,  13, 4,  15, 11, 18, 7,  17, 10, 13, 3,  14, 11, 19, 8,  21, 13, 18, 5,
  17, 12, 19, 7,  16, 9,  11, 2,  11, 9,  16, 7,  19, 12, 17, 5,  18, 13, 21,
  8,  19, 11, 14, 3,  13, 10, 17, 7,  18, 11, 15, 4,  13, 9,  14, 5,  11, 6,
  7,  1,  6,  5,  9,  4,  11, 7,  10, 3,  11, 8,  13, 5,  12, 7,  9,  2,  9,
  7,  12, 5,  13, 8,  11, 3,  10, 7,  11, 4,  9,  5,  6,  1,  5,  4,  7,  3,
  8,  5,  7,  2,  7,  5,  8,  3,  7,  4,  5,  1,  4,  3,  5,  2,  5,  3,  4,
  1,  3,  2,  3,  1,  2,  1,  1,  0
};
static const ulong fuscBA[] = {
  0,  1,  1,  2,  1,  3,  2,  3,  1,  4,  3,  5,  2,  5,  3,  4,  1,  5,  4,
  7,  3,  8,  5,  7,  2,  7,  5,  8,  3,  7,  4,  5,  1,  6,  5,  9,  4,  11,
  7,  10, 3,  11, 8,  13, 5,  12, 7,  9,  2,  9,  7,  12, 5,  13, 8,  11, 3,
  10, 7,  11, 4,  9,  5,  6,  1,  7,  6,  11, 5,  14, 9,  13, 4,  15, 11, 18,
  7,  17, 10, 13, 3,  14, 11, 19, 8,  21, 13, 18, 5,  17, 12, 19, 7,  16, 9,
  11, 2,  11, 9,  16, 7,  19, 12, 17, 5,  18, 13, 21, 8,  19, 11, 14, 3,  13,
  10, 17, 7,  18, 11, 15, 4,  13, 9,  14, 5,  11, 6,  7,  1,  8,  7,  13, 6,
  17, 11, 16, 5,  19, 14, 23, 9,  22, 13, 17, 4,  19, 15, 26, 11, 29, 18, 25,
  7,  24, 17, 27, 10, 23, 13, 16, 3,  17, 14, 25, 11, 30, 19, 27, 8,  29, 21,
  34, 13, 31, 18, 23, 5,  22, 17, 29, 12, 31, 19, 26, 7,  23, 16, 25, 9,  20,
  11, 13, 2,  13, 11, 20, 9,  25, 16, 23, 7,  26, 19, 31, 12, 29, 17, 22, 5,
  23, 18, 31, 13, 34, 21, 29, 8,  27, 19, 30, 11, 25, 14, 17, 3,  16, 13, 23,
  10, 27, 17, 24, 7,  25, 18, 29, 11, 26, 15, 19, 4,  17, 13, 22, 9,  23, 14,
  19, 5,  16, 11, 17, 6,  13, 7,  8
};
static const ulong fuscBB[] = {
  1,  1,  2,  1,  3,  2,  3,  1,  4,  3,  5,  2,  5,  3,  4,  1,  5,  4,  7,
  3,  8,  5,  7,  2,  7,  5,  8,  3,  7,  4,  5,  1,  6,  5,  9,  4,  11, 7,
  10, 3,  11, 8,  13, 5,  12, 7,  9,  2,  9,  7,  12, 5,  13, 8,  11, 3,  10,
  7,  11, 4,  9,  5,  6,  1,  7,  6,  11, 5,  14, 9,  13, 4,  15, 11, 18, 7,
  17, 10, 13, 3,  14, 11, 19, 8,  21, 13, 18, 5,  17, 12, 19, 7,  16, 9,  11,
  2,  11, 9,  16, 7,  19, 12, 17, 5,  18, 13, 21, 8,  19, 11, 14, 3,  13, 10,
  17, 7,  18, 11, 15, 4,  13, 9,  14, 5,  11, 6,  7,  1,  8,  7,  13, 6,  17,
  11, 16, 5,  19, 14, 23, 9,  22, 13, 17, 4,  19, 15, 26, 11, 29, 18, 25, 7,
  24, 17, 27, 10, 23, 13, 16, 3,  17, 14, 25, 11, 30, 19, 27, 8,  29, 21, 34,
  13, 31, 18, 23, 5,  22, 17, 29, 12, 31, 19, 26, 7,  23, 16, 25, 9,  20, 11,
  13, 2,  13, 11, 20, 9,  25, 16, 23, 7,  26, 19, 31, 12, 29, 17, 22, 5,  23,
  18, 31, 13, 34, 21, 29, 8,  27, 19, 30, 11, 25, 14, 17, 3,  16, 13, 23, 10,
  27, 17, 24, 7,  25, 18, 29, 11, 26, 15, 19, 4,  17, 13, 22, 9,  23, 14, 19,
  5,  16, 11, 17, 6,  13, 7,  8,  1
};

#define fusc8bits(a, b, idx)                                                   \
  {                                                                            \
    int i = (idx)&0xFF;                                                        \
    int newA = a * fuscAA[i] + b * fuscAB[i];                                  \
    b = a * fuscBA[i] + b * fuscBB[i];                                         \
    a = newA;                                                                  \
  }
static void
fusc_word(ulong u, ulong* a, ulong* b)
{
  *a = fuscAA[u & 0xFF];
  *b = 0;
  fusc8bits(*a, *b, u >> 8) fusc8bits(*a, *b, u >> 16)
    fusc8bits(*a, *b, u >> 24)
#ifdef LONG_IS_64BIT
      fusc8bits(*a, *b, u >> 32) fusc8bits(*a, *b, u >> 40)
        fusc8bits(*a, *b, u >> 48) fusc8bits(*a, *b, u >> 56)
#endif
}

Answer 1

this is a performance-critical function any replacements will need to keep this requirement in mind

Use restrict

fusc_word(ulong u, ulong* a, ulong* b) uses a, b and certainly the algorithm does not work if those pointers point to overlapping data.

A complier can not make that assumption though and so must emit code as if a and b potentially point to the same location.

Use restrict to allow the compiler to assume they do not overlap and potentially emit more efficient code.

// static void fusc_word(ulong u, ulong* a, ulong* b)
static void fusc_word(ulong u, ulong* restrict a, ulong* restrict b)

Smaller table

Depending on your platform, a smaller table may be faster - profile to find out.

// static const ulong fuscAA[] 
static const unsigned char fuscAA[] 
// or 
static const unsigned fuscAA[] 

Bug

int newA = a * fuscAA[i] + b * fuscAB[i]; may truncate. I'd expect newA to be ulong.

---

To return 2 int, the C library already has div_t, so code could use/mimic that.

Note: returning a struct, even a small one, tends to rankle some. For me I find returning a small struct fine for numeric coding such a making a struct fraction { int n, d; } a useful approach. The trick is that there are those who will employ ever larger and larger structs making very inefficient code. Use sane judgment.