SAGE Demo (SD6 Bristol, 10.11.2007)

1349 days ago by pub

You can use Sage to teach Calculus
2 + 3 
        
5
5
x = var('x') show(plot(sin(x^2), 0, pi)) 
       
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
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}
\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}}
X  =  \frac{-\left( \sqrt{ {b}^{2}  - {{4 a} c} } \right) - b}{{2 a}}

Finite Field Arithmetic
k.<a> = GF(2^8) k 
        
Finite Field in a of size 2^8
Finite Field in a of size 2^8
latex(k) 
        
\mathbf{F}_{2^{8}}
\mathbf{F}_{2^{8}}
show(k) 
        
\mathbf{F}_{2^{8}}
\mathbf{F}_{2^{8}}
P.<x> = GF(3)['x'] while True: p = P.random_element(degree=5) if p.is_irreducible() and p.degree() == 5: break p 
        
x^5 + 2*x^4 + 2*x^3 + 2*x^2 + x + 2
x^5 + 2*x^4 + 2*x^3 + 2*x^2 + x + 2
p = p/p.leading_coefficient() p 
        
x^5 + 2*x^4 + 2*x^3 + 2*x^2 + x + 2
x^5 + 2*x^4 + 2*x^3 + 2*x^2 + x + 2
k.<a> = GF(3^5, modulus=p) k.modulus() 
        
x^5 + 2*x^4 + 2*x^3 + 2*x^2 + x + 2
x^5 + 2*x^4 + 2*x^3 + 2*x^2 + x + 2

Dense Linear Algebra over the Rationals and Integers

The NTRUEncrypt Public Key Cryptosystem is based on the hard mathematical problem of finding very short vectors in lattices of very high dimension. Generate a ntru-like lattice of dimension (400 \times 400), with the coefficients h_i chosen as random 130 bits integers and parameter q=35:
\left( \begin{array}{cccccccc} 1 & 0 & \dots & 0 & h_0 & h_1 & \dots & h_{d-1}\\ 0 & 1 & \dots & 0 & h_1 & h_2 & \dots & h_0 \\ \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ 0 & 0 & \dots & 1 & h_{d-1} & h_0 & \dots & h_{d-1}\\ 0 & 0 & \dots & 0 & q & 0 & \dots & 0 \\ 0 & 0 & \dots & 0 & 0 & q & \dots & 0 \\ \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ 0 & 0 & \dots & 0 & 0 & 0 & \dots & q \\ \end{array} \right)
from sage.libs.fplll.fplll import gen_ntrulike A = gen_ntrulike(200,130,35) time B = A.LLL(algorithm='fpLLL:fast') 
        
Time: CPU 0.14 s, Wall: 0.24 s
Time: CPU 0.14 s, Wall: 0.24 s
n = 100 a = random_matrix(QQ,n, n+1, num_bound=2^64, den_bound=1) time a.echelonize() 
        
Time: CPU 0.14 s, Wall: 0.14 s
Time: CPU 0.14 s, Wall: 0.14 s
%magma n := 100; a := RMatrixSpace(RationalField(), n,n+1)![Random(1,2^64): i in [1..n*(n+1)]]; time e := EchelonForm(a); 
        
Time: 2.100
Time: 2.100

Dense Linear Algebra over Finite Fields
A = random_matrix(GF(127),2000,2000) B = random_matrix(GF(127),2000,2000) time C = A._multiply_linbox(B) 
        
Time: CPU 2.60 s, Wall: 3.13 s
Time: CPU 2.60 s, Wall: 3.13 s
time A.echelonize() 
        
CPU time: 1.50 s,  Wall time: 1.54 s
CPU time: 1.50 s,  Wall time: 1.54 s
n = 5000 A = random_matrix(GF(2),n,n) time E = A.echelon_form() 
        
Time: CPU 0.20 s, Wall: 0.21 s
Time: CPU 0.20 s, Wall: 0.21 s
%magma n := 5000; k := GF(2); A := Random(RMatrixSpace(k, n,n)); time E := EchelonForm(A); 
        
Time: 0.700
Time: 0.700
n = 20000 A = random_matrix(GF(2),n,n) time E = A.echelon_form() 
        
Time: CPU 12.70 s, Wall: 13.25 s
Time: CPU 12.70 s, Wall: 13.25 s
%magma n := 20000; k := GF(2); A := Random(RMatrixSpace(k, n,n)); time E := EchelonForm(A); 
        
Time: 14.510
Time: 14.510

Sparse Linear Algebra over Finite Fields
A = random_matrix(GF(127),3000,3000,density=1.0/3000, sparse=True) time A.echelonize() A.rank() 
        
