Faster high-precision logarithms

[fredrik-johansson]

Idea: with log(2), log(3), log(5) (and maybe log(7)) precomputed, quickly determine a 5-smooth (or 7-smooth) rational approximation of x.

Hypothesis: if 10-20 bits of argument reduction can be achieved quickly, it ought to make the bit-burst algorithm faster than the AGM for log.

Bonus: instant logarithms of smooth integers.

Can the same trick speed up the argument reduction for other elementary functions?

[rickyefarr]

I was looking at this issue, and perhaps I'm not understanding what you mean. What do you mean by "n-smooth rational approximation"? Do you have a reference where this is defined?

[simonpuchert]

I'm not @fredrik-johansson, but I think "n-smooth rational approximation" is supposed to mean "rational approximation that factors into primes <= n". If that's the case, it's quite a nasty problem (or I'm just blind to the "obvious" solution, as usual).

I found the following "algorithm", which might or might not be useful (in this example we use 2,3,5 and 7):

0. Precompute log(2), ..., log(7)
1. Use log(p/q) to reduce the argument x to the interval [1,p/q] with the following fractions. Actually, keep track of he coefficients a_1, a_2, a_3, a_4 and the current rational approximation
   P / Q = 2^a_1 * 3^a_2 * 5^a_3 * 7^a_4.
2. First, use p / q = 2. Here, it's probably easier to just multiply / shift the argument instead of manipulating the rational approximation.
3. Use e.g. 3/2, 5/4, 9/8, 16/15, 28/27, 50/49, 81/80, 126/125, 225/224 (at most once).
   If x > (P * p) / (Q * q), set P <- P * p and Q <- Q * q, adjust a_1, ..., a_4 accordingly, maybe eliminate common prime factors, proceed to the next fraction.
4. Freestyle, 2401/2400 and 4375/4374 might still be usable like in step 3.
5. log(x) = a_1 * log(2) + ... + a_4 * log(7) + log(x * Q/P)

It's not very fast, but still faster than using a square root to reduce the argument.

If you want more of these "esoteric" fractions, I wrote a program (not on Github, but it's trivial) that finds (probably all) fractions of the form (n+1)/n that factor into primes smaller than a given bound. It seems there are only finitely many of these (and that's very plausible) but I haven't looked for a proof yet.

[fredrik-johansson]

I had somehow forgotten about this.

@simonpuchert Your method seems promising. Do you have some code to experiment with (no need to actually compute the logs; just for the approximation problem)?

There is also a brute force solution: precompute a table of all 7-smooth numbers up to some bound, sorted. For each n in the table, try to use it as a denominator and use bisection to find a numerator. Find the best approximation this way. (There should be some tricks to speed up the search by avoiding common factors.) I'm quite sure that this will work, but I can't tell how it compares to your method.

[fredrik-johansson]

By the way, I suspect this will not actually be faster than the AGM. But it would be interesting to see how far it can be pushed.

[fredrik-johansson]

I think I see how to make this work now. Essentially using Simon's strategy, but a much larger table of smooth fractions for reduction can be precomputed using LLL. I will try some implementation experiments.

[fredrik-johansson]

This is now implemented for exp and log and does give roughly a 2x speedup (much better than I had hoped for).

Still working on a trigonometric version.

[fredrik-johansson]

The only problem now is that Appveyor tests are failing, for who knows what reason.

[fredrik-johansson]

Done.