UP | HOME

Séance 2 : Algo & Prog avec R

Table des matières

Un ordinateur c'est aussi fait pour calculer …

1 Calcul flottant   KEY

Définissez en R la variable a dont la valeur est \(\sqrt{3}\).

  1. Demandez si \(a^2\) est bien égal à 3.
  2. Expliquez. Indice : print(sprintf('%.17f',0.1)) et sur stackoverflow.

    a <- sqrt(3)
    b <-  a**2
    a == b
    print(paste('Attention à l affichage de a**2', b))
    print(sprintf('Voici la valeur de a**2 en interne : %.17f',b))
    
    [1] FALSE
    [1] "Attention à l affichage de a**2 3"
    [1] "Voici la valeur de a**2 en interne : 2.99999999999999956"
    
  1. Programmer \(100 \times 1.0\) par additions successives (à l'aide d'une boucle) : \(100 \times 1.0 = 1.0 + 1.0 + \dots + 1.0 + 1.0\). Que constatez-vous ?

    acc <- 0
    i <- 0
    while (i < 100) {
      i <- i + 1
      acc <- acc + 1.0
    }
    print(sprintf("100 * 1.0 = %.16f",acc))
    
    [1] "100 * 1.0 = 100.0000000000000000"
    
  2. Même question pour \(100 \times 0.1\).

    acc <- 0
    i <- 0
    while (i < 100) {
      i <- i + 1
      acc <- acc + 0.1
    }
    print(sprintf("100 * 0.1 = %.16f",acc))
    
    [1] "100 * 0.1 = 9.9999999999999805"
    

2 Factorielles

Les choix algorithmiques que vous ferez influent sur la qualité du programme en terme de complexité et donc de vitesse d'exécution.

2.1 Petites factorielles

Ecrivez un programme itératif faisant afficher toutes les factorielles des entiers de 0 à 20, une par ligne.

  1. En utilisant la fonction factorial, combien votre programme fera-t-il de multiplications ?
    print('Voici les factorielles des entiers de [1,20] :')
    i = 1                 
    while (i <= 20) {       
      ## la boucle d'iteration effectue en tout 1+2+3+...+20 = 210 multiplications
      print(paste(i,'! =', factorial(i)))
      i = i + 1
    }
    ## Si l'on remplace 20 par N, on obtiendrait N(N+1)/2 multiplications, polynome du
    ## second degre. On dit que l'algorithme a un cout QUADRATIQUE. Pas fameux...
    
    [1] "Voici les factorielles des entiers de [1,20] :"
    [1] "1 ! = 1"
    [1] "2 ! = 2"
    [1] "3 ! = 6"
    [1] "4 ! = 24"
    [1] "5 ! = 120"
    [1] "6 ! = 720"
    [1] "7 ! = 5040"
    [1] "8 ! = 40320"
    [1] "9 ! = 362880"
    [1] "10 ! = 3628800"
    [1] "11 ! = 39916800"
    [1] "12 ! = 479001600"
    [1] "13 ! = 6227020800"
    [1] "14 ! = 87178291200"
    [1] "15 ! = 1307674368000"
    [1] "16 ! = 20922789888000"
    [1] "17 ! = 355687428096000"
    [1] "18 ! = 6402373705728000"
    [1] "19 ! = 121645100408832000"
    [1] "20 ! = 2432902008176640000"
    
  2. Sans utiliser cette fonction. Tâchez de faire baisser le nombre de multiplications ! Combien votre programme fera-t-il de multiplications ?
    ## Si l'on veut faire tomber le nombre de multiplications, il ne FAUT PAS
    ## utiliser une fonction factorielle qui repart chaque fois du debut, mais maintenir
    ## dans une variable f la derniere factorielle calculee !
    print('Voici les factorielles (rapides) des entiers de [1,5] :')
    i = 1
    f = 1
    while (i <= 5) {
      ## la boucle effectue en tout 5 multiplications !!!
      print(paste(i,'! =', f))
      i = i + 1
      f = f * i
    }
    ## Le cout d'affichage des factorielles des entiers de [1,N] est donc maintenant de N, 
    ## le cout est devenu LINEAIRE. Difficile de faire mieux dans le cas present...
    
    [1] "Voici les factorielles (rapides) des entiers de [1,5] :"
    [1] "1 ! = 1"
    [1] "2 ! = 2"
    [1] "3 ! = 6"
    [1] "4 ! = 24"
    [1] "5 ! = 120"
    

