8
\$\begingroup\$

The natural logarithm can be approximated by the following formula:

$$ln(x) = 2\sum_{k = 0}^{\infty}\frac{(x-1)^{2k+1}}{(2k+1)(x+1)^{2k+1}}$$

I've implemented this formula using Java. Parameter n is the amount of iterations to approximate the natural logarithm.

    static double lnApproximation(double x, int n) {
        double ln = 0;
        for (int k = 0; k < n; k++) {
            double a = Math.pow(x - 1, 2*k + 1);
            double b = 2*k + 1;
            double c = Math.pow(x + 1, 2*k+1);
            ln += a/(b*c);
        }
        return 2*ln;
    }

I've noticed that my first version has potential to be optimized as calculated factors are not stored temporarily for the next iteration.

    /**
     * Faster than lnApproximation by factor ~5.
     * Still slower than Math.log that is calling the native system method.
     */
    static double lnApproximationOptimized(double x, int n) {
        double a = x - 1;
        double aIteration = Math.pow(x - 1, 2);
        double b = 1d;
        double bIteration = 2d;
        double c = x + 1;
        double cIteration = Math.pow(x + 1, 2);
        double ln = a/(b*c);

        for (int k = 1; k < n; k++) {
            a *= aIteration;
            b += bIteration;
            c *= cIteration;
            ln += a/(b*c);
        }
        return 2*ln;
    }

The optimized version is faster than the first version by the factor ~5. On my machine and the parameters x = 1024, n = 32 the first version takes around 50000 ns to finish. The second version takes around 10000 ns to finish the same thing.

I got curious and wonder if there are more micro optimizations possible here?


Two constraints: First, it has to be Java. Second, it has to be this formula.

\$\endgroup\$
  • \$\begingroup\$ Have you measured performance? From one sample, you can't tell that the second algorithm is 5 times faster. You should compare the two curves. \$\endgroup\$ – Thomas Weller Mar 28 at 18:57
  • \$\begingroup\$ @ThomasWeller I've ran each function 10 times and have written down times using System.nanoTime. The first version takes on average 50k ns and the second version takes on average 10k. Of course this is far from professional benchmarking but we can say that the optimized function is faster for sure. \$\endgroup\$ – akuzminykh Mar 28 at 19:06
  • \$\begingroup\$ Have you compared this approximation algorithm with actually invoking the functions? I am worried because the asker at stackoverflow.com/questions/45785705/… benchmarked C++'s std::log function and gets, if I calculate correctly, 6ns per logarithm... That makes me assume I do not calculate correctly, but, more importantly, unless I'm off by a factor of ten thousand, it is faster than 10000ns. There's also a random blog article that accuses the FYL2X assembly instruction of taking up to 700 cycles, which is definitely much less than 10000 nanoseconds. \$\endgroup\$ – my pronoun is monicareinstate Mar 29 at 13:22
  • \$\begingroup\$ @mypronounismonicareinstate The native functions in C++ are probably handwritten assembly. I think 6 ns sounds too fast to believe though. ^^ I've added an answer where you can find my final version. It takes 1000 - 2000 ns in average while the native Math.log function takes 400 - 600 ns in average. Edit: Also the native functions are probably using much more efficient algorithms. \$\endgroup\$ – akuzminykh Mar 29 at 17:13
6
\$\begingroup\$

To speed it up a fraction more add double dIteration = aIteration/cIteration; and double d = a/c; before the loop. In the loop remove a and c and insert

      d*= dIteration;
      ln += d/b;

You then have one multiplication, one division and an addition in the loop, cutting out two excess multiplications in your version.

Edit: you can get further optimization (though perhaps outside the scope of your question) if you look at maths of your series. For 1/sqrt(10) < x <= sqrt(10) you only need around 20 terms (because the ratio between successive terms is$<0.26). For x outside this range, find an n such th at x_0*10**n = x such that x_0 is in this range (by repeated division or multiplication by 10) and use ln(x) = ln(x_0) + n*ln(10).

| improve this answer | |
\$\endgroup\$
6
\$\begingroup\$

Since it's always \$n^2\$ in the second version, by not using java.lang.Math#pow (\$n * n\$), you can save computation time. It takes me ±3k nanoseconds instead of the ±20k nanoseconds with the java.lang.Math#pow.

    static double lnApproximationOptimized(double x, int n) {
        double a = x - 1;
        double aIteration = a * a;
        double b = 1d;
        double bIteration = 2d;
        double c = x + 1;
        double cIteration = c * c;
        double ln = a / (b * c);

        for (int k = 1; k < n; k++) {
            a *= aIteration;
            b += bIteration;
            c *= cIteration;
            ln += a / (b * c);
        }

        return 2 * ln;
    }
| improve this answer | |
\$\endgroup\$
1
\$\begingroup\$

After applying the suggestions the function has changed as follows.

    static final double ln2;
    static final int twoExpLimit;
    static final double precision;
    static {
        int k = 0;
        int n = Integer.MAX_VALUE;
        while (n > 1) {
            n /= 2;
            k++;
        }
        twoExpLimit = k;
        // WolframAlpha
        ln2 = 0.693147180559945309417232121458176568075500134360255254120;
        precision = 1.0e-16;
    }

    /**
     * No way to optimize that further.
     */
    static double lnApproximationPerfect(double x) {
        int exponent = 1;
        int twoExp = 2;
        while (x > twoExp && exponent < twoExpLimit) {
            twoExp *= 2;
            exponent++;
        }
        x /= twoExp;

        double a = x - 1;
        double aIteration = a*a;
        double b = 1d;
        double bIteration = 2d;
        double c = x + 1;
        double cIteration = c*c;
        double d = a/c;
        double dIteration = aIteration/cIteration;
        double iteration = d/b;
        double ln = iteration;

        while (iteration < -precision || iteration > precision) {
            b += bIteration;
            d *= dIteration;
            iteration = d/b;
            ln += iteration;
        }
        return 2*ln + exponent*ln2;
    }

I've added a nicer way to approach the result by using a while-loop. Another thing I've added is the extraction of \$2^k\$ from \$ln(x)\$, thus the approximation is simplified to \$ln(\frac{x}{2^k})+ln(2^k)\$. This resolved the problem where the previous versions were getting very slow with growing \$x\$.

I've tested it against Math.log.

        long startTime;
        double testNumber = 3.14;
        double lnNative, lnCustom;
        long nativeTook, customTook;

        // Run the methods and print something for initialization.
        Math.log(testNumber - 1);
        lnApproximationPerfect(testNumber - 1);
        System.nanoTime();
        System.out.println("Start measurement:");

        {
            startTime = System.nanoTime();
            lnNative = Math.log(testNumber);
            nativeTook = System.nanoTime() - startTime;

            startTime = System.nanoTime();
            lnCustom = lnApproximationPerfect(testNumber);
            customTook = System.nanoTime() - startTime;
        }

        System.out.println("Native: " + nativeTook);
        System.out.println("Custom: " + customTook);
        System.out.println("Native result: " + lnNative);
        System.out.println("Custom result: " + lnCustom);

It's still slower but not by much. I guess you cannot beat handwritten assembly code. ^^

| improve this answer | |
\$\endgroup\$
  • 1
    \$\begingroup\$ Re the code comment No way to optimize that further. -- There most certainly is, which is to use a completely different algorithm. The constraint in the question restricts solutions to the given series. Library implementations are not constrained in this way. \$\endgroup\$ – David Hammen Mar 29 at 18:13
  • \$\begingroup\$ @DavidHammen Please feel free to add an answer that is free from constraints. I've added the constraints just so the question is focused, rather than being too broad and off-topic. I'm well aware of the fact this can be programmed much more efficient with C/C++/Assembly and a different algorithm. There's probably even still potential in Java with bit operations. I'm programming this for learning reasons, so I'd love to see other solutions. \$\endgroup\$ – akuzminykh Mar 29 at 19:15
  • 1
    \$\begingroup\$ If you want no-constraints fast Log approximations, you're going to want to look at getting the FP bit pattern as an integer. The exponent part of the floating point format already uses log2() so you only need a polynomial approximation for log that works over a range like [0.5, 2) or [0.5, 1) to handle the linear mantissa. My answer on Efficient implementation of log2(__m256d) in AVX2 explains in more detail, with links to fast float implementations. The same algorithm that works on every element of a SIMD vector also works scalar if you want. \$\endgroup\$ – Peter Cordes Mar 29 at 21:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.