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\$
3
  • \$\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\$ Mar 28 '20 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 '20 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\$ Mar 29 '20 at 13:22
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).

\$\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;
    }
\$\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. ^^

\$\endgroup\$
2
  • 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\$ Mar 29 '20 at 18:13
  • 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\$ Mar 29 '20 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.