Newbie question about 'decimal' calculations

Could you be more precise in this statement, i.e. define “range and precision”, where this should hold?

Why not?

This seems to be a counterexample, if I correctly understood your hypothesis:

>>> from decimal import Decimal as D
>>> D(1e100) + D(-1e100) + D(1)
Decimal('1')
>>> D(1e100) + (D(-1e100) + D(1))
Decimal('-2.400819531639191436054718610E+72')

It’s with default settings, but you could adjust this example for the decimal64 context. BTW @tim.one, the gh-53032: support IEEE 754 contexts in the decimal module by skirpichev · Pull Request #122003 · python/cpython · GitHub — waits for you opinion;)

1 Like

Details aren’t just important in math, they’re pretty much everything :wink:

For example, you almost certainly mean you want arithmetic to satisfy the properties of a mathematical field. That’s precise, unambiguous, clear, and quite possibly exhaustive. “And the like” is too fuzzy to guide thought. Here’s a number of equivalent definitions.

Is that missing anything you want? Does it specify something you don’t actually care about?

And here I really have no clear idea of what you intend, starting with “try to!”. Tol? What is that? If you’re not comfortable with technical language, that’s OK, you can make up for it by giving carefully thought-out, concrete examples of things you insist on, would merely like, and can’t tolerate.

>>> from decimal import Decimal as D
>>> 1/(1/D(6))
Decimal('5.999999999999999999999999999')

That violates that, in a field, every non-zero element x satisfies 1/(1/x) == x, where “1” is the field’s multiplicative identity. It fails here. This is not a field.

But the input has only 1 digit. What possible “defined frames of range and precision” could worm around that? Yes, 1/6 has “a lot of digits” - in fact, as many non-zero digits as the working precision, no matter how small or large it may be. Will you give up on division? Or what?

BTW, in decimal you can enable the “inexact” trap, and then an exception will get raised if you ever do something that loses information (to rounding).

>>> import decimal
>>> decimal.getcontext().traps[decimal.Inexact] = True
>>> D(12) / D(6)
Decimal('2')
>>> 1 / D(6)
Traceback (most recent call last):
    ,,,
decimal.Inexact: [<class 'decimal.Inexact'>]
>>> D(6) ** 30
Decimal('221073919720733357899776')
>>> D(6) ** 40
Traceback (most recent call last):
    ...
decimal.Inexact: [<class 'decimal.Inexact'>]

So long as Inexact isn’t raised, results are 100% reliable. If it is raised, information would have been lost to rounding.

In general, “inexact” is so very common that enabling it is useless. However, in some contexts it is useful. Typically when a user isn’t doing “floating point” at all, and using decimal instead for its fast implementations of * and // for huge integers (asymptotically superior to CPython’s implementations of its own binary bigint * and //).

thanks for caring, but you misunderstood,

don’t overstretch 16 >> common! << decimals places in operands and result.

e.g.
1234567890123456 +
2345678901234567 will hold, while
1234.567890123452 +
9999.000000000003 will not, result requires rounding, rounding kills associativity and distributivity.
Besides: I suspect your counterexample to calculate with ‘converted floats’ rather than exact decimals. Try quoting the operands, or work with str( x ).
The 16 digits apply to IEEE decimal64, think python / mpdecimal have default 28 digit precision.

Yes, that’s the problem, that’s why I said ‘try to’. Can we say python decimal / mpdecimal datatypes behave like a field for operations which do not overstretch the defined precision?
It would be easier to explain to dummies: ‘while you stay in xx common digits without rounding your results are correct’. Much easier than ‘expect deviations everywhere, study numerical analysis, use numerically stable algorithms and hope for the best’. :wink:

