Algorithm optimization: Logarithm by calculation or lookup table

Question

I am trying to optimize an audio-algoritm which has to calculate two algorithms like the following in every step. Now I've read, that there is no algorithm for logarithms that runs in polynomial time. My question would be, if it would make sense to do all logarithms by a lookup table, since they are always the same, although the disadvantage of a large amount of memory access?

for(int f=1;f<11000;f++){
    for(int fk=1;fk<1000;fk++){
        int k = ceil(12 * log2f((f - 0.5) / fk));
    }
} 

I'm programming in C++.

Thanks alot!

Answer 1

If what you really need is

ceil(12 * log2(/* something */))

then there is a very simple O(1) computation which will work, using a table of only 12 values.

1. Use frexp() to split the value into exponent and mantissa. (This is just bit manipulation, so it just takes a couple of cycles.)

2. Look the mantissa up in a precomputed list of the powers of the 12th roots of 2.0 (divided by 2), which you can do with at most four comparisons.

3. The result is 12*(exponent - 1) + index of the smallest root >= mantissa.

Edited to add:

There's actually an even better solution, because the powers of the 12th root of two are reasonably evenly spread. If you divide [0.5, 1.0) (the range of the mantissa returned by frexp) into 17 evenly spaced subranges, each subrange will fall into one of two possible return values. So if you associate each subrange with an index into the vector of roots, you need only compare the target (within that range) with a single root.

It's too late for me to write coherent English, so you'll have to settle for code:

int Ceil12TimesLog2(double x) {
  int exp;
  double mantissa = std::frexp(x, &exp) * 2;
  int idx = indexes[int((mantissa - 1.0) * 17)];
  return 12 * (exp - 1) + (mantissa <= roots[idx] ? idx : idx + 1);
}

// Here are the reference tables.

double roots[12] = {
  0x1.0000000000000p+0,
  0x1.0f38f92d97963p+0,
  0x1.1f59ac3c7d6c0p+0,
  0x1.306fe0a31b715p+0,
  0x1.428a2f98d728bp+0,
  0x1.55b8108f0ec5ep+0,
  0x1.6a09e667f3bccp+0,
  0x1.7f910d768cfb0p+0,
  0x1.965fea53d6e3dp+0,
  0x1.ae89f995ad3adp+0,
  0x1.c823e074ec129p+0,
  0x1.e3437e7101344p+0
};
int indexes[17] = { 0, 1, 2, 3, 4, 5, 5, 6, 7, 8, 8, 9, 9, 10, 10, 11, 11 };

I tried this, and it brings the total computation time down from about 0.5 seconds (for log2f) to about 0.15 seconds, using the loop in the OP. The total size of the reference tables is 164 bytes.

Answer 2

Writing z for brevity and clarity instead of fk, your inner-loop computes ceil(12 * log2f((f - 0.5) / z)). Now 12 * log2f((f - 0.5) / z) = 12*log2f(f - 0.5) - 12*log2f(z). Beforehand, compute an array with 999 entries, allowing you to compute all the values with only 11998 logarithm calculations, instead of 10988001 of them:

for (int z=1; z<1000; ++z)
  z12[z] = 12 * log2f(z);

for (int f=1; f<11000; ++f) {
  w = 12 * log2f(f - 0.5);
  for (int z=1; z<1000; ++z) {
    int k = ceil(w - z12[z]);
  }
} 

Answer 3

I found:

• Though you have 11Kx1K = 11M combinations of (f, fk),
• the number of distinct value of k is only 294. (you just need 9 bits to represent it).

So definitely you can construct a 2d-array to store the mapping and load it in memory. It looks like

LOOKUP[f][fk] = v, f in 1..11000, fk in 1..1000 
--------------------
v1,1 v1,2 v1,3 ... v1,1000
v2,1 v2,2 v2,3 ... v2,1000
...    ...    ...  ...
v11000,1 , ...     v11000,1000

Since every v is two Bytes, you only need 11Kx1Kx2B = 22MB memory to load this table. That's nothing.

Answer 4

If the loop order is reversed, then the number of repeated values for k goes much higher. In the set there are then only 12977 cases, where k has a "run of one"; The longest run is 618.

This would suggest that a reverse approach would minimize the number log2f calls -- one has to calculate the index n, where k(z,f+n) != k(z,f) (and loop the n instances...)

Anyway, in the end product I doubt the benefit of a huge LUT. Even the approach of using 11000 + 1000 sized tables seems suboptimal to me. Instead I would guess that with only 11000 + 1000 integers there exists either a linear or max 2nd degree piecewise polynomial approximation to log2 that is composed of 8 to 16 pieces.

If such an approach is found, then the polynomial coefficients fit into NEON or XXM registers and allow parallel implementation without memory access.