Arithmétique avec CAML

Voici quelques fonctions Caml relatives au cours d'arithmétique:

OCaml
let rec quo(a,b) =
  if b = 0 then raise Division_by_zero
  else
    if b > 0 then
      if a >= 0 then
        if a<b then 0
        else 1 + quo(a-b,b)
      else
        if b*quo(-a,b) = -a then
          -quo(-a,b)
        else
          -1 -quo(-a,b)
    else
      - quo(a,-b);;
 
let rec rem(a,b) =
  if b = 0 then raise Division_by_zero
  else
    if b > 0 then
      if a >= 0 then
        if a<b then a
        else rem(a-b,b)
      else
        if rem(-a,b) = 0 then
          0
        else
          b - rem(-a,b)
    else
      rem(a,-b);;
 
 
let (%) a b = rem(a,b);;
 
 
 
(* exo 2-5 Nb chiffres d'un entier en décimal *)
 
let rec nb_chiffres(n) =
  if n / 10 = 0 then
    1
  else
    1 + nb_chiffres(n/10);;
 
(* 
exo 2-6 
Conversion de n (entier positif en décimal) 
en la liste de ses chiffres dans la base b (1 < b <= 10)
*)
 
exception Entier_negatif;;
 
let base(n,b) = 
let rec base_local(k,ba,acc) = 
  if k < 0 then 
    raise Entier_negatif
  else
    if k < ba then
      acc @ [k]
    else
      base_local(k/ba,ba,acc @ [k mod ba])
in base_local(n,b,[]);;
 
 
 
 
type 'a classe = {modulo:'a ; classe: 'a};;
 
exception  Classes_incompatibles;;
 
let iso op_int a b =
  {modulo =
      if a.modulo= b.modulo
      then a.modulo 
      else raise Classes_incompatibles ;
   classe =
      ((op_int) (a.classe) (b.classe)) % (a.modulo)
  };; 
 
let ( +% ) = iso  ( + ) ;;
let ( *% ) = iso  ( * ) ;;
 
{modulo =5; classe = 9%5} +% {modulo=5; classe =12%5};;
 
 
let cl n a = {modulo = n; classe = (a%n)};;
 
cl 5 12;;
 
let c_5 = cl 5;;
 
c_5 12;;
 
(c_5 12) +% (c_5 9);;
(c_5 12) *% (c_5 9);;
 
12*9 % 5;;
 
 
 
(*-----------  exo 2-7 -----------*)
 
(* étape 1 : n -> liste de ses chiffres en base 10 *)
let liste_chiffres = function
  |n -> base(n,10);; 
 
(* étape 2 : produit des éléments d'une liste *)
exception Liste_vide;;
 
let produit(liste) = 
let rec produit_local = function
  |[],acc -> raise Liste_vide
  |tete::[],acc -> tete*acc
  |tete::queue,acc -> produit_local(queue,tete*acc)
in produit_local(liste,1);;
 
(* étape 3 : calcul de la persistance d'un entier *)
let persistance(n) = 
let rec  pers_local(n,acc) =
  if n >= 0 then 
    if n < 10 then
      acc
    else
       pers_local(produit(liste_chiffres(n)),1+acc)
  else
    pers_local(-n,acc) 
in pers_local(n,0);;
 
(* étape 4 : recherche du + petit entier 
             de persistance donnée  (impératif)     *)
 
let petit_pers1(p) = 
  let k = ref 0 in
    while persistance(!k) < p do
      k := !k + 1;
    done;
  !k ;;
 
 
(* étape 4-bis : recherche du + petit entier 
             de persistance donnée  (récursif)     *)
 
let petit_pers2(p) =
  let rec  petit_local(n,debut) =
    if persistance(debut) = n then
      debut
    else
      petit_local(n,debut+1)
  in petit_local(p,0);;
 
 
 
 
exception Date_yennapas_exister;;
exception Date_avant_petit_JC;;
 
(*----------- Exo 2-9 -------------------------*)
 
(* prévoir aussi un vérificateur de bonne date *)
 
let zeller(j,m,a) = 
  if (a=1582) && (m=10) && ((j>4)&&(j<15)) then
    raise Date_yennapas_exister
  else
    if a <= 0 then
      raise Date_avant_petit_JC
    else
      let semaine = [|"Samedi";"Dimanche";"Lundi";"Mardi";"Mercredi";"Jeudi";"Vendredi"|] in
      let mm = if m < 3 then m+12 else m in
        (* décalage éventuel  du no de mois  *)
      let aa = if m < 3 then a-1 else a in
        (* décalage éve1ntuel  du no d'année  *)
      let na = aa mod 100 (* numéro d'année *) in 
      let s = aa / 100 (* numéro de siècle *) in
      let greg = if ((a<1582)||((a=1582)&&((m<10)||((m=10)&&(j<5))))) then (5-s) else (s/4 -2*s)  in 
        (* résidu à ajouter selon  gregorien/julien *)
      let nj = j + na + na/4 +(26*(mm+1))/10 +greg in
        semaine.(rem(nj,7)) ;;
 
    
 
 
 
 
 
 
(* pgcd en récursif  *)
let rec pgcd(a,b)=
  if a >= b then
    if b = 0 then
      a
    else
      pgcd(b,a mod b)
  else
    pgcd(b,a);;
 
    
 
(* pgcd en itératif *)
 
let pgcdit(a,b)=
  let aa =  if a >= b then ref a else ref b  in
  let bb =  if a >= b then ref b else ref a in
  let temp = ref 0 in
    while rem(!aa,!bb)  >  0 do
      temp := !aa ;
      aa := !bb ;
      bb := rem(!temp,!bb) ;
    done;
   !bb
;;
 
 
(* Euclide étendu : coeff de Bézout *)
 
 
let rec euclide2 a b =
  let rec euc(u,v,r,u',v',r') =
    if r' = 0 then
      [|u;v;r|]
    else
      let q = r/r' in
      euc(u',v',r',u - q * u', v - q * v', r - q * r')
  in euc(1,0,a,0,1,b);;
 
 
 
(* Nb de diviseurs d'un entier *)
 
(* récursif naïf :*)
 
let rec nb_div_naif(n) =
  let rec nb_div_naif_loc(nb_a_diviser,candidat) =
    if candidat > nb_a_diviser then
      0
    else
      if nb_a_diviser mod candidat = 0 then
        1 + nb_div_naif_loc(nb_a_diviser,candidat+1)
      else
        nb_div_naif_loc(nb_a_diviser,candidat+1)
  in nb_div_naif_loc(n,1);;
 
(* Impératif naïf *)
 
let nb_div_it(n) =
  let k = ref 0 in
    for i=1 to n do
      if n mod i = 0 then 
        k := !k + 1;
    done;
    !k ;;
 
(* Récursif amélioré *)
 
let rec nb_div(n) =
  let rec nb_div_loc(nb_a_diviser,candidat) =
    if candidat*candidat > nb_a_diviser then
      0
    else
      if candidat*candidat = nb_a_diviser then
        1 
      else
        if nb_a_diviser mod candidat = 0 then
          2 + nb_div_loc(nb_a_diviser,candidat+1)
        else
          nb_div_loc(nb_a_diviser,candidat+1)
  in nb_div_loc(n,1);;
 
 
 
 
 
(*----- exo 2-12 --------*)
 
(* liste des diviseurs *)
 
let rec liste_div(n) = 
  let rec liste_div_loc(m,candidat) = 
    if candidat*candidat > m then
      [] 
    else
      if (candidat*candidat = m) then
        (*pas besoin d'aller voir plus loin *)
        [candidat] 
      else
        if m mod candidat = 0 then
          (*on rajoute candidat et m/candidat à la liste *)
          (m/candidat)::candidat::liste_div_loc(m,candidat+1)
        else 
           (*on passe au suivant sans rien ajouter*)
          liste_div_loc(m,candidat+1)
  in liste_div_loc(n,1);; 
 
(* somme des éléments d'une liste *)
(* laisse une impression de déjà-vu... *)
 
exception Liste_vide;;
 
let rec som_liste = function
  |[] -> raise Liste_vide
  |tete::[] -> tete
  |tete::queue -> tete + som_liste(queue) ;;
 
let parfait(n) = 
  if som_liste(liste_div(n)) = n+n  then 
    true 
  else 
    false ;;
 
let cherche_parfait(max) =
  let rec cherche_local(candidat,borne) =
    if candidat > borne then
      []
    else if parfait(candidat) then
      candidat::cherche_local(candidat+1,borne)
    else
      cherche_local(candidat+1,borne)
  in cherche_local(1,max) ;;
 
 
 
(*----Plus efficace------*)
 
 
(* Calcule directement la somme des diviseurs de n *)
 
let rec som_div1(n) = 
  let rec som_div_loc(nb,candidat) = 
   if candidat*candidat > nb then
     0
   else
     if (candidat*candidat = nb) then
        candidat 
     else
      if nb mod candidat = 0 then
       candidat + (nb/candidat) + som_div_loc(nb,candidat+1)
      else
          som_div_loc(nb,candidat+1)
  in som_div_loc(n,1);;
 
 
(* Version terminale : plus efficace *)
 
let rec som_div(n) =  
  let rec som_div_loc(nb,candidat,acc) = 
   if candidat*candidat > nb then
     acc
   else
     if (candidat*candidat = nb) then
        candidat+acc 
     else
      if nb mod candidat = 0 then
        som_div_loc(nb,candidat+1,acc+candidat+(nb/candidat))
      else
          som_div_loc(nb,candidat+1,acc)
  in som_div_loc(n,1,0);;
 
 
 
 
(* test de perfection *)
 
let parfait(n) = 
  if som_div(n) = 2*n  then 
    true 
  else 
    false ;;
 
 
 
 
(* Recherche des parfaits <= max donné en argument *)
 
 
let cherche_parfait1(max) =
  let rec cherche_local(candidat,borne) =
    if candidat > borne then
      []
    else if parfait(candidat) then
      candidat::cherche_local(candidat+1,borne)
    else
      cherche_local(candidat+1,borne)
  in cherche_local(1,max) ;;
 
 
 
 
 
 
(* récursion terminale tout en 1  plus rapide *)
 
let cherche_parfait(max) =
let rec som_div(n) =  
  let rec som_div_loc(nb,candidat,acc) = 
   if candidat*candidat > nb then
     acc
   else
     if (candidat*candidat = nb) then
        candidat+acc 
     else
      if nb mod candidat = 0 then
        som_div_loc(nb,candidat+1,acc+candidat+(nb/candidat))
      else
          som_div_loc(nb,candidat+1,acc)
  in som_div_loc(n,1,0) in 
  let rec cherche_local(candidat,borne,acc) =
    if candidat > borne then
      acc
    else if som_div(candidat)=2*candidat then
      cherche_local(candidat+1,borne,candidat::acc)
    else
      cherche_local(candidat+1,borne,acc)
  in cherche_local(1,max,[]) ;;
 
 
 
 
(*------------------------PGCD & Co-------------------------------*)
 
(*------Algo d'Euclide I / calcul du PGCD version récursive-------*)
 
let rec pgcd(a,b) =
  if a >= b then
    if b = 0 then
      a
    else 
      pgcd(b,a mod b)
  else
    pgcd(b,a);;
 
(*----- Avec espion pour compter les récursions----*)
 
let pgcd_spy(a,b) =
  let rec pgcd(a,b,spy) =
    if a >= b then
      if b = 0 then
        (a,spy)
      else 
        pgcd(b,a mod b,spy+1)
    else
      pgcd(b,a,spy+1)
  in pgcd(a,b,0);;
 
 
(*------Fonction qui teste si 2 entiers sont 1ers entre eux--------*)
 
let premiers_entre_eux(a,b) =
  if pgcd(a,b) = 1 then
    true
  else
    false;;
 
 
(*-Euclide II / calcul du PGCD + coeff de Bézout version récursive-*)
 
let bezout(a,b) =
  let rec euc(u,v,r,u',v',r') =
    if r' = 0 then
      [|u;v;r|]
    else
      let q = r/r' in 
        euc(u',v',r',u-q*u',v-q*v',r-q*r')
  in  euc(1,0,a,0,1,b);;
 
(*---------Fonction qui recherche l'inverse de x modulo n---------*)
 
exception Non_inversible ;;
 
 
let inv_mod a =
  let b = euclide2 a.classe a.modulo in
  if b.(2)=1 then
    cl a.modulo b.(0)
  else
    raise Non_inversible;;
 
let inv_mod2 n a  =
  let b = euclide2 a n in
  if b.(2)=1  then
    cl n b.(0)
  else
    raise Non_inversible;;
 
 
 
 
 
      
 
(* ------- PRIME NUMBERS ----------*)
 
 
 
 
let rec est_premier1(n) =
  let rec test_local(t,k) =
    if t*t > k then
      true
    else
      if k mod t = 0 then
        false
      else
        test_local(t+1,k)
  in test_local(2,n);;
 
 
(* en enlevant les pairs *)
 
let rec est_premier2(n) =
  let rec test_local(t,k) =
    if (k=2) || (t*t > k) then
      true
    else
      if k mod t = 0 then
        false
      else
        test_local(t+2,k)
  in test_local(3,n);;
 
 
 
 
(* Ératosthène *)
 
 
(* On retire les multiple de n autres que n dans une liste*)
 
let rec retmultiples (l,n) =
  let rec ret_loc = function
    | [],n,acc -> acc
    | t::q,n,acc -> if (t mod n = 0) && (t > n) then 
                      ret_loc(q,n,acc)
                    else
                      ret_loc(q,n,t::acc)
  in ret_loc(l,n,[]);;
 
 
let rec retmultiples2 liste n =
  match liste with
  | [] -> []
  | tete::queue -> if tete mod n = 0 then 
                    retmultiples2 queue n
                  else
                    tete::(retmultiples2 queue n);;
 
 
let retmultiples3 liste n =
  let pred e = (e mod n != 0) in
  List.find_all pred liste;;
 
 
(*La liste des entiers de 2 à n*)
 
let  entiers(n) =
  let rec listentier(m,n,acc) =
    if m > n then
      acc
    else
      listentier(m+2,n,m::acc)
  in listentier(3,n,[2]);;
 
 
 
let retmultiples liste n =
  let pred e = (e mod n != 0) || (e=n) in
  List.find_all pred liste;;
 
 
(*Le crible : on barre les multiples de n sauf n jusqu'à atteindre
  la racine carrée de n *)
 
 
let crible(m)=
  let rec crible_rec n =
    if n>m then 
      []
    else
      if n*n>m then 
        n::(crible_rec (n+1))
      else 
        n::(retmultiples (crible_rec (n+1)) n)
  in crible_rec(2);;
 
let crible m=
  let rec crible_rec n acc =
    if n>m then 
      acc
    else
      if n*n>m then 
        crible_rec (n+1) (n::acc)
      else 
        retmultiples (crible_rec (n+1) (n::acc))  n
  in crible_rec 2 [];;
 
 
let crible m =
  let rec crible_rec n acc =
    match n with
      |n when n>m -> acc
      |n when n*n>m -> crible_rec (n+1) (n::acc)
      |_ -> retmultiples (crible_rec (n+1) (n::acc))  n
  in crible_rec 2 [];;
 
(*
crible 100;;
*)
 
 
 
 
 
 
 
 
 
(* lOngueur d'une liste pour compter les nb premiers inférieurs à n*)
 
let long liste =
  let rec long_loc liste acc =
    match liste with
    |[] -> acc
    |t::q -> long_loc q (1+acc)
  in long_loc liste 0;;
 
 
let pi n  = long (crible n);;
 
(* pi 1_000_000;; *)
 
(*
# long(crible(1_000_000));;
- : int = 78498
# long(crible(10_000_000));;
- : int = 664579
 
moi@moi-bur:~/IUT/INFO1/CAML$  ./a.out
  5.37  s pour trouver le nb de premiers 78498 
moi@moi-bur:~/IUT/INFO1/CAML$ ocamlopt arith_info1.ml 
moi@moi-bur:~/IUT/INFO1/CAML$  ./a.out
139.53  secondes pour trouver le nb de premiers < à  10000000 664579 
  
*)
 
 let time_sn fn n msg =
  let start = Sys.time() in
  let res = fn(n) in
  let stop = Sys.time() in
  Printf.printf "%6.2f  %s %n %n \n" (stop -. start) msg n res;; 
 
(*
  time_sn pi 1_000_000 "secondes pour trouver le nb de premiers < à ";;
*)
 
 
 
(*Décomposition en éléments simples (non terminale)*)
 
 
let decompo1(n)=
  let rec decompotemp = function
      | (0 | 1),l -> []
      | n,[] -> []
      | n,tete::queue ->
          if n mod tete = 0 
            then tete::decompotemp(n/tete,tete::queue)
            else decompotemp(n,queue)
  in 
  decompotemp(n,crible(n));; 
 
 
let decompo(n)=
  let rec decompotemp = function
      | (0 | 1),l,acc -> acc
      | n,[],acc -> acc
      | n,tete::queue,acc ->
          if n mod tete = 0 
            then decompotemp(n/tete,tete::queue,tete::acc)
            else decompotemp(n,queue,acc)
  in 
  decompotemp(n,crible(n),[]);; 
 
 
 
 
(* POur enlever les doublons et n'avoir que les diviseurs *)
 
let purge_liste(liste) =
  let rec loop  = function
  | [],acc -> acc
  | t::[],acc -> t::acc
  | a::b::queue,acc -> if a = b then
      loop(b::queue,acc)
    else
      loop(b::queue,a::acc)
  in loop(liste,[]);;
 
 
 
(*Indicatrice d'Euler : phi(k) = k * prod(p-1)/p *)
 
 
 
let euler(n) =
  let rec loop = function
    |[],acc1,acc2 -> n*acc1/acc2
    |t::q,acc1,acc2 -> loop(q,(t-1)*acc1,t*acc2)
  in loop(purge_liste(decompo(n)),1,1);;
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
(*La série d'exponentiations*)
 
 
 
let rec prap1 = function
  |(a,0) -> 1
  |(a,1) -> a
  |(a,n) -> if (n mod 2 = 0)
            then prap1(a*a,n/2)
            else a*prap1(a*a,(n-1)/2);;
 
 
(* version terminale *)
 
let prap(a,n) =
  let rec local = function
      |(a,0,acc) -> acc
      |(a,1,acc) -> a*acc
      |(a,n,acc) -> if (n mod 2 = 0)
            then local(a*a,n/2,acc)
            else local(a*a,(n-1)/2,a*acc)
  in local(a,n,1);;
    
 
let pmod(a,n,p) =
  let rec local = function
      |(a,0,p,acc) -> acc mod p
      |(a,1,p,acc) -> (a mod p) * (acc mod p)
      |(a,n,p,acc) -> if (n mod 2 = 0)
            then local((a mod p) * (a mod p),n/2,p,acc mod p)
            else local((a mod p) * (a mod p),(n-1)/2,p,(a mod p) * (acc mod p))
  in local(a,n,1,p);;
 
 
(*Utilisation des big zentiers*)
 
 
#load "nums.cma";;
 
open Big_int;;
 
let badd(a,b) =
    Big_int.add_big_int a b;;
 
(* la même en infixe  on fera a &+ b *)
 
let (&+) a b = Big_int.add_big_int a b;;
 
 
let bmul(a,b) =
    Big_int.mult_big_int a b;;
 
(* la même en infixe *)
 
let (&*) a b = Big_int.mult_big_int a b;;
 
let bmod(a,n) =
    Big_int.mod_big_int a n;;
 
(* la même en infixe *)
 
let (&%) a n = Big_int.mod_big_int a n;;
 
let bquo(a,b) =
    Big_int.div_big_int a b;;
 
(* la même en infixe *)
 
let (&/) a b = Big_int.div_big_int a b;;
 
let bsub(a,b) =
    Big_int.sub_big_int a b;;
 
(* la même en infixe *)
 
let (&-) a b = Big_int.sub_big_int a b;;
 
let begal(a,b) =
  Big_int.eq_big_int a b ;;
 
(* la même en infixe *)
 
let (&=) a b  =  Big_int.eq_big_int a b ;;
 
(* transforme un petit entier en grand entier *)
let big n =
  Big_int.big_int_of_int n;;
 
(* les petits entiers basiques *)
 
let zero = zero_big_int;;
 
let un = unit_big_int ;;
 
let deux = big(2) ;;
 
 
(* pour afficher un grand entier sous forme de chaîne *)
 
let montre(b) =
  string_of_big_int b;;
 
 
(* Le pgcd pour plus tard  *)
 
let bgcd a b =
  Big_int.gcd_big_int a b ;;
 
 
(* exponentiation rapide *)
 
let bprap(a,n) =
  let rec local(ba,bn,acc) = 
      if bn &= zero then acc
      else if bn &= un then ba &* acc
      else  if (bn &% deux) &= zero
            then local(ba &* ba, bn &/ deux , acc)
            else local(ba &* ba, (bn &- un) &/ deux ,ba &* acc)
  in local(a,n,un);;
 
 
(* let (&^) a n = bprap(a,n);; *)
 
let (&^) a n = power_big_int_positive_big_int a n;;
 
 
 
 
 
 
 
 
let big_iso op_big_int a b =
  {modulo =
      if a.modulo &= b.modulo
      then a.modulo 
      else raise Classes_incompatibles ;
   classe =
      ((op_big_int) (a.classe &% a.modulo) (b.classe &% a.modulo)) &% (a.modulo)
  };; 
 
let ( &&+ ) = big_iso  ( &+ ) ;;
let ( &&* ) = big_iso  ( &* ) ;;
let ( &&/ ) = big_iso  ( &/ ) ;;
let ( &&% ) = big_iso  ( &% ) ;;
 
 
 
let montre_classe b = montre b.classe;;
 
let bcl n a = {modulo = n; classe = (a &% n)};;
 
let bc_5 = bcl (big 5);;
 
montre_classe ( bc_5 (big 12));;
 
 
(* pour l'égalité *)
 
 
let (&&=)  a b =
  (a.modulo &= b.modulo) && ((a.classe &% a.modulo)  &= (b.classe &% a.modulo)) ;;
 
 
 
 
 
let mprap(a,n,p) =
  let rec local(ba,bn,acc) = 
      if bn &= zero then acc
      else if bn &= un then ba &&* acc
      else  if (bn &% deux) &= zero
            then local(ba &&* ba, bn &/ deux , acc)
            else local(ba &&* ba, (bn &- un) &/ deux ,ba &&* acc)
  in local({modulo = p; classe = a},n,{modulo = p ; classe = un});;
 
 
let (&&^) a (n,p) = (mprap(a,n,p)).classe;;
 
 
(deux &^ (big 1000)) &&^ ((big 1000),(big 10101));;
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
(* teste l'appartenance d'un elmt à une liste*)
 
let rec appartient = function
  | (el,[]) -> false
  | (el,tete::queue) -> if el = tete then true
                        else appartient(el,queue);;
 
 
(* teste si n est pseudo-premier de base a  
   avec a et b des entiers *)
let pseudo a n =
  let c = crible(n) in
  let rec local(a,n,i,acc) =
    if i >= n then acc
    else
      let bca = bcl (big i) (big a) in
      let bci = bcl (big i) (big (i-1)) in
      let p = bca.classe  &&^ (bci.classe,big i) in
      if (p &= un) && not(appartient(i,c)) then
        local(a,n,i+2,i::acc)
      else
        local(a,n,i+2,acc)
  in local(a,n,3,[]);;
 
 
 
pseudo 2 10000;;
 
 
(*-------------------CÉSAR---------------------- *)
 
 
 
let queue_chaine(chaine) = String.sub chaine 1 (String.length(chaine)-1);;
 
let tete_chaine(chaine) = chaine.[0];;
 
let string_of_char car =  String.make 1 car;;
 
let (%) a b = rem(a,b);;
 
let decalage(n,cle) = (n- 32 + cle) % 95 + 32;; 
 
let code(car,cle) = 
    string_of_char (char_of_int(decalage(int_of_char(car),cle)));;
 
let rec cesar(m,cle) =
    if m = "" then ""
    else
      (code(tete_chaine(m),cle))^cesar(queue_chaine(m),cle);;
 
let decalage_affine(n,a,b)= ((n-32)*a + b) % 95 + 32;;
 
let code_affine(car,a,b) = 
    string_of_char (char_of_int(decalage_affine(int_of_char(car),a,b)));;
 
let rec chiffre_affine(m,a,b) =
    if m = "" then ""
    else
      (code_affine(tete_chaine(m),a,b))^chiffre_affine(queue_chaine(m),a,b);;
 
let inv_decalage_affine(n,a,b)= (n - 32 - b)* (inv_mod {modulo = 95 ; classe =a}).classe % 95 + 32;;
 
 
let decode_affine(car,a,b) = 
    string_of_char (char_of_int(inv_decalage_affine(int_of_char(car),a,b)));;
 
 
let rec dechiffre_affine(m,a,b) =
    if m = "" then ""
    else
      (decode_affine(tete_chaine(m),a,b))^dechiffre_affine(queue_chaine(m),a,b);;
 
 
 
 
 
 
 
 
 
(*Rabin Miller *)
 
(*Grand entier aléatoire*)
exception Taille_negative;;
 
(* nombres aléatoires (d'après le module numerix de Michel Quercia) *)
(* land est le ET logique bit à bit : plus rapide que mod
   1 lsl p effectue un décalage de 1 bit : plus rapide que 2^p
   0xn est l'entier écrit n en base 16
   On travaille donc par paquet de 16^4
   renvoie un nombre < 2^n*)
let nrandom n =
  if n < 0 then raise Taille_negative
  else if n = 0 then big(0)
  else begin
    let p = if n land 15 = 0 then 16 else n land 15 in
    let mask = (1 lsl p) - 1 in
    let r = ref(big(((Random.int(0xfffe)) land (mask-2))+2)) in
    for i=1 to (n-p)/16 do
        r := add_int_big_int (Random.int(0x10000)) (mult_int_big_int 0x10000 !r)
    done;
    !r
  end
;;
 
 
 
 
 
 
 
 
 
  
(*Décomposition 2^t x s*)
 
let dec(n) =
  let rec loop(s,t) =
    if (s &% deux) &= zero then
      loop( s &/ deux , t &+ un )
    else
      [|s;t|]
  in loop(n &- un,zero);;
 
(* teste a^s mod n puis les v^(2^u)  *)
 
let carre(v,n,t) =
  if v &= un then
    true
  else
    let rec loop(v,u) =
      if v &= (n &- un)  then
        true
      else
        if u &= (t &- un) then
          false
        else
          loop( v &&^ (deux,n) , u &+ un )
    in loop(v,zero);;
 
 
 
 
(* renoie k, la partie entière de log_2(n) : 2^k <= n < 2^(k+1)
On devra considérer des nombres inférieurs à
   max_float ~ 10^(308) ~ 2^(1023)
*)
let log2(n) =
  int_of_float(log(float_of_big_int(n))/.log(2.));;
 
let blog2(n,k) =
  k*1024 + int_of_float(log(float_of_big_int( div_big_int n (power_int_positive_int 2 (k*1023))))/.log(2.));;
 
 
exception Nombre_pair;;
 
let rabin_miller(n) =
  if  (n &% deux) &= zero then
    raise Nombre_pair
  else
    let d = dec(n) in
    let t = d.(1) in
    let s = d.(0) in
    let rec loop(k) =
      let a = nrandom (log2(n)) in
      if k < 64 then
        let v = a &&^ (s,n) in
        if carre(v,n,t) then
          loop(k+1)
        else
          false
      else
        true
    in loop(0);;
     
let big_rabin_miller(n,p) =
  if  (n &% deux) &= zero  then
    raise Nombre_pair
  else
    let d = dec(n) in
    let t = d.(1) in
    let s = d.(0) in
    let rec loop(k) =
      let a = nrandom (blog2(n,p)) in
      if k < 128 then
        let v = a &&^ (s,n) in
        if carre(v,n,t) then
          loop(k+2)
        else
          false
      else
        true
    in loop(0);;
 
 
(*
 let start = Sys.time();;
 let res = big_rabin_miller(bsub(power_int_positive_int 2 19937,un),19);;
 let stop = Sys.time();;
 Printf.printf "%6.2f  %B  \n" (stop -. start)  res;;
 
La réponse en 20 minutes
  
  moi@moi-bur:~/IUT/INFO1/CAML$ ocamlopt nums.cmxa  arith_info1.ml
moi@moi-bur:~/IUT/INFO1/CAML$ ./a.out
1200.31  true  
*)
 
 
 
(* test de primalité pour de grands entiers
   On teste d'abord si n est divisible par les
   petits nombres premiers < 1000
*)
 
let (&>=) a b = compare_big_int a b >= 0;;
 
let est_premier(n) =
  let c = if  n &>=  (big(1000000)) then
      crible(1000)
    else
      crible(int_of_big_int(sqrt_big_int n))
  in
  let rec test_loop = function
    | [] -> rabin_miller(n)
    | tete::queue ->
      if (n &% big(tete)) &= zero then
        false
      else
        test_loop queue
  in test_loop(c);;
    
 
 
 
(*
Génère un nombre fortement probablement premier
  Dans l'intervalle [2^(n-1);2^n[
*)
 
 
 (*renvoie un nombre entre 2^(n-1) et 2^n*)
let nrandom1 n =
  if n < 0 then raise Taille_negative 
  else if n = 0 then zero
  else begin
    let p = if n land 15 = 0 then 16 else n land 15 in
    let mask = 1 lsl (p-1) in
    let r = ref(big_int_of_int((Random.int(0x10000) land (mask-1)) lor mask)) in
    for i=1 to (n-p)/16 do
        r := add_int_big_int (Random.int(0x10000)) (mult_int_big_int 0x10000 !r)
    done;
    !r
  end
;;
 
 
(* génère un fpp en faisant 10xn essais*)
exception Rate;;
 
let genere_premier n =
  let rec loop(k) =
    if k = 0 then raise Rate
      else
      let p = nrandom1 n in
      if est_premier(p) then
        (*-- montre(p),10*n-k pour voir le nb d'essais --*) 
        p
      else
        loop(k-1)
  in loop(10*n);;
 
 
(*
#- : string =
"9323965933044776889460796166024415293535391530550864193145919935327042258144900912163674609469446072081838554854780728105800118911495921310476453003669076137150734310674927109180878880209834705050148025111338966756468709557010173342077387401166199830900082541888253591201452076074749456950555323761317"
*)
 
 
(* RSA *)
 
(*-Big Euclide II / calcul du PGCD + coeff de Bézout version récursive-*)
 
let bbezout(a,b) =
  let rec euc(u,v,r,u',v',r') =
    if r' &= zero then
      [|u;r|]
    else
      let q = r &/ r' in 
        euc(u',v',r',u &- (q &* u'), v &- (q &* v'),r &- (q &* r'))
  in  euc(un,zero,a,zero,un,b);;
 
(*---------Fonction qui recherche l'inverse de a modulo n---------*)
 
let binv_mod(a,n) =
  let bb = bbezout(a,n) in
  if bb.(1) &= un then
    bb.(0) &% n
  else
    raise Non_inversible;;
 
 
 
(*--
let p = genere_premier 1000;;
let q = genere_premier 1000;;
let n = p &* q ;;
let phi_n = (p &- un) &* (q &- un) ;;
let e = genere_premier 1000;;
montre(bgcd e phi_n);;
let d = binv_mod(e,phi_n);;
let m = big(646566);;
let crypte = m &&^ (e,n);;
let decrypte = crypte &&^ (d,n);;
montre(decrypte);;
--*)
 
 
(*-------------------------POLYNÔMES--------------------------------*)
 
(* Polynômes à coeff entiers sans créer de module *)
 
type 'a polynome = {coeff : 'a list ; deg : int };;
 
let zero = 0;;
 
let trouve_degre l =
  let t = List.rev l in
  let d = List.length l - 1  in
  let rec loop = function
    |[],k -> k
    |tete::queue,k ->
      if tete = zero then loop(queue,k-1)
      else k
  in loop(t,d);;
 
let liste = [0;1;1;0;0;1;0;0;0;1;0;0;0];;
 
trouve_degre liste;;
 
let pol_of_list l = {coeff = l ; deg = trouve_degre l};;
 
let nul = pol_of_list [zero];;
 
let a = pol_of_list [0;1;1;0;0;1];;
let b = pol_of_list [0;1;1];;
 
let xor a b = if a <> b then 1 else 0;;

courtesy of webmatter.de