no, it has more, we fall into rounding. You have the same in schoolhouse decimal math users are used to. Precisely a system matching that or coming near to it isn’t ‘math’, but a subset users are used to. I’d consider users would consider it meaningful to aim for such.
That’s why I’m investigating how much of that can be fulfilled, with which datatypes, with which performance. And - I’m a nitpicker too - also if the results match ‘decimal’, or use binary under the hood, where I catched python in cheating.
No offense, but don’t hold it against me to want to develop and disseminate precise knowledge.

can you make that work for the initial example too?

import decimal
decimal.getcontext().traps[decimal.Inexact] = True
math.asin(decimal.Decimal("0.9999999999999995"))-math.asin(decimal.Decimal("0.9999999999999994"))

→ 0.0

1 Like

I don’t know exactly what that means to you, and can’t flesh it out for you, because I don’t see a sense in which it could be useful.

The “while you stay in xx common digits” part seems redundant: your results are correct if and only if no rounding was needed.

Although I expect that the more precise you make the details (which remain quite vague), “easier” will get more & more out of view.

So you’re saying you’re fine with not being able to compute 1/6 at all? It simply can’t be done, in binary or decimal fp, without rounding.

Then you need to crank up your precision :wink:

If decimal supported asin, that’s trivial: both asin() calls would raise Inexact. All fixed-precision floats (regardless of base) are rational numbers, so both asin() calls are of the form asin(rational), and the only rational argument with an arcsine that’s also rational (so could possibly be a representable float) is 0. That is, only asin(0) wouldn’t raise Inexact. asin(x) for any non-zero float x is not exactly representable.

This is a consequence of Niven’s theorem. If we were working in degrees instead of radians, there would be a few more, but still trivial to special-case. asin(x) is almost always inexact in fp regardless.

Of course none of these are issues in the constructive reals, which, again, is what you really want :smile:,

No, it doesn’t depend on conversion of floats. Similar example you can arrange with exact decimal floating-point numbers:

>>> Decimal('1.23E100') + Decimal('-1.23E100') + Decimal('1')
Decimal('1')
>>> Decimal('1.23E100') + (Decimal('-1.23E100') + Decimal('1'))
Decimal('0E+73')

I doubt there is some finite model of real numbers, where all field axioms holds.

That’s true, but you can run this example in the IEEE decimal64 context too.

So, as Tim, I don’t think you could clarify what you meant by “defined frames of range and precision”.

That corresponds to FE_INEXACT exception for IEC 60559-conforming implementations. Unfortunately, this is not exposed to the python level.

1 Like

It certainly is exposed in Python’s decimal module - in fact, I posted a verbatim REPL session showing decimal.Inexact in action. decimal is a highly conformant implementation.

It’s true that Python’s binary floats offer almost no support for any of 754’s special features.

There can’t be. To start with, it has to contain the multiplicative identity, which is 1. Then it has to be closed under addition. So 1+1+1+1+… all have to be present too.

One of my very favorite :wink: fields is GF(2), a smallest-possible field, with just 2 elements, 0 and 1. Bitwise & is its “multiplication”, and bitwise “^” is its “addition”. I think it sometimes surprises people that & “distributes over” ^:

a & (b ^ c) = (a & b) ^ (a & c)

define frame of range / precision …

I don’t think my idea is bad, I’m just weak to explain it in your words. I don’t think of a fixed range, but somewhat floating ‘while not overstretching available precision’, ‘while results don’t suffer from rounding’. With that understanding you can define / describe where calculations with decimal datatypes should be exact. More ‘which’ than ‘where’. And while they do I and most users will be pleased despite it’s not infinite precise infinite math.
You can do the same for binary datatypes, but not if you include values approximating decimal fractions.
And in next step I’d consider it meaningful to - try to - come as near as possible to this understanding with binary datatypes too. Despite it will not be possible ‘in total’ I’m fan of the idea that some improvements are possible.

Well, in an fp system in base B with prec digits of precision, a result is representable in that system if and only if the infinitely precise value can be expressed in the form:

I * B**J

