extend methods of decimal module

[Mark H. Harris]

Would it be possible to extend the methods of the decimal module just a bitto include atan(), sin(), cos(), and exp() ?

The module has methods for ln() and sqrt(); and that's great!

I have done some rudimentary searching of the pep history and I'm not finding any pep related to extending the decimal module with other scientific functions.

It is easy to write them in pure python, of course, but I am interested in having the same performance boost with atan(), sin(), cos(), and exp() as Isee with the rest of the decimal module on 3.3/ Is it possible anytime sooner than later?

By-the-by, my hat is off to the person(s) who did the 3.3 work on the decimal module --- the performance boost is astounding. My agm routine for pi100k (which runs on 3.2.3 2Ghz in 10.5 minutes) runs on the same processor in42 seconds on 3.3.4/ very nice.

 

[Terry Reedy]

Would it be possible to extend the methods of the decimal module just
a bit to include atan(), sin(), cos(), and exp() ?

The decimal module implements IEEE 854

The module has methods for ln() and sqrt(); and that's great!

that includes just these. I am not sure if we have an ironclad policy
against adding things not in the standard.

By-the-by, my hat is off to the person(s) who did the 3.3 work on the
decimal module --- the performance boost is astounding.

Stephen (?) Krah.

 

[Mark H. Harris]

The decimal module implements IEEE 854

Thanks Terry... ... long time.

I would like to find out if there is some iron-clad policy about extending the implementation of an IEEE standard... decimal module in this case; I'mjust thinking that this particular extension really fits the python "batteries included" philosophy.

I guess what I'm really asking for are the same routines found in "bc -l" math library. I've finally moved my number crunching stuff to python (from bc) because the performance of "decimal" is finally way better than bc for the moment, and wrapping python are the math routines for control and processing is so much better. Anyway, sure would be nice to have a very speedy atan() function built-in for decimal.

Thanks for the answer... I hope you enjoy the day!

 

[Zachary Ware]

Thanks Terry... ... long time.

I would like to find out if there is some iron-clad policy about extending
the implementation of an IEEE standard... decimal module in this case; I'm
just thinking that this particular extension really fits the python
"batteries included" philosophy.

I guess what I'm really asking for are the same routines found in "bc -l"
math library. I've finally moved my number crunching stuff to python (from bc)
because the performance of "decimal" is finally way better than bc for the
moment, and wrapping python are the math routines for control and processing
is so much better. Anyway, sure would be nice to have a very speedy atan()
function built-in for decimal.

You might consider suggesting a "decimal.math" module on python-ideas.

 

[Terry Reedy]

You might consider suggesting a "decimal.math" module on python-ideas.

Or just dmath. I think this is a better idea than suggesting additions
to decimal itself. For one thing, anything put in decimal would be
subject to change if the function were to be added to the standard. It
is worth noting in such a posting that Krah's speedup make such
functions really feasible. The algorithms could be similar, at least
initially, to the one used for floats.

 

[Zachary Ware]

Or just dmath.

The name (and location) is of course endlessly bikesheddable

I think this is a better idea than suggesting additions to
decimal itself. For one thing, anything put in decimal would be subject to
change if the function were to be added to the standard. It is worth noting
in such a posting that Krah's speedup make such functions really feasible.
The algorithms could be similar, at least initially, to the one used for
floats.

It might be good to even go a step further, and make it "amath" (or
"numbers.math", or ...), a type-agnostic math module. Such a module
would likely be slower than math, cmath, or "dmath", but should end
the proliferation of math modules.

 

[Oscar Benjamin]

Would it be possible to extend the methods of the decimal module just a bit to include atan(), sin(), cos(), and exp() ?

The module has methods for ln() and sqrt(); and that's great!

I have done some rudimentary searching of the pep history and I'm not finding any pep related to extending the decimal module with other scientific functions.

As Terry has pointed out, the decimal module follows IEEE 854 which
doesn't include those.

I think the module sort of has two big use-cases. On the one hand you
have people looking to do financial calculations etc. who are looking
for basic arithmetic with decimal type rounding. On the other hand you
have people like me who see it as a convenient multi-precision library
for when double precision just isn't enough. The module doesn't fully
live up to my use-case because as you say it lacks support for
transcendental functions.

