3
\$\begingroup\$

I'm looking for a simple way to calculate the logarithm without using java.lang.Math, just out of interest.

Here's my try (I adopted some code parts from this answer):

public class Log {
    public static double log10(int x) {
        if (x <= 0) {
            throw new IllegalArgumentException("x must be positive. log_10 not defined for non-positive numbers.");
        }

        int a = 0;
        while (x / 10 > 0) {
            a++;
            x /= 10;
        }

        double r = a;
        double frac = 0.5;
        double x1 = x;
        for (int i = 0; i < 25; i++) {
            x1 *= x1;
            if (x1 > 10) {
                x1 /= 10;
                r += frac;
            }
            frac /= 2;
        }
        return r;
    }

    public static void main(final String[] args) {
        System.out.println(log10(1));
        System.out.println(log10(5));
        System.out.println(log10(10));
        System.out.println(log10(20));
        System.out.println(log10(50));
        System.out.println(log10(1000));
    }
}

How to choose the number of iterations (here 25)? And in what order should the assignments in the while and if block be written? Both "arrangements" are possible without semantic change. Thanks for a review.

\$\endgroup\$
6
  • \$\begingroup\$ Why do this at all? For practice, other? \$\endgroup\$ Commented Dec 27, 2023 at 20:51
  • \$\begingroup\$ Have you read Sean Anderson's bithacks? \$\endgroup\$ Commented Dec 27, 2023 at 21:06
  • \$\begingroup\$ @Reinderien Part of the algorithm comes from this answer. But I came up with the rest myself. That's why I can't provide a citation. I only asked this question for my learning purposes. \$\endgroup\$ Commented Dec 27, 2023 at 23:20
  • \$\begingroup\$ "Have you read Sean Anderson's bithacks?" Looks good, but that seems to be a calculation for the base 2 logarithm. \$\endgroup\$ Commented Dec 27, 2023 at 23:30
  • \$\begingroup\$ He also shows base 10, though it seems to be in integer and not floating-point math, so some adaptation would need to be done \$\endgroup\$ Commented Dec 28, 2023 at 4:21

2 Answers 2

Reset to default
3
\$\begingroup\$

The unchecked IllegalArgumentException is lovely, thank you. Very nice diagnostic error. Referring to the function's name as log_10 seems like a typo. You might consider including the value of x in the message, for the benefit of some poor maintenance engineer chasing down a difficult-to-reproduce bug.