where I ia an integer with prec base-B digits, and J is an integer. It can’t be made simpler than that without “lying” about something. But that ignores the possibility for underflow/overflow too, so the full truth is actually more complicated. B=10 is no easier than B=2, although people can picture 10**J more easily than 2**J if they’re thinking in decimal.

Regardless, it’s a mass of picky details that nobody will “intuitively” grasp. There’s just no way to, e.g., explain that, in decimal fp, 3/75 is exact but 3/74 and 3/76 round without the tedium of trying to express them all in the constrained I*10**J form. People aren’t good at that in their heads.

For a start, in decimal fp, if they can write the infinitely precisely result as a rational P/Q, and reduce to lowest terms, there’s a chance it’s exactly representable if Q is 1 or its prime factorization contains powers of only 2 and/or 5. If not, then it can’t be expressed in the constrained I*10**J form.

So trying to make it easier piles on more details to try to take “shortcuts”.

I advise concentrating intensely on trying to explain just decimal division. Every detail needed. Then ask whether it’s the improvement you were hoping for. Then ask an actual newbie too :wink:.

Yes, yes, all of that is fine.

All I was referring to, is that making a built-in, in this case sum, implement Neumaier for everyone and not as an opt-in does:

and for a new user of Python, which can at the same time be a new user of floating point arithmetic, for which built-ins are some of the first features of the language that they see, will


I also do not think this is accurate, in the sense of the location, of which function implements the solution. The only honest problem that I see solved, by making sum do Neumaier, that wouldn’t be solved if implemented by other function, or even sum with a flag, is fewer new users of floating point arithmetic asking questions or complaining. This non-problem, I see as detrimental when partial solutions are tried. Of course, there is no full solution to it.

1 Like

Yes, I overstated it - mea culpa. But I think you undervalue it. There are a great many cases now where the builtin sum() delivers the best possible result, and in some common cases provably so. For example, pick 10 floats uniformly at random with -1e6 <= x <= 1e6. No matter how they’re permuted, sum(xs) returns the same as math.fsum(xs) now. That’s far from “some specific simple cases I bumped into”. The sum is almost never stable across the 10! == 3_628_800 permutations using left-to-right addition.

I just tried that 100 times, and the left-to-right sum was never stable across permutations. The number of different results it got across permutations ranged from a low of 3 (seen in 8 cases) to a high of 18 (in 1 case). The mode was 7 different results (in 22 cases). As expected, the current sum(xs) always returned the best possible result regardless of order.

Although, to be fair, the kind of person you describe wouldn’t have any idea at all of whether their sum was “right” or “wrong”. They start with vastly overblown confidence - it can’t be lulled into getting any worse :wink:

That’s not my personal experience, so I can’t accept it as universally true. I personally do use sum() now in some cases where I used to use math.fsum(). Very much faster, a little more convenient to type.

There’s a long history of computing aids trying to improve float sums, because floating + is so very far from meeting associativity (while floating * is very much tamer in that respect - I don’t think anything like a math.fprod() exists anywhere),

It’s why numpy.sum() uses tree addition, in a context where users are generally far more numerically sophisticated. Otherwise unacceptable error distributions when adding “a large” number of floats is about as common as adding a large number of floats.

Tree addition is more of a probabilistic gimmick, though. It still uses just N-1 float additions to add N floats, and isn’t in any way trying to emulate extra precision. Python’s sum() is slower, but easier to prove things about, because it is using twice the native number of float bits to store the running sum from the start. In a context where math.fsum() needs no more adjacent bits than that, it’s straightforward to see that it must return the same result.

Why do some functions get more attention by developers? Additions where relevant to early customers, telephone billing, and quite often used in common users common calculations.

Also here users have some capa to spot deviations and complain about.

Both very less for divisions and up.

Why are additions more often affected by deviations? They suffer from last bit rounding often directly behind or near their limit, for multiplications that mostly happens in bit ‘2 * precision’, and only harms in final if propagates through 00000 or 11111 chains.

