// Exercices de la série 1 du cours d'OS AM : boucles, sommes produits et récurrences.

funcprot(0); // pour ne pas afficher le message : "redefining function".
format(20);  // Pour afficher 20 chiffres par défaut

function ex01()
//=============  
// 
Connue = 0;
for kk = 1:100 do
  Connue = Connue + kk;
  end;

mprintf("La somme vaut : %8d\n", Connue);
endfunction
//ex01();

function ex02()
//=============  
// Sommes, curieuses conincidences.
// Les sommes suivantes sont presque identiques.
// Elles se rapprochent de plus en plus vers ln(10)
s1 = 0;
for kk = 10:99 do
  s1 = s1 + 1/kk;
  end;
mprintf("s1 = %12.8f\n", s1);

s2 = 0;
for kk = 100:999 do
  s2 = s2 + 1/kk;
  end;
mprintf("s2 = %12.8f\n", s2);

s3 = 0;
for kk = 1000:9999 do
  s3 = s3 + 1/kk;
  end;
mprintf("s3 = %12.8f\n", s3);

s4 = 0;
for kk = 10000:99999 do
  s4 = s4 + 1/kk;
  end;
mprintf("s4 = %12.8f\n", s4);

s5 = 0;
for kk = 100000:999999 do
  s5 = s5 + 1/kk;
  end;
mprintf("s5 = %12.8f\n", s5);

mprintf("ln(10) = %12.8f\n", log(10));

mprintf("10^1*(s2-s1) = %14.10f\n", (s2-s1)*10);
mprintf("10^2*(s3-s2) = %14.10f\n", (s3-s2)*100);
mprintf("10^3*(s4-s3) = %14.10f\n", (s4-s3)*1000);
mprintf("10^4*(s5-s4) = %14.10f\n", (s5-s4)*10000);
endfunction
//ex02();


function ex03()
//=============  
// Une somme convergente
s1 = 0;
for kk = 1:50 do
  s1 = s1 + 1/(2*kk-1) - 1/(2*kk);
  end;
mprintf("s1 = %12.8f\n", s1);

s2 = 0;
for kk = 1:500 do
  s2 = s2 + 1/(2*kk-1) - 1/(2*kk);
  end;
mprintf("s2 = %12.8f\n", s2);

s3 = 0;
for kk = 1:5000 do
  s3 = s3 + 1/(2*kk-1) - 1/(2*kk);
  end;
mprintf("s3 = %12.8f\n", s3);

mprintf("ln(2) = %12.8f\n", log(2));
endfunction
//ex03();


function ex04()
//=============  
// Une somme calculable à la main
simple = 0;
for kk = 1:99 do
  simple = simple + 1/(kk * (kk+1));
  end;
mprintf("simple = %12.8f\n", simple);

endfunction
//ex04();


function ex05()
//=============  
// Curieux produits
p1 = 1;
for kk = 2:100000 do
  p1 = p1 * (1 - 1/(kk * kk));
  end;
mprintf("p1 = %12.8f\n", p1);

p2 = 1;
for kk = 1:100 do
  p2 = p2 * (1 + 1/(2*kk-1)) * (1 - 1/(2*kk));
  end;
mprintf("p2 = %12.8f\n", p2);
// Quelle que soit la limite supérieure, le produit donne toujours 1 !

p3 = 1;
for kk = 1:100000 do
  p3 = p3 * (1 - 1/(2*kk+1)^2);
  end;
mprintf("p3 = %12.8f\n", p3);
mprintf("pi/4 = %12.8f\n", %pi/4);

p4 = 1;
for kk = 2:100 do
  p4 = p4 * (1 - 2/(kk*(kk+1)));
  end;
mprintf("p4 = %12.8f\n", p4);
// Après des tests en quadruple précision (Big_numbers/DD_QD_0.2/dd_0020_sommes_produits.sce)
// le résultat est exactement 0.3334,  avec juste le nombre de "3" qui change
// en fonction de la puissance de 10.
// En mettant au dénominateur commun chaque facteur, puis en les factorisant,
// on remarque que presque tout se simplifie, il reste :  (Max + 2) / (3*Max). Max = 10^n

p5 = 1;
for kk = 2:10000 do
  p5 = p5 * (1 - 2/(kk^3+1));
  end;
mprintf("p5 = %12.8f\n", p5);