Time: CPU 0.11 s, Wall: 0.16 s
1635
Time: CPU 0.11 s, Wall: 0.16 s
1635
A = random_matrix(GF(127),1500,1500,density=10/1500, sparse=True) b = random_matrix(GF(127),1500,1) time c = A\b (A*c) == b 
        
Time: CPU 3.27 s, Wall: 3.72 s
True
Time: CPU 3.27 s, Wall: 3.72 s
True
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 *
203703597633448608626844568840937816105146839366593625063614044935438129\
9763336706183397533
CPU time: 3.99 s,  Wall time: 4.82 s
1099511627791 * 2037035976334486086268445688409378161051468393665936250636140449354381299763336706183397533
CPU time: 3.99 s,  Wall time: 4.82 s
time ecm.factor(next_prime(2^40) * next_prime(2^300)) 
        
[1099511627791,
203703597633448608626844568840937816105146839366593625063614044935438129\
9763336706183397533]
CPU time: 0.18 s,  Wall time: 1.40 s
[1099511627791, 2037035976334486086268445688409378161051468393665936250636140449354381299763336706183397533]
CPU time: 0.18 s,  Wall time: 1.40 s
v,t = qsieve(next_prime(2^90)*next_prime(2^91),time=True) print v, t[:4] 
        
[1237940039285380274899124357, 2475880078570760549798248507] 3.99
[1237940039285380274899124357, 2475880078570760549798248507] 3.99

Elliptic Curves
e = EllipticCurve("37a") # Cremona Label show(e) 
        
y^2 + y = x^3 - x 
y^2 + y = x^3 - x
show(plot(e)) 
       
#e.aplist e.a_invariants() 
        
[0, 0, 1, -1, 0]
[0, 0, 1, -1, 0]
print e.aplist(20) time n = e.aplist(10^6) 
        
[-2, -3, -2, -1, -5, -2, 0, 0]
Time: CPU 5.22 s, Wall: 6.00 s
[-2, -3, -2, -1, -5, -2, 0, 0]
Time: CPU 5.22 s, Wall: 6.00 s
%magma print GetVersion(); E := EllipticCurve([0,0,1,-1,0]); print TracesOfFrobenius(E, 20); time n := TracesOfFrobenius(E,10^6); 
        
2 14 14

[ -2, -3, -2, -1, -5, -2, 0, 0 ]
Time: 5.980
2 14 14

[ -2, -3, -2, -1, -5, -2, 0, 0 ]
Time: 5.980
E = e.change_ring(GF(997)) show(E.plot()) 
       
k = GF(next_prime(10^7)) E = EllipticCurve(k, [k.random_element(),k.random_element()]) E 
        
Elliptic Curve defined by y^2  = x^3 + 3004046*x + 7315639 over Finite
Field of size 10000019
Elliptic Curve defined by y^2  = x^3 + 3004046*x + 7315639 over Finite Field of size 10000019
P = E.random_element() P.order() 
        
9997377
9997377
E.cardinality() 
        
9997377
9997377
2*P + P 
        
(6921336 : 2053572 : 1)
(6921336 : 2053572 : 1)

p-Adic Numbers
R = Zp(5) e = R(3125346) f = R(5*34234234) t = cputime() for i in xrange(10^6): z = e*f cputime(t) 
        
0.45693099999999731
0.45693099999999731
%magma R := pAdicRing(5); e := R!3125346; f := R!5*34234234; time for i in [1..10^6] do z := e*f; end for 
        
Time: 0.750
Time: 0.750
gp.eval("e = %s"%(e._gp_().name())) gp.eval("f = %s"%(f._gp_().name())) 
        
'4*5 + 5^2 + 4*5^3 + 3*5^4 + 4*5^5 + 4*5^6 + 3*5^8 + 2*5^9 + 2*5^10 +
3*5^11 + O(5^21)'
'4*5 + 5^2 + 4*5^3 + 3*5^4 + 4*5^5 + 4*5^6 + 3*5^8 + 2*5^9 + 2*5^10 + 3*5^11 + O(5^21)'
%gp gettime; for(i=1,10^6,z=e*f); gettime/1000.0 
        
0.76500000000000000000000000000000000000
0.76500000000000000000000000000000000000

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
True
D.radius() 
        
5
5

Univariate Polynomials over \mathbf{Z}
P.<x> = ZZ[] deg = 32; coeff=64 f = P.random_element(degree=deg,x=0,y=2^coeff) g = P.random_element(degree=deg,x=0,y=2^coeff) 
       
type(f) 
        
<type
'sage.rings.polynomial.polynomial_integer_dense_flint.Polynomial_integer\
_dense_flint'>
<type 'sage.rings.polynomial.polynomial_integer_dense_flint.Polynomial_integer_dense_flint'>
%time for _ in xrange(10^5): w = f*g 
        