Subtractions suffer even more because the deviations / rounding effects are pulled into the bits left for the result.

However having x * x * x often different to pow( x, 3 ) reg. two fit-in roundings leaves air for improvements. ‘FMM’?
And the effort invested for Kahan - Babushka - Neumaier summations and similar justifies a demand to do similar for other functions, for good calculations helping in a technical world all computations should try to be reliable and easy for users.

And to deliver good service to users we ought to provide good calculations, and! qualified info about limitations and alternatives. In that respect his thread - IMHO - outperforms most former attempts including ‘Goldberg’.

1 Like

Please try to explain it in any well defined words.

You are correct that (as we assume a finite model) it should be isomorphic to some finite field, GF(p^n) . And neither resembles real numbers in any way (starting from “wrapping around” all sums a+a+a... due to nonzero characteristic).

Why do you think so? Why fewer, not more? Did you try to search bugtracker for bugs like “sum depends on the order of elements in the list”?

1 Like

I’ll address that a bit since it’s an interesting sociological question that isn’t often asked.

The glib answer is “it depends on the developer”. In general, Python is an open-source project, and in most such most developers “scratch their own itches”. It’s all done in “spare time”, without pay, and most often too with no real recognition of time and effort.

For numerics in particular, Python isn’t at all aimed at numeric computation, and has attracted relatively few core devs with “the right stuff” to do this kind of work.

For the most part, at the start we only strived to write thin wrappers around the platform C’s libm functions, and to implement - in simple and portable ways - the “big int” arithmetic functions needed to support Python’s unbounded ints (for which the platform libm has no support).

There has been no “grand plan” beyond that. People work on what they want to work on, as and when they can make time for it. They’re certainly motivated by legitimate bug reports against code they’re familiar with, but “new features” have a very hard time getting adopted.

Even a world-class expert like Mark Dickenson got some flak for, e.g., adding an innovative (& formally proved correct) algorithm for computing the mathematical floor(sqrt(n)) function exactly for unbounded ints (math.isqrt()).

The original decimal.py was written by a high-school student who got paid a little by a “feel-good” internship program. He did an amazingly good job on it! But it didn’t become suitable (too slow in Python) for “real work” until Stefan Krah released his own libmpdec and contributed careful Python bindings around it.

Beyond that it’s been “a bit here, a bit there”. The last concerted effort I recall wasn’t all that long ago: implementing asymptotically faster ways to do conversions between decimal strings and giant binary ints. All the “number geeks” contributed to that.

Ironically, while that relies in part on the decimal module for its asymptotically superior multiplication algorithms, converting between giant binary ints and decimal.Decimal objects remains quadratic-time. That’s on my own list of “scratches to itch some day”.

Then there are unplanned things. For example, I answered a question on StackOverflow about poor results from numpy’s log1p() function applied to complex arguments. In some ways that was hilarious :wink: Digging into it, I didn’t find an implementation anywhere that actually gave reasonably good results. I thought at first that mpmath was doing a decent job on it, but eventually found cases where it didn’t get any correct bits unless boosting precision to absurdly high values.

It’s hilarious because I thought I could fix it in an afternoon, but ended up crawling on my belly for over a pretty busy week. Looks look everyone just assumed the error analysis done for the real-valued log1p() would just carry over to the complex-valued case, but it doesn’t at all. The complex-valued case was open to entirely new kinds of numerical problems, including catastrophic cancellation so extreme as to destroy all the significant bits.

But non-obvious, applying only in certain cases.

Anyway, I eventually contributed a different implementation to mpmath, which I believe is “almost always” correctly rounded everywhere, and generally using as little “extra” precision as actually need in various cases.

But everything we do carries “opportunity cost” too: what didn’t I do in Python because I was off doing that instead in mpmath? Which is the other hilarious part: log1p(complex) is basically a goofy function to begin with. The real-valued log1p() serves a real purpose, but not the complex-valued one that I can see. numpy’s version, for a complex argument, just adds 1 to the argument and passes the result to its complex-valued log(). A completely mindless “check the box” implementation. So that’s the funny part: all that effort went into something that likely wasn’t worth doing at all :wink:

