#19964 closed enhancement (fixed)
tight complex interval inverse
Reported by:  mmarco  Owned by:  

Priority:  major  Milestone:  sage7.2 
Component:  numerical  Keywords:  interval, root 
Cc:  tscrim, vdelecroix, dkrenn, jdemeyer, tmonteil, cheuberg  Merged in:  
Authors:  Miguel Marco  Reviewers:  Vincent Delecroix, Marc Mezzarobba 
Report Upstream:  N/A  Work issues:  
Branch:  39558b6 (Commits)  Commit:  
Dependencies:  Stopgaps: 
Description (last modified by )
The current method to compute the multiplicative inverse of a complex interval uses schoolbok formula that produces intervals much bigger than the actual result.
This ticket implements a method that produces a tight interval enclosure of the result.
As a result, we get better precission in root isolating methods (and hence, less steps needed).
Here are some examples that illustrate the impact of the changes:
Old behaviour:
sage: a = CIF((1, 2), (3, 4)) sage: a.real().lower(), a.real().upper() (1.00000000000000, 2.00000000000000) sage: a = CIF((1, 2), (3, 4)) sage: b = a.__invert__() sage: b.real().lower(), b.real().upper() (0.0499999999999999, 0.200000000000001) sage: b.imag().lower(), b.imag().upper() (0.400000000000001, 0.149999999999999) sage: b.diameter() 1.20000000000000 sage: %time a.__invert__() CPU times: user 0 ns, sys: 0 ns, total: 0 ns Wall time: 23.1 µs 1.?  1.?*I
new behaviour:
sage: a = CIF((1, 2), (3, 4)) sage: a.real().lower(), a.real().upper() (1.00000000000000, 2.00000000000000) sage: b = a.__invert__() sage: b.real().lower(), b.real().upper() (0.0588235294117647, 0.153846153846154) sage: b.imag().lower(), b.imag().upper() (0.300000000000001, 0.200000000000000) sage: b.diameter() 0.893617021276596 sage: %time a.__invert__() CPU times: user 0 ns, sys: 0 ns, total: 0 ns Wall time: 19.1 µs 0.1?  0.3?*I
Another example with a smaller interval,
old:
sage: a = CIF((5, 5.01), (7, 7.01)) sage: b = a.__invert__() sage: b.diameter() 0.00523867991752868 sage: b.real().lower(), b.real().upper() (0.0673489564952680, 0.0677027027027028) sage: b.imag().lower(), b.imag().upper() (0.0947297297297298, 0.0942885390933752) sage: %time a.__invert__() CPU times: user 0 ns, sys: 0 ns, total: 0 ns Wall time: 23.1 µs 0.068?  0.095?*I
new:
sage: a = CIF((5, 5.01), (7, 7.01)) sage: b = a.__invert__() sage: b.diameter() 0.00253766599329326 sage: b.real().lower(), b.real().upper() (0.0674398874563158, 0.0676112447891434) sage: b.imag().lower(), b.imag().upper() (0.0945945945945946, 0.0944232370063658) sage: a = CIF((5.01, 5), (7, 7.01)) sage: %time a.__invert__() CPU times: user 0 ns, sys: 0 ns, total: 0 ns Wall time: 49.1 µs 0.068?  0.0945?*I
As you can see, the diameter of the result can easily get cut to half or even smaller.
The timings of the inversion can vary, but they remain in the same order of magnitude as before. It has a noticeable effect in root isolation:
old:
sage: R.<x> = QQ[] sage: f = 7/6*x^7  x^6 + 1/2*x^4 + 2/17*x^3 + 2*x^2  6*x  9 sage: %time f.roots(CC) CPU times: user 6 ms, sys: 0 ns, total: 6 ms Wall time: 5.34 ms [(1.03157268729039, 1), (1.18964736821330  0.850802715570186*I, 1), (1.18964736821330 + 0.850802715570186*I, 1), (0.0780792982605128  1.39288589462175*I, 1), (0.0780792982605128 + 1.39288589462175*I, 1), (1.19878298502655  0.599304323571307*I, 1), (1.19878298502655 + 0.599304323571307*I, 1)] sage: %time f.roots(QQbar) CPU times: user 17 ms, sys: 0 ns, total: 17 ms Wall time: 15.7 ms [(1.031572687290387?, 1), (1.189647368213303?  0.8508027155701860?*I, 1), (1.189647368213303? + 0.8508027155701860?*I, 1), (0.07807929826051275?  1.392885894621755?*I, 1), (0.07807929826051275? + 1.392885894621755?*I, 1), (1.198782985026555?  0.5993043235713073?*I, 1), (1.198782985026555? + 0.5993043235713073?*I, 1)]
new:
sage: R.<x>=QQ[] sage: f = 7/6*x^7  x^6 + 1/2*x^4 + 2/17*x^3 + 2*x^2  6*x  9 sage: %time f.roots(CC) CPU times: user 4 ms, sys: 0 ns, total: 4 ms Wall time: 3.28 ms [(1.03157268729039, 1), (1.18964736821330  0.850802715570186*I, 1), (1.18964736821330 + 0.850802715570186*I, 1), (0.0780792982605128  1.39288589462175*I, 1), (0.0780792982605128 + 1.39288589462175*I, 1), (1.19878298502655  0.599304323571307*I, 1), (1.19878298502655 + 0.599304323571307*I, 1)] sage: %time f.roots(QQbar) CPU times: user 9 ms, sys: 1 ms, total: 10 ms Wall time: 9.3 ms [(1.031572687290387?, 1), (1.189647368213303?  0.8508027155701860?*I, 1), (1.189647368213303? + 0.8508027155701860?*I, 1), (0.07807929826051275?  1.392885894621755?*I, 1), (0.07807929826051275? + 1.392885894621755?*I, 1), (1.198782985026555?  0.5993043235713073?*I, 1), (1.198782985026555? + 0.5993043235713073?*I, 1)]
Change History (19)
comment:1 Changed 5 years ago by
 Branch set to u/mmarco/tight_complex_interval_inverse
