gmpy support

[jensgrabner]

gmpy support

[GoogleCodeExporter]

A brilliant patch, thanks! I've committed your changes. It works beautifully; 
the
tests are slightly faster overall and for example mp.dps=10000; exp(pi) is 5.5 
times
faster.

A few comments:

Looks like you diffed in the wrong direction? In either case, I applied the 
patch
manually (which wasn't much work, and gave me a chance to review all the 
changes).

In e.g. sin_taylor, you changed

    k = 3
    ...
        a = ((a * x2) >> prec) // (k*(1-k))
        k += 2

to use k = THREE. Since k is always small, and there are three arithmetic 
operations
done to k for each conversion to MPBASE in the division, I would expect that 
keeping
k as a Python int would be slightly faster. The difference is probably 
marginal, but
some benchmarking would be useful just to be sure.

It seems that pickling mpz doesn't work, at least not on Windows:

Python 2.5.2 (r252:60911, Feb 21 2008, 13:11:45) [MSC v.1310 32 bit (Intel)] on
win32
Type "help", "copyright", "credits" or "license" for more information.
>>> from gmpy import mpz
>>> from pickle import dumps
>>> dumps(mpz(3))
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\Python25\lib\pickle.py", line 1366, in dumps
    Pickler(file, protocol).dump(obj)
  File "C:\Python25\lib\pickle.py", line 224, in dump
    self.save(obj)
  File "C:\Python25\lib\pickle.py", line 313, in save
    (t.__name__, obj))
pickle.PicklingError: Can't pickle 'mpz' object: mpz(3)

This broke the pickling tests in mpmath for me. I added a fix, changing MPBASE 
to
long and long back to MPBASE in __getstate__ and __setstate__. This needed doing
anyway, to make sure __setstate__ works with mpf instances pickled in different
sessions (a test to check that this works would be useful).

In any case, this seems like a bug in gmpy that should be fixed.

It would be nice to add an optional environment variable, e.g. MPMATH_USE_GMPY, 
and
allowing it to be toggled on/off with a command line parameter to runtests.py.

Original comment by [email protected] on 26 Jun 2008 at 3:28

[GoogleCodeExporter]

I just improved the number-to-string conversion functions to make better use of 
gmpy
at high precision, and fixed pidigits.py to work with gmpy.

Calculating 1 million digits of pi now takes just 30.8 seconds (string 
conversion
included), rather than 30 minutes.  Beautiful! This is essentially on par with 
really
fast pi calculating software.

pidigits.py is actually slightly faster than gmpy's pi function;
s=str(gmpy.pi(3321940)) takes 33.0 seconds for me. Strange, because apparently 
gmpy
uses the same algorithm...?

Original comment by [email protected] on 26 Jun 2008 at 3:59

[GoogleCodeExporter]

Thanks for accepting the patch. I'll update my svn copy and do some more 
testing. 
I'll look into optimizing the sqrt functions and bitcount usage next.

Original comment by casevh on 26 Jun 2008 at 5:25

[GoogleCodeExporter]

If you like, I can give you svn write access (I will need your Google account 
email
address).

Original comment by [email protected] on 26 Jun 2008 at 5:36

[GoogleCodeExporter]

Great job casevh!

How is the mpmath speed compared to mpfr now? Do you have some benchmarks 
Fredrik?
What is the bottleneck now?

Original comment by [email protected] on 27 Jun 2008 at 6:37

[GoogleCodeExporter]

Much better than before, but I don't have any benchmarks run on the same system.

The speed of some functions could be improved much further by taking advantage 
of
gmpy's asymptotically fast multiplication and square root.

For example, it becomes feasible to use arithmetic-geometric mean iteration to
compute log (and maybe exp as well). In python-mode, this is slightly slower 
than the
current method up to at least 10,000 digits, but in gmpy-mode it is faster 
already
slightly over 1000 digits (maybe earlier with some optimization) and more than 
twice
as fast at 10,000.

I'll add this code as soon as I've tested it a little more.

Original comment by [email protected] on 27 Jun 2008 at 7:42

[GoogleCodeExporter]

Write access would be great.

I've found a couple of minor changes I'd like to make. I noticed that the 
variable
"one" is used in the code. I think that is too close to "ONE" so I think is 
should
rename all the MPBASE constants to MPONE or MP_ONE.

