#4446: New module complex_mpc using lib mpc for complex multiprecision
arithmetic
-----------------------------------------------------------------+----------
Reporter: thevenyp | Owner:
mabshoff
Type: enhancement | Status:
positive_review
Priority: major | Milestone:
sage-4.6
Component: basic arithmetic | Keywords:
Author: Philippe Theveny, Alex Ghitza, Yann Laigle-Chapuy | Upstream:
N/A
Reviewer: William Stein, David Kirkby, Paul Zimmermann | Merged:
Work_issues: |
-----------------------------------------------------------------+----------
Changes (by zimmerma):
* status: needs_review => positive_review
* reviewer: William Stein, David Kirkby => William Stein, David Kirkby,
Paul Zimmermann
Comment:
Replying to [comment:52 was]:
> I'm "out there". Is it necessary and not slower than what is in Sage
already? I'm willing to take your (=Paul's) word for it, assuming you say
you've done some benchmarks. Also, rounding might be (vastly) better in
mpc, which would be another point in its favor.
William, here are some results of benchmarks I've done on my Core 2 Duo
laptop, with Sage 4.5.3,
with Yann latest patch and the associate spkg. Efficiency:
{{{
sage: K = MPComplexField(53)
sage: a = K(3.3902384+1.23456789*i)
sage: b = CC(3.3902384+1.23456789*i)
sage: timeit('a*a')
625 loops, best of 3: 1.03 µs per loop
sage: timeit('b*b')
625 loops, best of 3: 1.08 µs per loop # CC is slower
sage: timeit('a+a')
625 loops, best of 3: 566 ns per loop
sage: timeit('b+b')
625 loops, best of 3: 571 ns per loop # CC is slower
sage: K=MPComplexField(1000)
sage: a = K(3.3902384,9203483)
sage: b = ComplexField(1000)(3.3902384,9203483)
sage: timeit('a*a')
625 loops, best of 3: 3.26 µs per loop
sage: timeit('b*b')
625 loops, best of 3: 2.99 µs per loop # CC is slightly faster
sage: timeit('a+a')
625 loops, best of 3: 680 ns per loop
sage: timeit('b+b')
625 loops, best of 3: 694 ns per loop # CC is slower
sage: timeit('a.sin()')
625 loops, best of 3: 294 µs per loop
sage: timeit('b.sin()')
625 loops, best of 3: 343 µs per loop # CC is slower
sage: timeit('a+17')
625 loops, best of 3: 2.56 µs per loop
sage: timeit('b+17')
625 loops, best of 3: 3.02 µs per loop # CC is slower
sage: a2=a+17
sage: b2=b+17
sage: timeit('a/a2')
625 loops, best of 3: 30.3 µs per loop
sage: timeit('b/b2')
625 loops, best of 3: 14.6 µs per loop
sage: timeit('a.real()')
625 loops, best of 3: 854 ns per loop
sage: timeit('b.real()')
625 loops, best of 3: 5.13 µs per loop
sage: timeit('a.conjugate()')
625 loops, best of 3: 736 ns per loop
sage: timeit('b.conjugate()')
625 loops, best of 3: 750 ns per loop
sage: timeit('a.argument()')
625 loops, best of 3: 146 µs per loop
sage: timeit('b.argument()')
625 loops, best of 3: 155 µs per loop
sage: timeit('a**12345')
625 loops, best of 3: 126 µs per loop
sage: timeit('b**12345')
625 loops, best of 3: 176 µs per loop
sage: timeit('a.log()')
625 loops, best of 3: 315 µs per loop
sage: timeit('b.log()')
625 loops, best of 3: 400 µs per loop
sage: timeit('a.exp()')
625 loops, best of 3: 221 µs per loop
sage: timeit('b.exp()')
625 loops, best of 3: 211 µs per loop
}}}
For the accuracy, I've computed the maximal relative error on the real or
imaginary part after
1000 random tries for several operations. My program also prints the
inputs that give the
corresponding maximal error:
{{{
multiplication:
sage: foo(53) # 53 bits
MPC: 1.09250136812570e-16 (-0.199676992382071 - 0.194223361415138*I,
0.9413020\
06592817 + 0.645200777306369*I)
CC: 1.32263143752823e-14 (0.834168924383086 - 0.525629438090125*I,
-0.4766533\
42074290 + 0.752454804563813*I)
sage: foo(100) # 100 bits
MPC: 7.83034842148009e-31 (-0.43034701147206399980669837638 -
0.81582963965013\
399110093432939*I, -0.83969050263610625199834532314 +
0.17390362688794167984134\
423048*I)
CC: 9.28509698731955e-29 (-0.36479911591428496460512307920 +
0.22548928494028\
018164534780632*I, -0.49018610148599010061638163259 -
0.30229330441199619133786\
980680*I)
division:
sage: foo(53)
MPC: 1.08485488717363e-16 (0.906903106647522 + 0.991328300040510*I,
0.26885465\
0418765 - 0.188551687646414*I)
CC: 4.27558190140949e-13 (-0.947451505768202 - 0.424071772860430*I,
-0.499313\
865619141 - 0.223444662259920*I)
sage: foo(100)
MPC: 7.66357970552281e-31 (0.92755016191057837385980632404 -
0.004266231128213\
8322794515257916*I, 0.14072961243734751540073515895 -
0.21101836467335963334809\
014965*I)
CC: 1.06937889944401e-27 (-0.41008840630275384495897016614 +
0.43420284734725\
368743056683764*I, -0.61829368333881933443455346326 -
0.58354470182146315365872\
353144*I)
sqrt:
sage: foo(53)
MPC: 1.09880191814582e-16 (-0.621894659027438 - 0.0246731753089378*I,
0.823805\
504837539 - 0.981310087843700*I)
CC: 2.08688816342048e-16 (0.968316682679636 + 0.637686286592402*I,
-0.2999297\
58140749 + 0.680399103485303*I)
sage: foo(100)
MPC: 7.83762497864303e-31 (0.87484855451536280543412317489 +
0.744936094843967\
94278935970162*I, -0.93716227305934196720704783450 -
0.765294741452728993051615\
69068*I)
CC: 1.59191232901277e-30 (0.22724304905610094307701469163 +
0.290884560229311\
05826353325335*I, -0.44728415428244253958062250934 +
0.959404399673905010209671\
68949*I)
exp:
sage: foo(53)
MPC: 1.08513702552684e-16 (0.801188951612012 - 0.440834327392402*I,
-0.6487549\
15945509 - 0.558408459317570*I)
CC: 2.56306731927053e-16 (-0.534146492930164 - 0.126305922104917*I,
-0.525552\
199671560 + 0.0746660767007601*I)
sage: foo(100)
MPC: 7.85000220770852e-31 (-0.25871533881495460870167038331 +
0.70802010475507\
807939080332669*I, 0.99833000148645085071979336755 +
0.103046934118108408058895\
43882*I)
CC: 2.02731881365962e-30 (0.033507794108965175610485564405 +
0.53119476356929\
639125934355035*I, 0.46545588987697950055582059300 +
0.896373882394705924442797\
08331*I)
log:
sage: foo(53)
MPC: 1.09778133774750e-16 (-0.391618732071762 + 0.824073966701766*I,
0.3486864\
59442427 - 0.191201122946134*I)
CC: 6.10819034602284e-14 (0.306684274970498 + 0.953229884394732*I,
0.69862405\
1151183 + 0.276699733597418*I)
sage: foo(100)
MPC: 7.78374256155462e-31 (0.50154056740914466911692379172 +
0.796641256036122\
24316946320184*I, -0.76918900642484651661867461393 -
0.286514343536794358035934\
91438*I)
CC: 6.82825438431515e-28 (-0.13721125119514449776575251187 +
0.99010163436560\
678238821778034*I, 0.92169332492985590698251634723 +
0.433132422276759879534575\
89642*I)
sin:
sage: foo(53)
MPC: 1.08720457020742e-16 (-0.369315245123802 + 0.854572675318644*I,
0.8169569\
42853907 - 0.383919250130617*I)
CC: 2.59601650874087e-16 (-0.546197856444743 - 0.534628610514192*I,
0.3036961\
63087668 - 0.231790426712382*I)
sage: foo(100)
MPC: 7.55871247471999e-31 (0.20947476591320905500539204843 -
0.626588978693616\
08926271268655*I, -0.55364821811403565597840750754 +
0.743341403117384889712851\
11372*I)
CC: 1.72458459418951e-30 (0.42758001678856166974028750611 -
0.923412131000390\
51795327859750*I, -0.99149234221470907346167026703 -
0.690427397584128136419461\
37603*I)
cos:
sage: foo(53)
MPC: 1.08584899379431e-16 (0.0646782540849349 - 0.111085748624151*I,
0.4701880\
86702821 + 0.252843262283133*I)
CC: 2.85959441825441e-16 (0.0850486257569383 - 0.330860565935531*I,
0.0564130\
612472886 + 0.0847664237260943*I)
sage: foo(100)
MPC: 7.78394156732906e-31 (0.93946166432029710121550615978 -
0.592139906105641\
30050851821134*I, 0.53601265027477876871715914306 -
0.1322166340317310632649551\
5769*I)
CC: 1.94260689254542e-30 (-0.55385187621901359297440670364 +
0.50415610417338\
788296133563973*I, -0.13371850315082336799522666801 -
0.96241946806142646000989\
112726*I)
}}}
I used the following programs:
{{{
def check(res,ref):
if ref.real() == 0.0:
if res.real() <> 0.0:
err_real = +Infinity
else:
err_real = 0.0
else:
err_real = RR((res.real()-ref.real()).abs()/ref.real().abs())
if ref.imag() == 0.0:
if res.imag() <> 0.0:
err_imag = +Infinity
else:
err_imag = 0.0
else:
err_imag = RR((res.imag()-ref.imag()).abs()/ref.imag().abs())
return max(err_real,err_imag)
def foo(p):
K=MPComplexField(p)
L=ComplexField(p)
K2=MPComplexField(p+20)
L2=ComplexField(p+20)
max_erra = 0
max_errb = 0
ina = 0
inb = 0
for i in range(10^3):
b = L.random_element()
br = b.real()
bi = b.imag()
a = K(b)
c = L.random_element()
cr = b.real()
ci = b.imag()
d = K(c)
aa=cos(a) # change to a*d for multiplication
bb=cos(b) # change to b*c for multiplication
ref = cos(K2(a)) # change to K2(a)*K2(d) for multiplication
erra = check(K2(aa),ref)
errb = check(K2(bb),ref)
if erra > max_erra:
max_erra = erra
ina = a, d
if errb > max_errb:
max_errb = errb
inb = b, c
print "MPC: ", max_erra, ina
print "CC: ", max_errb, inb
}}}
The conclusion is that MPC is in the majority of the cases faster than
ComplexField (the only
exception being the division) and always more accurate, with a ratio of a
factor up to
3941 (i.e., about 12 bits) for the 53-bit division.
For MPC we observe that the maximal relative error is --- as guaranteed by
MPC --- always less than 1/2 ulp (unit in last place) which corresponds to
a relative error of 2**(-p), i.e., 1.12e-16 for p=53 and 7.89e-31 for
p=100.
I thus give a positive review to Yann's patch.
Paul
--
Ticket URL: <http://trac.sagemath.org/sage_trac/ticket/4446#comment:53>
Sage <http://www.sagemath.org>
Sage: Creating a Viable Open Source Alternative to Magma, Maple, Mathematica,
and MATLAB
--
You received this message because you are subscribed to the Google Groups
"sage-trac" group.
To post to this group, send email to [email protected].
To unsubscribe from this group, send email to
[email protected].
For more options, visit this group at
http://groups.google.com/group/sage-trac?hl=en.
