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:  sage5.13 
Component:  basic arithmetic  Keywords:  
Cc:  Paul Zimmermann  Merged in:  sage5.13.beta3 
Authors:  Jeroen Demeyer  Reviewers:  Paul Zimmermann 
Report Upstream:  N/A  Work issues:  
Branch:  Commit:  
Dependencies:  Stopgaps: 
Description (last modified by )
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)
Change History (22)
comment:1 Changed 9 years ago by
Description:  modified (diff) 

comment:2 Changed 9 years ago by
Cc:  Paul Zimmermann added 

Status:  new → needs_review 
comment:3 Changed 9 years ago by
comment:4 followup: 6 Changed 9 years ago by
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
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 Changed 9 years ago by
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
I'd like to know which command you used to get those results.
in devel/sage15337/sage
I ran ../../../sage tp 4 long *
, which seems equivalent to the 3rd form above.
Paul
comment:8 Changed 9 years ago by
All doctests pass with ./sage tp 4 long devel/sage15337/sage
.
I will now start my review.
Paul
comment:9 Changed 9 years ago by
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
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 followup: 12 Changed 9 years ago by
Status:  needs_review → needs_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) <ipythoninput222b55db8a739a> in <module>() > 1 a.exact_rational() /localdisk/tmp/sage5.12/local/lib/python2.7/sitepackages/sage/rings/real_mpfr.so in sage.rings.real_mpfr.RealNumber.exact_rational (sage/rings/real_mpfr.c:18776)() /localdisk/tmp/sage5.12/local/lib/python2.7/sitepackages/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 nonzero number, then the next number should be 2*a
. In other words, the difference between two floatingpoint numbers is always a multiple of a
(and thus always a floatingpoint number), which is not the case here:
sage: ba 0.000000000000000
Paul
comment:12 Changed 9 years ago by
Work issues:  segfault in fp_rank, error in nextabove → documentation 

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 nonzero number, then the next number should be2*a
.
In MPFR, this isn't true since MPFR doesn't have denormal numbers.
comment:13 Changed 9 years ago by
Status:  needs_work → needs_review 

Work issues:  documentation 
comment:14 Changed 9 years ago by
Description:  modified (diff) 

comment:16 Changed 9 years ago by
Description:  modified (diff) 

comment:17 Changed 9 years ago by
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.94065645841e324) sage: b, b.ulp() (4.94065645841e324, 0.0)
The same holds with RR:
sage: a=RR(0) sage: a, a.ulp() (0.000000000000000, 8.50969131174084e1388255822130839284) sage: b=a.ulp() sage: b, b.ulp() (8.50969131174084e1388255822130839284, 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
Reviewers:  → Paul Zimmermann 

Status:  needs_review → positive_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.94065645841e324) sage: b=a.ulp() sage: b, b.ulp() (4.94065645841e324, 4.94065645841e324)
and with RR:
sage: a=RR(0) sage: a, a.ulp() (0.000000000000000, 8.50969131174084e1388255822130839284) sage: b=a.ulp() sage: b, b.ulp() (8.50969131174084e1388255822130839284, 8.50969131174084e1388255822130839284)
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
Attachment:  15337_ulp.patch added 

comment:19 Changed 9 years ago by
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 692 692 +infinity 693 693 sage: (a).ulp() 694 694 +infinity 695 sage: a = R R('nan')695 sage: a = RDF('nan') 696 696 sage: a.ulp() is a 697 697 True 698 698
in the doctest for the method ulp()
of RDF
.
comment:20 Changed 9 years ago by
I made the additional change [...]
correct, sorry for not having seen this.
Paul
comment:21 Changed 9 years ago by
Merged in:  → sage5.13.beta3 

Resolution:  → fixed 
Status:  positive_review → closed 
now in vacation, I will look at it next week, unless someone beats me.
Paul