2.2 Fonctions factorielle

Programmez plusieurs fonctions prenant un entier positif n et retournant la factorielle n! de n. Par exemple, 5! est égal à 120.

  1. Par récurrence, sous la forme d'une fonction FactRec(n), en supposant le problème résolu pour n-1. Vérifier FactRec(5) grâce à la fonction factorial.
    FactRec <- function(n) {
      if(n > 1) return(n*FactRec(n-1))
      else return (1)
    }
    paste('FactRec(-1) =',FactRec(-1))
    paste('FactRec(0) =',FactRec(0))
    paste('FactRec(5) =',FactRec(5))
    
    [1] "FactRec(-1) = 1"
    [1] "FactRec(0) = 1"
    [1] "FactRec(5) = 120"
    
  2. Par une itération (boucle while) sous la forme d'une fonction Fact(n).
    Fact <- function(n) {
      acc=1;
      while(n > 1) {
        acc = acc * n
          n = n - 1
      }
      return(acc)
    }
    paste('Fact(-1) =',Fact(-1))
    paste('Fact(0) =',Fact(0))
    paste('Fact(5) =',Fact(5))
    
    [1] "Fact(-1) = 1"
    [1] "Fact(0) = 1"
    [1] "Fact(5) = 120"
    
  3. R est-il capable de calculer 1000! ?
    factorial(1000)
    
    [1] Inf
    
  4. Pour chronométrer un calcul en R, il suffit d'importer la fonction system.time : system.time(replicate(1000,factorial(100))). Qui est le plus rapide entre Fact(100) et FactRec(100) ?
    FactVec <- function(n) {
      return(prod(1:n))
    }
    for(i in c(50, 100)) {
      print(paste('fact(', i, ') =', factorial(i)))
      for(f in c(factorial, Fact, FactRec, FactVec)) {
        print(system.time(replicate(2000,do.call(f,list(i)))))
      }  
    }
    
    [1] "fact( 50 ) = 3.0414093201713e+64"
    utilisateur     système      écoulé 
          0.004       0.000       0.004 
    utilisateur     système      écoulé 
          0.057       0.000       0.057 
    utilisateur     système      écoulé 
          0.080       0.000       0.081 
    utilisateur     système      écoulé 
          0.006       0.000       0.006 
    [1] "fact( 100 ) = 9.33262154439422e+157"
    utilisateur     système      écoulé 
          0.005       0.001       0.005 
    utilisateur     système      écoulé 
          0.108       0.000       0.109 
    utilisateur     système      écoulé 
          0.159       0.000       0.160 
    utilisateur     système      écoulé 
          0.008       0.000       0.007
    

3 Formule d'Euler

Le cours de maths de seconde année vous démontrera peut-être la formule suivante : \(\sum_{n=1}^{+\infty} \frac{1}{n^2}=\frac{\pi^2}{6}\). Comment écririez-vous un programme R permettant d'en déduire une valeur approchée de \(\pi\) ?

PiEuler <- function(n) {
  ## calcul de pi avec n termes dans la somme
  i = 1
  ## acc represente le resultat en construction
  acc = 1           
  while (i < n) {
    i = i + 1
    acc = acc + 1/(i*i)
  }
  return(sqrt(6 * acc))
}
cat(sprintf('Calcul approximatif de pi=%.16f :\n',pi))
k <- 1
while (k <= 7){
  n <- 10**k
  cat(sprintf("PiEuler(%d) --> %.10f\n",n, PiEuler(n)))
  k <- k+1
}
Calcul approximatif de pi=3.1415926535897931 :
PiEuler(10) --> 3.0493616360
PiEuler(100) --> 3.1320765318
PiEuler(1000) --> 3.1406380562
PiEuler(10000) --> 3.1414971639
PiEuler(100000) --> 3.1415831043
PiEuler(1000000) --> 3.1415916987
PiEuler(10000000) --> 3.1415925581

4 Nombres premiers   KEY

