On Dec 3, 2007 5:40 AM, Bill Hart <[EMAIL PROTECTED]> wrote:
>
> I've just been looking at SAGE ticket number 173:
>
> http://www.sagemath.org:9002/sage_trac/ticket/173
>
> The idea is that Mathematica raises a 3 dimensional matrix M over QQ
> to the power 20,000 much faster than either SAGE or Magma.
>
> I don't know any algorithm for doing this efficiently. I only know one
> algorithm:
>
> 1) Compute the characteristic polynomial p(x) of M (time 0.00s)
> 2) Compute x^20000 mod p (time 0.22s)
> 3) Substitute M into the result (time 0.00s)
>
> It's pretty obvious where the time is going here - polynomial
> arithmetic. I guess this is the algorithm being used. Is Pari or NTL
> being used for the polynomial expmod?
Where are you getting the above timings from?
I tried the following just now on sage.math:
sage: m = matrix(QQ,3,range(9)); m[0,0] =20; n=m^(-1)
sage: time k = n^20000 # Sage is now 4 times faster than when
#173 was posted!
CPU times: user 0.16 s, sys: 0.00 s, total: 0.16 s
Wall time: 0.16
sage: f = n.charpoly()
sage: K.<a> = NumberField(f)
sage: time b=a^20000
CPU times: user 0.16 s, sys: 0.00 s, total: 0.16 s
It would be better to rescale f, do the arithmetic there, then scale back.
sage: x = f.parent().gen()
sage: g = f(x/(2*3*5))*27000
sage: K.<a> = NumberField(g)
sage: (a/30).minpoly()
x^3 + 37/10*x^2 - 8/15*x + 1/60
sage: n.charpoly()
x^3 + 37/10*x^2 - 8/15*x + 1/60
sage: c = a/30
sage: time b = c^20000
CPU times: user 0.04 s, sys: 0.00 s, total: 0.04 s
Wall time: 0.04
I think this is using NTL for underlying arithmetic. I think the exponentiation
is just using a completely generic binary powering algorithm (i.e.,
object creation
overhead, not using only NTL, etc.) So the power could be sped up a lot.
Here's Mathematica again (I'm having tons of trouble using its Timing command,
so just wall time):
sage: nn = mathematica(n)
sage: time k=nn^20000
CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
Wall time: 0.01
sage: time k=nn^200000
CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
Wall time: 0.27
sage: time k=nn^2000000
CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
Wall time: 4.53
I think the nn^200000 and n^2000000 are probably good benchmarks to shoot for.
For comparison, in Sage:
sage: time k = n^200000
CPU times: user 4.36 s, sys: 0.05 s, total: 4.41 s
Wall time: 4.41
...
Too long to do 2000000.
Regarding using number field arithmetic:
sage: sage: time b = c^200000
CPU times: user 1.11 s, sys: 0.01 s, total: 1.12 s
Wall time: 1.13
sage: sage: time b = c^2000000
CPU times: user 26.85 s, sys: 0.91 s, total: 27.75 s
Wall time: 27.75
...
By the way, PARI is *not* competitive here:
sage: nn = gp(n)
sage: gp.eval('gettime; a = %s^20000; gettime/1000.0'%nn.name())
'0.42800000000000000000000000000000000000'
sage: gp.eval('gettime; a = %s^200000; gettime/1000.0'%nn.name())
'14.084000000000000000000000000000000000'
Magma is not competitive at all here either (Sage is already nearly
3 times faster):
sage: mm = magma(n)
sage: magma.eval('time a := %s^20000;'%mm.name())
'Time: 0.390'
sage: magma.eval('time a := %s^200000;'%mm.name())
'Time: 18.090'
I want to time maple, but I just spent 10 minutes and couldn't
even figure out how to raise a matrix to a power!
So if we could get the number field powering to be about 7 times
faster in this case and implement the algorithm Bill is suggesting,
we could beat Mathematica, and make Sage the fastest
program in the world at this problem. Sage already is the second
fastest for this particular benchmark
>
> I reckon we can speed this up. What do people think?
>
> Bill.
> >
>
--
William Stein
Associate Professor of Mathematics
University of Washington
http://wstein.org
--~--~---------~--~----~------------~-------~--~----~
To post to this group, send email to sage-devel@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/
-~----------~----~----~----~------~----~------~--~---
