15-Sage Demo.sagews

Sage Demo

Calculus

x = var('x')
@interact
def _(f=sin(x)+1/(x+1),xmax=(pi,(0,4*pi,pi/8))):
  show(f)
  show(plot(f, 0, xmax))

WARNING: Output truncated!

full_output.txt

f  xmax  [removed] [removed] [removed] [removed]

CPU time: 1.24 s, Wall time: 2.98 s

a = integrate(sin(x^2),x); a

sqrt(pi)*((sqrt(2)*I + sqrt(2))*erf((sqrt(2)*I + sqrt(2))*x/2) + (sqrt(2)*I - sqrt(2))*erf((sqrt(2)*I - sqrt(2))*x/2))/8

show(a)

\frac{{\sqrt{ \pi } \left( {\left( {\sqrt{ 2 } i} + \sqrt{ 2 } \right) \text{erf} \left( \frac{{\left( {\sqrt{ 2 } i} + \sqrt{ 2 } \right) x}}{2} \right)} + {\left( {\sqrt{ 2 } i} - \sqrt{ 2 } \right) \text{erf} \left( \frac{{\left( {\sqrt{ 2 } i} - \sqrt{ 2 } \right) x}}{2} \right)} \right)}}{8}

latex(a)

\frac{{\sqrt{ \pi } \left( {\left( {\sqrt{ 2 } i} + \sqrt{ 2 } \right) \text{erf} \left( \frac{{\left( {\sqrt{ 2 } i} + \sqrt{ 2 } \right) x}}{2} \right)} + {\left( {\sqrt{ 2 } i} - \sqrt{ 2 } \right) \text{erf} \left( \frac{{\left( {\sqrt{ 2 } i} - \sqrt{ 2 } \right) x}}{2} \right)} \right)}}{8}

var('a,b,c,X')
s = solve(a*X^2 + b*X + c == 0, X)
show(s[0])

X = \frac{-\left( \sqrt{ {b}^{2} - {{4 a} c} } \right) - b}{{2 a}}

Linear Algebra

A = random_matrix(GF(127),3000,3000,density=1.0/3000,sparse=True)
time A.echelonize()
A.rank()

Time: CPU 0.10 s, Wall: 0.10 s 1626

n = 300
A = random_matrix(GF(127),n,n+100,sparse=True,density=2/n)
A.echelonize()
A.visualize_structure()

sr = mq.SR(4,2,2,4,gf2=True, allow_zero_inversions=True)
F,s = sr.polynomial_system() 
A,v = F.coefficient_matrix()
A.visualize_structure(maxsize=600)

Factoring

time factor(next_prime(2^40) * next_prime(2^300),verbose=0)

1099511627791 * 2037035976334486086268445688409378161051468393665936250636140449354381299763336706183397533 CPU time: 3.78 s, Wall time: 3.86 s

time ecm.factor(next_prime(2^40) * next_prime(2^300))

[1099511627791, 2037035976334486086268445688409378161051468393665936250636140449354381299763336706183397533] CPU time: 0.15 s, Wall time: 0.57 s

v,t = qsieve(next_prime(2^90)*next_prime(2^91),time=True)
print v, t[:4]

[1237940039285380274899124357, 2475880078570760549798248507] 3.55

Elliptic Curves

e = EllipticCurve("37a") # Cremona Label
show(e)

y^2 + y = x^3 - x

show(plot(e))

E = e.change_ring(GF(997))
show(E.plot())

Graph Theory

D = graphs.DodecahedralGraph()
D.show()

D.show3d()

gamma = SymmetricGroup(20).random_element()
E = D.copy()
E.relabel(gamma)
D.is_isomorphic(E)

True

D.radius()

5

Multivariate Polynomial Rings

P.<x,y,z> = PolynomialRing(QQ,3)
p = (x + y + z + 1)^20 # the Fateman fastmult benchmark
q = p + 1
time r = p*q

Time: CPU 1.03 s, Wall: 1.06 s

P.<x,y,z> = PolynomialRing(GF(32003),3)
p = (x + y + z + 1)^20 # the Fateman fastmult benchmark
q = p + 1
time r = p*q

Time: CPU 0.12 s, Wall: 0.12 s

P = PolynomialRing(GF(32003),10,'x')
I = sage.rings.ideal.Cyclic(P,7)
t = cputime()
gb = I.groebner_basis('libsingular:std')
print 'Sage/Singular', cputime(t)

I = sage.rings.ideal.Cyclic(P,7)
t = magma.cputime()
gb = I.groebner_basis('magma:GroebnerBasis')
print 'MAGMA', magma.cputime(t)

Sage/Singular 2.365641 MAGMA 0.38

Combinatorics

C = Combinations(range(5)); C

Combinations of [0, 1, 2, 3, 4]

C.list()

[[], [0], [1], [2], [3], [4], [0, 1], [0, 2], [0, 3], [0, 4], [1, 2], [1, 3], [1, 4], [2, 3], [2, 4], [3, 4], [0, 1, 2], [0, 1, 3], [0, 1, 4], [0, 2, 3], [0, 2, 4], [0, 3, 4], [1, 2, 3], [1, 2, 4], [1, 3, 4], [2, 3, 4], [0, 1, 2, 3], [0, 1, 2, 4], [0, 1, 3, 4], [0, 2, 3, 4], [1, 2, 3, 4], [0, 1, 2, 3, 4]]

C.unrank(10)

[1, 2]

Interfaces

%gap
a := 1;
for i in [1..100] do if IsPrime(i) then a:=a+1; else a:=a+2; fi;  od; a;

1 176

%singular
int a = 1; int i = 1;
for(i=1; i<=100; i=i+1) { if(prime(i) == i){  a=a+1; } else {a=a+2;} };
a;

176

%gp
a=1; for(X=1,100,if(isprime(X),a+=1,a+=2)); a

176

a= gap(1) # or gp(1), magma(1), singular(1)
for i in range(100):
    if gap(i+1).IsPrime():
        a+=1
    else:
        a+=2
a

176

a = gap(176)
type(a)

<class 'sage.interfaces.gap.GapElement'>

pari(a).factor()

[2, 4; 11, 1]

Curve Fitting

3D Plotting

# Ms strip:
u,v = var("u,v")
parametric_plot3d([cos(u)*(1+v*cos(u/2)), sin(u)*(1+v*cos(u/2)), 0.2*v*sin(u/2)], (u,0, 4*pi+0.5), (v,0, 0.3),plot_points=[50,50])

t = Tachyon(xres=500, yres=500, camera_center=(2,7,4), look_at=(2,0,0), raydepth=4)
t.light((10,3,2), 1, (1,1,1))
t.light((10,-3,2), 1, (1,1,1))
t.texture('black', color=(0,0,0))
t.texture('red', color=(1,0,0))
t.texture('grey', color=(.9,.9,.9))
t.plane((0,0,0),(0,0,1),'grey')
t.cylinder((0,0,0),(1,0,0),.01,'black')
t.cylinder((0,0,0),(0,1,0),.01,'black')
E = EllipticCurve('37a')
P = E([0,0])
Q = P
n = 100
for i in range(n):  
   Q = Q + P
   c = i/n + .1
   t.texture('r%s'%i,color=(float(i/n),0,0))
   t.sphere((Q[0], -Q[1], .01), .04, 'r%s'%i)
t.show()

Sequences of Integers

for sq in sloane_find([2,3,5,7], 2):
  print sq[0], sq[1]

Searching Sloane's online database... 40 The prime numbers. 41 a(n) = number of partitions of n (the partition numbers).