Reprenons le cours pages 11-12 sur les nombres premiers. La fonction estPremier(n) programmée en page 12 retourne False si elle trouve un diviseur non trivial (dans [2,n-1]) de n. Nous voudrions connaître ce diviseur !

  1. Programmez avec une boucle une fonction PlusPetitDiv(n) prenant un entier n ≥ 2 et retournant son plus petit diviseur \(\geq 2\). Par exemple PlusPetitDiv(35) renvoie la valeur 5. Notez au passage que le résultat de PlusPetitDiv(n) est le plus petit facteur premier de n (pourquoi ?).
    PlusPetitDiv <- function(n) {
      ## le plus petit diviseur de n dans [2,n-1], avec n >= 2
      if (n %% 2 == 0) return(2)  # je regle le cas pair en premier
      ## donc je peux maintenant monter de 3 vers n, sur les impairs uniquement !
      d <- 3
      while(n %% d != 0) {
        ## deux fois plus rapide que le d=d+1 du cours. Bof...
        d <- d + 2       
      }      
      ## d == n dans le pire des cas, coût = au plus n operations
      return(d)                  
    }
    paste('PlusPetitDiv(77) =',PlusPetitDiv(77))
    paste('PlusPetitDiv(29) =',PlusPetitDiv(29))
    
    p <- 1000003
    paste('PlusPetitDiv(',p,') Version O(n)')
    system.time(replicate(5,PlusPetitDiv(p)))
    
    [1] "PlusPetitDiv(77) = 7"
    [1] "PlusPetitDiv(29) = 29"
    [1] "PlusPetitDiv( 1000003 ) Version O(n)"
    utilisateur     système      écoulé 
          1.231       0.001       1.234
    
  2. En déduire une nouvelle version en une ligne de la fonction EstPremier(n).
    EstPremier <- function(n) (n >= 2) && (PlusPetitDiv(n) == n)      
    
    paste('EstPremier(77) =',EstPremier(77))
    paste('EstPremier(29) =',EstPremier(29))
    
    [1] "EstPremier(77) = FALSE"
    [1] "EstPremier(29) = TRUE"
    
  3. Faites afficher sur une ligne tous les nombres premiers jusqu'à 100.
    cat("Voici les nombres premiers compris entre 3 et 100 :\n")
    n <- 3
    while (n <= 100) {
      if (EstPremier(n)) {
        cat(n)
        cat(" ")
      }
      n <- n + 2
    }
    cat("\n")
    
    Voici les nombres premiers compris entre 3 et 100 :
    3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 
    
  4. Il est souvent intéressant d'accélérer un algorithme. Accélérons donc PlusPetitDiv. Montrez par un raisonnement mathématique très simple que s'il n'existe aucun diviseur de n dans \([2,\sqrt{n}]\) alors PlusPetitDiv(n) == n, ce qui évite de chercher dans le gros intervalle \([\sqrt{n},n]\). Implémentez cette optimisation. Au passage, accélérez encore en évitant de considérer les nombres pairs qui à part 2, ne sont pas premiers.
    PlusPetitDiv <- function(n) {
      ## avec au plus sqrt(n) operations
      if (n %% 2 == 0 ) return(2)
      d <- 3
      ## acc to compute the d*d
      acc <- 9 
      while (n %% d != 0 ) {
        ## <==> d > sqrt(n) mais plus rapide !
        if(acc > n) return(n)    
        d <- d + 2
        acc <- acc + 4*(d-1)
      }
      return(d)
    }
    paste('PlusPetitDiv(77) =',PlusPetitDiv(77))
    paste('PlusPetitDiv(84201097) =',PlusPetitDiv(84201097),'donc il est premier !')
    
    paste('PlusPetitDiv(',p,') Version O(sqrt(n))')
    system.time(replicate(5,PlusPetitDiv(p)))
    
    [1] "PlusPetitDiv(77) = 7"
    [1] "PlusPetitDiv(84201097) = 84201097 donc il est premier !"
    [1] "PlusPetitDiv( 1000003 ) Version O(sqrt(n))"
    utilisateur     système      écoulé 
          0.003       0.000       0.003
    

5 Méthode des tangentes de Newton   HARD

5.1 Formule de Newton

La formule de Newton [cours 2 page 16] pour améliorer une approximation \(a\) de \(\sqrt{r}\) est obtenue de la manière suivante. On trace la courbe d'équation \(y = f(x) = x^2-r\) qui coupe Ox précisément en \(\sqrt{r}\). formule_newton.png

  1. Quelle est l'équation de la tangente (T) à la courbe au point d'abscisse \(a\) ?
  2. La tangente n'étant pas horizontale, calculez l'abscisse b du point d'intersection de (T) avec Ox. Vous devez retrouver la formule de Newton !

5.2 Généralisation