I think, though, that there's a case for having say a "dmath" module
that would export a similar interface to the math and cmath modules
but with functions that work with Decimals in full precision. Someone
has already created a pure Python version of this idea here:
https://code.google.com/p/dmath/

It is easy to write them in pure python, of course, but I am interested in having the same performance boost with atan(), sin(), cos(), and exp() as I see with the rest of the decimal module on 3.3/ Is it possible anytime sooner than later?

Actually the performance difference isn't as big as you might think.
So given the following exponential function:

from decimal import Decimal, localcontext

e = Decimal('2.7182818284590452353602874713527')

def exp(x):
'''7.389056098930650227230427461
'''
# Work in slightly higher precision
with localcontext() as ctx:
ctx.prec += 2
xi, xf = divmod(x, 1)
# Use integer exponentiation
yi = e ** xi
# Now use the Maclaurin series for the fractional part
lastyf = -1
yf = 1
n = 0
fact = 1
xfn = 1
while yf != lastyf:
lastyf = yf
n += 1
fact *= n
xfn *= xf
yf += xfn / fact
# Downgrade precision
return yi * yf

import doctest
doctest.testmod()


I get the following timings:
$ python3.3 -m timeit -s 'from decimal import Decimal as D;
d=D('0.123'); from tmp import exp' 'exp(d)'
10000 loops, best of 3: 32.3 usec per loop
$ python3.3 -m timeit -s 'from decimal import Decimal as D;
d=D('0.123'); from tmp import exp' 'd.exp()'
10000 loops, best of 3: 26.5 usec per loop

So the pure Python exponential function (using the C-accelerated
decimal module) weighs in at 32usec and the pure C version at 26usec.
The intensity of software arithmetic even in C is still dominating the
performance here. The difference becomes more noticeable as you
approach an integer value from below e.g. something like 24.9 but not
more than a factor of 2.

For comparison here's how it looks on the FPU (from Python's perspective):

$ python3.3 -m timeit -s 'd=0.123; from math import exp' 'exp(d)'
10000000 loops, best of 3: 0.149 usec per loop

So that's 2 orders of magnitude faster. It makes sense if you think
about it in terms of CPU instructions: on x87 there are about 5
instructions to call exp in (extended) double precision and takes
about 20-30 cycles (with the time dominated by the F2XM1 instruction).
The number of integer instructions required to compute the above with
the decimal module is massive in comparison.


Oscar

 

[casevh]

I guess what I'm really asking for are the same routines found in "bc -l"
math library. I've finally moved my number crunching stuff to python (from
bc) because the performance of "decimal" is finally way better than bc for
the moment, and wrapping python are the math routines for control and
processing is so much better. Anyway, sure would be nice to have a very
speedy atan() function built-in for decimal.

Have you looked at the gmpy2 ( https://code.google.com/p/gmpy/ ) module?

It supports all the transcendental function available in the MPFR library. I
did a quick performance test of sqrt() and ln() at around 1000 decimal digits.
gmpy2 was about ~200 times faster than the corresponding functions in decimal.

casevh

 

[Mark H. Harris]

Have you looked at the gmpy2 ( https://code.google.com/p/gmpy/ ) module?

casevh

No... was not aware of gmpy2... looks like a great project! I am wondering why it would be sooo much faster? I was hoping that Stefan Krah's C accelerator was using FFT fast fourier transforms for multiplication at least.... but, as Oscar pointed out a native python series algorithm runs right along with the built-in (at least for exp() and ln() ) I have not looked at Krah's code, so not sure what he did to speed things up... (something more than just writing it in C I would suppose).

Thanks for the heads up on gmpy

 

[Mark H. Harris]

Actually the performance difference isn't as big as you might think.

Oscar

You're right. At least my own benchmark on my native exp() vs the built-in was about ~same ~same.

I was hoping that Stefan had used FFT... but I need to talk with him...

Thanks for the reply.

 

[Mark H. Harris]

Or just dmath. I think this is a better idea than suggesting additions
to decimal itself. For one thing, anything put in decimal would be
subject to change if the function were to be added to the standard. It
is worth noting in such a posting that Krah's speedup make such
functions really feasible. The algorithms could be similar, at least
initially, to the one used for floats

Terry Jan Reedy

I have created a project here:

https://code.google.com/p/pythondecimallibrary/

I wrote a dmath.py library module for use with the C accelerated decimal module, that I would like to see merged into the C Python distribution so that folks have it by default... without having to pull it down with GIT, or whatever.

Anyway, thanks for responding (everyone) and for your consideration.

Oh, and one more thing... whoever is doing the work on IDLE these days, nice job! It is stable, reliable, and just works/ appreciate it!

Kind regards,
Mark H Harris

 

[Oscar Benjamin]

I have created a project here:

https://code.google.com/p/pythondecimallibrary/

I wrote a dmath.py library module for use with the C accelerated decimal module, that I would like to see merged into the C Python distribution so that folks have it by default... without having to pull it down with GIT, or whatever.

Hi Mark,

Some points:

1) Why have you committed the code as a .tar.gz file?

2) This function is not such a good idea:
def D(numform):
return Decimal(str(numform))
The Decimal constructor already accepts strings and many types of
numbers. Going via str like this reduces accuracy if e.g. someone
passes a float in.