CPU time: 1.39 s,  Wall time: 1.42 s
CPU time: 1.39 s,  Wall time: 1.42 s
gp.eval('f = Pol(%s)'%f.list()) _ = gp.eval('g = Pol(%s)'%g.list()) 
       
%gp gettime; for(i=1,10^5,w=f*g); gettime/1000.0 
        
9.1200000000000000000000000000000000000
9.1200000000000000000000000000000000000
%magma R<x> := PolynomialRing(IntegerRing()); 
       
magma.eval('f := R!%s'%f.list()) _ = magma.eval('g := R!%s'%g.list()) 
       
%magma time for i in [1..10^5] do w := f*g; end for; 
        
Time: 14.260
Time: 14.260
# NTL ff = ntl.ZZX(f.list()) gg = ntl.ZZX(g.list()) 
       
time for i in xrange(10^5): w = ff*gg 
        
CPU time: 8.19 s,  Wall time: 8.51 s
CPU time: 8.19 s,  Wall time: 8.51 s

Multivariate Polynomial Rings over Finite Fields
#In SAGE 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.09 s, Wall: 1.11 s
Time: CPU 1.09 s, Wall: 1.11 s
%magma P<x,y,z> := PolynomialRing(RationalField(),3); p:= (x + y + z + 1)^20; q:= p + 1; time r:= p*q; 
        
Time: 0.450
Time: 0.450
#In SAGE 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
Time: CPU 0.12 s, Wall: 0.12 s
%magma P<x,y,z> := PolynomialRing(GF(32003),3); p := (x + y + z + 1)^20; q := p + 1; time r := p*q; 
        
Time: 0.260
Time: 0.260
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.254657
MAGMA 0.38
Sage/Singular 2.254657
MAGMA 0.38

Algebraic Cryptanalysis
sr = mq.SR(2,1,1,4,gf2=True) sr 
        
SR(2,1,1,4)
SR(2,1,1,4)
F,s = sr.polynomial_system() F 
        
Polynomial System with 104 Polynomials in 36 Variables
Polynomial System with 104 Polynomials in 36 Variables
gb = F.groebner_basis() Ideal(gb).variety() 
        
[{s001: 1, s103: 1, s101: 1, x103: 0, s000: 1, x101: 1, k003: 0, k100:
1, k001: 0, k200: 1, x200: 1, k202: 1, x202: 0, w102: 1, w100: 0, w201:
0, s002: 1, w203: 0, k101: 0, s102: 1, s100: 1, x102: 1, x100: 0, k002:
0, k000: 1, x201: 0, k201: 0, x203: 1, k203: 1, k103: 0, w103: 1, k102:
0, w101: 1, w200: 0, s003: 1, w202: 1}, {s001: 1, s103: 0, s101: 1,
x103: 1, s000: 1, x101: 1, k003: 0, k100: 0, k001: 1, k200: 0, x200: 1,
k202: 1, x202: 0, w102: 1, w100: 1, w201: 1, s002: 0, w203: 1, k101: 0,
s102: 1, s100: 1, x102: 0, x100: 0, k002: 0, k000: 0, x201: 0, k201: 1,
x203: 0, k203: 0, k103: 1, w103: 1, k102: 1, w101: 0, w200: 1, s003: 1,
w202: 1}, {s001: 0, s103: 1, s101: 1, x103: 0, s000: 1, x101: 1, k003:
0, k100: 0, k001: 0, k200: 0, x200: 1, k202: 0, x202: 1, w102: 0, w100:
1, w201: 1, s002: 0, w203: 1, k101: 1, s102: 0, s100: 1, x102: 0, x100:
0, k002: 1, k000: 0, x201: 0, k201: 0, x203: 1, k203: 0, k103: 0, w103:
1, k102: 0, w101: 1, w200: 0, s003: 1, w202: 0}]
[{s001: 1, s103: 1, s101: 1, x103: 0, s000: 1, x101: 1, k003: 0, k100: 1, k001: 0, k200: 1, x200: 1, k202: 1, x202: 0, w102: 1, w100: 0, w201: 0, s002: 1, w203: 0, k101: 0, s102: 1, s100: 1, x102: 1, x100: 0, k002: 0, k000: 1, x201: 0, k201: 0, x203: 1, k203: 1, k103: 0, w103: 1, k102: 0, w101: 1, w200: 0, s003: 1, w202: 1}, {s001: 1, s103: 0, s101: 1, x103: 1, s000: 1, x101: 1, k003: 0, k100: 0, k001: 1, k200: 0, x200: 1, k202: 1, x202: 0, w102: 1, w100: 1, w201: 1, s002: 0, w203: 1, k101: 0, s102: 1, s100: 1, x102: 0, x100: 0, k002: 0, k000: 0, x201: 0, k201: 1, x203: 0, k203: 0, k103: 1, w103: 1, k102: 1, w101: 0, w200: 1, s003: 1, w202: 1}, {s001: 0, s103: 1, s101: 1, x103: 0, s000: 1, x101: 1, k003: 0, k100: 0, k001: 0, k200: 0, x200: 1, k202: 0, x202: 1, w102: 0, w100: 1, w201: 1, s002: 0, w203: 1, k101: 1, s102: 0, s100: 1, x102: 0, x100: 0, k002: 1, k000: 0, x201: 0, k201: 0, x203: 1, k203: 0, k103: 0, w103: 1, k102: 0, w101: 1, w200: 0, s003: 1, w202: 0}]
s 
        
