math module for Decimals

[jerry.carl.mi]

Hi,

I have been looking for a Python module with math functions that would
both eat and spit Decimals. The standard math module eats Decimals
allright but spits floats... herefore exp(sin(Decimal())) produces exp
() of a float :-(

So far, i found:

-AJDecimalMathAdditions (http://www.ajgs.com/programming/
PythonForDownload.html)
-decimalfuncs (http://pypi.python.org/pypi/decimalfuncs/1.4)
-and dmath (http://pypi.python.org/pypi/dmath/0.9)

I tried using the AJDecimalMathAdditions, but ran into issues like dSin
(1E4) would take forever to calculate and would result in sin() >
1 ... If i understand it correctly, the program is using maclaurin
series to calculate the value and since it does not chop off all the
multiples of 2*pi, the maclaurin approximation becomes useless when
its too far from x=0.

I also ran into some issues with the decimalfuncs, but i have never
tried the dmath thing.

Now, i can go and start modifying these modules to behave more like
the standard math module but since i am neither mathematician or
programer, i really hesitate.

So my questions are:

(1) what do folks use when they need to calculate something like exp
(sin(Decimal())) or even more complex things? Any recommendations? Or
am I completely missing something?

(2) Is there any plan to provide a standard python module that would
do that? (Python 4.0? ;-)

Thanks for your comments and help,

Jerry

 

[Gabriel Genellina]

I have been looking for a Python module with math functions that would
both eat and spit Decimals. The standard math module eats Decimals
allright but spits floats... herefore exp(sin(Decimal())) produces exp
() of a float :-(

Which math functions? ln, log10, exp, sqrt already exist as methods of
Decimal instances. At the end of the Decimal docs there are a few
examples, including computing sin and cos (but apparently they naïvely use
a McLaurin series like you noticed in other module).
There is a way to convert a float into a Decimal object with no loss in
precision, see the FAQ (using float.as_integer_ratio(), I think such
function was implemented for exactly this use case)
So you might convert to float, operate, and go back -- if the precision of
"double" operands is enough for you.

 

[jerry.carl.mi]

Which math functions? ln, log10, exp, sqrt already exist as methods of  
Decimal instances. At the end of the Decimal docs there are a few  
examples, including computing sin and cos (but apparently they naïvely use  
a McLaurin series like you noticed in other module).

Hi Gabriel - thanks! For example all goniometric functions are
missing. Or the log(x, base). Or rand(). Sure I can spend time trying
to put it all together but I thought somebody would have done that
already. It seems though that the codes that are out there are not
ready - every one of the modules i mentioned above has some issues.
Maybe I can put bits and pieces together, but if anyone knows of a
well proven module (as is), I would feel much safer using that (again
I am not a mathematician and poking into these algorithms makes me
feel like trying to fix an automatic transmission).

 

[Steven D'Aprano]

Hi,

I have been looking for a Python module with math functions that would
both eat and spit Decimals. The standard math module eats Decimals
allright but spits floats... herefore exp(sin(Decimal())) produces exp
() of a float :-(

You can write your own wrappers. Assuming you always want to return
decimals instead of floats:

from decimal import Decimal
import math

def exp(x):
if isinstance(x, Decimal):
return x.exp()
f = math.exp(x)
return make_decimal(f)

def sin(x):
f = math.sin(x)
return make_decimal(f)

You can get fancier if you like:

def cos(x):
try:
return x.cos()
except AttributeError:
f = math.cos(x)
return make_decimal(f)

That version will automatically use Decimal.cos() in the future, assuming
Decimals grow a cos() method. If they don't, it will continue to work
(although a tad slower than the sin() above.


Converting floats to decimals is tricky, because you have to decide how
much precision you care about. You can keep all the precision of floats.
Depending on your needs, this is either a feature or a gotcha:
(6004799503160661L, 18014398509481984L)


Anyway, here are at least three ways to write make_decimal:


# Python 2.6 only
def make_decimal(f): # exact, no loss of precision
# however, will fail for f = NAN or INF
a, b = f.as_integer_ratio()
return Decimal(a)/Decimal(b)

def make_decimal(f):
return Decimal("%r" % f)

def make_decimal(f, precision=16):
# choose how many decimal places you want to keep
return Decimal.from_float(f, precision)


You choose which suits you best.

So far, i found:

-AJDecimalMathAdditions (http://www.ajgs.com/programming/
PythonForDownload.html)
-decimalfuncs (http://pypi.python.org/pypi/decimalfuncs/1.4) -and dmath
(http://pypi.python.org/pypi/dmath/0.9)

I tried using the AJDecimalMathAdditions, but ran into issues like dSin
(1E4) would take forever to calculate and would result in sin() > 1 ...
If i understand it correctly, the program is using maclaurin series to
calculate the value and since it does not chop off all the multiples of
2*pi, the maclaurin approximation becomes useless when its too far from
x=0.

You might try something like this:

pi = Decimal("%r" % math.pi)
twopi = 2*pi
pi2 = pi/2


def sin(x):
# *** UNTESTED ***
if x < 0:
return -sin(-x)
elif x > twopi:
x = remainder(x, twopi)
return sin(x)
elif x > pi:
return -sin(x - pi)
elif x > pi2:
return sin(pi - x)
else:
return AJDecimalMathAdditions.dSin(x)


You can treat cos and tan similarly but not identically.

However, you might find that it quicker and simpler to just do your maths
calculations as floats, then convert to decimals, as in my earlier
examples.

I also ran into some issues with the decimalfuncs, but i have never
tried the dmath thing.

Now, i can go and start modifying these modules to behave more like the
standard math module but since i am neither mathematician or programer,
i really hesitate.

And so you should. Getting numerical functions correct is damn hard.

However, if all you need is "close enough", then don't be scared -- the
sort of tricks I've shown are probably "close enough", although you may
want a disclaimer about not using them for running nuclear reactors or
flying aircraft. The basic tools Python supplies are pretty good, and you
can probably get a 95% solution with them without being an expert in
computational mathematics.

So my questions are:

(1) what do folks use when they need to calculate something like exp
(sin(Decimal())) or even more complex things? Any recommendations? Or am
I completely missing something?

(2) Is there any plan to provide a standard python module that would do
that? (Python 4.0? ;-)

Guido is very reluctant to put numeric libraries in the standard library,
precisely because it is so hard to get it right. Don't expect Python to
become Mathematica any time soon!

 

[Steven D'Aprano]

def make_decimal(f, precision=16):
# choose how many decimal places you want to keep return
Decimal.from_float(f, precision)

Ah crap, I forgot that from_float() has been left out of the decimal API.
That's very annoying.

Sorry about that.

 

[Steven D'Aprano]

Hi Gabriel - thanks! For example all goniometric functions are missing.

As my earlier email suggested, the easiest way to fix that is to simply
write some wrappers that convert floats to decimal.

Of course, that assumes that you don't need more precision than floats
offer -- it's easy to change the result to get less precision, but if you
really do need more, then you will need to bite the bullet and find (or
write) your own trigonometric functions.

Or the log(x, base).

Easy peasey.
.... return x.log10()/Decimal(base).log10()
....Decimal('5.999999999999999999999999999')


If that's not accurate enough, you will need to do some additional work:

.... power = 1
.... while base**power <= x:
.... power += 1
.... power -= 1
.... x = x/Decimal(base**power)
.... return x.log10()/Decimal(base).log10() + power
....

Decimal('6')


Or rand().

Generating random numbers is an art in of itself. Again, unless you
really need the extra precision, just convert the random number to a
string, then the string to a Decimal.

An alternative is to find one of the many, many published algorithms for
generating random numbers, and re-write that to use Decimals. However,
most such algorithms aren't very random; they only give you a little bit
of randomness, and increasing the precision of the numbers they return
doesn't make them any more random.

Sure I can spend time trying to put it
all together but I thought somebody would have done that already.

I believe they're all waiting for you to do it *wink*

Development of Python is driven by need. If you need something, you can
make a feature request, or you can build it yourself. But if everyone
sits around hoping somebody else will develop the feature, nothing will
ever happen.

It seems though that the codes that are out there are not ready

Not ready for what? Perhaps they are perfectly ready for the needs of the
person who wrote them.

- every one
of the modules i mentioned above has some issues. Maybe I can put bits
and pieces together, but if anyone knows of a well proven module (as
is), I would feel much safer using that (again I am not a mathematician
and poking into these algorithms makes me feel like trying to fix an
automatic transmission).

Decimal is a relatively new module for Python, so it's not surprising
that the functionality is relatively light. But all the hard parts are
done. If you need extra functionality, you have a number of choices:

(1) Find another way to solve your problem, or lower your expectations.
(e.g. do you need arbitrary precision?)

(2) Build the extra functionality yourself.

(3) Pay somebody who is an expert to build it.

(4) Ask nicely and hope somebody eventually decides to build it.

(5) Use another language that already provides exactly all the features
your problem needs.


All of these are valid strategies. Unfortunately, a very popular strategy
is counter-productive:

(6) Complain loudly that Python is rubbish because it doesn't have the
functionality you need, and Somebody Better Fix That RIGHT NOW or else
there will be trouble you betcha.

I recommend against strategy number six

 

[jerry.carl.mi]

Hi Steven... thanks for your kind and extensive reply. Lots of good
food for thought. I know it's easy to complain about lack of
functionality, but it really was not my intention. Python is very cool
as it is and I still know too little about it to even suggest
anything. I just thought maybe I was missing something.

Sure it would make my life much easier if i had a full set of Decimal
(c)math functions and I suppose there are quite a few folks who would
like it as well. *** It's not that it's 100% necessary, it just makes
things easier. *** Just the fact that at least 3 people already tried
to write a decimal math module speaks for how desired and how
difficult it is.

It's really just the goniometric functions that I am missing most at
the moment, so maybe I can figure it out with help of what you said
plus the already existing imperfect modules. Meantime maybe this
discussion will caught Guido's eye... ;-) And btw I do expect that
Python becomes better than Mathematica one day because it's free and
open Maybe when Wolfram retires ;-) Thanks again!

 

[Mark Dickinson]

I have been looking for a Python module with math functions that would
both eat and spit Decimals. The standard math module eats Decimals
allright but spits floats.

Yes: it just converts the input (whether float, int, Fraction or
Decimal) to a float using the __float__ method, and then
applies the usual float->float maths function.

I tried using the AJDecimalMathAdditions, but ran into issues like dSin
(1E4) would take forever to calculate and would result in sin() >
1 ... If i understand it correctly, the program is using maclaurin
series to calculate the value and since it does not chop off all the
multiples of 2*pi, the maclaurin approximation becomes useless when
its too far from x=0.

Implementing these functions *is* kinda tricky, especially if you
want them to be correctly rounded, efficient, and to deal correctly
with all the special cases (infinities, nans, ...). But it's far
from impossible. sin and cos present particular difficulties
for large arguments; probably the range of accepted inputs
would have to be restricted for those functions.

(1) what do folks use when they need to calculate something like exp
(sin(Decimal())) or even more complex things? Any recommendations? Or
am I completely missing something?

Generally, something other than Python. (For example, GP/Pari
or MPFR.) I'd like to be able to use Python for this, though.

A couple of questions for you (I'm curious what additions would
suit your particular use case best):

- are you using Decimal for the base-10-ness or the
extra precision Decimal provides? Or significant zeros?
Or compatibility with existing Decimal code, or what?

- what 3 functions would you most like to see added?
For me, I think it would be something like sin, cos
and atan (or possibly atan2). Once you've got those
three, everything else is fairly easy. In particular,
atan/atan2 at least gives you access to pi.

(2) Is there any plan to provide a standard python module that would
do that? (Python 4.0? ;-)

No, there's no such plan. It's highly unlikely that the Decimal
module will grow any extra math functions: it's designed to
stick very closely to the IBM standard. It's not impossible
that extra Decimal functions could be added in some other
module, but it seems fairly unlikely.

FWIW, I'm the author of the current Decimal log, log10, exp
and pow functions, so I'm probably in a fairly good position
to try to implement reasonably high-quality versions of some
other elementary functions (again, just as an external
addition to the decimal module, not as part of the decimal
module itself). This is an itch I've often wanted
scratched, as well. I might just have a go....
(Help in the form of code, tests, suggestions, etc.
would be welcome!)

Mark

 

[Mark Dickinson]

Ah crap, I forgot that from_float() has been left out of the decimal API.
That's very annoying.

Agreed. It's maybe even annoying enough that a feature request
at bugs.python.org might be honoured. (Hint, hint!)

It's fairly easy to emulate in Python 2.6 and above, using the
as_integer_ratio float method:
Decimal('3.141592653589793115997963469')

Mark

 

[jerry.carl.mi]

- are you using Decimal for the base-10-ness or the
  extra precision Decimal provides?  Or significant zeros?
  Or compatibility with existing Decimal code, or what?

Oh boy, now I will finally prove myself illiterate... well, so be it.
But i am after the extra precision:
0.0

....in this trivial example above I will lose the 1e-16... which may be
an issue if you code something that evaluates slightly more complex
expressions. I would feel much more comfortable if I lose 1e-60. But
in physics, one can get parts of an expression equal to 1e-16 while
(by mistake or not) other parts are > 1. Hence it becomes a greater
puzzle to debug the calculation. Having the possibility to increase
precision would help.

Sure, you can say, there is such a small market for this application,
and maybe I should use other tools. Well, I found Python so much
easier to use for other reasons. And, again, it seems like there is a
desire for it outside of my own office.

- what 3 functions would you most like to see added?
  For me, I think it would be something like sin, cos
  and atan (or possibly atan2).  Once you've got those
  three, everything else is fairly easy.  In particular,
  atan/atan2 at least gives you access to pi.

Agree: sin, cos and atan would do it.

FWIW, I'm the author of the current Decimal log, log10, exp
and pow functions, so I'm probably in a fairly good position
to try to implement reasonably high-quality versions of some
other elementary functions (again, just as an external
addition to the decimal module, not as part of the decimal
module itself).  This is an itch I've often wanted
scratched, as well.  I might just have a go....
(Help in the form of code, tests, suggestions, etc.
would be welcome!)

Mark

Wow, i would never think my posting would go that high in the Python
world. I can't wait to tell my colleagues after these holidays ;-)

If I improve (in my view that is) the existing modules (dmath) etc. i
will keep you posted. For now I am reducing large arguments of
goniometric functions by adding the following into the dmath's sin(x)
and cos(x):

x=Decimal.__mod__(x,Decimal('2')*pi())

Works fine for what i need, but i am sure it's not the right way to do
it.

Thanks Mark!

 

[Mark Dickinson]

But i am after the extra precision:


0.0

Sounds like you don't care too much about the base-10 part,
so there may be other solutions out there.

Have you tried:

1. mpmath?
2. sympy?
3. Sage?

Any of these is likely to be faster than a decimal
solution.

One thing that it would be *really* nice to have is a
set of Python bindings to MPFR---to me, MPFR is the
one obvious way to do high-precision math. Maybe
bindings already exist somewhere. Sage includes them, but
also includes many hundreds of megabytes of other stuff
that you probably don't need or want.

Maybe you could get access to MPFR via ctypes; I don't
have any idea how feasible or complicated this might be.

Sure, you can say, there is such a small market for this application,
and maybe I should use other tools. Well, I found Python so much
easier to use for other reasons. And, again, it seems like there is a
desire for it outside of my own office.

Well, I'd certainly find high-precision sin, cos, ...
within Python useful. I'd guess the market isn't *so* small.

Agree: sin, cos and atan would do it.

I'll see if I can find time. I've just discovered a half-finished
version of atan for Decimal sitting on my computer, complete with
error analysis.

Wow, i would never think my posting would go that high in the Python
world. I can't wait to tell my colleagues after these holidays ;-)

LOL---hardly! I'm just one of hundreds of Python core devs, and I've
only been that for around a year...

If I improve (in my view that is) the existing modules (dmath) etc. i
will keep you posted. For now I am reducing large arguments of
goniometric functions by adding the following into the dmath's sin(x)
and cos(x):

x=Decimal.__mod__(x,Decimal('2')*pi())

Works fine for what i need, but i am sure it's not the right way to do
it.

I don't know of any better way to deal with large arguments.
The main problem is that the reduction step can introduce fairly
large errors: for example, if you're using a value of pi
that's accurate to 10**-20, say, then reducing something of
magnitude 10**5*pi will give a result with error of around
10**-15. As far as I know, this problem is essentially
unavoidable, and it's the reason why implementing sin for inputs
like 10**999999999 isn't feasible.

Mark

 

[jerry.carl.mi]

1. mpmath?
2. sympy?
3. Sage?

Haven't tried those, i guess i have some studying to do.

I don't know of any better way to deal with large arguments.
The main problem is that the reduction step can introduce fairly
large errors:  for example, if you're using a value of pi
that's accurate to 10**-20, say, then reducing something of
magnitude 10**5*pi will give a result with error of around
10**-15.  As far as I know, this problem is essentially
unavoidable, and it's the reason why implementing sin for inputs
like 10**999999999 isn't feasible.

Good point. No tool will work in all parts of the universe (which is
especially true for the universal ski wax).

Let me check the 3 modules you listed above!

 

[jerry.carl.mi]

mpmath... wow... just did what i needed

Thanks, Mark! Hopefully i did not waste too much of your time... and
perhaps this discussion will send other lost sheeps in the right
direction.

(Still, it would make sense to have the goniometric functions in
decimal.)

 

[Steven D'Aprano]

Agreed. It's maybe even annoying enough that a feature request at
bugs.python.org might be honoured. (Hint, hint!)

I can take a hint

http://bugs.python.org/issue4796#


Thank you for looking at this.

 

[Mark Dickinson]

I can take a hint

http://bugs.python.org/issue4796#

Decimal.from_float() implemented by Raymond Hettinger for Python 2.7
and Python 3.1, within 72 hours of Steven submitting the feature
request. If only all issues could be resolved this quickly.

Mark

 

[Mark Dickinson]

Decimal.from_float() implemented by Raymond Hettinger for Python 2.7
and Python 3.1, within 72 hours of Steven submitting the feature
request.  If only all issues could be resolved this quickly.

Rats. I left out the crucial line of that post, namely:

Thank you, Steven!

 

[Steven D'Aprano]

Rats. I left out the crucial line of that post, namely:

Thank you, Steven!

And many thanks to both you and Raymond, firstly for doing the actual
work, and secondly for your patience with my questions on the bug tracker.

 

[Mark Dickinson]

Also, another reason why I'm posting to this thread. I noticed some
error/typo in line 683 of decimal.py located on public svn repo. This is
what is looks like.

        sign = 0 if _math.copysign(1.0, f) == 1.0 else 1

This line looks okay to me; can you say why you think there's a typo
here?

If it's the 'sign = 0' part that's bothering you, that's just
a peculiarity of Decimal: the stored sign value is 0 for
positive numbers, 1 for negative numbers, which I agree
is a little counterintuitive. (Think of the sign
as analogous to the sign *bit* in an IEEE 754 floating-point
number.)

Mark

 

[Mark Dickinson]

haha python-svn # python Lib/decimal.py
  File "Lib/decimal.py", line 683
    sign = 0 if _math.copysign(1.0, f) == 1.0 else 1
              ^
SyntaxError: invalid syntax

Although, It may be only because I ran it through python 2.4.3

Ah yes, that's it: the 'x if b else y' syntax wasn't
introduced until Python 2.5; see

http://docs.python.org/dev/whatsnew/2.5.html#pep-308-conditional-expressions

for more.

Hmm. Maybe we shouldn't be using this syntax in from_float, if it's
the only thing that prevents the trunk version of decimal.py from
being used with Python 2.4. On the other hand, from_float isn't
going to work until 2.7 anyway, since it uses a whole bunch of
new stuff: as_integer_ratio and copysign (both introduced in 2.6),
and bit_length (introduced in 2.7).

Mark

 

[Mark Dickinson]

I get this when importing decimal:

Python 2.7a0 (trunk:68339M, Jan 5 2009, 05:18:41)
[GCC 3.4.6] on linux2
Type "help", "copyright", "credits" or "license" for more information.Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "/usr/local/lib/python2.7/decimal.py", line 138, in <module>
import math as _math
ImportError: No module named math

It's working fine for me. I assume that a simple "import math"
also fails? A couple of possibilities:

(1) the math module built, but Python is looking in the wrong
place for it, for some reason.

(2) the math module failed to build. If you look at the end of
the make output, you'll probably see some lines that look like
the following (but with different modules listed).

Python build finished, but the necessary bits to build these modules
were not found:
bsddb185 gdbm linuxaudiodev
ossaudiodev spwd sunaudiodev

If math is included in the modules listed, then I'd very much
like to know about it.

Mark