3) In many places you've written code like this:
prec=dscale(getcontext().prec +7)
sqr = (D(x).sqrt()).__round__(prec)
retscale=dscale(prec)

The preferred way is:
with localcontext() as ctx:
ctx.prec += 7
sqr = round(D(x).sqrt(), prec)

i.e. use a context manager to restore the context even if an error
occurs and use the round builtin rather than the dunder method.

4) You shouldn't be using round/__round__ for precision rounding: it
uses decimal places rather than significant figures. If you want to
round to context precision just use unary +. i.e.:
return +sqr

5) The Decimal.sqrt method already rounds to context precision.
There's no need to compute in higher precision and then round it
yourself (in fact you're invalidating the correctness of the rounding
by double-rounding like this). So really it's just:
def sqrt(x):
return Decimal(x).sqrt()

6) You should organise it in such a way that you're not progressively
increasing the precision each time one function calls another.


Oscar

 

[Mark H. Harris]

Some points:

Thanks so much... you have clarified some things I was struggling with....

1) Why have you committed the code as a .tar.gz file?

um, to save space... well, I know its tiny, but its just a habit I have.... 5kb instead of 25kb...

2) This function is not such a good idea:

def D(numform):
return Decimal(str(numform))

The reason for this is not for strings, but for float literals. All ofmy functions may take a float literal and things still work, because the D(float) function converts the float literal to a string, which is then passed to the Decimal constructor. so... this works sqrt(0.1) or, sqrt(2.01) Without the D() function float literals may not be passed to the Decimal because it really cannot handle them... 0.1 and others cause it a fit.... is there another way to do what I'm looking for here..?

3) In many places you've written code like this:

prec=dscale(getcontext().prec +7)
sqr = (D(x).sqrt()).__round__(prec)
retscale=dscale(prec)

The preferred way is:
with localcontext() as ctx:
ctx.prec += 7
sqr = round(D(x).sqrt(), prec)

Thanks...

i.e. use a context manager to restore the context even if an error
occurs and use the round builtin rather than the dunder method.

This has confused me, because when I look into the module help text I don't see roun() all I see is
__round__(...) how do I know when the round() method exists vs. __round__(...) ??

4) You shouldn't be using round/__round__ for precision rounding: it
uses decimal places rather than significant figures. If you want to
round to context precision just use unary +. i.e.:

return +sqr

whoohoo! thank you.

5) The Decimal.sqrt method already rounds to context precision.
There's no need to compute in higher precision and then round it
yourself (in fact you're invalidating the correctness of the rounding
by double-rounding like this). So really it's just:

def sqrt(x):
return Decimal(x).sqrt()
ok...

6) You should organise it in such a way that you're not progressively
increasing the precision each time one function calls another.

right... well, and if you look at the calls to my __atan__ funcs I am assuming that they are still in the context of the calling func atan()..... right...? so no need to even have the bump up in prec there.