comment:2 Changed 5 years ago by
 Commit set to 97f4504c20b5bfbf315f5e00013dab27df25fbd3
comment:3 Changed 5 years ago by
 Cc tscrim vdelecroix dkren jdemeyern added
 Component changed from PLEASE CHANGE to numerical
 Description modified (diff)
 Keywords interval root added
 Status changed from new to needs_review
 Type changed from PLEASE CHANGE to enhancement
comment:4 Changed 5 years ago by
 Cc jdemeyer added; jdemeyern removed
comment:5 Changed 5 years ago by
 Cc dkrenn added; dkren removed
Can you give an example which clearly illustrates the better bounds? Also, have you run timings versus the schoolbook version? (I don't have the technical abilities to review this ticket though.)
comment:6 Changed 5 years ago by
 Description modified (diff)
Added some examples to the description.
comment:7 Changed 5 years ago by
Is this better in all cases?
 0.?e182  0.7598602580415435?*I,  0.?e249 + 0.7598602580415435?*I, + 0.?e215  0.7598602580415435?*I, + 0.?e229 + 0.7598602580415435?*I,
comment:8 Changed 5 years ago by
It should be, up to floating point artifacts. I don't know much about the subtleties of floating point numbers though.
comment:9 Changed 5 years ago by
I was wondering about the e249
becoming a e229
in the above example...
comment:10 Changed 5 years ago by
Yes, that is what i imagined, but i would guess that in those two cases the discrepancy is due to those pesky numerical artifacts that can appear just by, for example, performing the same operations in different order.
comment:11 Changed 5 years ago by
 Cc tmonteil cheuberg added
comment:12 Changed 5 years ago by
IMHO the tests should cover all branches of the new implementation (ideally with test cases where the rounding modes matter, if that's not too hard to achieve).
comment:13 Changed 5 years ago by
 Branch changed from u/mmarco/tight_complex_interval_inverse to u/vdelecroix/19964
 Commit changed from 97f4504c20b5bfbf315f5e00013dab27df25fbd3 to d57c62632a612d8f098515ff604ba153a3566ec2
 Milestone changed from sage7.1 to sage7.2
 Reviewers set to Vincent Delecroix
Following comment:12 I added a doctest. And I also removed the trailing whitespaces in a second commit.
This is good for me to go. Any remarks?
New commits:
1e5fe26  Trac 19964: add doctest for all branches

d57c626  Trac 19964: get rid of trailing whitespaces

comment:14 Changed 5 years ago by
 Commit changed from d57c62632a612d8f098515ff604ba153a3566ec2 to d57fee2259bf69d0f0471f1bc89c3a81f1e12e27
comment:15 Changed 5 years ago by
lgtm except that I'd suggest the following change in order to really test all possible cases:
@@ 1130,12 +1130,15 @@ cdef class ComplexIntervalFieldElement(sage.structure.element.FieldElement): Check that the code is valid in all regions of the complex plane:: sage: rats = [15423/906, 337/59976, 145151/145112]  sage: for a in rats:  ....: for b in rats:  ....: x = CIF(a, b)  ....: assert (x * (~x)  1).contains_zero()  ....: x = CIF(a,b)  ....: assert (x * (~x)  1).contains_zero() + sage: ivs = [RIF(aa, bb) for a in rats for b in rats + ....: for aa in [a, a] for bb in [b, b] + ....: if aa <= bb] + sage: for x in ivs: + ....: for y in ivs: + ....: z = CIF(x, y) + ....: inv = ~z + ....: assert ((z*inv  1).contains_zero() + ....: or z.contains_zero() and inv.real().is_NaN())
comment:16 Changed 5 years ago by
 Commit changed from d57fee2259bf69d0f0471f1bc89c3a81f1e12e27 to 39558b6299664e10bd008e2f97bb9c64e62096b3
Branch pushed to git repo; I updated commit sha1. New commits:
39558b6  Trac 19964: extend doctest

comment:17 Changed 5 years ago by
 Reviewers changed from Vincent Delecroix to Vincent Delecroix, Marc Mezzarobba
 Status changed from needs_review to positive_review
You really don't want to test the cases where the interval contains zero ;)
. But so be it.
comment:18 Changed 5 years ago by
 Branch changed from u/vdelecroix/19964 to 39558b6299664e10bd008e2f97bb9c64e62096b3
 Resolution set to fixed
 Status changed from positive_review to closed
comment:19 Changed 4 years ago by
 Commit 39558b6299664e10bd008e2f97bb9c64e62096b3 deleted
Breakage: #22834
Branch pushed to git repo; I updated commit sha1. New commits:
compute tightly the inverse of complex interval.
adapted division to new inverse
fixed doctests
Merge branch 'inversecomplexintervals' into t/19964/tight_complex_interval_inverse