p6 = 1;
for kk = 1:100000 do
  p6 = p6 * (1 - 0.25/(kk * kk));
  end;
mprintf("p6 = %12.8f\n", p6);
mprintf("2/pi = %12.8f\n", 2/%pi);

endfunction
//ex05();


function ex06()
//=============  
// Suite de Fibonacci, étude
// On constate que le rapport de deux termes consécutifs de la suite
// tend vers le nombre d'or

gold = (1 + sqrt(5))/ 2; // nombre d'Or

fibnm1 = 1;
fibn = 1;
for nn=2:30
  fibnp1 = fibn + fibnm1;
  rap = fibnp1 / fibn;
  fibnm1 = fibn; // passage de nn à nn+1
  fibn   = fibnp1;
  mprintf("fib(%3d) = %9d,  rapport = %18.15f,   fib-Or = %7.3e\n", nn, fibnp1, rap, rap - gold);
  end;

mprintf("                          gold = %18.15f\n", gold);
endfunction
//ex06();


function ex07()
//=============  
// Calcul de racine carrée d'un nombre par itération

aa = 2; // nombre dont on veut calculer la racine carrée.
racine_aa = sqrt(aa);

approx = aa/2;
approx_last = 0;
nn = 0;
while (abs(approx - approx_last) > 1e-15) do
  nn = nn + 1;
  approx_last = approx;
  approx = (approx + aa/approx) / 2;
  mprintf("nn = %2d,  approx = %18.15f,   différence = %7.3e\n", nn, approx, approx - racine_aa);
  end;

mprintf("approx^2 = %18.15f\n", approx^2);
endfunction
//ex07();


function ex07b()
//=============  
// Calcul de racine cubique d'un nombre par itération

aa = 2; // nombre dont on veut calculer la racine cubique.
racine_aa = aa^(1/3);

approx = aa/2;
approx_last = 0;
nn = 0;
while (abs(approx - approx_last) > 1e-15) do
  nn = nn + 1;
  approx_last = approx;
  approx = (2*approx + aa/approx^2) / 3;
  mprintf("nn = %2d,  approx = %18.15f,   différence = %7.3e\n", nn, approx, approx - racine_aa);
  end;

mprintf("approx^3 = %18.15f\n", approx^3);
endfunction
//ex07b();


function ex08()
//=============  
// Un curieux test, qui montre la différence entre les nombres en 
// informatique et en mathématique.

ee = 1;

nn = 0;
while (1 + ee > 1) do
  nn = nn + 1;
  ee = ee / 2;
  mprintf("nn = %2d,  ee = %9.3e\n", nn, ee);
  end;

endfunction
//ex08();


function ex08b()
//==============  
// Un curieux test, qui montre la différence entre les nombres en 
// informatique et en mathématique.

ee = 1;

nn = 0;
while (ee > 0) do
  nn = nn + 1;
  ee = ee / 2;
  mprintf("nn = %2d,  ee = %9.3e\n", nn, ee);
  end;

endfunction
//ex08b();


function ex09()
//=============  
// Étude du dernier chiffre de la suite de Fibonacci
  
  function y=fib(n)
  //===============  
  // Retourne les trois derniers chiffres du nombre  f_n  de la suite de Fibonacci.
  
  fibnm1 = 1;
  fibn = 1;
  for kk=2:n
    fibnp1 = modulo(fibn + fibnm1, 1000);
    fibnm1 = fibn; // passage de nn à nn+1
    fibn   = fibnp1;
    end;
  
  y = fibn;
  endfunction

for kk=1:100 do
  mprintf("fib(%3d) =%4d\n", kk, fib(kk));
  end;

kk = 1000;  mprintf("fib(%3d) =%4d\n", kk, fib(kk));
kk = 10000;  mprintf("fib(%3d) =%4d\n", kk, fib(kk));
kk = 100000;  mprintf("fib(%3d) =%4d\n", kk, fib(kk));
// On constate que le dernier chiffre est toujours  1.
// La périodicité du dernier chiffre du nombre de Fibonacci est de 60
// Donc  fib(kk)  et  fib(kk - 60*entier)  ont le même chiffre des unités.
// fib_unite(100) = fib_unite(40) = 1
// fib_unite(1000) = fib_unite(1000- 16*60) = fib_unite(40) = 1
// fib_unite(10000) = fib_unite(10000- 166*60) = fib_unite(40) = 1
// etc.
endfunction
//ex09();