No idea what “FMM” might mean. Regardless, that’s not going to change in Python. The result of x * x * x is defined by the IEEE-754 standard, and it doesn’t matter to the standard that it may not (and in fact sometimes does not) return the correctly rounded value of the cube of x.

When they differ, pow(x ,3) is typically the better result. But Python defers to the platform libm to do float powers, and has no say in what they return. Python passes on whatever libmreturns, to be compatible with C/C++ code on the same platform.

Both of those (follow the std; play nice with C/C++) are highly desirable for their own sakes.

An experienced numeric programmer who values accuracy above all will stick with the pow(x, 3) spelling, because something that “looks like” a single operation is open to a correctly-rounded implementation. x**3 similarly for the same reason. x*x*x appears guaranteed to suffer at least two roundings in most languages.

Although libm’s pow() need not come with accuracy guarantees either.

2 Likes

A - not funny - point aside … as german news TV was dominated by the catastrophes of Los Angeles Fires and Donald Trump it went unnoticed for me that Starship 7 failed, wasted Millions of Dollars pesting the atmosphere and throwing debris on innocent people.

I don’t have evidence what happened exactly, but long time ago I’d name rocket science as critical in my - yet unpublished - papers. Why? because thousands of people and millions of calculations ‘at the limit’ contribute.

And somewhere between october and december '24 I’d note the following argument that we are not done with IEEE 754:

“For example the IEEE workgroup is discussing demands for “correctly rounded” functions and a new rounding mode “Round to Prepare for Shorter Precision”, and space-X rockets still fail. According to some critizizing videos on the internet also the patriot system still has more fails than hits, but I don’t have evidence on that.”

1 Like

Yes, numeric errors have cost money and lives in the real world. Then again, lots of things have cost money and lives in the real world. Perfect numerics won’t help a bit if, e.g., program logic doesn’t account for a clock “jumping” at a “daylight saving time” transition. I forget the details now, but in some Eurpoean country it’s illegal for a train to have 256 axles, because some cheap automated counter had an 8-bit register that flipped over to 0 at 256, causing “central control” to believe the track was empty instead of carrying a monstrously long chain of cars. It’s always something :wink:

They’re no longer pioneers, though - they’re followers now. A trend toward better rounding for more functions has been in place for decades, and has been making slow progress all along.

"Round to Prepare for Shorter Precision” is already implemented by Python’s decimal module (the otherwise baffling decimal.ROUND_05UP rounding mode).

What they have little chance of getting anymore is new hardware support. Post-754 “fused multiply-add” hardware was driven by industry demand for faster code.

comments on some hints I got in private, pinpointing shortcomings in my argumentation:

don’t want to talk abstractions …

IMHO we / computer math needs understanding of an abstract model to orientate to, and then can discuss which steps come near to it, e.g. in abstract: do we consider binary deputy values as exact ( IEEE ), or can an understanding as deputies for decimal values get further ( ‘normal users’ implicitely think in deputy until seeing deviations ).

what do you want 1/6 to do?

I thought I had mentioned, we have matroshkas of irrationals, rationals / fractions, decimal fractions, binary fractions. Fractions are a superpower to decimals, and I’m grateful that Tim showed how easy to use. If they - not yet extensively tested - work well they are ‘better math’ than decimal, however users seek for decimal and are afraid for fractions since school. So 1/6 is a fraction, is a topic for another day, and having same functionality as in decimal, translating 1/6 + 1/6 into 0.1666667 + 0.1666667 with a result of 0.3333334 instead of 0.3333333 isn’t great, however easier to understand for people than 0.2 + 0.1 != 0.3. One can say one approach is to get nearer to fractions - current understanding of IEEE ‘exact’, but not full success, it’s better than decimal only ‘by chance’, or alternatively try to get nearer to decimal, where IMHO is an option for a limited frame with math logic fulfilled … not full, but as good as possible with limited decimals.

