7 ms·
How does a computer/calculator compute logarithms?
- constantcrying 2y ago>for if our calculators and computers calculated logarithms inaccurately, as well as exponentials, trig functions, and square roots, to name but a few, a lot of scientific and engineering work would be broken and end in catastrophe. I think I have some news for the writer of the article.
- bongodongobob 2y agoWait until they learn about floating point.
- stephencanon 2y agoWhat are they going to learn? That’s it’s by far the best fixed-size approximation to real number arithmetic we have, sufficient for essentially all of our scientific and engineering needs?
- bongodongobob 2y agoI know that, but there are assumptions you cannot make about equivalence etc. So if he's commenting on inaccuracies of numbers, he may find that surprising. Tough crowd, sheesh.
- account42 2y ago> That’s it’s by far the best fixed-size approximation to real number arithmetic we have (*) * for some use cases (**) ** conditions apply
- constantcrying 2y agoBesides what you said (which I mostly agree with), that floating point calculations can at times be extremely counter intuitive, weird and totally inaccurate. While a single floating point operation gives you the best possible floating point result, the error for two consecutive floating point operations is unbounded. They need to be treated carefully and people need to be aware of their shortcomings. You would be surprised how often I have seen even people here getting hung up over the fact that floating point numbers aren't behaving identically to real numbers.
- deleted 2y ago[deleted]
- kps 2y ago“It would be better to never make a dime of profit than to have a product out there with a problem.” (Today, the generator attached to David Packard powers half of California.)
- 082349872349872 2y agosee also https://www.feynmanlectures.caltech.edu/I_22.html#footnote_source_1 https://www.feynmanlectures.caltech.edu/I_22.html#footnote_s...
- rhelz 2y agoIts polynomials all the way down.
- 082349872349872 2y agoThe VAX POLY instruction is generally held by s/w folk to be the poster child for CISC over-complication. Evidently they were mostly correct, as it's disappeared, but one thing that is missing from many criticisms: POLY as part of the ISA should have enabled the h/w folk to play games with representations of intermediate results (higher precision than the visible register set, or maybe even carry save?) that would not have been possible with RISC algorithmic implementations. [I once worked with people who provided additional crunch for your POLY-heavy problems via extra VAXBI boards. Their bugaboo was that their approach (like GPU offload today?) often succeeded in turning what had been compute-bound problems into i/o-bound ones]
- adgjlsfhk1 2y agoit turns out that fma is the right framework here. it gets you the extra accuracy you want, and can be really fast.
- rhelz 2y agoDamn good point. Sure, there’s no point in microcoding something if it is just doing the same steps as a hand-coded version could do. But that doesn’t mean we just, say, stop implementing floating point operations!! We just make the hardware do it better than hand-coding can.
- myth2018 2y agoBtw: I'm looking for efficient ways to represent numbers of arbitrary precision, and performing mathematical operations on the them. It's for a toy financial calculator I'm writing in my spare time. My knowledge is really lacking in this field and I'm not even sure about what I should be asking/searching for. Could someone suggest me some sources? Preferably online, but books would also be very welcome.
- deleted 2y ago[deleted]
- kolinko 2y agoYou can check how python does it, and also ask chatgpt - I found it to be very helpful in cases like these.
- pavpanchekha 2y agoThe typical library folks use is called MPFR: https://mpfr.org https://mpfr.org It's fast, high quality, and very complete.
- vlovich123 2y agoIt might be also important to note that mpfr is really for scientific purposes that need to deal with such large precision - there are significant memory and performance costs. For financial accounting applications, a fixed precision floating point library (aka fixed point arithmetic) is probably a better choice. Also mpfr is c/c++ but the author didn’t really specify their language and there may be other better options although bindings likely exist.
- denton-scratch 2y agoI once wrote some accounting software for a tiny team in an investment bank. The bosses decided that we should use our brand-new platform, which was cheap to develop for, because most (normal) things could be done without code. All went well for a couple of months; then the customer reported that the system failed for currency conversions to and from Yen. There wasn't enough precision. We had to write an arbitrary-precision maths library, resulting in the project overrunning budget by a factor of three.
- owlbite 2y agoIf I remember correctly, real libm implementations rarely use an analytically derived approximation, they just fit a polynomial rational approximation to minimize the error on the required range and precision.
- teleforce 2y agoPlease check RLIBM project, it's building a collection of correctly rounded elementary functions for multiple representations ( 32-bit float, posits, bfloat16, tensorfloat32) for multiple rounding modes. RLIBM: Rutgers Architecture and Programming Languages Lab's Correctly Rounded Libm: https://people.cs.rutgers.edu/~sn349/rlibm/ https://people.cs.rutgers.edu/~sn349/rlibm/
- pavpanchekha 2y agoThis is a well done post, but please note that computers do not in fact use Taylor series for evaluating polynomials. They do use polynomials! (Or, sometimes, rational functions, but I think the recent trend has been toward polynomials only.) They also use transformations like the log(1+x) - log(1-x) one described in this article. (That one is used in fdlibm, a very influential math library, though more modern libraries seem to not use it.) But the polynomials themselves are usually generated through some other technique, especially the Remez algorithm (https://en.wikipedia.org/wiki/Remez_algorithm https://en.wikipedia.org/wiki/Remez_algorithm) and modern improvements to it like Sollya's (https://sollya.org/ https://sollya.org/) LLL-based floating-point polynomial fitting. It's still an active area of research; the RLibM project (https://people.cs.rutgers.edu/~sn349/rlibm/ https://people.cs.rutgers.edu/~sn349/rlibm/) from Rutgers has introduced a totally new way to fit low-precision polynomials (say, up to 32 bits) using massive linear programming problems (roughly, 4 constraints on 10 variables). Source: am researcher in this field.
- pclmulqdq 2y agoHi Pavel. For people who don't know, this is the person whose lab produced Herbie (https://herbie.uwplse.org/ https://herbie.uwplse.org/). Tacking on to this, I have seen rational approximants in "fast approximate" math libraries, while the ones targeting the highest numerical accuracy are often still using polynomials. Setting up a solver for 0.5 ULP when you have a rational function is definitely a lot harder. In my own work, it also empirically seems that the division creates some problems for precision of the last bit, so your polynomials are longer than you might expect. One interesting approach for bit-accurate rational approximants is to have a rational approximant get most of the way there while fitting the error of your rational approximant to a polynomial, but I don't think there's a good way to figure out where the computationally-efficient split of the problem is (ie how big to make the rational side and the polynomial side).
- pavpanchekha 2y agoYeah, you're right on all of that. I think there have been recent breakthroughs on quadratic (I think at ARITH 23? Yep, just looked, it's the best paper: https://arith2023.arithsymposium.org/program.html https://arith2023.arithsymposium.org/program.html) so hopefully this will become more accessible with time. Though OTOH RLibM-style techniques can get you very short polynomials, to the point that it's hard to imagine beating them by much given the high cost of division operations (like 3x a multiplication).
- tasty_freeze 2y agoWang Labs made an amazing (for its time) calculator, the LOCI-1, quickly replaced by the much more capable LOCI-2. They were from 1964/1965 and were built with core memory and discrete transistors, 1200 of them, no ICs. It operated with 10 digits of precision, and its native hardware could add and subtract. Rather than doing multiplication and division and square root directly, instead the calculator was able to compute logs and antilogs via successive addition and subtraction. Multiplying two numbers might return 199.9999999 instead of 200.0, but since this was built for engineers who were mathematically literate, it was acceptable. There were other calculators that could be programmed to compute logs and antilogs, but they were slow. On the LOCI machines, the results were instant (on a human timescale). Description of the LOCI-2: https://www.oldcalculatormuseum.com/wangloci.html https://www.oldcalculatormuseum.com/wangloci.html Just found this description of the log algorithm -- it used only six constants in its iterative add/subtract/shift algorithm to compute logs and antilogs: https://osgalleries.org/journal/pdf_files/20.2/V20.2P51.pdf https://osgalleries.org/journal/pdf_files/20.2/V20.2P51.pdf
- martin293 2y ago> antilog So exp? e to the power of x?
- anamexis 2y agoYes. I think antilog name used to be more common in engineering, and commonly referred to 10^x rather than e^x. Some calculators even had a log^-1 button.
- eesmith 2y agoNeat Google Ngram for it at https://books.google.com/ngrams/graph?content=antilog&year_start=1800&year_end=2019&corpus=en-2019&smoothing=3 https://books.google.com/ngrams/graph?content=antilog&year_s... showing a peak for "antilog" in 1976, decreasing now to prewar usage of 25% that of peak.
- Thoreandan 2y agoIsn't https://en.wikipedia.org/wiki/CORDIC https://en.wikipedia.org/wiki/CORDIC used a lot?
- teleforce 2y agoMost of calculator and hardware do use CORDIC for real-time implementation. Why the CORDIC algorithm lives rent-free in my head: https://news.ycombinator.com/item?id=40326563 https://news.ycombinator.com/item?id=40326563
- r00tanon 2y agoGood article. I, too, am fascinated by calculators and how they work, especially the HP scientific RPN calculators. So much so, that I ended up developing this: https://play.google.com/store/apps/details?id=com.limpidfox.hp67_2&hl=en_US&pli=1 https://play.google.com/store/apps/details?id=com.limpidfox....
- agumonkey 2y agoIs there a recursive computation for logarithms just like subtraction or division ?
- IIAOPSW 2y agoe^x = sum_n x^n/n! = 1 + x + (x^2)/2 + (x/3)/6... = 1 + x(1 + (x/2)(1 + (x/3)(1 +... exp(x,n) = 1 + (x/n)*exp(x,n+1) this gives you the exponential function as a recursion formula. Just set some condition for when to truncate. From there, a recursive function that does the exact opposite is easy enough. I'm sure there are optimizations you can do on top in the case that you don't get convergence within the first few terms.
- agumonkey 2y agoThanks that's the kind of process I was aiming at but was stuck early.
- xg15 2y ago> While these equations of polynomials contain a finite number of terms, we can have polynomials with an infinite number of terms. These are called series, and one of the simplest types of series is called a geometric series. Slightly OT, but what made infinite sums a lot more sane to me was the understanding that the sums are in fact not infinite - it's just syntactic sugar for a perfectly ordinary limit over a series of finite sums, each with one more term than the last. E.g., sum(1 / 2^n) for n from 1 to +infinity "really" means: lim(sum(1 / 2^n) for n from 1 to m) for m -> +infinity So no infinite loops of additions are actually involved.
- bluenose69 2y agoFor anyone who remembers the Sinclair Scientific calculator, the site http://files.righto.com/calculator/sinclair_scientific_simulator.html http://files.righto.com/calculator/sinclair_scientific_simul... has a good explanation of its way of calculating. It is unlike other machines, overcoming serious electronic limitations with some simple but clever approximate methods. This little machine was my first calculator, a present from my Dad. Man, I loved that thing. It was so small that I could easily hold it, and strike the keys, with one hand. That single-hand aspect proved helpful in my undergraduate physics and chemistry labs. The RPN methodology also made it easy to work through complicated formulae without getting lost in parentheses. Over time, the keys stopped working well. And then stopped working at all. I kept it for a year or two, but then got rid of it, having moved on to another calculator that was more powerful but much less fun. I don't even remember the brand. Looking back four decades, I wish I had kept the Sinclair Scientific, for a memento. But I did not appreciate that a dead little machine like that, in its thin plastic case, would bring back so many fond memories. There's a life lesson in that.
- Razengan 2y agoRIP Sir Clive :')
- cb321 2y agoThe article's impl of its own series [1] centers on the octave [2/3, 4/3] but the error is actually more balanced across [1/sqrt(2), sqrt(2)]. (A series expansion of the error probably explains this - some of the experts here like pavpanchekha | pclmulqdq might be able to explain pithily.) So, here is one way to evaluate the discussed series evaluated across the centered interval and I think it only needs 9 terms for "near IEEE double" NOT the author's 15 terms: import sys, math # arg = number of terms of series n = int(sys.argv[1]) if len(sys.argv) > 1 else 3 def lnATH(x): # Ln via ArcTanH form/series u = (x - 1)/(x + 1) # Can expand loop to expr return sum(2/(2*j+1)*u**(2*j+1) for j in range(n)) x = math.sqrt(0.5) # sqrt(1/2) centers error x1 = x*2 # IEEE FP exponent gets within 1 octave m = 3840 # 4K monitor; adjust to whatever dx = (x1 - x)/m for i in range(m): # Plot error in favorite tool print(x, lnATH(x) - math.log(x)) x += dx If you run it with 8 you begin to see round-off effects. With 9 terms, it becomes dominated by such effects. Anyway 15/9 = 1.66X speed cost which seemed enough to be noteworthy. I mean, he does call it "log_fast.py" after all. (4 seems adequate for IEEE single precision - though you might be 1-bit shy of 24-bits of relative precision at the very edges of octaves if that matters to you.) [1] https://github.com/zachartrand/SoME-3-Living/blob/main/scripts/log_fast.py https://github.com/zachartrand/SoME-3-Living/blob/main/scrip... (and yeah, the code says "16 terms" & article "15 terms", but any way you slice it, I think it's a lot of extra terms).
- pclmulqdq 2y agoI'm just guessing here. What stands out is that the extrema of that octave have the same value. I would guess that this means that things are probably especially well-behaved in that region. If you look at the Chebyshev series of log in that range, you may see some very small coefficients.
- njuffa 2y agoRe: centers on the octave [2/3, 4/3] but the error is actually more balanced across [1/sqrt(2), sqrt(2)] For the computation of the natural logarithm using binary floating-point arithmetic there is no practical difference, i.e. the overall accuracy that can be achieved is basically the same. Here is a worked example based on minimax core approximations extracted from my personal code collection (design and implementation are my own; no guarantee of fitness for any particular purpose). /* natural logarithm */ float my_logf (float a) { float m, r, s, t, i, f; int32_t e; #if LOG_VARIANT == 1 // max ulp err = 0.86280 const float cutoff_f = 0.707106781f; #elif LOG_VARIANT == 3 // max ulp err = 0.85089 const float cutoff_f = 0.666666667f; #endif // LOG_VARIANT if ((a > 0.0f) && (a <= 0x1.fffffep+127f)) { // 3.40282347e+38f m = frexpf (a, &e); if (m < cutoff_f) { m = m + m; e = e - 1; } i = (float)e; f = m - 1.0f; s = f * f; #if LOG_VARIANT == 1 /* Compute log1p(m) for m in [sqrt(0.5)-1, sqrt(2.0)-1] */ r = -0x1.384000p-4f; // -0.076232910 t = 0x1.084000p-3f; // 0.129028320 r = fmaf (r, s, -0x1.0f218cp-3f); // -0.132388204 t = fmaf (t, s, 0x1.226b04p-3f); // 0.141805679 r = fmaf (r, s, -0x1.5427a4p-3f); // -0.166091233 t = fmaf (t, s, 0x1.99a35ap-3f); // 0.200018600 r = fmaf (r, s, -0x1.000416p-2f); // -0.250015587 r = fmaf (t, f, r); r = fmaf (r, f, 0x1.55555cp-2f); // 0.333333433 r = fmaf (r, f, -0x1.fffffap-2f); // -0.499999911 #elif LOG_VARIANT == 3 /* Compute log1p(f) for f in [-1/3, 1/3] */ r = -0x1.0ae000p-3f; // -0.130310059 t = 0x1.208000p-3f; // 0.140869141 r = fmaf (r, s, -0x1.f1988ap-4f); // -0.121483363 t = fmaf (t, s, 0x1.1e5740p-3f); // 0.139814854 r = fmaf (r, s, -0x1.55b36ep-3f); // -0.166846141 t = fmaf (t, s, 0x1.99d8b2p-3f); // 0.200120345 r = fmaf (r, s, -0x1.fffe02p-3f); // -0.249996200 r = fmaf (t, f, r); r = fmaf (r, f, 0x1.5554fap-2f); // 0.333331972 r = fmaf (r, f, -0x1.000000p-1f); // -0.500000000 #endif // LOG_VARIANT r = fmaf (r, s, f); r = fmaf (i, 0x1.62e430p-01f, r); // 0.693147182 // log(2) } else { r = a + a; // silence NaNs if necessary if (a < 0.0f) r = 0.0f/0.0f; // QNaN INDEFINITE if (a == 0.0f) r = -INFINITY; // -INF } return r; }
- bvan 2y agoGreat post and fascinating series of comments.