function ex10()
//=============  
// Somme des inverses des entiers positifs, ne coportant pas le chiffre "9"
// La somme vaut :
// s90 + s91 + s92 + s93 + s94 + s95 + ...  qui vaut environ
// s90 + s91 + s92 + s93 + s94 + s95 + s96 + s96*0.9 + s96*0.9^2 + s96*0.9^3 + s96*0.9^4 + ...  =
// s90 + s91 + s92 + s93 + s94 + s95 + s96 * (1 + 0.9 + 0.9^2 + 0.9^3 + 0.9^4 + ...) =
// s90 + s91 + s92 + s93 + s94 + s95 + s96 * 10 =         22.920678 549999998
// s90 + s91 + s92 + s93 + s94 + s95 + s96 + s97 * 10 =   22.920676 794415
// s90 + s91 + s92 + s93 + s94 + ... + s97 + s98 * 10 =   22.920676 635118
// s90 + s91 + s92 + s93 + s94 + ... + s98 + s99 * 10 =   22.920676 620689
// s90 + s91 + s92 + s93 + s94 + ... + s99 + s910* 10 =   22.920676 61939
// Des calculs plus poussés donnent :             somme = 22.920676 619264 150 ...

s90 = 2.717857142857143;
s91 = 2.053991622224917;
s92 = 1.817871425200977;
s93 = 1.633364212583175;
s94 = 1.469783389239973;
s95 = 1.322783057766752;
s96 = 1.190502772693381;
s97 = 1.071452317287522;

disp(s91 / s90); // 0.75573936163546851
disp(s92 / s91); // 0.88504325543450857
disp(s93 / s92); // 0.89850370221176756
disp(s94 / s93); // 0.89985036889935399
disp(s95 / s94); // 0.89998503687728093
disp(s96 / s95); // 0.89999850368759693
disp(s97 / s96); // 0.89999985036866370
disp(s90 + s91 + s92 + s93 + s94 + s95 + s96*10);
disp(s90 + s91 + s92 + s93 + s94 + s95 + s96 + s97*10);

// Probablement que :
s98 = s97 * 0.89999998503686;
disp(s90 + s91 + s92 + s93 + s94 + s95 +  s96 + s97 + s98*10);

// Probablement que :
s99 = s98 * 0.899999998503686;
disp(s90 + s91 + s92 + s93 + s94 + s95 +  s96 + s97 + s98 + s99*10);

// Probablement que :
s910 = s99 * 0.8999999998503686;
disp(s90 + s91 + s92 + s93 + s94 + s95 +  s96 + s97 + s98 + s99 + s910*10);

s91 = 0;
//  Erreur après plusiques dizaines de minutes, si nn = 10000000;
nn = 100;
nMax = 10*nn-1
while (nn <= nMax) do
  s91 = s91 + 1/nn;
  nn = nn + 1;
  if (modulo(nn, 10) == 9) then // saute les nombres terminant par 9
     nn = nn + 1;
     if (modulo(nn/10, 10) == 9) then // saute les nombres dont le chiffre des dizaines vaut 9
       nn = nn + 10;
        if (modulo(nn/100, 10) == 9) then // saute les nombres dont le chiffre des centaine vaut 9
          nn = nn + 100;
          if (modulo(nn/1000, 10) == 9) then // saute les nombres dont le chiffre des milliers vaut 9
            nn = nn + 1000;
            if (modulo(nn/10000, 10) == 9) then // saute les nombres dont le chiffre des 10^4 vaut 9
              nn = nn + 10000;
              if (modulo(nn/100000, 10) == 9) then // saute les nombres dont le chiffre des 10^5 vaut 9
                nn = nn + 100000;
                if (modulo(nn/1000000, 10) == 9) then // saute les nombres dont le chiffre des 10^6 vaut 9
                  nn = nn + 1000000;
                  end;
                  if (modulo(nn/10000000, 10) == 9) then // saute les nombres dont le chiffre des 10^7 vaut 9
                    nn = nn + 10000000;
                    end;
                end;
              end;
            end;
          end;
       end;
     end;
   end;

mprintf("nn = %2d,  s9 = %19.15f\n", nn, s91);
endfunction
ex10();


