5-You can use Sage to teach Calculus.sagews

You can use Sage to teach Calculus

2 + 3

5

x = var('x')
show(plot(sin(x^2), 0, pi))

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

1/16*sqrt(pi)*((I + 1)*sqrt(2)*erf((1/2*I + 1/2)*sqrt(2)*x) + (I - 1)*sqrt(2)*erf((1/2*I - 1/2)*sqrt(2)*x) - (I - 1)*sqrt(2)*erf(sqrt(-I)*x) + (I + 1)*sqrt(2)*erf((-1)^(1/4)*x))

show(a)

$\displaystyle \frac{1}{16} \, \sqrt{\pi} {\left(\left(i + 1\right) \, \sqrt{2} \operatorname{erf}\left(\left(\frac{1}{2} i + \frac{1}{2}\right) \, \sqrt{2} x\right) + \left(i - 1\right) \, \sqrt{2} \operatorname{erf}\left(\left(\frac{1}{2} i - \frac{1}{2}\right) \, \sqrt{2} x\right) - \left(i - 1\right) \, \sqrt{2} \operatorname{erf}\left(\sqrt{-i} x\right) + \left(i + 1\right) \, \sqrt{2} \operatorname{erf}\left(\left(-1\right)^{\frac{1}{4}} x\right)\right)}$

search_src("integrate")

Interact: please open in CoCalc

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

(a, b, c, X)

$\displaystyle X = -\frac{b + \sqrt{b^{2} - 4 \, a c}}{2 \, a}$

Finite Field Arithmetic

k.<a> = GF(2^8)
k

Finite Field in a of size 2^8

latex(k)

\Bold{F}_{2^{8}}

show(k)

$\displaystyle \Bold{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^2 + x + 2

p = p/p.leading_coefficient()
p

x^5 + 2*x^4 + 2*x^2 + x + 2

k.<a> = GF(3^5, modulus=p)
k.modulus()

x^5 + 2*x^4 + 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')

Error in lines 1-1 Traceback (most recent call last): File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 1013, in execute exec compile(block+'\n', '', 'single') in namespace, locals File "", line 1, in <module> ImportError: No module named fplll.fplll

n = 100
a = random_matrix(QQ,n, n+1, num_bound=2^64, den_bound=1)
time a.echelonize()

Error in lines 3-3 Traceback (most recent call last): File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 1013, in execute exec compile(block+'\n', '', 'single') in namespace, locals File "<string>", line 1 time a.echelonize() ^ SyntaxError: invalid syntax

%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

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 A.echelonize()

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

%magma
n := 5000;
k := GF(2);
A := Random(RMatrixSpace(k, n,n));
time E := EchelonForm(A);

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

%magma
n := 20000;
k := GF(2);
A := Random(RMatrixSpace(k, n,n));
time E := EchelonForm(A);

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

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

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.99 s, Wall time: 4.82 s

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

[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

Elliptic Curves

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

y^2 + y = x^3 - x

show(plot(e))

#e.aplist
e.a_invariants()

[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

%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

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

P = E.random_element()
P.order()

9997377

E.cardinality()

9997377

2*P + P

(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

%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

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)'

%gp
gettime;
for(i=1,10^6,z=e*f);
gettime/1000.0

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

D.radius()

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'>

%time
for _ in xrange(10^5):
   w = f*g

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

%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

# 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

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

%magma
P<x,y,z> := PolynomialRing(RationalField(),3);
p:= (x + y + z + 1)^20;
q:= p + 1;
time r:= p*q;

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

%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

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

Algebraic Cryptanalysis

sr = mq.SR(2,1,1,4,gf2=True)
sr

SR(2,1,1,4)

F,s = sr.polynomial_system()
F

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}]

s

{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]