Opened 9 years ago

Closed 9 years ago

#15337 closed enhancement (fixed)

Speed up ulp() method of real_mpfr.pyx

Reported by: Jeroen Demeyer Owned by:
Priority: minor Milestone: sage-5.13
Component: basic arithmetic Keywords:
Cc: Paul Zimmermann Merged in: sage-5.13.beta3
Authors: Jeroen Demeyer Reviewers: Paul Zimmermann
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Status badges

Description (last modified by Jeroen Demeyer)

Before:

sage: a=RealField(10**6).one(); timeit('a.ulp()')
625 loops, best of 3: 42.9 µs per loop

After:

sage: a=RealField(10**6).one(); timeit('a.ulp()')
625 loops, best of 3: 6.13 µs per loop

The patch also adds a field argument, which can be used to result in a different RealField. It also adds a related epsilon() method, defined as a.epsilon() = a/2^a.precision().

Attachments (1)

15337_ulp.patch (14.3 KB) - added by Jeroen Demeyer 9 years ago.

Download all attachments as: .zip

Change History (22)

comment:1 Changed 9 years ago by Jeroen Demeyer

Description: modified (diff)

comment:2 Changed 9 years ago by Jeroen Demeyer

Cc: Paul Zimmermann added
Status: newneeds_review

comment:3 Changed 9 years ago by Paul Zimmermann

now in vacation, I will look at it next week, unless someone beats me.

Paul

comment:4 Changed 9 years ago by Paul Zimmermann

I tested the patch against Sage 5.12, and got the following failures

sage -t --long __init__.pyc  # AttributeError in doctesting framework
sage -t --long env.pyc  # AttributeError in doctesting framework
sage -t --long version.pyc  # AttributeError in doctesting framework
sage -t --long tests/interrupt.pyx  # Time out
sage -t --long tests/cmdline.py  # 11 doctests failed
sage -t --long calculus/desolvers.py  # 8 doctests failed
sage -t --long misc/interpreter.py  # 1 doctest failed
sage -t --long misc/trace.py  # 2 doctests failed

This is under Ubuntu 12.04. I can give more detail if needed.

Paul

comment:5 Changed 9 years ago by Jeroen Demeyer

I don't know how you ran the doctests, but it must have been in an invalid way because the errors clearly have nothing to do with the patch and you're testing .pyc files which are normally not tested (I bet you would get the same results without my patch).

So please run, from SAGE_ROOT something like

$ make ptestlong

or

$ ./sage -t --long --all

or

$ ./sage -t --long devel/sage/sage  # Skipping notebook and doc

comment:6 in reply to:  4 Changed 9 years ago by Jeroen Demeyer

Replying to zimmerma:

I can give more detail if needed.

I'd like to know which command you used to get those results.

comment:7 Changed 9 years ago by Paul Zimmermann

I'd like to know which command you used to get those results.

in devel/sage-15337/sage I ran ../../../sage -tp 4 --long *, which seems equivalent to the 3rd form above.

Paul

comment:8 Changed 9 years ago by Paul Zimmermann

All doctests pass with ./sage -tp 4 --long devel/sage-15337/sage. I will now start my review.

Paul

comment:9 Changed 9 years ago by Paul Zimmermann

I confirm the speedup for normal numbers. Before:

sage: a=RealField(10**6).one(); timeit('a.ulp()')
625 loops, best of 3: 118 µs per loop

After:

sage: a=RealField(10**6).one(); timeit('a.ulp()')
625 loops, best of 3: 7.74 µs per loop

However for +Inf I notice a slowdown. Before:

sage: a=RealField(10**6)("Inf"); timeit('a.ulp()')
625 loops, best of 3: 101 ns per loop

After:

sage: a=RealField(10**6)("Inf"); timeit('a.ulp()')
625 loops, best of 3: 330 ns per loop

For -Inf it is faster now, and I see no difference for NaN. Paul

comment:10 Changed 9 years ago by Jeroen Demeyer

The slowdown for +Inf was because we used to simply return self in that case. My code instead always creates a new number, so the slowdown you see is the overhead of creating a new number. However, I personally prefer the current code, which is simpler but not optimal for +Inf (which should be a rare case anyway). I don't think it is justified to make the code more complicated to deal with the case of +Inf in a special way.

comment:11 Changed 9 years ago by Paul Zimmermann

Status: needs_reviewneeds_work
Work issues: segfault in fp_rank, error in nextabove

"Adding or subtracting something less than half an ulp always gives the same number (unless the number is exactly a power of 2...)": the case of the opposite of a power of 2 is missing here, also this assumes the rounding mode is to nearest.

I don't see the logic of returning self for NaN but not for infinities.

Since the patch also adds tests for fp_rank (which I discovered) I tested it too:

sage: a=RR(0).ulp()
sage: a.exact_rational()
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-22-2b55db8a739a> in <module>()
----> 1 a.exact_rational()

