Closed GoogleCodeExporter closed 9 years ago
I will work on adding these functions. Just to confirm the behavior, add(xman,
xexp,
yman, yexp, prec, rounding) should return (xman*2^xexp + yman*2^yexp) rounded
to
prec digits with the specified rounding mode.
Should I also include mult?
Original comment by casevh
on 31 Aug 2009 at 10:56
Correct. I think mult is unnecessary because it can be written trim(xman*yman,
xexp+yexp) without loss of either speed or clarity (but benchmarking could
prove me
wrong).
Original comment by fredrik....@gmail.com
on 1 Sep 2009 at 4:22
Can I make any assumptions about the maximum size of an exponent or should I
allow
them to be unbounded?
Should I raise a ZeroDivisionError when dividing by 0?
Should I raise a ValueError when taking the square root of a negative number?
Original comment by casevh
on 2 Sep 2009 at 4:50
Should I require a value for prec?
Original comment by casevh
on 2 Sep 2009 at 4:53
> Can I make any assumptions about the maximum size of an exponent or should I
allow
them to be unbounded?
Mpmath needs them to be unbounded (Python int and long).
> Should I raise a ZeroDivisionError when dividing by 0?
> Should I raise a ValueError when taking the square root of a negative number?
Yes and yes.
> Should I require a value for prec?
Mpmath requires that prec=0 is supported and perform an exact operation for
trim and
add. (The way to specify this could be modified.)
Original comment by fredrik....@gmail.com
on 2 Sep 2009 at 7:44
I have a few more questions.
Can prec be assumed to be a C long? (That constraint is present in the existing
_mpmath_normalize.)
Is there default rounding mode that is most frequently used?
Does the requested prec value change frequently?
I'm trying different approaches to parsing the arguments in an attempt to
minimize
calling overhead. If I can decrease the number of objects that are passed, I
can
decrease the overhead. BTW, the overhead for calling the prototype for trim
(just
parsing and returning values) with Python 3.1 is significantly less than for
Python
2.6: ~200 uSec vs. ~250 uSec.
Original comment by casevh
on 4 Sep 2009 at 7:28
> Can prec be assumed to be a C long?
Yes, definitely.
> Is there default rounding mode that is most frequently used?
I wrote 'f' above but I think it should probably be 'd'. While 'n' is probably
most
common, it will always be passed explicitly, at least the way the code is
written
right now. This might be possible to change.
Original comment by fredrik....@gmail.com
on 4 Sep 2009 at 7:45
> Does the requested prec value change frequently?
Yes; even when the user-set working precision is fixed, many functions change
the
temporary precision by all kinds of amounts.
Original comment by fredrik....@gmail.com
on 4 Sep 2009 at 7:47
I've added _mpmath_trim and _mpmath_mult to svn. Using _mpmath_mult(...) is
faster
than _mpmath_trim(xman*yman, xexp+yexp) so I'll leave it in. I'm working on
_mpmath_add next.
Original comment by casevh
on 6 Sep 2009 at 5:05
Nice, thanks! The trim function recovers the performance lost to having to wrap
_mpmath_normalize. There is a 'fprintf(stderr, "zero\n");' left in there though
which
I had to comment out.
I don't think _mpmath_mult will have any effect. In the route through
mpf.__mul__,
the difference is dwarfed by overheads. Besides trim, the helper function for
addition is probably the most important one.
Obtaining a much more substantial speedup will involve implementing the number
type
for mpmath entirely in C or Cython. Overall, the changes I'm experimenting with
right
now are intended to facilitate such a switch by substantially reducing the
amount of
interfacing between code at different levels of abstraction. In the short term,
the
only significant speedup (even with the helper functions) will be for complex
arithmetic.
Original comment by fredrik....@gmail.com
on 6 Sep 2009 at 2:53
The trim function causes a segfault somewhere... I'll try to hunt it down.
Original comment by fredrik....@gmail.com
on 6 Sep 2009 at 3:26
Ah, the answer is simple: it returns an mpz as exponent but the code in mpmath
expects it to be int or long.
Original comment by fredrik....@gmail.com
on 6 Sep 2009 at 5:53
Although the exponent type should be fixed, the direct cause of the segfault is
comparison of an mpz with a Python float infinity:
mpz(3) > 1e1000
Original comment by fredrik....@gmail.com
on 6 Sep 2009 at 6:01
Thanks for finding the mpz(3) > 1000 bug. I'll get that fixed today.
I'll also convert the exponent back to a Python integer.
After I get gmpy 1.10 released, I want to replace the mpf layer with something
better
(i.e. correctly rounded and basic transcendental functions). One option is MPFR
but
that introduces additional dependencies and possible licensing issues. I'd
prefer to
provide better support to mpmath.
Original comment by casevh
on 6 Sep 2009 at 6:27
With the upcoming changes in place, it should be feasible for mpmath to
directly use
a gmpy-based floating-point type, assuming that it implements compatible
functionality.
Certainly the biggest source of friction is mpmath's use of a global (or rather,
bound to a context object), rather than instance-controlled, precision. What
would
you think of using something similar in gmpy? One alternative that could work
is to
bind the precision to the mpf type. Then gmpy.mpf(3.14, 100) would be
equivalent to
gmpy._create_mpf_type(prec=100)(3.14). If the precision of the type objects were
mutable, mpmath could subclass gmpy.mpf and modify the precision of that type
in-place when changing precision.
Also, mpmath obviously needs both real and complex numbers, so gmpy mpfs should
ideally support complex values natively. Ditto for infinities.
Original comment by fredrik....@gmail.com
on 7 Sep 2009 at 10:38
A question on the mpf_add implementation in libmpf.py: What is purpose of the
"offset
> 100" portion of the line "if offset > 100 and prec:" ? (My best guess is that
it
keeps values close together in the fastest code path or am I missing something?)
I'm still thinking about how best to control precision.
Original comment by casevh
on 10 Sep 2009 at 6:13
Yes; it assumes that shifting a couple of extra bits is cheaper than comparing
the
exact sizes. This is not necessarily true in C. (Though maybe comparing limb
sizes
first can save a little bit of time.)
For reference, I'm attaching the naive Cython implementation I wrote for Sage
(not
finished). It uses mpz for exponents too (as I think is just as fast or faster
in a
pure C(ython) implementation. I'm not sure if it's correct; I haven't tested it
extensively.
Original comment by fredrik....@gmail.com
on 10 Sep 2009 at 9:39
Attachments:
I've committed _mpmath_add and fixed a bug in _mpmath_trim. I've done limited
testing. The functions still return the exponent as an mpz. I will get that
fixed.
Original comment by casevh
on 16 Sep 2009 at 4:54
I think the versions of _mpmath_add, _mpmath_div, and _mpmath_mult are correct.
They
agree with libmpf over a wide range of arguments. Can you check if there is a
significant performance advantage, especially for _mpmath_add?
I still need to change _mpmath_trim to return a Python integer for the exponent
and
then I'll do _mpmath_sqrt.
Original comment by casevh
on 12 Oct 2009 at 7:23
Hmm, while there is a performance advantage, I think now that it would be
better to
not bother changing the low-level function interface in mpmath (although we
could
still use these implementations to speed them up slightly), and instead aim for
implementing mpf and mpc base classes in gmpy.
Do you think this sounds realistic? We'd need a context class (basically just
holding
the precision), mpf and mpc classes implementing basic arithmetic. They'd
basically
need to be subclassable, and then coercions and everything else messy could
fall back
to Python code. (Though coercions mpf <-> mpc and with int/long and
float/complex
should be implemented directly in C for speed.)
(Side note: with classes done in C, it's probably better to just use mpz for
both
mantissa and exponent.)
The mp4 branch in mpmath trunk implements two separate backends -- one based on
the
proposed new format, and one based on Python floats / complexes. I should fix
it so
it also supports an old-format backend. That code should provide a basis for
implementing a gmpy-based backend.
Unfortunately, I'm a little preoccupied at the moment (mostly due to school)...
I'll
see when I get time to work more on this. I greatly appreciate your work so far.
BTW, I've been experimenting some with fixed-point implementations of
transcendental
functions on top of GMP/MPIR and it seems that it's possible to beat MPFR by a
factor
2-4 for elementary functions at low precisions. There will of course be a little
additional overhead when adding rounding and floating-point conversions.
I might set up a separate project for this. If it goes anywhere, gmpy could
include
it (it'd be vanilla C depending only on GMP/MPIR, and hopefully relatively
small) and
use it as a base library for floating-point transcendental functions.
Original comment by fredrik....@gmail.com
on 12 Oct 2009 at 3:01
I'm still trying determine what I want to do next (besides getting 1.10
officially
released.) I've had a couple of requests to support MPFR. I also want to support
mpmath. But I don't know how much time I'll have.
My primary goal for 1.20 is to replace the GMP/MPIR mpf layer with a better
behaved
floating point type. These routines are a good start and I'll look at using
your new
format. I want to maintain compatibility with prior gmpy releases if I can.
If I do support MPFR, I'll probably call it gmpy2 so I can break backwards
compatibility. I've even thought of making it Python 3.x only. But I am worried
about
increasing the dependencies and I don't know if license for MPFR will change to
LGPL
3 in the future.
I like the idea of your fixed-point library. What would be the maximum practical
precision?
Original comment by casevh
on 13 Oct 2009 at 5:16
> But I am worried about increasing the dependencies and I don't know if
license for
MPFR will change to LGPL 3 in the future.
They already switched to LGPL 3 for the next release, afaik.
> I like the idea of your fixed-point library. What would be the maximum
practical
precision?
Should work fine for any precision as long as shift offsets fit in ints, like
the
current mpmath fixed-point code. Though I'd probably concentrate on low
precision
optimization first.
Original comment by fredrik....@gmail.com
on 13 Oct 2009 at 9:37
_mpmath_trim now converts the exponent back to a Python integer. I'm planning to
release gmpy 1.10 in the next week or two. If you have time to test, that would
be great.
Original comment by casevh
on 18 Oct 2009 at 7:50
Original issue reported on code.google.com by
fredrik....@gmail.com
on 28 Aug 2009 at 12:20