NameError: name 'y' is not defined
}}}
{{{id=163|
show(a)
///
{x}^{y} + {\pi}^{e} + \frac{26}{3}
}}}
{{{id=20|
show(expand(a^2))
///
{x}^{{2 y}} + {{2 {\pi}^{e} } {x}^{y} } + \frac{{52 {x}^{y} }}{3} + {\pi}^{{2 e}} + \frac{{52 {\pi}^{e} }}{3} + \frac{676}{9}
}}}
Solve equations
{{{id=101|
var('a,b,c,x')
show(solve(x^2 + sqrt(17)*a*x + b == 0, x))
///
\begin{array}{l}[x = \frac{-\left( \sqrt{ {17 {a}^{2} } - {4 b} } \right) - {\sqrt{ 17 } a}}{2},\\
x = \frac{\sqrt{ {17 {a}^{2} } - {4 b} } - {\sqrt{ 17 } a}}{2}]\end{array}
}}}
{{{id=181|
var('a,b,c,x')
show(solve(a*x^3 + b*x + c == 0, x)[0])
///
x = {\left( \frac{{-\sqrt{ 3 } i}}{2} - \frac{1}{2} \right) {\left( \frac{\sqrt{ \frac{{{27 a} {c}^{2} } + {4 {b}^{3} }}{a} }}{{{6 \sqrt{ 3 }} a}} - \frac{c}{{2 a}} \right)}^{\frac{1}{3}} } - \frac{{\left( \frac{{\sqrt{ 3 } i}}{2} - \frac{1}{2} \right) b}}{{{3 a} {\left( \frac{\sqrt{ \frac{{{27 a} {c}^{2} } + {4 {b}^{3} }}{a} }}{{{6 \sqrt{ 3 }} a}} - \frac{c}{{2 a}} \right)}^{\frac{1}{3}} }}
}}}
{{{id=200|
QQ
///
Rational Field
}}}
{{{id=204|
a = 290823098423/2349023094823084
///
}}}
{{{id=205|
a
///
290823098423/2349023094823084
}}}
{{{id=206|
a^20
///
18731732272438885579309471548372824428588245150340221923545430485496458533735980578867523409620449920918256575231891509136504267736928258649691464628040600331999325542157848840014655340456935577675516272236483489312137742807911201/26166505842766494283561401292089146973182567348625157751659094334594689380717152836075367631811803015423294643035334638273762419270762402702701922882086138228898233623197512253298025231805031827112521943486787630994426538549855684872362768419151385105459598073048535928789031985181870396985106351589924274176
}}}
{{{id=207|
QQ[sqrt(2)]
///
Number Field in sqrt2 with defining polynomial x^2 - 2
}}}
{{{id=201|
RDF
///
Real Double Field
}}}
{{{id=202|
R = RealIntervalField(100); R
///
Real Interval Field with 100 bits of precision
}}}
{{{id=203|
R((1.01030,1.010332))
///
1.0103?
}}}
{{{id=164|
A = random_matrix(QQ, 500); v = random_matrix(QQ,500,1)
time x = A \ v
///
Time: CPU 1.56 s, Wall: 1.54 s
}}}
{{{id=3|
len(str(x[0]))
///
1484
}}}
{{{id=32|
///
}}}
Example: A Huge Integer Determinant
{{{id=210|
///
}}}
{{{id=31|
a = random_matrix(ZZ,200,200,x=-2^127,y=2^127)
time d = a.determinant()
len(str(d))
///
Time: CPU 3.78 s, Wall: 3.85 s
7786
}}}
We can also copy this matrix over to Maple and compute the same determinant there...
{{{id=167|
a[0,0]
///
-91532941219663566551217034654788002094
}}}
{{{id=27|
maple.with_package('LinearAlgebra')
B = maple(a)
t = maple.cputime()
time c = B.Determinant()
maple.cputime(t)
///
Time: CPU 0.00 s, Wall: 23.18 s
22.613
}}}
{{{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: 10.74 s
10.51
}}}
{{{id=208|
e == c
///
True
}}}
{{{id=112|
bernoulli
///
}}}
{{{id=25|
///
}}}
Example: A Symbolic Expression
{{{id=50|
x = var('x')
f(x) = sin(3*x)*x+log(x) + 1/(x+1)^2
show(f)
///
x \ {\mapsto}\ {x \sin \left( {3 x} \right)} + \log \left( x \right) + \frac{1}{{\left( x + 1 \right)}^{2} }
}}}
Plotting functions has similar syntax to Mathematica:
{{{id=171|
show(f.integrate(x))
///
x \ {\mapsto}\ \frac{\sin \left( {3 x} \right) - {{3 x} \cos \left( {3 x} \right)}}{9} + {x \log \left( x \right)} - \frac{1}{x + 1} - x
}}}
{{{id=109|
plot(f,(0.01,2), thickness=4) + text("Mathematica-style plotting in Sage", (1,-2), rgbcolor='black')
///
}}}
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) = sin(3*x)*x + log(z) + 1/(1+y)^2
///
}}}
{{{id=74|
g = f._fast_float_(x,y,z)
timeit('g(4.5r,3.2r,5.7r)')
///
625 loops, best of 3: 557 ns per loop
}}}
{{{id=75|
%python
import math
def g(x): return math.sin(3*x)*x + log(z) + 1/(1+y)**2
///
}}}
{{{id=77|
timeit('g(4.5r)')
///
625 loops, best of 3: 163 µs per loop
}}}
{{{id=48|
plot3d(sin(3*x*y) + x*y + sin(y^3), (x,0,3), (y,0,3))
///
}}}
{{{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))
///
\frac{\sin \left( {3 x} \right) - {{3 x} \cos \left( {3 x} \right)}}{9} + {x \log \left( x \right)} - \frac{1}{x + 1} - x
}}}
The command maple(...) fires up Maple (if you have it!), and creates a reference to a live object:
{{{id=0|
m = maple(f); m
///
sin(3*x)*x+ln(x)+1/(x+1)^2
}}}
{{{id=89|
type(m)
///
}}}
{{{id=90|
m.parent()
///
Maple
}}}
{{{id=91|
m.parent().pid()
///
24038
}}}
{{{id=92|
os.system('ps ax |grep 24038')
///
24038 s007 Ss+ 0:00.01 /bin/sh /Users/wstein/bin/maple -t
24233 s015 S+ 0:00.00 sh -c ps ax |grep 24038
24235 s015 R+ 0:00.00 grep 24038
0
}}}
Use Maple objects via a Pythonic notation:
{{{id=88|
show(m.integrate('x'))
///
1/9\,\sin \left( 3\,x \right) -1/3\,\cos \left( 3\,x \right) x+\ln
\left( x \right) x-x- \left( x+1 \right) ^{-1}
}}}
{{{id=42|
mathematica(f).Integrate(x)
///
-x - (1 + x)^(-1) - (x*Cos[3*x])/3 + x*Log[x] + Sin[3*x]/9
}}}
{{{id=43|
///
}}}
Example: Interactive Image Compression
This illustrates pylab (matplotlib + numpy), Sage plotting, html output, and @interact.
{{{id=97|
# first just play
import pylab
A = pylab.imread(DATA + 'emoryimage.png')
graphics_array([matrix_plot(A), matrix_plot(1-A[0:,0:,2])]).show(figsize=[10,4])
///
}}}
{{{id=174|
A[0,0,]
///
array([ 0.76862746, 0.89803922, 0.96470588, 1. ], dtype=float32)
}}}
{{{id=44|
import pylab
A_image = pylab.mean(pylab.imread(DATA + 'emoryimage.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
}}}
Rewrite using CYTHON
{{{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
}}}
A speedup by a factor of 238!
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=192|
a = random_matrix(RDF, 3); show(a)
///
\left(\begin{array}{rrr}
-0.898388350554 & -0.245114654439 & 0.630541996552 \\
0.654081247919 & 0.696119569249 & -0.957945963168 \\
0.00662634420862 & -0.704817828713 & -0.517608714138
\end{array}\right)
}}}
{{{id=193|
show(a.eigenvalues())
///
\left[0.961673429437,
-0.4808854176,
-1.20066550728\right]
}}}
{{{id=195|
a = random_matrix(RDF, 1000)
time v = a.eigenvalues()
///
Time: CPU 5.81 s, Wall: 6.21 s
}}}
{{{id=196|
show(points(v), axes=False, aspect_ratio=1, figsize=4)
///
}}}
{{{id=184|
a = random_matrix(RDF, 200) # read doubles 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=187|
@interact
def f(n=(10..400), field=[RDF,CDF]):
G = points(random_matrix(field,n).eigenvalues())
show(G, axes=False, frame=True)
///
}}}
{{{id=186|
///
}}}
{{{id=144|
///
}}}
{{{id=133|
///
}}}
{{{id=61|
///
}}}
{{{id=108|
///
}}}