obscure comparison

        while (x / 10 > 0) {

That seems a slightly cryptic way to communicate Author's Intent. Consider rephrasing as (x >= 10), which I think is what you meant. Alternatively, mention a loop variant in a comment, to emphasize we're predicting the effect of the x /= 10 assignment.

meaningful identifiers

        double r = a;

Neither one of these is a great name. And it's unclear why we abandoned a in favor of r. For a local variable to be single letter is not necessarily a problem. But the meaning, and usually the "mental pronunciation", should be clear, e.g. a comment hinting that "the Radix starts out as the Approximation". (No, I don't believe Author intended either of those nouns. IDK, perhaps the Root of some unspecified equation we're solving?) If you had cited a reference, perhaps a wikipedia URL, then the Gentle Reader could come up to speed by tracing that reference's variable names, which you preserved, down into your source code.

Absent comments or a cited reference, we likely want more than single-character names for these. OTOH names like x or i are perfect as-is. The frac identifier is well-chosen.

extract helper

The point of that while loop is to compute very special values of x and of a. Tell us the meaning of those values. You could use a comment, but it would be much better to break out a tiny helper function. It will have an informative name. We can /** document */ its return value(s). We can separately unit test it.

magic number

    ...; i < 25; ...

I imagine there's some special relationship between a power of two and a power of ten. But I'm not seeing it. You have to spell it out for us, minimally with a comment, preferably with a helpfully named MANIFEST_CONSTANT. Here's what I see, but it appears to be off by a factor of about a million:

$ python -c 'print(f"{2 ** 64}\n{int(1e19)}")'
18446744073709551616
10000000000000000000

width of x

            frac /= 2;

Oh, wait! We're cycling through 25 bits?!? How is that relevant?

Don't we want at least 53 bits of significand? Tell us the meaning of what you're computing.

Minimally you have to document that we don't correctly process all input integers, as we only offer limited precision. Or put another way, we don't compute what the Math library computes, so we can't recycle its javadoc documentation for this function. (Well, plus a signature difference. And we substitute an exception for a NaN return, a change which I feel is perfectly reasonable.)

Wow, that's a surprise I didn't see coming!

Note that Math.log10 computes a different function from what java's StrictMath.log10 computes. It's an engineering tradeoff which supports returning an answer quickly, even if it isn't completely accurate. You should decide which specification you choose to embrace.

automated test suite

        System.out.println(log10(1000));

This is nice enough, I guess. Certainly it exercises the target code, and could reveal e.g. embarrassing divide-by-zero errors.

But it is not self-evaluating code. It does not "know the answer". It cannot display a Green bar, nor break the build with a Red bar.

You don't have to use JUnit for automated unit tests. But it's easy, so you may as well take advantage of what it offers.

I note in passing that there's no "negative x" test call.

You have the advantage of being able to call the well-tested Math.log10() function in your tests, so it's easy to draw input values from a wide range of numbers and verify your function returns the correct result, or in any event returns what the Math library returns, which is pretty much the same thing.

algorithm

The Integer.toString() function contains some well-tested looping code, and it offers a convenient .length() getter. You could choose to use that as the core of a rather different log10 implementation.

in what order should the assignments in the while and if block be written?

(Since both arrangements are possible without semantic change.)

I do not especially care; it doesn't have a huge effect on readability.

        while (x / 10 > 0) {
            a++;
            x /= 10;
        }

I have a weak preference for reordering that so the a increment comes last. The idea is to visually highlight that x divided by 10 appears twice. But as mentioned above, a less cryptic comparison would be nice, in which case order wouldn't matter.

            if (x1 > 10) {
                x1 /= 10;
                r += frac;
            }

I really like that we mess with r last, here. We focus on value of x1, comparing it and then altering it. Nice.

It is unclear whether this code achieves its design goals. This is due to the lack of several things:

  • meaningful variable names
  • meaningful helper function names
  • cited reference
  • /** javadoc */ spelling out log10's contract
  • automated unit tests
  • absent the above, // comments could be a poor substitute that partly makes up for a lack of above items

I would not be willing to delegate or accept maintenance tasks on this codebase in its current form.

\$\endgroup\$
3
  • \$\begingroup\$ "You have the advantage of being able to call the well-tested Math.log10() function in your tests, so it's easy to draw input values from a wide range of numbers and verify your function returns the correct result" Do you think a pure comparison, rather than an "is approximately equal" comparison would work? I'd be surprised if the final bit isn't ever different. \$\endgroup\$ Commented Dec 28, 2023 at 5:38
  • \$\begingroup\$ @Acccumulation, recall the "You should decide which specification..." remark. The OP needs to write down a spec, perhaps mentioning ULP, or Absolute error, or Relative error. There can be no code defect absent a spec. We cannot begin to author an automated test suite if we haven't been told what the target code promises to compute. As it happens, OP code has trouble getting even a couple of sig figs right, so that's the least of our worries ATM. \$\endgroup\$ Commented Dec 28, 2023 at 6:03
  • \$\begingroup\$ @J_H I figured out that my code wasn't correct for all input values, for example for log10(1111) the result isn't correct. Perhaps I will create a follow-up question. Thanks for your review. \$\endgroup\$ Commented Jan 2, 2024 at 6:31
0
\$\begingroup\$

The posted codebase would benefit from the addition of an automated test suite, as noted before. I initially took OP at their word, believing the code runs correctly as-is. When I took a moment to write a self.assertAlmostEqual( ... ) check, I quickly saw that it falls apart when we use multiple sig figs. Essentially, the OP while loop is broken.

Here is an implementation that better conveys the intent of the original @Meower68 author.

# from https://stackoverflow.com/questions/739532/logarithm-of-a-bigdecimal
from functools import partial
import math
import unittest
from beartype import beartype
from hypothesis import given
from hypothesis import strategies as st

@beartype
def log_meow(base: int, x: float) -> float:
    assert x >= 1  # This is a pretty significant restriction.
    input_ = x
    result = 0.0

    while input_ > base:
        input_ /= base
        result += 1
    assert 0 < input_ <= base

    fraction = 0.5
    input_ *= input_

    # 1st conjunct is verifying that fraction is still above machine FP precision.
    while fraction + result > result and input_ > 1:
        if input_ > base:
            input_ /= base
            result += fraction
        input_ *= input_
        fraction /= 2.0

    return result

Exercise it with:

log2 = partial(log_meow, 2)
log10 = partial(log_meow, 10)

class LogTest(unittest.TestCase):
    def test_log2(self) -> None:
        self.assertAlmostEqual(0.0, log2(1.0))
        self.assertAlmostEqual(1.0, log2(2.0))
        self.assertAlmostEqual(2.0, log2(4.0))
        self.assertAlmostEqual(3.0, log2(8.0))
        self.assertAlmostEqual(3.0, math.log2(8.0))
        self.assertAlmostEqual(3.16992500, math.log2(9.0))
        self.assertAlmostEqual(3.16992500, log2(9.0))
        self.assertAlmostEqual(3.32192809, math.log2(10.0))
        self.assertAlmostEqual(3.32192809, log2(10.0))
        self.assertAlmostEqual(3.45943161, math.log2(11.0))
        self.assertAlmostEqual(3.45943161, log2(11.0))

        self.assertAlmostEqual(-1.0, math.log2(0.5))
        self.assertAlmostEqual(-1.0, log2(0.5))
        self.assertAlmostEqual(-2.0, math.log2(0.25))
        self.assertAlmostEqual(-2.0, log2(0.25))
        self.assertAlmostEqual(-2.32192809, math.log2(0.2))
        self.assertAlmostEqual(-2.32192809, log2(0.2))

    def test_log10(self) -> None:
        self.assertAlmostEqual(0.0, log10(1.0))
        self.assertAlmostEqual(1.0, log10(10.0))
        self.assertAlmostEqual(2.0, log10(100.0))
        self.assertAlmostEqual(3.0, math.log10(1000.0))
        self.assertAlmostEqual(3.0, log10(1000.0))
        self.assertAlmostEqual(0.47712125, log10(3.0))
        self.assertAlmostEqual(1.11394335, log10(13.0))
        self.assertAlmostEqual(1.47712125, log10(30.0))
        self.assertAlmostEqual(2.47712125, log10(300.0))
        self.assertAlmostEqual(3.30102999, log10(2000.0))
        self.assertAlmostEqual(3.47712125, log10(3000.0))
        self.assertAlmostEqual(3.60205999, log10(4000.0))
        self.assertAlmostEqual(3.5385737, log10(3456.0))

    def test_both(self, limit: int = 200) -> None:
        for actual, log in [
            (math.log2, log2),
            (math.log10, log10),
        ]:
            j = 1.0
            for i in map(float, range(1, limit)):
                # sequential args
                self.assertAlmostEqual(actual(i), log(i))
                self.assertAlmostEqual(actual(1 / i), log(1 / i))

                # power-of-two args
                self.assertAlmostEqual(actual(j), log(j))
                self.assertAlmostEqual(actual(1 / j), log(1 / j))
                j *= 2

Those are just some arbitrary values I found convenient for the edit-debug test cycle when trying to get it to work. Then I tortured it more systematically with a few hundred example {in, out} pairs.

I am very fond of the hypothesis library, which is good at finding edge cases. Haven't seen similar technology on offer in the JVM environment.

    @given(st.floats(min_value=1 + 1e-15, max_value=3.4e38))
    def test_log_hypothesis(self, x):
        self.assertAlmostEqual(math.log2(x), log2(x))
        self.assertAlmostEqual(math.log10(x), log10(x))

Before tackling that dodgy target code I found it convenient to work with an algorithm I was very confident would work -- Newton's method.

def log_newton(base: int, x: float, epsilon: float = 1e-14) -> float:
    """Calculates the root r = log_{base} x, to within epsilon."""
    # https://en.wikipedia.org/wiki/Natural_logarithm#High_precision
    assert x > 0
    pow_ = partial(math.pow, base)
    r = 1.0  # initial guess
    velocity = 1.0

    while abs(relative_error(pow_(r), x)) > epsilon:
        r += 2 * velocity * (x - pow_(r)) / (x + pow_(r))
        velocity *= 0.9999  # avoids stable oscillation around the true root
    return r


def relative_error(x: float, y: float) -> float:
    assert y != 0
    return abs(x - y) / y

Unlike the first routine, Newton will additionally accept positive inputs on the (0, 1] unit interval. In the OP routine I took a stab at using sqrt() convergence for that interval, but no joy. I would be grateful for math or code advice on what would make sense there for positive inputs less than one.

Perhaps these routines will be of help to some future implementor.

\$\endgroup\$
1
  • 2
    \$\begingroup\$ This isn't Java...? \$\endgroup\$ Commented Dec 29, 2023 at 0:04

You must log in to answer this question.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.