Oscar

thanks Oscar, I really appreciate your comments

marcus

 

[Chris Angelico]

um, to save space... well, I know its tiny, but its just a habit I have... 5kb instead of 25kb...

When you commit changes, though, it has to treat it as a completely
changed binary file - utterly opaque, and it destroys your space
savings too. Normal source code diffs are tiny, and they're readable.
Check out the history of this file; you can easily see every change
I've made:

https://github.com/Rosuav/ExceptExpr/commits/master/find_except_expr.py

Click on any commit to see what happened there. It's hard to do that
with .tar.gz blobs.

ChrisA

 

[Oscar Benjamin]

Thanks so much... you have clarified some things I was struggling with....


um, to save space... well, I know its tiny, but its just a habit I have... 5kb instead of 25kb...

It made it awkward for me to see the code you were referring to. In a
separate thread you asked about how to share code samples. These days
it's common practice to host code somewhere it can be viewed from a
browser. This means that I can do things like link to a line in Chris'
code:
https://github.com/Rosuav/ExceptExpr/blob/master/examples.py#L169

There are many other reasons (as mentioned by Chris) why committing a
binary archive instead of just the plain old files is a generally bad
idea with standard version control systems. A relevant one is that
every change you make to your code, even if you only change 2 lines
will require an additional 5kb of space to store the new .tar.gz file
alongside the old rather than 200 bytes or whatever it is to store the
diff of the code files.

The reason for this is not for strings, but for float literals. All of my functions may take a float literal and things still work, because theD(float) function converts the float literal to a string, which is then passed to the Decimal constructor. so... this works sqrt(0.1) or, sqrt(2.01) Without the D() function float literals may not be passed to the Decimal because it really cannot handle them... 0.1 and others cause it a fit... is there another way to do what I'm looking for here..?

I think you're misunderstanding what's going on:
Decimal('0.1000000000000000055511151231257827021181583404541015625')

There is no binary floating point value that is exactly equal to the
mathematical number 0.1. So when you write 0.1 in Python the number is
rounded to the nearest binary floating point number. When you create a
Decimal from a float it will create a Decimal with the *exact* same
value that the float had. That's exactly what you're seeing above.

When you print a float normally it pretends to have the value 0.1 but
that's really a lie:False

When you convert a float to a string it does not show the exact value
of the float but rather a rounded, shortened decimal form. By
converting float->str->Decimal you are introducing unnecessary
rounding. The Decimal constructor deliberately never rounds its inputs
and every float can be represented as a Decimal with a finite number
of digits. So Decimal(float) is exact but your D function is not.

It is precisely the fact that binary floating point cannot represent
seemingly simple non-integer decimal numbers like 0.1 that leads to
the creation of the decimal module. But why would anyone pass a float
into one of your dmath functions since they're surely much slower than
the math module? The main reason that I can imagine is that they would
like a really accurate result in which case rounding all floats on
input is unacceptable.

This has confused me, because when I look into the module help text I don't see roun() all I see is
__round__(...) how do I know when the round() method exists vs.. __round__(...) ??

The round builtin function calls __round__ on any object. Dunder
methods are not supposed to be called directly.


Oscar

 

[Terry Reedy]

Oh, and one more thing... whoever is doing the work on IDLE these
days, nice job! It is stable, reliable, and just works/
appreciate it!

As one of 'them', thank you for the feedback. There are still some bugs,
but I hit them seldom enough that I am now looking at enhancement issues.

 

[Chris Angelico]

As one of 'them', thank you for the feedback. There are still some bugs, but
I hit them seldom enough that I am now looking at enhancement issues.

Aside from having to be careful not to paste in a non-BMP character (a
Tk limitation, I believe, and currently being tracked at issue 13153),
I'm pretty happy with IDLE too. It's my standard way to fiddle with
Python code. Easier to use than the inbuilt interactive mode, as Alt-P
will recall an entire block command, rather than one line at a time.

So, I echo that sentiment of thanks.

ChrisA

 

[Mark H. Harris]

Decimal('0.1000000000000000055511151231257827021181583404541015625')

