Creating a viable free open source alternative to Maple, Mathematica, Matlab and Magma
William Stein, Associate Professor, University of Washington
BTW: I really, REALLY like sage. I'm just surprised i haven't heard of it before 5 days ago.
-- Clinton Bowen on the sage-support mailing list Mon, Feb 9, 2009 at 11:03 PM.
{{{id=2|
///
}}}
History
I started Sage at Harvard in January 2005.
Existing math software not good enough(free or commercial).
Sage-1.0 released February 2006 at Sage Days 1 (San Diego).
Sage Days Workshops 1, 2, ..., 11, at UCLA, UW, Cambridge, Bristol, Austin, France, San Diego, Seattle, etc.
Sage won first prize in Trophees du Libre (November 2007)
Funding from Microsoft, UW, NSF, DoD, Google, Sun, private donors, etc.
{{{id=117|
///
}}}
Sun's Interest in Sage
Sage will very soon (weeks) fully support Solaris. This porting work was fully funded by the Department of Defense, and will continue.
Sage is the only project whose goal is to create a viable open source alternative to all of Mathematica, Maple, Matlab and Magma. It is the analogue of Star Office or Firefox, as alternatives to Microsoft Office or Internet Explorer.
Free webapp -- sagenb.org -- has about 3000 users (and there are other servers at universities around the world..)
sagemath.org
Each day about 1500-2000 different real human visitors each day.
Sage Devel Headquarters: Four 24-core Sun X4450's with 128GB RAM each + 1 Sun X4540 with 24TB disk
sagenb.org
{{{id=161|
A. = QuaternionAlgebra(QQ, -1,-1)
A
///
Quaternion algebra with generators (i, j, k) over Rational Field
}}}
{{{id=162|
i*j + k
///
2*k
}}}
{{{id=62|
///
}}}
Many Advantages of Sage over the Ma*'s
Python-- a general purpose language at core of Sage; huge user base compared to Matlab, Mathematica, Maple and Magma
Cython -- write blazingly fast compiled code in Sage
Free and Open source
Excellent Peer review of Code: "I do really and truly believe that the Sage refereeing model results in better code." -- John Cremona
{{{id=119|
email('wstein@gmail.com', 'The calculation finished!')
///
Child process 50883 is sending email to wstein@gmail.com...
}}}
{{{id=118|
///
}}}
Sage Is...
Over 300,000 lines of new Python/Cython code
Distribution of mathematical software; easy to build from source (over 5 million lines of code).
}}}
{{{id=164|
A = random_matrix(QQ, 500); v = random_matrix(QQ,500,1)
time x = A \ v
///
Time: CPU 1.69 s, Wall: 2.20 s
}}}
{{{id=3|
len(str(x[0]))
///
1486
}}}
{{{id=32|
///
}}}
Example: A Huge Integer Determinant
{{{id=31|
a = random_matrix(ZZ,200,x=-2^127,y=2^127)
time d = a.determinant()
len(str(d))
///
Time: CPU 3.14 s, Wall: 4.25 s
7786
}}}
We can also copy this matrixover to Maple and compute the same determinant there...
{{{id=167|
a[0,0]
///
-10495856349494057617178227541504239545
}}}
{{{id=27|
maple.with_package('LinearAlgebra')
B = maple(a)
t = maple.cputime()
time c = B.Determinant()
maple.cputime(t)
///
21.969999999999999
}}}
{{{id=168|
c == d
///
True
}}}
This ability to easily move objects between math software is unique to Sage.
{{{id=82|
B = magma(a)
t = magma.cputime()
time e = B.Determinant()
magma.cputime(t)
///
Time: CPU 0.00 s, Wall: 11.36 s
10.69
}}}
{{{id=112|
///
}}}
IRC a few days ago...
22:41 < linuxgus> I can show people at work tomorrow that they don't have to abandon Matlab [if they switch to Sage]. Can Sage convert Matlab data types to its own? That would be awesome!
Sage also has 2d plotting that is almost identical to MATLAB:
{{{id=159|
import pylab as p
p.figure()
t = p.arange(0.01, 2.0, 0.01)
s = p.sin(2 * p.pi * t)
s = p.array([float(f(x)) for x in t])
P = p.plot(t, s, linewidth=4)
p.xlabel('time (s)'); p.ylabel('voltage (mV)')
p.title('Matlab-style plotting in Sage')
p.grid(True)
p.savefig('sage.png')
///
}}}
{{{id=110|
///
}}}
{{{id=73|
///
}}}
_fast_float_ yields super-fast evaluation of Sage symbolic expressions -- e.g., here it is 10 times faster than native Python!
{{{id=172|
f(x,y,z) = x^3 * sin(x^2) + cos(x*y) - 1/(x+y+z)
///
}}}
{{{id=74|
g = f._fast_float_(x,y,z)
timeit('g(4.5r,3.2r,5.7r)')
///
625 loops, best of 3: 709 ns per loop
}}}
{{{id=75|
%python
import math
def g(x): return math.sin(3*x)*x + log(x) + 1/(1+x)**2
///
}}}
{{{id=77|
timeit('g(4.5r)')
///
625 loops, best of 3: 6.89 µs per loop
}}}
{{{id=48|
///
}}}
{{{id=76|
///
}}}
Example: Compare Answers with Maple
{{{id=83|
var('x')
f = sin(3*x)*x+log(x) + 1/(x+1)^2
show(integrate(f))
///
This illustrates pylab (matplotlib + numpy), Sage plotting, html output, and @interact.
{{{id=97|
# first just play
import pylab
A = pylab.imread(DATA + 'seattle.png')
graphics_array([matrix_plot(A^(-1)), matrix_plot(1-A[0:,0:,2])]).show(figsize=[10,4])
///
}}}
{{{id=174|
A[0,0,]
///
array([ 0.37254903, 0.42745098, 0.80784315, 1. ], dtype=float32)
}}}
{{{id=44|
import pylab
A_image = pylab.mean(pylab.imread(DATA + 'seattle.png'), 2)
@interact
def svd_image(i=(20,(1..100)), display_axes=True):
u,s,v = pylab.linalg.svd(A_image)
A = sum(s[j]*pylab.outer(u[0:,j], v[j,0:]) for j in range(i))
g = graphics_array([matrix_plot(A),matrix_plot(A_image)])
show(g, axes=display_axes, figsize=(8,3))
html('
Compressed using %s eigenvalues
'%i)
///
}}}
{{{id=94|
///
}}}
{{{id=55|
///
}}}
Example: 3d Plots
{{{id=84|
var('x y')
plot3d( 4*x*exp(-x^2-y^2), (x,-2,2), (y,-2,2) )
///
}}}
{{{id=124|
( sphere((0,0,0), opacity=0.7) + sphere((0,1,0), color='red', opacity=0.5)
+ icosahedron((1,1,0), color='green') )
///
}}}
{{{id=116|
L = []
@interact
def random_list(number_of_points=(10..50), dots=True):
n = normalvariate
global L
if len(L) != number_of_points:
L = [(n(0,1), n(0,1), n(0,1)) for i in range(number_of_points)]
G = list_plot3d(L,interpolation_type='nn', texture=Color('#ff7500'),num_points=120)
if dots: G += point3d(L)
G.show()
///
}}}
{{{id=54|
///
}}}
3d plotting (using jmol) is fast even though it does not use Java3d or OpenGL or require any special signed code or drivers.
{{{id=53|
# Yoda! -- over 50,000 triangles.
from scipy import io
X = io.loadmat(DATA + 'yodapose.mat')
from sage.plot.plot3d.index_face_set import IndexFaceSet
V = X['V']; F3=X['F3']-1; F4=X['F4']-1
Y = IndexFaceSet(F3,V,color='green') + IndexFaceSet(F4,V,color='green')
Y = Y.rotateX(-1)
Y.show(aspect_ratio=[1,1,1], frame=False, figsize=4)
html('"Use the source, Luke..."')
///
"Use the source, Luke..."
}}}
{{{id=145|
///
}}}
Real World Example: Super Fast Code (using Cython)
to sage-support
date Sat, Jan 31, 2009 at 11:15 AM
Hi,
I received first a MemoryError, and later on Sage reported:
uitkomst1=[]
uitkomst2=[]
eind=int((10^9+2)/(2*sqrt(3)))
print eind
for y in srange(1,eind):
test1=is_square(3*y^2+1,True)
test2=is_square(48*y^2+1,True)
if test1[0] and test1[1]%3==2: uitkomst1.append((y,(2*test1[1]-1)/3))
if test2[0] and test2[1]%3==1: uitkomst2.append((y,(2*test2[1]+1)/3))
print uitkomst1
een=sum([3*x-1 for (y,x) in uitkomst1 if 3*x-1<10^9])
print uitkomst2
twee=sum([3*x+1 for (y,x) in uitkomst2 if 3*x+1<10^9])
print een+twee
If you replace 10^9 with 10^6, the above listing works properly.
Maybe I made a mistake?
Rolandb
I rewrote Roland's code slightly so it wouldn't waste memory by constructing big lists... but the result was slow.
{{{id=150|
def f_python(n):
uitkomst1=[]
uitkomst2=[]
eind=int((n+2)/(2*sqrt(3)))
print eind
for y in (1..eind):
test1=is_square(3*y^2+1,True)
test2=is_square(48*y^2+1,True)
if test1[0] and test1[1]%3==2:
uitkomst1.append((y,(2*test1[1]-1)/3))
if test2[0] and test2[1]%3==1:
uitkomst2.append((y,(2*test2[1]+1)/3))
print uitkomst1
een=sum(3*x-1 for (y,x) in uitkomst1 if 3*x-1<10^9)
print uitkomst2
twee=sum(3*x+1 for (y,x) in uitkomst2 if 3*x+1<10^9)
print een+twee
///
}}}
{{{id=182|
time f_python(10^5)
///
28868
[(1, 1), (15, 17), (209, 241), (2911, 3361)]
[(1, 5), (14, 65), (195, 901), (2716, 12545)]
51408
Time: CPU 0.72 s, Wall: 0.77 s
}}}
{{{id=8|
time f_python(10^6)
///
288675
[(1, 1), (15, 17), (209, 241), (2911, 3361), (40545, 46817)]
[(1, 5), (14, 65), (195, 901), (2716, 12545), (37829, 174725)]
716034
Time: CPU 7.14 s, Wall: 7.65 s
}}}
3d Plots
{{{id=1|
%cython
from sage.all import is_square
cdef extern from "math.h":
long double sqrtl(long double)
def f(n):
uitkomst1=[]
uitkomst2=[]
cdef long long eind=int((n+2)/(2*sqrt(3)))
cdef long long y, t
print eind
for y in range(1,eind):
t = sqrtl( (3*y*y + 1))
if t * t == 3*y*y + 1:
uitkomst1.append((y, (2*t-1)/3))
t = sqrtl( (48*y*y + 1))
if t * t == 48*y*y + 1:
uitkomst2.append((y, (2*t+1)/3))
print uitkomst1
een=sum([3*x-1 for (y,x) in uitkomst1 if 3*x-1<10^9])
print uitkomst2
twee=sum([3*x+1 for (y,x) in uitkomst2 if 3*x+1<10^9])
print een+twee
///
}}}
{{{id=45|
time f(10^5)
///
28868
[(1L, 1L), (4L, 4L), (15L, 17L), (56L, 64L), (209L, 241L), (780L, 900L), (2911L, 3361L), (10864L, 12544L)]
[(1L, 5L), (14L, 65L), (195L, 901L), (2716L, 12545L)]
2
Time: CPU 0.00 s, Wall: 0.00 s
}}}
{{{id=4|
time f(10^6)
///
288675
[(1L, 1L), (4L, 4L), (15L, 17L), (56L, 64L), (209L, 241L), (780L, 900L), (2911L, 3361L), (10864L, 12544L), (40545L, 46817L), (151316L, 174724L)]
[(1L, 5L), (14L, 65L), (195L, 901L), (2716L, 12545L), (37829L, 174725L)]
2
Time: CPU 0.03 s, Wall: 0.03 s
}}}
{{{id=5|
time f(10^9)
///
288675135
[(1L, 1L), (4L, 4L), (15L, 17L), (56L, 64L), (209L, 241L), (780L, 900L), (2911L, 3361L), (10864L, 12544L), (40545L, 46817L), (151316L, 174724L), (564719L, 652081L), (2107560L, 2433600L), (7865521L, 9082321L), (29354524L, 33895684L), (109552575L, 126500417L)]
[(1L, 5L), (14L, 65L), (195L, 901L), (2716L, 12545L), (37829L, 174725L), (526890L, 2433601L), (7338631L, 33895685L), (102213944L, 472105985L)]
2
Time: CPU 25.60 s, Wall: 26.50 s
}}}
{{{id=6|
7.14/0.03
///
238.000000000000
}}}
This is not a contrived example. This is a real world example that came up a this weekend. For C-style computations, Sage (via Cython) is as fast as C.
{{{id=151|
///
}}}
Numerical Matrix Algebra
{{{id=184|
a = random_matrix(RDF, 500) # entries in [-1,1]
G = graphics_array([matrix_plot(a^i) for i in [-3,-2,-1,1,2,3]], 2, 3)
show(G, axes=False)
///
}}}
{{{id=148|
@interact
def f(i=slider([-10..-1]+[1..10])):
show(matrix_plot(a^i),axes=False)
///
to sage-support date Sat, Jan 31, 2009 at 11:15 AM Hi,
I received first a MemoryError, and later on Sage reported:
uitkomst1=[] uitkomst2=[] eind=int((10^9+2)/(2*sqrt(3))) print eind for y in srange(1,eind): test1=is_square(3*y^2+1,True) test2=is_square(48*y^2+1,True) if test1[0] and test1[1]%3==2: uitkomst1.append((y,(2*test1[1]-1)/3)) if test2[0] and test2[1]%3==1: uitkomst2.append((y,(2*test2[1]+1)/3)) print uitkomst1 een=sum([3*x-1 for (y,x) in uitkomst1 if 3*x-1<10^9]) print uitkomst2 twee=sum([3*x+1 for (y,x) in uitkomst2 if 3*x+1<10^9]) print een+twee
If you replace 10^9 with 10^6, the above listing works properly.
Maybe I made a mistake?
Rolandb
I rewrote Roland's code slightly so it wouldn't waste memory by constructing big lists... but the result was slow.
{{{id=128|
///
}}}
Plotting an elliptic curve over a finite field
{{{id=127|
p = 389
E = EllipticCurve(GF(p).random_element())
E.plot()
///
}}}
This is not a contrived example. This is areal world example that came up last weekend. For C-style computations, Sage (via Cython) is as fast as C.
{{{id=137|
E.abelian_group()
///
(Multiplicative Abelian Group isomorphic to C180 x C2, ((262 : 190 : 1), (228 : 0 : 1)))
}}}
Creating another random elliptic curve over a bigger finite field and compute its cardinality in seconds:
{{{id=129|
p = next_prime(10^60)
E = EllipticCurve(GF(p).random_element())
show(E)
///