I showed too few interest in mpmath and python-flint.

Yes, sorry, it’s just a plethora of options, and my focus is what can be achieved as ‘standard’, with little special care, learn, import etc. by users. And have a focus more on general than python. I did appreciate and look into the modules Tim explained in short crash courses and learned a lot in that, if mpmath, flint or others can be presented similar easy I’ll like to look.

this forum is more for python than general considerations,

accepted, sorry, my initial question was about something in python, but after that the discussion esp. with Tim was so nice pinpointing and boiling down issues, ideas and POV’s that I didn’t want to stop it. Normally stackoverflow would be an appropriate forum, alas overwhelmed by downvoters less interested in ‘matters’ than formalisms :frowning: .

think less and code more,

already did that, long time of little by little ping pong step up in understanding and realisation, I’m talking about things that can be done. As I’m not experienced in presentations, papers and the like I’d like to find someone qualified for some pre-release reviews …

And comments here:

Perfect numerics won’t help a bit if, e.g., …

yes, but that’s messing apples with pears, or pears with bears, IMHO it’s not justified ‘there are other issues, thus we stop working on ours’ …

They’re no longer pioneers, though - they’re followers now.

:slight_smile: that helps a lot in understanding and dealing with, thank you!

and has been making slow progress all along.

yes, there is a - abstract - blocker while ignoring the special shortcomings for decimal fractions and insisting on ‘exact’ for their binary deputies.

"Round to Prepare for Shorter Precision” is already implemented by Python’s decimal module …

there’s evtl. a misunderstanding, ‘ties away’ is in decimal, alas not ‘full’ e.g. not for conversion to binaries as we discussed some time ago, but ‘RtPSP’ or ‘round to odd’ is likely something different, trying to avoid accumulation of devia in double rounding, don’t ask me ( now ) about details, evtl. later.

What they have little chance of getting anymore is new hardware support.

one point: “FMM” mentioned earlier was an idea of fused-multiply-multiply to avoid one rounding and reduce differences between e.g. pow( x, 3 ) vs. x * x * x what we discussed there, and a question, evtl. a hardware pro is reading here … I think to recall that recent processors have some space of ‘microcode’ reserved for repairs of production bugs or similar. Does someone know how much code can be put there or evtl. in bios? At the moment I’m coding my ideas on application level, once that’s all well I’d like to let it go down into languages, libraries and evtl. hardware.

Unfortunately, size is a problem — sometimes it can be a high price for field axioms. Decimals are of fixed size (precision, exponent limits) — rational numbers might have numerators/denominators of an arbitrary length.

And you can notice, that in long calculations (try e.g. harmonic sums) usually both components of the rational answer tend to grow. Not a surprise: ~60% of integers are relatively prime — not easy to cancel! Try to implement a polynomial gcd for univariate polynomials with rational coefficients, using ordinary Euclid algorithm.

That’s why I expect(ed) there might be issues … a spontaneous idea … if one wants a qualified substitute for decimal math he’d normally do with limited digit length … desk calculator compatible … calculate in fractions, and have some point where fractions become to long and a defined / wanted amount of decimal digits is safely covered … convert to nearest decimal fraction. Could work?
If one starts with decimal fractions ( user friendly UI ) additions, subtractions and multiplications will stay decimal, overtaxing digits can be rounded decimal / schoolhouse compatible, divisions will get uneven denominators, but their rounding not harming more than would occur with pencil and paper and limited to x significant digits.
Tasks like 1/6 + 1/6 could be shown as fractions or 0.16…7 + 0.16…7, but in background kept as fraction so that the sum goes to 2/6 → 0.33…3, on first glance meaningful.
Handele nominator and denominator as product of an array of prime factors could evtl. make some calculations easier?