hi Oscar, well, that's not what I'm doing with my D()... I'm not just making D() mimic Decimal... look inside it... there's a str() call.... consider the following experiment and you'll see what I'm talking about...

Decimal does not keep 0.1 as a floating point format (regardless of size) which is why banking can use Decimal without having to worry about the floating formatting issue... in other words, 0.0 is not stored in Decimal as any kind of floating value... its not rounded.... it really is Decimal('0.1').

Ok, so for the experiment: consider this exchange from IDLE:

You will notice that Decimal(.1) and my D(.1) work against PI differently... in fact, Decimal(.1) decimates PI.... <sorry for the PUN>

The reason is that Decimal(.1) stores the erroneous float in the Decimal object including the float error for .1 and D(.1) works correctly because the D(.1) function in my dmath.py first converts the .1 to a string value before handing it to Decimal's constructor(s)

Now, try this experiment: again from IDLE and dmath.py

You see? All of my functions make certain that when the Decimal objects are created that the string constructor gets called.... so that, in this case, .1 really is 0.1 precisely, not floating format.

and take a look at this:

With my dmath D() function the object holds 0.1 precisely... its not the float 0.1 rounded ...


shifting gears a bit....

The real reason I pushed a commit of a gzip tarball to code.google is because I've never used code.google before (thought I'd give it a try) and I didn't realize that its only pull / download is a ZIP... so I am going to take your advice and remove the tarballs and just put the code up there... if I can get my GIT to work ... its GIVing me fits (probably because I don't have it configured right on this system.

Thanks for the discussion... its helping me get this stuff into my brain, and its giving me a chance to discuss what's at issue between floats and Decimal.

kind regards,
marcus

 

[Wolfgang]

hi Oscar, well, that's not what I'm doing with my D()... I'm not justmaking D() mimic Decimal... look inside it... there's a str() call.... consider the following experiment and you'll see what I'm talking about....
You will notice that Decimal(.1) and my D(.1) work against PI differently... in fact, Decimal(.1) decimates PI.... <sorry for the PUN>
The reason is that Decimal(.1) stores the erroneous float in the Decimal object including the float error for .1 and D(.1) works correctly because the D(.1) function in my dmath.py first converts the .1 to a string value before handing it to Decimal's constructor(s)

What Oscar tried to make clear is that Decimal(.1) does not store an "erroneous" float in the Decimal object, but an "honest" one.
Maybe this becomes clearer with calculations instead of literal values:
0.1

now try:
D(1000000000000000055511151231257827021181583404541015625/10**55)
vs
Decimal(1000000000000000055511151231257827021181583404541015625/10**55)
for an example that it is not better, in general, to go through the str()-mediated rounding you are doing.

Of course, if the user *intended* that float(.1) to mean exactly 0.1 then your D class does a good job at guessing, but only then.
The Decimal class, on the other hand, leaves the choice to the user as you are illustrating perfectly in your own example below:

Now, try this experiment: again from IDLE and dmath.py

Decimal('0.31415926535897932384626433832795')

see? Decimal(.1) gives the exact representation of the float, but Decimal('..1') is the way to get your D()'s default behavior.

You see? All of my functions make certain that when the Decimal objectsare created that the string constructor gets called.... so that, in thiscase, .1 really is 0.1 precisely, not floating format.

Yes, just that this is almost never what your users will want to happen ...
The Decimal constructor is really well designed, so stick to it !

Best,
Wolfgang

 

[Oscar Benjamin]

hi Oscar, well, that's not what I'm doing with my D()... I'm not just making D() mimic Decimal... look inside it... there's a str() call.... consider the following experiment and you'll see what I'm talking about...

I understood what your code is doing but I'm not sure if you do.
Calling str on a float performs an inexact binary to decimal
conversion. Calling Decimal on a float performs an exact binary to
decimal conversion. Your reasoning essentially assumes that every
float should be interpreted as an approximate representation for a
nearby decimal value. This is probably true if the user wrote "a =
0.1" but is generally not true in the kind of numeric code that is
likely to be using the transcendental functions defined in your dmath
module.

Calling Decimal(str(float)) introduces entirely avoidable inaccuracy
in your code when the primary purpose of your code as accuracy!


Oscar