Ce document dans d'autres formats: feuille de travail, source RST.
Moyens:
Premiers calculs avec Sage :
1 + 1
( 1 + 2 * (3 + 5) ) * 2
2^3
2**3
20/6
2^10
2^100
2^1000
20.0 / 14
numerical_approx(20/14)
numerical_approx(2^1000)
numerical_approx(20/14, digits=60)
Premier exemple d'instabilité numérique:
(1 + 10^50) - 10^50
(1.0 + 10^50) - 10^50
En exigeant suffisamment de précision:
x = n(1, digits=49)
(x+10^50)-10^50
x = n(1, digits=48)
(x+10^50)-10^50
Quelques exemples supplémentaires:
factorial(100)
factor(2^(2^5)+1)
Calcul formel avec des fonctions et constantes usuelles:
arccos(sin(pi/3))
sqrt(2)
exp(I*pi/6)
simplify(arccos(sin(pi/3)))
simplify(exp(i*pi/6))
numerical_approx( 6*arccos( sin(pi/3)), digits=60 )
numerical_approx( sqrt(2), digits=60 )
m = 7 % 4; m
3 * m + 1
Et si l'on veut faire tout les calculs suivants modulo $4$ :
Z4 = IntegerModRing(4); Z4
m = Z4(7); m
Par la suite, tous les calculs faisant intervenir m sont fait modulo $4$. Ainsi, dans l'exemple suivants, $3$ et $1$ sont automatiquement convertis dans $\ZZ/n\ZZ$ :
3 * m + 1
Corps finis:
Z3 = GF(3); Z3
a = matrix(QQ, [[1,2,3],[2,4,8],[3,9,27]])
(a^2 + 1) * a^(-1)
P = QQ['x']; P
F = P.fraction_field(); F
p = P(x+1) * P(x); p
p + 1/p
parent(p + 1/p)
k.<a> = NumberField(x^3 + x + 1)
a^3
a^4+3*a
demo-symbolics
Calcul formel =
Calcul mathématique (computational mathematics) =
Systèmes généralistes:
Systèmes spécialisés:
var('a,b,c,d,e,f,g')
F = a + b * c + d * e * sin(f)^g
F.operands()
Que se passe-t'il lorsque l'on fait:
F = 0
(comptage de références ou glaneur de cellule)
Listes, ensembles et tables d'association:
liste = [sin(1+x), 3, sin(1+x)]; liste
ensemble = { sin(1+x), 3, sin(1+x) }; ensemble
tableAssociative = { sin(1+x) : 1, 3 : 2 }
tableAssociative[3]
tableAssociative[sin(1+x)]
Exécution conditionnelle, boucles, fonctions, ...
Début de l'utilisation de l'ordinateur comme outil pour la recherche en mathématique:
SageSage est développé par 300 enseignants et chercheurs dans le monde entierSage est orienté objet¶Python et Sage utilisent fortement la programmation orientée objet. Même si cela reste relativement transparent pour une utilisation occasionnelle, il est utile d’en savoir un minimum, d’autant que ce fait est très naturel dans un contexte mathématique.
Le paradigme de la programmation orientée objet s’appuie sur un principe: «toute entité du monde physique ou mathématique que l’on souhaite manipuler avec l’ordinateur est modélisé par un objet»; de plus cette objet est une instance d’une classe. Ainsi, le nombre rationnel $o=12/35$ est modélisé par un objet qui est une instance de la classe Rational :
o = 12/35
type(o)
Noter que cette classe est vraiment associée à l’objet $12/35$, et non seulement à la variable o qui le contient:
type(12/35)
Précisons les définitions. Un objet est une portion de la mémoire de l’ordinateur qui contient l’information nécessaire pour représenter l’entité qu’il modélise. La classe quant à elle définit deux choses:
Rational stipule qu’un nombre rationel comme $12/35$ est représenté, en gros, par deux nombres entiers: son numérateur et son dénominateur.numer, abs, {__mult__}, {__add__}, ...).Pour factoriser un nombre entier $o$, on va donc appeller la méthode factor avec la syntaxe suivante:
o = 720
o.factor()
que l’on peut lire comme: «prendre la valeur de o et lui appliquer la méthode factor sans autre argument». Sous le capot, effectue le calcul suivant:
type(o).factor(o)
De gauche à droite: «demander à la classe de (la valeur de) o (type(o)) la méthode appropriée de factorisation (type(o).factor), et l’appliquer à o».
Notons au passage que l’on peut appliquer une opération à une valeur sans passer par une variable:
720.factor()
et donc en particulier enchaîner les opérations, de la gauche vers la droite. Ici, on prend le numérateur d’un nombre rationnel, que l’on factorise ensuite:
o = 720 / 133
o.numerator().factor()
En quoi cela nous concerne-t-il? Tout d’abord, l’orientation objet permet le polymorphisme: quelque soit l’objet o que l’on veut factoriser, on peut toujours utiliser la notation o.factor() (ou son raccourci factor(o)). De même, calquant les notations mathématiques usuelles, le produit de deux objets a et b peut toujours être noté a*b même si l’algorithme utilisé dans chaque cas est différent (Pour une opération arithmétique binaire comme le produit, la procédure de sélection de la méthode appropriée est un peu plus compliquée que ce qui a été décrit précédemment. En effet elle doit gérer des opérations mixtes comme la somme $2 + 3/4$ d’un entier et d’un nombre rationnel. En l’occurence, $2$ sera converti en nombre rationnel $2/1$ avant l’addition. C’est le modèle de coercion de Sage qui est en charge de cela.). Voici un produit de deux nombres entiers:
3 * 7
un produit de deux nombres rationnels, obtenu par produit des numérateurs et dénominateurs puis réduction:
(2/3) * (6/5)
Un produit de deux nombres complexes, utilisant $I^2=-1$ :
(1 + I) * (1 - I)
des produits commutatifs formels de deux expressions:
(x + 2) * (x + 1)
(x + 1) * (x + 2)
Outre la simplicité de notation, cela permet d’écrire des programmes génériques comme:
def puissance_quatre(a):
qui s’appliquent à tout objet admettant les opérations utilisées (ici la multiplication):
puissance_quatre(2)
puissance_quatre(3/2)
puissance_quatre(I)
puissance_quatre(x+1)
M = matrix([[0,-1],[1,0]]); M
puissance_quatre(M)
Plus prosaïquement, l’orientation objet permet l’introspection: on peut ainsi accéder à l’aide en ligne spécifique à la factorisation des nombres entiers avec:
o = 720
x.factor?
voire à l’implantation de cette fonction, précédée de son aide en ligne:
o.factor?
En passant au dessus des détails techniques, on distingue bien que Sage délègue le calcul à d’autres logiciels (Pari, Kash, ...).
Enfin, on peut utiliser la complétion automatique pour demander interactivement à un objet o quelles sont toutes les opérations que l’on peut lui appliquer:
o.<tab>
Ici, il y en a beaucoup; voici celles qui commencent par n :
o.n<tab>
Dans la section précédente, nous avons vu la notion technique de classe d’un objet. Dans la pratique, il est suffisant de savoir que cette notion existe; on a rarement besoin de regarder explicitement le type d’un objet. En revanche Sage introduit une contrepartie plus conceptuelle de cette notion que nous allons aborder maintenant: celle du parent d’un objet.
Supposons par exemple que l’on veuille déterminer si un élément $a$ est inversible. La réponse ne va pas seulement dépendre de l’élément lui-même, mais de l’ensemble $A$ auquel il est considéré appartenir. Par exemple, le nombre $5$ n’est pas inversible dans l’ensemble $\ZZ$ des entiers, son inverse $1/5$ n’étant pas un entier:
a = 5; a
a.is_unit()
En revanche, il est inversible dans l’ensemble des rationnels:
a = 5/1; a
a.is_unit()
Comme nous l’avons vu dans la section précédente, Sage répond différemment à ces deux questions car les éléments $5$ et $5/1$ sont dans des classes différentes:
type(5)
type(5/1)
Dans certains systèmes de calcul formel orientés objet, tels que MuPAD ou Axiom l’ensemble $X$ auquel $x$ est considéré appartenir (ici $\ZZ$ ou $\QQ$) est simplement modélisé par la classe de $x$. Sage suit l’approche de Magma, et modélise $X$ par un objet supplémentaire associé à $x$, et appelé son parent:
parent(5)
parent(5/1)
On peut retrouver ces deux ensembles avec les raccourcis:
ZZ
QQ
et les utiliser pour convertir aisément un élément de l’un à l’autre lorsque cela a un sens:
QQ(5).parent()
ZZ(5/1).parent()
ZZ(1/5)
Voici $1$ en tant qu’entier $1\in\ZZ$, nombre rationnel $1\in\QQ$, et approximations flottantes $1{,}0\in\RR$ et $1{,}0+0{,}0i \in\CC$ :
ZZ(1), QQ(1), RR(1), CC(1)
Selon le même principe, lorsque l'on demande toutes les partitions de l'entier 5, le résultat est un objet qui modélise cet ensemble:
P = Partitions(5); P
Pour obtenir la liste de ces objets, il faut le demander explicitement:
P.list()
Cela permet de manipuler formellement des grands ensembles:
Partitions(100000).cardinality()
Et de calculer paresseusement avec. Ici, on tire au hasard une main de cinq cartes à jouer:
Symboles = Set(["Coeur", "Carreau", "Pique", "Trefle"])
Valeurs = Set([2, 3, 4, 5, 6, 7, 8, 9, 10, "Valet", "Dame", "Roi", "As"])
Cartes = CartesianProduct(Valeurs, Symboles).map(tuple)
Mains = Subsets(Cartes, 5)
Mains.cardinality()
Mains.random_element() # random
et là on manipule un mot infini défini comme point fixe d'un morphisme:
m = WordMorphism('a->acabb,b->bcacacbb,c->baba')
m.fixed_point('a')
Les parents étant eux-même des objets, on peut leur appliquer des opérations. Ainsi, on peut construire le produit cartésien $\QQ^2$ :
cartesian_product([QQ, QQ])
retrouver $\QQ$ comme corps des fractions de $\ZZ$ :
ZZ.fraction_field()
construire l’anneau des polynômes en $x$ à coefficients dans $\ZZ$ :
ZZ['x']
Par empilement successif, on peut construire des structure algébriques avancées comme l’espace des matrices $3\times 3$ à coefficients polynomiaux sur un corps fini:
Z5 = GF(5); Z5
P = Z5['x']; P
M = MatrixSpace(P, 3, 3); M
dont voici un élément:
M.random_element() # random
M.det()
Dans les exemples ci-dessous, nous ferons de l'algèbre linéaire sur le corps fini $\ZZ/7\ZZ$ :
K = GF(7); K
list(K)
Nous avons choisi ce corps à titre d'illustration pour avoir des résultats lisibles. On aurait pu prendre des coefficients entiers, rationnels, ou numériques à plus ou moins haute précision. Les aspects numériques seront abordés plus en détail dans l'exposé suivant. Notons au passage que, même en calcul exact, il est possible de manipuler de relativement grosses matrices:
n = 500
M = random_matrix(K, n, sparse=True, density=3/n)
M.visualize_structure() # not tested
n = 10000
M = random_matrix(K, n, sparse=True, density=3/n)
M.rank() # random
Définissons donc une matrice à coefficients dans $\ZZ/7\ZZ$ :
A = matrix(K, 4, [5,5,4,3,0,3,3,4,0,1,5,4,6,0,6,3]); A
Calculons le polynôme caractéristique de cette matrice:
P = A.characteristic_polynomial(); P
On vérifie le théorème de Cayley-Hamilton sur cet exemple:
P(A)
Notons que l'information sur le corps de base est préservée:
P.parent()
ce qui influe directement sur la factorisation de ce polynôme:
factor(P)
factor(x^4 + 5*x^3 + 6*x + 2)
Le calcul ci-dessus nous donne les valeurs propres: -3=4,-6=1 et -5=2. Quels sont les espaces propres?
A.eigenspaces_left()
Récupérons ces espaces propres:
E = dict(A.eigenspaces_left())
E[2]
E[2] n'est pas une liste de vecteurs ni une matrice, mais un objet qui modélise l'espace propre $E_2$, comme le sous-espace de $(\ZZ/7\ZZ)^4$ décrit par sa base échelon réduite. On peut donc lui poser des questions:
E[2].dimension()
E[2].basis()
V = E[2].ambient_vector_space(); V
Voire faire des calculs avec:
E[2] + E[4]
v = V([1,2,0,3])
v in E[2]
E[2].echelon_coordinates(v)
E[2].is_subspace(E[4])
E[2].is_subspace(V)
Q = V/E[2]; Q
Q( V([0,0,0,1]) )
On veut maintenant manipuler $A$ comme un morphisme sur $V$ :
phi = End(V)(A); phi
v = V.an_element()
v
phi(v)
(phi^-1)(v)
phi^4 + 5*phi^3 + 6*phi + 2
(phi - 1).image()
(phi - 1).kernel() == E[1]
phi.restrict(E[2])
Un parent n’a, en général, pas lui-même un parent, mais une catégorie qui indique ses propriétés:
C = QQ.category(); C
De fait Sage sait que $\QQ$ est un corps:
QQ in Fields()
et donc, par exemple, un groupe additif commutatif:
QQ in CommutativeAdditiveGroups()
Voici tous les axiomes satisfaits par $\QQ$ :
C.axioms()
et les catégories de $\QQ$ :
G = C.category_graph()
G.set_latex_options(format="dot2tex")
view(G, tightpage=True, viewer="pdf")
Sage en déduit que $\QQ[x]$ est un anneau euclidien:
QQ['x'].category()
En général, il peut combiner des axiomes et des structures:
Magmas().Associative() & Magmas().Unital().Inverse() & Sets().Finite()
Et appliquer par exemple le théorème de Wedderburn:
Rings().Division() & Sets().Finite()
Toutes ces propriétés sont utilisées pour calculer rigoureusement et plus efficacement sur les éléments de ces ensembles.
Dans cette section, nous donnons quelques exemples typiques pour lesquels il est important de contrôler le domaine de calcul. En première lecture, on peut passer rapidement sur les exemples plus avancés pour arriver directement à la synthèse de fin de section.
Soit $c$ une expression un tout petit peu compliquée:
a = var('a')
c = (a+1)^2 - (a^2+2*a+1)
et cherchons à résoudre l’équation en $x$ donnée par $cx=0$ :
eq = c * x == 0
L’utilisateur imprudent pourrait être tenté de simplifier cette équation par $c$ avant de la résoudre:
eq2 = eq / c; eq2
solve(eq2, x)
Heureusement, Sage ne fait pas cette erreur:
solve(eq, x)
Ici, Sage a pu résoudre correctement le système car le coefficient $c$ est une expression polynomiale. Il est donc facile de tester si $c$ est nul; il suffit de le développer:
expand(c)
Et d’utiliser le fait que deux polynômes sous forme développée identiques sont égaux. On dit que la forme développée d’un polynôme est une forme normale.
En revanche, sur un exemple à peine plus compliqué, Sage commet une erreur:
c = cos(a)^2 + sin(a)^2 - 1
eq = c*x == 0
solve(eq, x)
alors même qu’il sait faire la simplification et même le test à zéro correctement:
c.simplify_trig()
c.is_zero()
Cet exemple illustre l’importance du test de nullité, et plus généralement des formes normales, dans un domaine de calcul. Sans lui, tout calcul faisant intervenir une division devient hasardeux. Les algorithmes comme le pivot de Gauss en algèbre linéaire sont particulièrement sensibles à ces considérations.
Construisons l’anneau $\QQ[x_1,x_2,x_3,x_4]$ des polynômes en $4$ variables:
R = QQ['x1,x2,x3,x4']; R
x1, x2, x3, x4 = R.gens()
Les éléments de $R$ sont automatiquement représentés sous forme développée:
x1 * (x2 - x3)
qui comme nous l’avons vu est une forme normale. On dit alors que $R$ est à représentation normale. En particulier le test à zéro y est immédiat:
(x1+x2)*(x1-x2) - (x1^2 -x2^2)
Mais ce n’est pas toujours un avantage. Par exemple, si l’on construit le déterminant de Vandermonde $\prod_{1\leq i < j \leq n} (x_i-x_j)$ :
prod( (a-b) for (a,b) in Subsets([x1,x2,x3,x4],2) )
on obtient $4!=24$ termes. Alors que la même construction avec une expression reste sous forme factorisée qui est ici beaucoup plus compacte et lisible:
x1, x2, x3, x4 = var('x1, x2, x3, x4')
prod( (a-b) for (a,b) in Subsets([x1,x2,x3,x4],2) )
De même, une représentation factorisée ou partiellement factorisée permet des calculs de { pgcd} bien plus rapides. Réciproquement, il ne serait pas judicieux de mettre automatiquement tout polynôme sous forme factorisée, même s’il s’agit aussi d’une forme normale, car la factorisation est coûteuse et non compatible avec l’addition.
De manière générale, selon le type de calcul voulu, la représentation idéale d’un élément n’est pas toujours sa forme normale. Cela amène les systèmes de calcul formel à un compromis avec les expressions. Un certain nombre de simplifications basiques, comme la réduction des rationnels ou la multiplication par zéro, y sont effectuées automatiquement; les autres transformations sont laissées à l’initiative de l’utilisateur auquel des commandes spécialisées sont proposées.
Considérons la factorisation de l’expression polynomiale suivante:
x = var('x')
p = 54*x^4+36*x^3-102*x^2-72*x-12
factor(p)
Cette réponse est-elle satisfaisante? Il s’agit bien d’une factorisation de $p$, mais son optimalité dépend fortement du contexte! Pour le moment Sage considère p comme une expression symbolique, qui se trouve être polynomiale. Il ne peut pas savoir si l’on souhaite factoriser $p$ en tant que produit de polynômes à coefficients entiers ou à coefficients rationnels (par exemple). Pour prendre le contrôle, nous allons préciser dans quel ensemble (domaine de calcul?) nous souhaitons considérer $p$. Pour commencer, nous allons considérer $p$ comme un polynôme à coefficient entiers. Nous définissons donc l’anneau $R=\ZZ[x]$ de ces polynômes:
R = ZZ['x']; R
Puis nous convertissons $p$ dans cet anneau:
q = R(p); q
À l’affichage on ne voit pas de différence, mais $q$ sait qu’il est un élément de $R$ :
parent(q)
Du coup, sa factorisation est sans ambiguïté:
factor(q)
On procède de même sur le corps des rationels:
R = QQ['x']; R
q = R(p); q
factor(R(p))
Dans ce nouveau contexte, la factorisation est encore non ambiguë; mais différente de précédemment. Notons au passage que Sage sait que $R$ est un anneau euclidien:
R.category()
et donc en particulier un anneau où la factorisation est unique (voir Figure {fig:premierspas:catégories}).
Cherchons maintenant une factorisation complète sur les nombres complexes. Une première option est de s’autoriser une approximation numérique des nombres complexes avec 16 bits de précision:
R = ComplexField(16)['x']; R
q = R(p); q
factor(R(p))
Une autre est d’agrandir un peu le corps des rationnels; ici, on va rajouter $\sqrt{2}$.
R = QQ[sqrt(2)]['x']; R
q = R(p); q
factor(R(p))
Enfin, peut-être souhaite-t’on que les coefficients soient considérés modulo $5$?
R = GF(5)['x']; R
q = R(p); q
factor(R(p))
Dans les exemples précédents, nous avons illustré comment l’utilisateur peut contrôler le niveau de rigueur dans ses calculs. D’un côté il peut utiliser les expressions symboliques. Ces expressions vivent dans l’anneau des expressions symboliques:
parent(sin(x))
que l’on peut aussi obtenir avec:
SR
Les propriétés de cet anneau sont assez floues; il est commutatif:
SR.category()
et les règles de calcul font en gros l’hypothèse que toutes les variables symboliques sont à valeur dans $\CC$. Le domaine de calcul (expressions polynomiale? rationnelles? trigonométriques?) n’étant pas spécifié explicitement, le résultat d’un calcul nécessite le plus souvent des transformations manuelles pour être mis sous la forme désirée (voir {sec:calculus:simplifications}), en utilisant par exemple expand, combine, collect et simplify. Pour bien utiliser ces fonctions, il faut savoir quel type de transformations elles effectuent et à quel domaine de calcul ces transformations s’appliquent. Ainsi, l’usage aveugle de la fonction simplify peut conduire à des résultats faux. Des variantes de simplify permettent alors de préciser la simplification à effectuer.
D’un autre côté, l’utilisateur peut construire un parent qui va spécifier explicitement le domaine de calcul. Cela est particulièrement intéressant lorsque ce parent est à forme normale: c’est-à-dire que deux objets éléments sont mathématiquement égaux si et seulement si ils ont la même représentation.
Pour résumer, la souplesse est l’avantage principal des expressions:
Les avantages de la déclaration explicite du domaine de calcul sont:
Sage n’est pas un système de calcul certifié; il peut donc toujours y avoir un bogue informatique; mais il n’y aura pas d’utilisation d’hypothèse implicite);L'objectif de la séance est de prendre en main le logiciel Sage. À vous d'explorer ses fonctionnalités selon vos goûts et connaissances préalables, et de préparer une mini-illustration de trois minutes que vous présenterez en fin de séance.
Voici quelques pistes:
Instructions pour les mini-illustrations: