whow! … rarely created a thread with that intense activity …
I am interested in all aspects of computer math from the nickname of Prof. Kahan, rather ‘Velvel’ than ‘Welvel’, over accuracy, speed, space, usability, universality …
and to have and spread qualified info about it, which is often not present in discussions about shortcoming of ‘as it is’.
For that this thread, besides answering my original question(s), gave lot! of info I’ll need some time to go through.
In pausing this thread … a little awareness test … about neglectible small deviations … pls. guess before trying to code!: setting ‘a’ to ‘0.2 + 0.1’, how many times can you calculate ‘a = 2 * a - 0.3’ in binary64 ( python floats ), pls. post your guesses while not spoilering the result. Then estimate how many calculations are used in building a spaceship, and after that take a look at How Not to Land an Orbital Rocket Booster
The correct rounding is useful at a lower level. The fact that the decimal module has correctly rounded arbitrary precision arithmetic provides a good basis for building intervals and other things over the top and means that a lot of the hard work is already done.
Within arb the arf type reimplements mpmath’s mpf type and provides correct rounding for elementary arithmetic. The arb type is then built on top. Functions like sin are implemented for arb but not for arf and provide bounds rather than correct rounding. While correct rounding is useful it is generally only guaranteed for a single operation so at a higher level it is usually more useful to have something like arb that gives hard bounds on the result of a multistep calculation.
There are many aspects to this but one of them is that to make the intervals fully useful you really want to have multiprecision interval arithmetic so that there is always a general way to make the intervals smaller. In fact this is really the reason I want multiprecision arithmetic in the first place. I don’t actually want thousands of digits but just a few correct ones but I am sometimes prepared to use thousands of digits if necessary to get there.
With arb this is catered for because you have two properties that are very useful as a foundation for rigorous numerics:
Computing at any fixed “precision” p will use reasonably bounded computational expense to give an interval that guarantees to enclose the true result.
If p is made large enough then the interval can be made arbitrarily small at the expense of increased computational cost.
These two properties together make it possible to have controlled computational cost while still being able to solve hard problems without error when needed. Most of the time you won’t need to increase p from the default 53 bits but when you do need to you can.
There is a third property that arb provides loosely and is not guaranteed in general:
For simple operations the relative width of the interval will be something like ~1/2^p for a given precision p.
This is a much weaker guarantee than correct rounding. Strengthening this to something that resembles correct rounding would compromise the “reasonably bounded computational expense” property though. In many applications that is not a good trade off and it is better to compute a wide interval quickly and then afterwards consider whether or not there is any reason to spend more resources making it smaller.
(Oscar did work better in explaining arb approach, so I’m giving up in improving my post on that side. Hopefully, this might be at least partially useful.)
In fact, they (arbitrary-precision floating-point numbers) are used. The arf_t (used for the ball midpoint) is more or less same as the mpfr_t (but no signed zero and “unlimited” exponent). On another hand, the ball radius has limited mantissa (30bit) and fmpz_t as exponent, just as arf_t.
The rest is more or less similar to the interval arithmetic: the result of an operation over balls is a ball (not minimal!) which contains the result of the (mathematically exact) operation applied to any choice of points in the input balls.
No. Signed integers are implemented (fmpz_t, might be better than mpz_t for small values, but same or worse in “big enough” integers) and used, see above.
I don’t see a big difference. Hardly we can do anything here:
>>> x = flint.arb(mid=0, rad=3)
>>> x
[+/- 3.01]
>>> x - x
[+/- 6.01]
It’s possible we’re talking about different things. I’m going solely by Fredrik Johansson’a paper (which Oscar linked to). While Fredrik wrote the arb library, that’s not what the paper is about. The paper is about how some transcendental functions were implemented insidearb. The paper couldn’t be clearer about some of the things I said about it, like:
We use the mpn layer of GMP, since the mpz layer has unnecessary overhead.
and
The main contribution of the paper, described in section 4, is an improved version of Smith’s rectangular splitting algorithm for evaluating Taylor series, in which we use fixed-point arithmetic efficiently …
Whether or when “magical” interval arithmetic is helpful isn’t really relevant here. It’s not used in the function implementations Fredrik sketched. They’re using fixed-point arithmetic (conceptually scaled non-negative integers, and all with the same scale factor so “we can add numbers without any rounding or shift adjustments”). Software help isn’t needed to keep track of rigorous error bounds for the very constrained kinds of arithmetic these specific function implementations use. “Head arguments” are (believed to be) enough in this context, and don’t cost machine cycles
Yes, but it’s not so much “correctly rounded” on its own that’s helpful here but that there are multiple rounding modes supported “as if to infinite precision”. The “directed” rounding modes (to plus or minus infinity) were added to 754 specifically to enable more-efficient computation of narrowest-possible rigorous interval bounds. Intervals do “blow up” in practice, and narrowest-possible is a very worthy goal.
But this is all far afield (to my eyes) to the question of whether decimal would benefit users by adding some support for more common transcendental functions. I’m a “good enough” guy in general, but like the possibility for extreme “purity” when that’s needed. However, I’m essentially never in favor of letting “extreme purity” block “good enough” forever. Particularly in CPython, where other libraries with solid Python bindings already exist, and attract world-class experts to expend the sometimes extreme effort required to achieve extreme purity.
Usually the same here, but as a pragmatic matter if I “just” double the precision and run again, getting essentially the same result to single precision, “good enough” usually rules for me.
Way back when I worked on compiler and library development for Cray Research, selling very expensive early “supercomputers” with novel (at the time) 64-bit HW .floats. People complained they were getting wrong results. So, as a hack, we added a compiler option to chop each float result back to N bits. Jaws dropped around the world Plots of results against the number of bits retained often dramatically showed that their code was numerically unstable. The more bits were retained, the more consistent the results - but generally matched the results they claimed were “right” if N was picked to match the HW they were porting from.
“Extra” precision covers a multitude of (but not all) numeric sins.
Correction gratefully accepted! It’s a mistake I’ve made before, the one and only blind spot I have
Yup, what “serious numeric programmers want” is “all of the above, but stressing different aspects at different times”. “The fast drives out the good” is a general rule of reality, although 754 and Python both worked some against it.
Using interval arithmetic is revealing here. Here’s a start:
>>> from mpmath import iv
>>> p1 = iv.mpf(.1)
>>> p2 = iv.mpf(.2)
>>> p3 = iv.mpf(.3)
Those intervals are zero-width: the input floats are taken as being exact(*):
The higher upper bound on p1+p2 then relentlessly grows across iterations.
>>> a = p1 + p2
>>> for i in range(50):
... a = 2.0 * a - p3
... print(a)
...
[0.2999999999999999889, 0.30000000000000009992]
[0.2999999999999999889, 0.30000000000000021094]
[0.2999999999999999889, 0.30000000000000043299]
[0.2999999999999999889, 0.30000000000000087708]
[0.2999999999999999889, 0.30000000000000176525]
[0.2999999999999999889, 0.30000000000000354161]
[0.2999999999999999889, 0.30000000000000709433]
[0.2999999999999999889, 0.30000000000001419975]
[0.2999999999999999889, 0.30000000000002841061]
[0.2999999999999999889, 0.30000000000005683232]
[0.2999999999999999889, 0.30000000000011367574]
[0.2999999999999999889, 0.30000000000022736257]
[0.2999999999999999889, 0.30000000000045473625]
[0.2999999999999999889, 0.3000000000009094836]
[0.2999999999999999889, 0.3000000000018189783]
...
About every 3 iterations, the uncertainty spreads to cover another decimal digit, and there’s no end to that. While mpmath emulates binary floats with 53 bits of precision by default, its exponent range is essentially unbounded, so this continues far, far beyond the range of Python’s floats.
(*) Although “the input floats” are binary floats, and 0.1 does not mean “one tenth” here - it means “whatever binary float Python created in response to the 0.1 literal”). If you mean “one tenth”, the argument has to be a string instead:
>>> from mpmath import iv
>>> iv.mpf(.1)
mpi('0.10000000000000001', '0.10000000000000001')
>>> iv.mpf(".1")
mpi('0.099999999999999992', '0.10000000000000001')
Your description here suggests that you are just eyeballing the results. If you want to build a robust algorithm over the top of these numerics then you want more than this. This is the situation that motivates having something like arb and where I am quite happy for the intervals to grow as long as I know that it is happening and never end up with a falsely accurate-seeming result. I can deal with a result like x \in [-\infty,\infty] but I don’t want x = 0.1 without even the knowledge that x is definitely a well defined number.
It’s not a property of any specific implementation of interval arithmetic, it’s that the approach itself is inherently blind to the higher-level reasoning humans see at once: “if I subtract x from itself, it doesn’t matter what values x may span: the result is exactly 0”.
Correlations among inputs are common in real life, but IV is blind to them. Much the same holds for any univariate polynomial. For example, x*(x-1). Suppose we know x is in the interval [-3, 3]. IV thinks "OK, x may be as close to +inf as 3, and x-1 may be as close to -inf as -4, so the product may be as close to -inf as -12`.
But that can’t actually happen, because x has the same value in both subexpressions, whatever it may be. The actual smallest value the expression can have is -0.25, if x happens to be 0.5.
Johansson heavily relies on univariate Taylor-series polynomials, and can “out think” these things by head.
Right - but note the “usually”. It means what it says
Which I usually don’t want. “One & done”, speculation, playing around, … plus the obvious one: when I suspect a numeric problem and try this, it’s usually the case that the doubled-precision result shows “yup! numeric problem - the result is waaaay off from the single-precision result”.
Sure. Me too. That’s just not usually the business I’m in, though, and “try it with twice the precision” is good enough for most of my actual purposes.
If I need more than that (which I usually don’t), I’ll usually think about proof next. “Oops! the result isn’t plausible - start over with more precision” (lather, rinse, repeat), would be a last-resort approach to my sensibilities. Although it may well be the most pragmatic approach to solving Johansson’s problem (implementing “correct rounding” in a context where the amount of precision needed apparently cannot be known in advance - “table maker’s dilemma”).
In any case, what does any of this suggest to you about Python’s future? I don’t really see any connection to the OP’s complaint (decimal is missing a bunch of transcendental functions). Perhaps that was the point of your much earlier:
You think Python should support robust interval arithmetic first? If so, that’s pretty much the same to me as saying “no way”. As you also said:
People don’t “find time” for such big undertakings. They have to ruthlessly force time for it, by dropping most other interests for the duration. Somebody has to want it a whole lot, and fight for it. “Interval arithmetic for Python” would predictably be an uphill battle, starting with that there are already capable packages implementing it for binary floats.
So I expect (but don’t know) that the only realistic hope for this path would be for Stefan Krah to implement it inside libmpdec.
You are right, I was talking more about data structures in the current flint (not sure how much they differ wrt ~10 year ago). Indeed, Fredrik uses mpn_* functions directly in code for some elementary transcendental functions, while e.g. in MPFR usually use own primitives (like mpfr_mul, etc, which in turn use mpn layer). Not sure how big gain is, I doubt it’s too much.
That depends on what is “x” and what is “subtraction”.
Just as the ball arithmetic, which is one form of IV. You should expect here some magic:
>>> x = flint.arb(0, 3)
>>> x
[+/- 3.01]
>>> x*(x-1)
[+/- 12.1]
>>> y = mpmath.iv.mpf([-3, 3]);y
mpi('-3.0', '3.0')
>>> y*(y-1)
mpi('-12.0', '12.0')
My mind model for x - x is that I want to subtract two quantities a - b. Each span some interval, maybe infinite, maybe with some nontrivial probability distribution. To account correlations I should use a correlation function, right? Suggested above delta-function sounds too strong for me in real world.
That’s a fair point, but - alas - is typically overlooked entirely when enthusiasts try to “sell” interval arithmetic to the public. The typical user:
has some numerical code they’re queasy about
is told “redeclare all your floats to be of an interval type - then you’ll get guaranteed rigorous bounds on their final values! all by magic”
they do that
they do get rigorous bounds
which are so very wide they’re most often useless
In their original code, the meaning of, e.g., x*(x-1) was utterly clear. After converting to an interval type, its actual meaning is something of a mental nightmare. They weren’t told about this in advance, and there’s very little accessible info about how they “should” go on to rewrite their code to get useful results from intervals.
Kahan has written extensively elsewhere about that there are essentially only two gimmicks he’s seen that “working programmers” (including scientists - not limited to those “naïve” about floating point) can actually use successfully to determine when their code is numerically rotten:
The ability to use extra precision, and preferably significantly more than they “believe they should need”.
The ability to change the default rounding mode.
Those match my experience. Those can very quickly expose problems. The second may be surprising, but it’s really not in hindsight: numerically naïve code is often so very delicate that the tiniest changes can “blow up” very quickly. For example, catastrophic cancellation can change a difference in the last bit into a difference in one of the most-significant bits via a single + or -. In modern HW, changing the rounding mode doesn’t even cost more cycles (although very few programming languages make it easy to change it).
Python’s decimal module makes both very easy for users (except for the initial one-time burden of learning how to spell them). Those are powerful tools just about anyone can use without significant training.
Fixing flawed code remains expert-level work, although it’s often possible to get away with just using more precision.
hi,
being short in time and lacking skills to understand all the points this discussion has branched to … just let me throw in 6 conceptual points and one technical question:
let’s accept: after nearly 100 years of computing, and 40 years of the standard, we are near to but not fully in a position to provide users with what they need, ‘numerical mathematics’.
let’s recall Prof. Kahan: ‘Redesigning computers costs less than retraining people, so it behooves us to adapt computers to the way people have evolved rather than try to adapt people to the way computers have evolved.’
in that sense Tim is one of few people who repeatedly mentions this goal, in other words but in meaning.
reading this thread for a normal user is … chinese … but to safely use computers he needs to know about interna once they deviate from what he’s expecting, in this respect deviate from ‘math’. And they do.
even all our simple hints ‘round your results’, ‘calculate scaled to integers’, ‘use decimal datatypes’ do help often, however not always.
we require users - ‘users’ … you know what type of people that is? - to ‘know about their use case’, ‘choose appropriate datatypes’ and ‘use numerically stable algorithms’ … WHAT!!!
And the technical question:
astonished, if I recall correctly I’d read quite some warnings that switching rounding modes severely harms performance in modern systems because of the discarding of the cache lines … wrong?
I am saying that I would find it useful if there was a library based on the decimal module that could provide real and complex interval arithmetic and some functions like sin and cos. I have written basic prototypes of this but not to the level that is even worth showing to anyone else. My interest in this comes from SymPy but I think it would also be useful elsewhere.
There are a few reasons why it would be useful for SymPy to have this. Currently SymPy has a hard dependency on mpmath and uses it for most numerical evaluation. It is also possible for both SymPy and mpmath to make use of gmpy2 to speed up high precision calculations. It is also now possible for SymPy to make use of FLINT via the python-flint wrapper although currently this is only used in a limited way. In particular SymPy does not yet make use of Arb even though it would be obviously beneficial.
There are two difficulties in SymPy making use of Arb stemming from the fact that FLINT is an optional dependency and is released under LGPL as compared to SymPy and mpmath’s BSD license. Also SymPy and mpmath can both be used as pure Python only without any C extensions and it would be problematic for many downstream users to add hard dependencies on things that are not pure Python. Probably if things were being redesigned today then requiring everything to be pure Python would seem less of a consideration but I think that the license issue would still seem important.
What this means is that without upending lots of things for SymPy to use Arb it would be useful to be able to have a fallback to something implemented in Python. It is acceptable that this fallback be a lot slower than arb as long as it can provide the 2-3 guarantees I outlined above. It is also acceptable that the fallback does not have all the special functions that mpmath/arb have. People who want the maximum speed should be expected to use SymPy with python-flint or just use python-flint or FLINT directly depending on what they are doing. Many people want to use SymPy without needing that speed or even without needing the particular features that might use Arb and so don’t need to bring in python-flint as a hard dependency.
An obvious question then is why that fallback cannot just be mpmath. Interval/ball arithmetic really needs to be implemented carefully from the ground up though and mpmath’s interval arithmetic is more of an incomplete afterthought that would likely need to be rewritten anyway to be sure that it always satisfies the primary guarantee of definitely enclosing the true result. Also many mpmath functions allow unbounded computational expense and building in those bounds as in Arb needs to be done from the ground up as well.
There are some places in the internals of SymPy where real and complex interval arithmetic with exact rational numbers are used. This can be slow even if GMP is used and this code should be rewritten to use adaptive multiprecision interval arithmetic instead. In this application functions like sin and cos are not necessarily needed but it would be good to extend the relevant features to be able to accommodate more functions like that. Much of this code in mpmath, SymPy etc predates libmpdec but if starting afresh today then Python’s decimal module looks a lot more attractive as a foundation for these things.
There are also some situations in which it is better to use interval arithmetic rather than ball arithmetic. The advantages of ball arithmetic apply when you are computing a high precision approximation of something that is in principle known exactly. When you actually want to do true interval calculations with potentially wide intervals it is better to use actual interval arithmetic rather than ball arithmetic especially if you can keep track of the end points exactly to avoid this sort of thing:
>>> a = flint.arb(1, 1) # [0, 2]
>>> a.sqrt()
nan
I have therefore mused about the possibility of implementing real and complex interval arithmetic over the top of the decimal module and implementing some functions like sin and cos. At some point Sergey suggested something like this to me and I think I said that I didn’t see the value specifically of having transcendental functions when working with decimal floating point. Having thought a lot about this since though I can see why this would be useful but it is at least partly just because of the peculiarities of the particular libraries that happen to exist in the Python stdlib and ecosystem rather than because it particularly makes sense to have transcendental functions in a decimal library.
There are also situations where decimal arithmetic specifically is just “better” in some sense primarily because users insist on using numbers like 0.1. Often Arb ends up giving less useful results just because the inputs are not representable in binary floating point. Of course this situation does not apply to calculations that involve transcendental functions that cannot be represented exactly in any floating point representation. Arb is very careful to keep track of exactness for as long as possible in a calculation but with 0.1 in the input it is defeated immediately. If libmpdec had existed at the time then there could have been a good case to say that all of SymPy’s floating point should have been based on the decimal module so that this problem could be “solved” for users.
This is an example of where this interval arithmetic would be used in SymPy:
In [1]: from sympy import *
In [4]: e = cos(pi/7) + sin(2*pi/7)
In [5]: minpoly(e)
Out[5]:
6 5 4 3 2
64⋅x - 64⋅x - 160⋅x + 160⋅x - 20⋅x - 8⋅x + 1
This polynomial is the minimal polynomial of the algebraic number e. Knowing the minpoly is very important/useful in many exact/symbolic calculations. The way that this is computed is that first exact symbolic methods are used to prove that e is a root of this polynomial:
We know now that e is a root of exactly one of these polynomial factors p1(x), p2(x) and p3(x) but we need to know which one. What we can do is compute p1(e), p2(e) and p3(e) and see whether or not they are equal to zero. If we could prove that p3(e) = 0 then we would know the minpoly but no version of floating point can provie that p3(e) = 0. With interval/ball arithmetic though we can prove that p1(e) != 0 and p2(e) != 0 and therefore by a process of elimination that p3(e) = 0.
When we do this though we will find ourselves attempting to compute p3(e) numerically which is an expression that is exactly equal to zero. This is the worst possible case for SymPy’s evalf built over mpmath which goes into overdrive using more and more precision until it hits the limit and returns a result “without precision”:
In [20]: minpoly(e, x).subs(x, e).evalf()
Out[20]: -0.e-124
This sort of thing is exactly the situation where interval/ball arithmetic is really what is needed. We don’t want correct rounding and we don’t want to compute p3(e) with any level of precision (it is zero and cannot be computed with any precision). The precision that we need for p1(e) and p2(e) is basically 1 digit: we just want to know that they are not zero. What we want is something that guarantees not to use massive computational resource and quickly returns an interval/ball containing or not containing zero. We can do this with Arb:
In [39]: p = flint.fmpz_poly([64, -64, -160, 160, -20, -8, 1][::-1])
In [40]: e = (flint.arb.pi()/7).cos() + (2*flint.arb.pi()/7).sin()
In [41]: e
Out[41]: [1.68280035037045 +/- 1.56e-15]
In [42]: p
Out[42]: 64*x^6 + (-64)*x^5 + (-160)*x^4 + 160*x^3 + (-20)*x^2 + (-8)*x + 1
In [43]: p(e)
Out[43]: [+/- 8.12e-13]
In [44]: p(e).mid()
Out[44]: [-8.30446822419617e-14 +/- 9.24e-30]
In [45]: p(e).rad()
Out[45]: [7.28571638918896e-13 +/- 5.42e-30]
In [46]: 0 in p(e)
Out[46]: True
(Note that arbs display strangely: the mid and rad are exact arbs but the +/- refers to the way they are rounded in the decimal string.)
With Arb the algorithm to choose the polynomial factor here is very simple and robust:
Compute p1(e), p2(e), … at some precision p.
Discard any factor that is proven nonzero.
If there are two or more factors remaining then double the precision and try again.
The vast majority of the time we will not need to double the precision but when we do need to we can. The properties that we need from Arb are as I said above:
Currently SymPy will make use of FLINT for things like factorising the polynomial but does not use Arb for these numeric calculations. Wanting to use Arb leaves us wanting something that can be used as a fallback but in general I don’t think that mpmath/evalf work in the way that we really want. I could imagine making a library that is much simpler than mpmath or arb but provides these sorts of calculations. The goal would be to have something small and simple and ideally faster but at least not too much slower than mpmath at low to mid precision while providing the style of computational and accuracy bounds that Arb provides.
This is what I meant when I said that I have thought about this but not found the time to do it. I don’t mean something that would be as big and complex as either Arb or mpmath. I think it is not too hard to make something useful over the top of the decimal module because it provides a good foundation for this. It would not be as fast as Arb but when we want the speed of Arb we will just use Arb. A web search reveals many libraries providing interval arithmetic in Python but I don’t believe any satisfy these requirements:
Real and complex arithmetic.
Multiprecision.
Satisfies the boundedness guarantees.
Pure Python with liberal license.
Just having that for arithmetic without even any transcendental functions like sin would already be useful. I could also imagine situations where we would prefer to use this rather than Arb to get true interval arithmetic and exact decimals but then combine it with Arb for transcendental functions when possible. Realistically SymPy’s fallback in the no-Arb case would always have to use mpmath somehow for many special functions but ideally we could have something Arb-like to handle the more common cases.
Implementing everything with intervals is slower than just making a multiprecision sin function but all of this stuff is slow anyway (why not just use hardware floats?). If you have the multiprecision interval version of the functions then at least for these transcendental functions it is very easy to implement a slow version of correct rounding over the top: exact ties are impossible so just increase p until the interval can be rounded unambiguously.
Would it be better to implement any of these things in C as part of libmpdec? Maybe but ultimately I want Arb-like things with upper and lower bounds, complex arithmetic and so on. That all seems way out of scope for stdlib or libmpdec to me.
The only way that adding functions like sin etc to the decimal module would be useful for the purposes I am describing is if we can use them to implement the interval versions of those functions. That would only work if these functions can provide something like correct rounding or some other proven accuracy bound though: otherwise how would we turn the output into guaranteed upper and lower bounds?
I see one point, avoiding the representation imprecision in binary datatypes. An extreme example: the next representable to pi()/2 in binary64 weights something near 1.57079632679489655800 in decimal. The tangent of that is ( Wolfram|alpha ) ~1.63315E16, the float value is - mostly - presented to the user as 1.5707963267948966, and the tangent for that more than 3 times bigger: ~5.19985E16. The input difference is neglectible small, the output deviation isn’t. That won’t happen with decimal datatypes where some constants and results are also not exact, however ‘What You See Is What Is Computed’. One less step in irritations.
yes, I know that, and then it calculates a tangent for it which is based on the ‘binary deputy’ which is used internally, and is factor 3 apart from the tangent of what is displayed
My quote was in reference to my earlier statement:
where “default” is key. To use rounding mode as a rot-detector, for the first attempt you change rounding mode once, at the start. For example, instead of running under the default nearest/even rounding, change to make round-to-0 (“chop”) the default.
The rounding mode in effect doesn’t change the speed of +-*/ in modern FPUs, and it’s in that sense “changing the rounding mode doesn’t even cost more cycles”. The other strategy (changing precision) can cost a whole lot more in runtime.
As to the speed of changing rounding mode once at the start, it will vary across processor types, but is trivial overall. I can’t imagine why it would affect cache lines (nothing about data representation or addresses changes), but it may well require internal computational pipelines to drain, so that computations that started before the change instruction was reached can complete before the rounding mode (which is typically stored in a “global” FPU control register) changes.
There are ways they could optimize that, but I doubt (but don’t know) that any FPU design bothers.
Which matches Python’s result (at least on my box - it’s really up to the platform’s libm):
>>> math.tan(math.pi/2)
1.633123935319537e+16
The input difference isn’t negligible in context: tan() is zooming off to an infinity near pi/2 (-inf if approaching from above, +inf if approaching from below), and the first derivative’s magnitude is massive for inputs close to pi/2 (on the order of the square of the reciprocal of the input’s distance to pi/2). That’s inherent to tan(), having nothing to do with how floats are represented.
It depends on how much of the internal representation you see. As above, there’s no surprise here if you use the actual decimal value of the binary float approximation to pi/2. If you were working in decimal instead, you’d see the same kinds of “surprises” if values rounded back the actual, full decimal value to something more convenient for human-readable output.
That is, if Python changed to display its binary floats to their exact decimal equivalents, there would be nothing here to complain about. But instead people would complain about “hey! I entered 0.1, but Python displayed 0.1000000000000000055511151231257827021181583404541015625 - how messed up is that?!” .
The problems working in base 10 solves are actually shallow from a numerical analysis view - indeed, rigorous error analysis is actually a bit easier in base 2 (look up “wobble”).
But those shallow problems do bite non-“experts” hard and often.
I don’t say decimals are a universal solution, just was pointed to them, investigated a little and see all implementations I have seen yet as weak, best still mpdecimal. I see them stronger in handling decimal fractions, and weak in performance. And would expect ‘higher’ functions in them with less irritations for users than binaries provide, and with less - not none but less - risks of ‘runaways’.
I as well see binaries better in handling of fractions, ‘free fractions’, 1 / 3 * 3 comes back to 1, and better in speed, partly technical and partly just better matured.
So there is no universal solution, for ‘simple math’ for ‘normal users’ I’d consider decimals preferable. Some of the decimal problems they are used to, they learned in school that free fractions often become difficult in decimals, at my time we did not! learn that decimal fractions become difficult in computers. Evtl. that changed since them, but from comments of irritated users I would deduce that they did not understand it.
About the challenge to estimate ‘* 2 - 0.3’ calculation options, only 2 persons picked up, Tim: sitdown failed, I asked not to spoiler and to calculate in floats. But many many thanks for the crash course in interval arithmetic!!!
Franklinvp: also failed, after less than 1080 iterations you blast away the roof of representable floats, and that’s a very very big value.
My idea is we need to get rid of the small deviations in the first steps like 0.2 + 0.1, else things will - rarely but they will - run out of control at points where we couldn’t need it.