Original comment by casevh on 27 Jun 2008 at 12:34

[GoogleCodeExporter]

Duh. My google email address is casevh at gmail dot com.

Original comment by casevh on 27 Jun 2008 at 12:35

[GoogleCodeExporter]

Or maybe INT, INT_ONE etc to avoid confusion with mp-floats.

Original comment by [email protected] on 27 Jun 2008 at 1:03

[GoogleCodeExporter]

I implemented a faster algorithm for e as well (binary splitting). One million 
digits
now takes 5 seconds (of which half is computation and half is conversion to 
string).
Surprisingly, this algorithm appears to be slightly faster in Python mode too, 
at
least up to 100,000 digits (I didn't test any further).

Binary splitting could be used for many other constants as well; the others are 
all a
bit more difficult to implement though.

Original comment by [email protected] on 28 Jun 2008 at 10:47

[GoogleCodeExporter]

Thanks for giving me svn access.

I'm running some tests regarding the extra conversions in sin_taylor and 
others. It
looks like it is a little faster with 3 instead of MPBASE(3). I want to work 
through
all the conversions and see which ones are needed so it will be a day or two 
before I
submit the changes.

Regarding the names for the ONE and MPBASE, what do you think of MP_ONE and 
MP_BASE?
I think it is a little more consistent.

I looked at optimizing normalize() using scan1() to count the number of 
trailing zero
bits. There was no significant improvement so I won't include that change. I do 
have
one question about normalize(). It looks like normalize(-0.0) will return 0.0. 
Is
that the expected behavior or should it respect the sign bit when it returns 
fzero?

The binary splitting algorithm for e is nice. My first real high-precision
calculation was calculating e to 3110 decimal places on an HP-41CV calculator 
over 25
years ago. The running time was around 5 days. :)

Original comment by casevh on 30 Jun 2008 at 7:00

[GoogleCodeExporter]

I'm fine with the MP_ convention.

Adding support for negative zero would require changes to a few places in the 
code
where zero is special-cased. Basically, x == fzero would have to be changed to 
check
for either fzero and fnzero (and handle fnzero properly). Negative zero is one 
of
those features of floating-point arithmetic that I've personally never actually
needed, but I don't mind it being added.

Actually, it might also be useful to support scaled zeros (i.e. with an 
arbitrary
exponent). The representation of special numbers would have to be modified 
slightly.

Original comment by [email protected] on 30 Jun 2008 at 10:53

[GoogleCodeExporter]

Blogged: 
http://fredrik-j.blogspot.com/2008/06/gmpy-makes-mpmath-more-speedy.html

Original comment by [email protected] on 1 Jul 2008 at 6:37

[GoogleCodeExporter]

Nice write-up. 

I just added support for enabling/disabling gmpy and the extra assert() tests. 
The
command line options to runtests.py are -nogmpy and -strict. The corresponding
environment variables are MPMATH_NOGMPY and MPMATH_STRICT.

Original comment by casevh on 2 Jul 2008 at 6:05

[GoogleCodeExporter]

Good work. Another use for -strict might be to check that bc == bitcount(man).

Let's also add all this to the documentation.

Original comment by [email protected] on 2 Jul 2008 at 7:29

[GoogleCodeExporter]

The new pidigits is impressive. 1 million digits in 1.91 seconds plus 0.43 
seconds to
convert to decimal. But 10 million digits triggered a segmentation fault in 
gmpy.digits. 

I committed a numeral_gmpy based on numeral_python that can handle extremely 
large
values. 10 million digits is 33.4 seconds to calculate and 10.5 seconds to 
convert.

Original comment by casevh on 3 Jul 2008 at 6:53

[GoogleCodeExporter]

You have a superfast computer then. It takes 7.5s on my laptop:

In [1]: from mpmath import mpf, pi, mp

In [2]: mp.dps = 10**6

In [3]: %time a = mpf(pi)
CPU times: user 7.47 s, sys: 0.06 s, total: 7.53 s
Wall time: 7.57 s

Original comment by [email protected] on 3 Jul 2008 at 9:11

[GoogleCodeExporter]

Is there anything left to be done?

Original comment by [email protected] on 20 Oct 2008 at 12:24

[GoogleCodeExporter]

Probably not. Everything has been working fine for months now.

Original comment by [email protected] on 20 Oct 2008 at 12:30

• Changed state: Fixed