/localdisk/tmp/sage-5.12/local/lib/python2.7/site-packages/sage/rings/real_mpfr.so in sage.rings.real_mpfr.RealNumber.exact_rational (sage/rings/real_mpfr.c:18776)()

/localdisk/tmp/sage-5.12/local/lib/python2.7/site-packages/sage/rings/integer.so in sage.rings.integer.Integer.__pow__ (sage/rings/integer.c:14022)()

RuntimeError: Segmentation fault

Finally this looks strange:

sage: a=RR(0).ulp()
sage: a.str(16)
'1.0000000000000@-1152921504606846976'
sage: b=a.nextabove()
sage: b.str(16)
'1.0000000000001@-1152921504606846976'

If a is the smallest positive non-zero number, then the next number should be 2*a. In other words, the difference between two floating-point numbers is always a multiple of a (and thus always a floating-point number), which is not the case here:

sage: b-a
0.000000000000000

Paul

comment:12 in reply to:  11 Changed 9 years ago by Jeroen Demeyer

Work issues: segfault in fp_rank, error in nextabovedocumentation

Replying to zimmerma:

"Adding or subtracting something less than half an ulp always gives the same number (unless the number is exactly a power of 2...)": the case of the opposite of a power of 2 is missing here, also this assumes the rounding mode is to nearest.

Fair enough, this documentation should be fixed.

I don't see the logic of returning self for NaN but not for infinities.

So what do you prefer: never returning self? returning self for NaN and +Inf but not for -Inf?

Since the patch also adds tests for fp_rank (which I discovered)

That's an indepedent bug deserving a new ticket: #15363.

If a is the smallest positive non-zero number, then the next number should be 2*a.

In MPFR, this isn't true since MPFR doesn't have denormal numbers.

comment:13 Changed 9 years ago by Jeroen Demeyer

Status: needs_workneeds_review
Work issues: documentation

comment:14 Changed 9 years ago by Jeroen Demeyer

Description: modified (diff)

comment:15 Changed 9 years ago by Paul Zimmermann

I will look at that one once I complete my review of #15363.

Paul

comment:16 Changed 9 years ago by Jeroen Demeyer

Description: modified (diff)

comment:17 Changed 9 years ago by Paul Zimmermann

So what do you prefer: never returning self? returning self for NaN and +Inf but not for -Inf?

I would prefer consistency: either always returning self, or never returning self (unless there is a good reason for not being coherent).

While playing with ulp I noticed something strange:

sage: a=RDF(0)
sage: a, a.ulp()
(0.0, 4.94065645841e-324)
sage: b, b.ulp()
(4.94065645841e-324, 0.0)

The same holds with RR:

sage: a=RR(0)
sage: a, a.ulp()
(0.000000000000000, 8.50969131174084e-1388255822130839284)
sage: b=a.ulp()
sage: b, b.ulp()
(8.50969131174084e-1388255822130839284, 0.000000000000000)

Finally:

In MPFR, this isn't true since MPFR doesn't have denormal numbers.

agreed.

Paul

comment:18 Changed 9 years ago by Paul Zimmermann

Reviewers: Paul Zimmermann
Status: needs_reviewpositive_review

I retract part of my comment 17 since it works correctly with the patch:

sage: a=RDF(0)
sage: a, a.ulp()
(0.0, 4.94065645841e-324)
sage: b=a.ulp()
sage: b, b.ulp()
(4.94065645841e-324, 4.94065645841e-324)

and with RR:

sage: a=RR(0)
sage: a, a.ulp()
(0.000000000000000, 8.50969131174084e-1388255822130839284)
sage: b=a.ulp()
sage: b, b.ulp()
(8.50969131174084e-1388255822130839284, 8.50969131174084e-1388255822130839284)

Thus the only remaining issue is whether we return self or a new object. Since this is a minor issue, I give a positive review.

Paul

Changed 9 years ago by Jeroen Demeyer

Attachment: 15337_ulp.patch added

comment:19 Changed 9 years ago by Jeroen Demeyer

I made the additional change

  • sage/rings/real_double.pyx

    diff --git a/sage/rings/real_double.pyx b/sage/rings/real_double.pyx
    a b  
    692692            +infinity
    693693            sage: (-a).ulp()
    694694            +infinity
    695             sage: a = RR('nan')
     695            sage: a = RDF('nan')
    696696            sage: a.ulp() is a
    697697            True
    698698

in the doctest for the method ulp() of RDF.

comment:20 Changed 9 years ago by Paul Zimmermann

I made the additional change [...]

correct, sorry for not having seen this.

Paul

comment:21 Changed 9 years ago by Jeroen Demeyer

Merged in: sage-5.13.beta3
Resolution: fixed
Status: positive_reviewclosed
Note: See TracTickets for help on using tickets.