New module complex_mpc using lib mpc for complex multiprecision arithmetic

[thevenyp]

Herewith and there ( http://www.loria.fr/~thevenyp/complex_mpc.patch 38K) is a patch with new classes using the MPC library for complex multi-precision arithmetic (see ticket #4308 for the associated spackage).

This is an adaptation of the module real_mpfr and of ComplexField? and ComplexNumber? classes. It adds a class MPComplexField with precision (common to both part) and rounding modes (specific to each part) and a class MPComplexNumber.

This first attempt implements only the complex arithmetic.

The test suite does fail due to coercion problems I can't solve.

[AlexGhitza]

replace the previous patch (rebase against 3.2)

[AlexGhitza]

There seems to have been some bitrot due to recent changes in Sage. I have uploaded a very slightly modified version of Philippe's patch, which applies cleanly against 3.2. A number of doctests fail for various reasons. I will try to look into this soon.

[thevenyp]

coercion problems resolved, more functions

[thevenyp]

The coercion problems have been solved, all doctests now succeed in new version. The whole set of mpc functions involving complex numbers is now in the interface.

[AlexGhitza]

apply after complex_mpc.p0.patch

[AlexGhitza]

Philippe,

Great work! I'll do my best to review this for 3.2.2. Right now I notice that the patch doesn't work in 3.2.1 because of #4580. I'm attaching a tiny patch that fixes that (should be applied after complex_mpc.p0.patch.

[AlexGhitza]

There are a few small issues that I ran into, but since this is going to be (for now) an optional part of Sage and nothing depends on it, I think we can fix these later (I'm keeping a list of them so they don't get lost).

However, there is a policy of 100% doctest for all new code. The coverage is now:

[ghitza@artin sage]$ sage -coverage rings/complex_mpc.pyx
----------------------------------------------------------------------
rings/complex_mpc.pyx
SCORE rings/complex_mpc.pyx: 68% (46 of 67)

Missing documentation:
	 * _repr_ (self):
	 * _latex_(self):
	 * _an_element_(self):
	 * random_element(self, bound):
	 * __hash__(self):
	 * prec(self):
	 * _set(self, z, y=None, base=10):
	 * _repr_(self):
	 * _latex_(self):
	 * ModuleElement _add_(self, ModuleElement right):
	 * ModuleElement _sub_(self, ModuleElement right):
	 * RingElement _mul_(self, RingElement right):
	 * RingElement _div_(self, RingElement right):
	 * ModuleElement _neg_(self):
	 * __create__MPComplexField_version0 (prec, rnd):
	 * __create_MPComplexNumber_version0 (parent, s, base=10):


Missing doctests:
	 * _rnd(self):
	 * _rnd_re(self):
	 * _rnd_im(self):
	 * _real_field(self):
	 * _imag_field(self):


Possibly wrong (function name doesn't occur in doctests):
	 * _element_constructor_(self, z):
	 * _coerce_map_from_(self, S):
	 * bint is_exact(self) except -2: return False def is_finite(self):
	 * Element _call_(self, z):
	 * Element _call_(self, x):

----------------------------------------------------------------------

So the missing docstrings and doctests will have to be added before this patch can be merged. I would happily do this but the earliest I can get to it is Thursday or so (and that's being optimistic).

To end on a positive note: this looks good. It will take a bit more work, but if we can show that (a) the performance is improved, or doesn't suffer too much, and (b) the switch can be made seamlessly, then it will be easy to convince people that we should switch the core of our complex numbers functionality over to MPC.

[AlexGhitza]

apply after complex_mpc.p0.patch and trac4446_fix.patch

[AlexGhitza]

I added a patch trac4446_doctests.patch, which does a number of things:

• adds doctests for all functions except three internal use only functions
• changes _repr_ of complex numbers so that it agrees with the way complex numbers are currently printed in Sage
• makes MPComplexField inherit from ParentWithGens?, being generated over its real field by the square root of -1 (just as it is now); this required adding a few functions. So now one can do

sage: from sage.rings.complex_mpc import MPComplexField
sage: MPC.<j> = MPComplexField()
sage: j^2
-1.00000000000000 + 0.000000000000000*I

Note that only the last three patches should be applied, in order: complex_mpc.p0.patch, trac4446_fix.patch, and trac4446_doctests.patch

[thevenyp]

One more patch trac4446_norm.patch with:

* Alex Ghitza listed in authors list and copyright notice

* complex_mpc uses a copy of mpfr rounding mode list instead of its private one

* the abs and norm documentation has been improved.

[was]

REFEREE REPORT:

1. How can this be optional? If the patch is applied then this is included:

 	625	    Extension('sage.rings.complex_mpc', 
 	626	              sources = ['sage/rings/complex_mpc.pyx'], 
 	627	              libraries = ['mpc', 'mpfr', 'gmp']), \ 
 	628	

So the only way I can see this going in like this is if it is standard. Am I missing something (probably)?

2. Timings on Xeon 2.6Ghz (sage.math). Multiplication is faster with mpc, but addition is significantly *slower*:

sage: K = MPComplexField(53)
sage: a = K(3.3902384)
sage: b = CC(3.3902384)

sage: timeit('a*a')
625 loops, best of 3: 376 ns per loop
sage: timeit('b*b')
625 loops, best of 3: 466 ns per loop

sage: timeit('a+a')
625 loops, best of 3: 368 ns per loop
sage: timeit('b+b')
625 loops, best of 3: 304 ns per loop

3. Powering doesn't work for mpc but it does for the current CC:

sage: a**a 
boom

4. Default rounding on coercion (or something?!) is different between MPC and CC:

sage: MPComplexField(100)(3.3902384)
3.3902383999999998742680418218 + 0.00000000000000000000000000000*I
sage: ComplexField(100)(3.3902384)
3.3902384000000000000000000000

5. With higher precision and nontrivial imaginary part, the timings are *all* in favor of the existing ComplexField? already in Sage:

sage: K = MPComplexField(1000)
sage: a = K(3.3902384,9203483)
sage: b = ComplexField(1000)(3.3902384,9203483)
sage: a
3.39023839999999987426804182177875190973281860351562500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 + 9.20348300000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e6*I
sage: b
3.39023840000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 + 9.20348300000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e6*I
sage: timeit('a*a')
625 loops, best of 3: 2.27 µs per loop
sage: timeit('b*b')
625 loops, best of 3: 2.2 µs per loop
sage: timeit('a+a')
625 loops, best of 3: 515 ns per loop
sage: timeit('b+b')
625 loops, best of 3: 445 ns per loop
sage: timeit('a.sin()')
^P625 loops, best of 3: 221 µs per loop
sage: timeit('b.sin()')
625 loops, best of 3: 130 µs per loop

Thus until the above issues are addressed, I see no point in mpc going into sage.

[zimmerma]

Philippe Theveny, who had a 2-year contract with us, just left, thus he won't be able to continue working on this. In case somebody wants to pursue his work, he has left the attached files. You need mpc-0.6 available at  http://www.loria.fr/~thevenyp/mpc-0.6.spkg. According to Philippe, the performance was improved, but still not as good as the current CC (but more reliable due to correct rounding).

[ylchapuy]

With the updated patch, and the corresponding spkg available at  http://yann.laiglechapuy.net/spkg/mpc-0.8.1.p0.spkg almost every methods for ComplexNumbers? are defined for MPComplexNumbers to.

The performance is quite good: with a precision of 500bits,

==========
add
ComplexField: a+b 40000 loops, best of 5: 700 ns per loop
MPComplexField: m+n 40000 loops, best of 5: 755 ns per loop
==========
add int
ComplexField: a+17 40000 loops, best of 5: 2.75 µs per loop
MPComplexField: m+17 40000 loops, best of 5: 2.43 µs per loop
==========
mul
ComplexField: a*b 40000 loops, best of 5: 2.66 µs per loop
MPComplexField: m*n 40000 loops, best of 5: 2.68 µs per loop
==========
div
ComplexField: a/b 40000 loops, best of 5: 5.76 µs per loop
MPComplexField: m/n 40000 loops, best of 5: 10.7 µs per loop
==========
real
ComplexField: a.real() 40000 loops, best of 5: 1.7 µs per loop
MPComplexField: m.real() 40000 loops, best of 5: 764 ns per loop
==========
conj
ComplexField: a.conjugate() 40000 loops, best of 5: 756 ns per loop
MPComplexField: m.conjugate() 40000 loops, best of 5: 633 ns per loop
==========
arg
ComplexField: a.argument() 2000 loops, best of 5: 104 µs per loop
MPComplexField: m.argument() 2000 loops, best of 5: 102 µs per loop
==========
cos
ComplexField: a.cos() 2000 loops, best of 5: 57.2 µs per loop
MPComplexField: m.cos() 2000 loops, best of 5: 86.4 µs per loop
==========
pow
ComplexField: a**b 2000 loops, best of 5: 240 µs per loop
MPComplexField: m**n 2000 loops, best of 5: 515 µs per loop
==========
pow int
ComplexField: a**12345 2000 loops, best of 5: 55.1 µs per loop
MPComplexField: m**12345 2000 loops, best of 5: 549 µs per loop
==========
log
ComplexField: a.log() 4000 loops, best of 5: 177 µs per loop
MPComplexField: m.log() 4000 loops, best of 5: 153 µs per loop
==========
exp
ComplexField: a.exp() 10000 loops, best of 5: 53.7 µs per loop
MPComplexField: m.exp() 10000 loops, best of 5: 57.5 µs per loop

The only big slow down is for pow.

It also allows to compute various trigonometric functions without convertion to pari.

[ylchapuy]

Erratum: The only big slow down are for pow and div

[ylchapuy]

apply after mpc-0.8.1.patch

[ylchapuy]

I put it as needs review to know what has eventually to be done.

[drkirkby]

Has this been checked on Solaris? There have been issues with complex maths recently on Solaris, which were solved, but this would need testing.

There's general information about building on Solaris at

 http://wiki.sagemath.org/solaris

Information specifically for 't2' at

 http://wiki.sagemath.org/devel/Building-Sage-on-the-T5240-t2

Both the source (4.3.0.1 is the latest to build on Solaris) and a binary which will run on any SPARC can be found at http://www.sagemath.org/download-source.html

Dave

[enge]

The poor performance of "pow int" is due to the fact that binary exponentiation is not currently implemented in mpc; "pow int" just calls "pow", and is present as a convenience function. Binary exponentiation is on the todo list.

Except for special cases, "pow" uses log and exp, so the timing should essentially be the sum for these two functions. Is it possible that you hit a corner case with your experiments, or are the timings consistent over several arguments?

Similarly, what are the arguments for division?

Concerning solaris, mpc itself builds without problems, see

 http://www.multiprecision.org/index.php?prog=mpc&page=platforms

[zimmerma]

Except for special cases, "pow" uses log and exp, so the timing should essentially be the sum for these two functions. Is it possible that you hit a corner case with your experiments, or are the timings consistent over several arguments?

then we should get about 210us for powint instead of 549us. Indeed, it would be good to know the exact arguments, so that we can try to reproduce with vanilla MPC (maybe the problem is in the interface between Sage and MPC).

For division, GMP 5 has a much faster division, which might benefit to MPC too. But again, we could investigate more with a vanilla MPC, if we know the exact inputs used (with exact_rational).

[ylchapuy]

First, I added a small patch to rebase this on sage 4.3.3.

One need to apply:

• mpc-0.8.1.patch
• bug_literal.patch
• trac4446-rebase_sage4.3.3.patch

in this order; and of course install the spkg:

•  http://yann.laiglechapuy.net/spkg/mpc-0.8.1.p0.spkg

Then, here is an example with the exact input (I only put the last part concerning pow)

a.real().exact_rational()
-1558022726154386835635734102276225872461248356241845468520145243680772756232147166398716395960573602794592840676040928991781722240298082039889553874139/1636695303948070935006594848413799576108321023021532394741645684048066898202337277441635046162952078575443342063780035504608628272942696526664263794688
a.imag().exact_rational()
38743143824869548998159480051850490241205026500878413377278464034858416981434423320408166976481329585622042476431430168493334123089110626626927265425/102293456496754433437912178025862473506770063938845774671352855253004181137646079840102190385184504910965208878986252219038039267058918532916516487168
------------------------
a.log()
CC 5000 loops, best of 3: 170 µs per loop
MPC 5000 loops, best of 3: 143 µs per loop
------------------------
a.exp()
CC 5000 loops, best of 3: 45.6 µs per loop
MPC 5000 loops, best of 3: 49.1 µs per loop
------------------------
a**12345
CC 5000 loops, best of 3: 48.7 µs per loop
MPC 5000 loops, best of 3: 501 µs per loop

[ylchapuy]

and for div (same value for a)

b.real().exact_rational()
-1335661139895223805455191750131391965080983427340135233479028132544481051566385358748782698880424407717119644643310969153530607585226701160632690491359/1636695303948070935006594848413799576108321023021532394741645684048066898202337277441635046162952078575443342063780035504608628272942696526664263794688
b.imag().exact_rational()
-191023642554886415830676606502854964270751830339806573746974716999454958912475563080101457087388705931975286670167330679505998875813037449880206143447/409173825987017733751648712103449894027080255755383098685411421012016724550584319360408761540738019643860835515945008876152157068235674131666065948672
------------------------
a/b
CC 100000 loops, best of 3: 4.84 µs per loop
MPC 100000 loops, best of 3: 9.57 µs per loop

[zimmerma]

with the attached MPC program I get on a 2.83Ghz Core 2 under Fedora Core 2, with GMP 5.0.1 and MPFR 2.4.2

tarte% ./sage
mpc_log took 74us per loop
mpc_exp took 24us per loop
mpc_pow(12345) took 260us per loop
mpc_div took 3.4us per loop

and within Sage with ComplexField? I get (where I took the best of 3 %timeit runs):

sage: %timeit c=a.log()
625 loops, best of 3: 173 µs per loop
sage: %timeit c=a.exp()
625 loops, best of 3: 37.1 µs per loop
sage: %timeit c=a**12345
625 loops, best of 3: 26.6 µs per loop
sage: %timeit c=a/b
625 loops, best of 3: 2.71 µs per loop

thus we indeed need to work on powint, but the division is only 25% slower.

While doing this test I also noticed this strange thing with ComplexField?:

sage: smallb = CC(b)
sage: %timeit c=a/b
625 loops, best of 3: 2.69 µs per loop
sage: %timeit c=a/smallb
625 loops, best of 3: 290 µs per loop

thus dividing by a 53-bit complex is 100 times slower than dividing by a 500-bit complex! Should I open a separate ticket for that?

Did somebody check the 500-bit outputs of CC and MPC agree?

[zimmerma]

By studying what happens with the powint example, there was a failure in Ziv's strategy, thus we did compute two approximations, one with precision 511 bits (which was not enough), and one with a larger precision. By adding one more guard bit, the first approximation with 512 bits is enough to get correct rounding, and we now get 131us per loop instead of 260us (in vanilla MPC). A further gain could be obtained if we get rid of the mpc_log call to detect overflow or underflow: we would then get 113us. But of course the definite solution will be to implement binary exponentiation.