{k000: 1, k001: 0, k002: 0, k003: 0}
{k000: 1, k001: 0, k002: 0, k003: 0}
A,v = F.coefficient_matrix() 
       
A.visualize_structure() 
       
F.round(2) 
        
[w200 + k100 + x100 + x102 + x103, w201 + k101 + x100 + x101 + x103 + 1,
w202 + k102 + x100 + x101 + x102 + 1, w203 + k103 + x101 + x102 + x103,
x100*w100 + x103*w100 + x102*w101 + x101*w102 + x100*w103, x100*w100 +
x101*w100 + x100*w101 + x103*w101 + x102*w102 + x101*w103, x101*w100 +
x102*w100 + x100*w101 + x101*w101 + x100*w102 + x103*w102 + x102*w103,
x102*w100 + x101*w101 + x100*w102 + x103*w103 + 1, x100*w100 + x101*w100
+ x103*w100 + x101*w101 + x100*w102 + x102*w102 + x100*w103 + x100,
x102*w100 + x100*w101 + x101*w101 + x103*w101 + x101*w102 + x100*w103 +
x102*w103 + x101, x100*w100 + x101*w100 + x102*w100 + x102*w101 +
x100*w102 + x101*w102 + x103*w102 + x101*w103 + x102, x101*w100 +
x100*w101 + x102*w101 + x100*w102 + x101*w103 + x103*w103 + x103,
x100*w100 + x102*w100 + x103*w100 + x100*w101 + x101*w101 + x102*w102 +
x100*w103 + w100, x101*w100 + x103*w100 + x101*w101 + x102*w101 +
x100*w102 + x103*w102 + x101*w103 + w101, x100*w100 + x102*w100 +
x100*w101 + x102*w101 + x103*w101 + x100*w102 + x101*w102 + x102*w103 +
w102, x101*w100 + x102*w100 + x100*w101 + x103*w101 + x101*w102 +
x103*w103 + w103, x100^2 + x100, x101^2 + x101, x102^2 + x102, x103^2 +
x103, w100^2 + w100, w101^2 + w101, w102^2 + w102, w103^2 + w103]
[w200 + k100 + x100 + x102 + x103, w201 + k101 + x100 + x101 + x103 + 1, w202 + k102 + x100 + x101 + x102 + 1, w203 + k103 + x101 + x102 + x103, x100*w100 + x103*w100 + x102*w101 + x101*w102 + x100*w103, x100*w100 + x101*w100 + x100*w101 + x103*w101 + x102*w102 + x101*w103, x101*w100 + x102*w100 + x100*w101 + x101*w101 + x100*w102 + x103*w102 + x102*w103, x102*w100 + x101*w101 + x100*w102 + x103*w103 + 1, x100*w100 + x101*w100 + x103*w100 + x101*w101 + x100*w102 + x102*w102 + x100*w103 + x100, x102*w100 + x100*w101 + x101*w101 + x103*w101 + x101*w102 + x100*w103 + x102*w103 + x101, x100*w100 + x101*w100 + x102*w100 + x102*w101 + x100*w102 + x101*w102 + x103*w102 + x101*w103 + x102, x101*w100 + x100*w101 + x102*w101 + x100*w102 + x101*w103 + x103*w103 + x103, x100*w100 + x102*w100 + x103*w100 + x100*w101 + x101*w101 + x102*w102 + x100*w103 + w100, x101*w100 + x103*w100 + x101*w101 + x102*w101 + x100*w102 + x103*w102 + x101*w103 + w101, x100*w100 + x102*w100 + x100*w101 + x102*w101 + x103*w101 + x100*w102 + x101*w102 + x102*w103 + w102, x101*w100 + x102*w100 + x100*w101 + x103*w101 + x101*w102 + x103*w103 + w103, x100^2 + x100, x101^2 + x101, x102^2 + x102, x103^2 + x103, w100^2 + w100, w101^2 + w101, w102^2 + w102, w103^2 + w103]