Les calculettes modernes possèdent souvent une touche Solve permettant de calculer une racine d'une équation \(f(x)=0\). Par exemple, il est difficile sans machine de trouver une solution réelle à l'équation \(x^5-3x+1=0\).

  1. Pourquoi sommes-nous sûrs qu'il y en a au moins une ?
  2. Nous allons faire abstraction de la fonction , la supposer dérivable et à dérivée non nulle en tout point (de sorte que la tangente à la courbe existe et n'est jamais horizontale) pour appliquer la méthode des tangentes de Newton vue ci-dessus dans un cas particulier …

On suppose que l'approximation courante est \(a > 0\). Calculez l'équation de la droite tangente à la courbe de \(f\) au point \((a,f(a))\).

  1. Calculez la valeur de \(b\) qui est une amélioration de \(a\).
  2. Ecrivez la condition pour que l'approximation courante \(a\) soit suffisamment proche de la solution. On nommera \(h\) la constante > 0 de précision.
  3. Ecrivez une expression mathématique utilisant \(f\), \(a\) et \(h\), qui approche la valeur de la dérivée \(f^{\prime}(a)\) de \(f\) au point \(a\).
  4. Programmez en R la fonction Solve(f,a,h) retournant une approximation d'une solution de f en partant de l'approximation \(a\). L'argument \(h\) gouvernera la précision.
Solve <- function(f,a,h) {    
  ## une solution de f(x) = 0 en partant de a, et h gouverne la precision
  ## tant que la precision n'est pas atteinte...
  while (abs(f(a)) > h)   {      
    dfa = (f(a+h)-f(a))/h     # approximation de f'(a)
    a = a - f(a) / dfa        # amelioration de Newton...
  }
  return(a)                      # et hop !
}
g <-function(x) x**5 - 3 * x + 1
sol = Solve(g,1,0.0001)
print('x**5-3x+1 = 0 admet au moins une racine reelle car le degre est impair, ok ?')
paste('Une solution de x**5-3x+1 = 0 :',sol)
paste('Verification : g(sol) =',g(sol),'où e-05 signifie *10**(-5)')
[1] "x**5-3x+1 = 0 admet au moins une racine reelle car le degre est impair, ok ?"
[1] "Une solution de x**5-3x+1 = 0 : 1.2146493830778"
[1] "Verification : g(sol) = 1.05670215906351e-05 où e-05 signifie *10**(-5)"

5.3 Applications numériques à faire sur ordinateur

Utilisez la fonction Solve(f,a,h) pour faire afficher une valeur approchée :

  • de \(\sqrt{2}\),
  • puis de \(\sqrt[3]{2}\),
  • puis d'une solution de l'équation \(x^5-3x+1=0\),
  • puis de l'équation \(cos(x)=x\),
  • et enfin du nombre \(\pi\).
paste('Approximation de la racine cubique de 2 :',Solve( function(x) x**3-2, 1, 0.0001))
paste('(la "vraie" valeur est ',2**(1/3),')')
paste('Approximation de pi comme solution de sin(x) = 0 :',Solve(sin,3,0.0001))
paste('(la "vraie" valeur est ',pi,')')
sol = Solve((function(x) x-cos(x)),1,0.0001)
paste('Solution de cos(x) = x :',sol)
paste('Verification de cos(',sol,') == ',cos(sol), ':', sol == cos(sol))
[1] "Approximation de la racine cubique de 2 : 1.25993381738363"
[1] "(la \"vraie\" valeur est  1.25992104989487 )"
[1] "Approximation de pi comme solution de sin(x) = 0 : 3.14159265325442"
[1] "(la \"vraie\" valeur est  3.14159265358979 )"
[1] "Solution de cos(x) = x : 0.739113153543173"
[1] "Verification de cos( 0.739113153543173 ) ==  0.73906625809501 : FALSE"

N.B. Au moment d'utiliser la fonction Solve(f,a,h), il n'est pas nécessaire que la fonction f soit déjà définie. On peut passer une "fonction anonyme" en paramètre de Solve.

Par exemple, la fonction QUI N'A AUCUN NOM s'écrit en R :

function(x) x**2 – 1

On pourra donc demander par exemple

Solve((function(x) x**2 – 1), 3, 0.01)

6 Produit de nombres premiers   HOME

D'après le cours Python 6.00 du MIT. La théorie des nombres sait démontrer que si \(P_n\) dénote le produit de tous les nombres premiers jusqu'à \(n\), alors : \[ \lim_{n\rightarrow +\infty} \frac{P_n}{e^n}=1 \quad (\text{en maths } P_n \sim e^n) \] Mais calculer le produit d'un grand nombre d'entiers premiers va produire un très grand nombre ce qui peut nous occasionner des problèmes de calcul. Une manière d'y remédier est de remplacer les multiplications par des additions. Les logarithmes existent précisément dans ce but. Nous allons donc remplacer le produit des nombres premiers par la somme \(S_n\) de leur logarithme.

  1. Montrez que le théorème précédent s'exprime alors sous la forme : \(S_n \sim n\). Attention, ne croyez pas que \(f \sim g\) entraîne toujours \(\log(f) \sim \log(g)\). Mais c'est vrai ici car elles tendent vers \(+\infty\). Tâchez de faire une petite démonstration si possible, ou parlez-en à votre prof de maths.
  2. La fonction logarithme (népérien) de R est log. Comment utiliseriez-vous l'ordinateur et le langage R pour vous convaincre du bien-fondé de ce théorème ?
S <-function(n) {
  ## somme des logarithmes des nombres premiers <= n
  p = 2 ; acc = 0
  while (p <= n){
    acc = acc + log(p)
    p = Suivant(p)      
    ## je saute de p au nombre premier suivant. Voir en-dessous...
  }
  return(acc)
}

Suivant <- function(p) {
  ## retourne le plus petit nombre premier apres p
    essai=ifelse(p %% 2 == 0, p+1, p+2)
    while ( ! EstPremier(essai)) {
      ## je suis certain d'en trouver, ok ?
      essai = essai + 2
    }
    return(essai)
}
paste('Suivant(13) =',Suivant(13))

print('Je regarde si S(n)/n tend bien vers 1 :')
 for (i in 1:5) {
      n = 10**i
      somme = S(n)
      print(paste('n =',n,':',somme/n))
  }
  print('Mmmm. Oui, on dirait bien que la limite vaut 1.')
[1] "Suivant(13) = 17"
[1] "Je regarde si S(n)/n tend bien vers 1 :"
[1] "n = 10 : 0.534710753071747"
[1] "n = 100 : 0.837283903990639"
[1] "n = 1000 : 0.956245265120059"
[1] "n = 10000 : 0.989599137915698"
[1] "n = 1e+05 : 0.996853892686125"
[1] "Mmmm. Oui, on dirait bien que la limite vaut 1."

7 Fonction "Lambda"   HARD HOME

  1. Définissez en R la fonction \(f(x)=\frac{sin(x)}{\sqrt{x^4+1}}\).
  2. Calculez une valeur approchée de la dérivée seconde \(f^{\prime\prime}(\sqrt{2})\). Réponse : 0.036705…
## Prenons une fonction f particuliere :
f <- function(x) sin(x)/sqrt(x**4 + 1)

## Calcul de la derivee premiere de f quelconque en un point x avec une precision de h
Deriv <- function(f,x,h) (f(x + h) - f(x)) / h          
## formule vue au lycee...

paste('f\'(sqrt(2)) =',Deriv(f,sqrt(2),0.0001))

## Probleme : pour exprimer que la derivee seconde est la derivee de la derivee, j'ai
## besoin d'obtenir la FONCTION DERIVEE f' et pas seulement sa valeur en un point f'(x).
## Pour cela on peut utilise la notation lambda en python. 
## En R, on définit une sous-fonction a l'interieur de deriv et en rendant en resultat cette sous-fonction.
## C'est ce que les specialistes nomment une "fermeture" :
## http://fr.wikipedia.org/wiki/Fermeture_(informatique)

Deriv <- function(f,h) return (function(x) (f(x + h) - f(x)) / h)            
## Du coup, je peux exprimer mes maths sans etat d'ame :
df = Deriv(f,0.0001)
d2f = Deriv(df,0.0001)
paste('f"(sqrt(2)) =',d2f(sqrt(2)))    # et zou ! Intellectuel non ?...
[1] "f'(sqrt(2)) = -0.430032458566543"
[1] "f\"(sqrt(2)) = 0.0367057773065227"

Si cet exercice vous passe par-dessus la tete, ce n'est pas si grave que cela. Vous y reviendrez plus tard ! simple question de maturité…

Created: 2017-10-04 mer. 10:02