Previous Up
Version pdf - Version archive

4.6  Polynômes

fichier source

[(-1.,3);(1., 1); (2.,0)];; (*******************************************************************) (* 4X^4 + 2X^3 -3X^2 + 3 *) (*******************************************************************) let polynomeUnite = [(1., 0)];; let polynomeNul = [(0., 0)];; (*******************************************************************) let strOfMonome (coeff, exp) = string_of_float coeff ^ "*X^" ^ string_of_int exp;; (*******************************************************************) let rec strOfPolynome = function | [] -> strOfPolynome polynomeNul | h::[] -> strOfMonome h | h::t -> strOfMonome h ^ " + " ^ strOfPolynome t;; (*******************************************************************) let rec expOfPolynome = function [] -> Const 0. | (c, e)::t -> Somme (Produit(Const(c), Puiss(X, e)), expOfPolynome t);; (*******************************************************************) let absFloat x = if x > 0. then x else (-.x);; (*******************************************************************) let estNul x = absFloat x < 1.e-15;; (*******************************************************************) let rec addPolynomes p q = match (p, q) with | [] , _ -> q | _, [] -> p | (coeffp, expp)::tailp, (coeffq, expq)::tailq -> if (expp = expq) then ( if (estNul (coeffp +. coeffq)) then (addPolynomes tailp tailq) else (coeffp +. coeffq, expp)::(addPolynomes tailp tailq) ) else if (expp > expq) then (coeffp, expp)::(addPolynomes tailp q) else (coeffq, expq)::(addPolynomes p tailq);; (*******************************************************************) let rec multConst p k = match p with [] -> [] | (c, e)::t -> (if (absFloat (k *. c) > 1.e-15) then [(k *. c, e)] else [] )@(multConst t k);; (*******************************************************************) let subPolynomes a b = addPolynomes a (multConst b (-1.));; (*******************************************************************) let rec multXk p k = match p with [] -> [] | (c, e)::t -> (c, e+k)::(multXk t k);; (*******************************************************************) let multMonome p (c, e) = multXk (multConst p c) e;; (*******************************************************************) let rec multPolynomes p q = match p with [] -> [] | h::t -> addPolynomes (multMonome q h) (multPolynomes t q);; (*******************************************************************) let rec expPolynome p = function 0 -> polynomeUnite | n -> multPolynomes p (expPolynome p (n-1));; (*******************************************************************) let rec polynomeOfExp = function Const(x) -> [(x,0)] | X -> [(1.,1)] | Somme(left, right) -> addPolynomes (polynomeOfExp left) (polynomeOfExp right) | Produit(left, right) -> multPolynomes (polynomeOfExp left) (polynomeOfExp right) | Puiss(left, n) -> expPolynome (polynomeOfExp left) n;; let k = polynomeOfExp (Puiss(Somme(Const(1.), X), 4));; strOfPolynome k;; (*******************************************************************) let rec floatExp x = function 0 -> 1. | n -> (if n mod 2 = 0 then 1. else x ) *. floatExp (x *. x) (n/2);; (*******************************************************************) let rec evalueMoche p x = match p with [] -> 0. | (c, e)::t -> c *. (floatExp x e) +. evalueMoche t x;; (*******************************************************************) (* un monome de degre k s'evalue en O(log k), qui est le temps necessaire pour mettre X a la puissance k. On determine le temps que met un polynome de degre n pour s'evaluer en calculant u_n = log 1 + log 2 + ... + log n qui est la somme des temps d'exponentiation de chaqe monome. On a u_n <= log n + ... + log n <= n log n evalueMoche s'execute en O(n log n). *) (*******************************************************************) (* ((x - 2)x + 7)x -1 (((4x^4 + 2)x - 2)x - 1)x + 1 *) (*******************************************************************) let rec evalHornerAcc p x acc = match p with [] -> acc | (c, e)::[] -> (acc +. c) *. (floatExp x e) | (c, e1)::tail -> (match tail with | (_, e2)::_ -> evalHornerAcc tail x ((acc +. c) *. (floatExp x (e1 - e2))) | _ -> failwith "Oups...");; let evalHorner p x = evalHornerAcc p x 0. ;; (*******************************************************************) (* Supposons que les degres decrive la suite n, n-1, n-2, ..., 2, 1, 0 Alors chaque chacune des (n + 1) iterations s'execute en temps constant. Dans ce cas l'algorithme s'execute en O(n). Dans le cas il existe au moins un entier k < n tel qu'il n'y a pas de monome de degre k, alors il suffit de considerer dans le calcul que le polynome contient un monome 0.X^k, ce qui ne change rien a la complexite qui a ete calculee. *) (*******************************************************************) let rec derivePolynome = function [] -> [] | (coeff, 0)::[] -> [] | (coeff, exp)::tail -> (coeff *. (float_of_int exp), exp - 1)::(derivePolynome tail);; (*******************************************************************) let rec primitivePolynome = function [] -> [] | (coeff, exp)::tail -> (coeff /. (float_of_int exp +. 1.), exp + 1)::(primitivePolynome tail);; (*******************************************************************) let rec integrePolynome p a b = let q = primitivePolynome p in evalHorner q b -. (evalHorner q a);; (*******************************************************************) let rectangles p a b n = let nFloat = float_of_int n in let rec aux p a b m = (b -. a) *. (evalHorner p a) /. nFloat +. ( if (m = 1) then 0. else aux p ((b -. a)/. nFloat +. a) b (m-1) ) in aux p a b n;; (*******************************************************************) let trapezes p a b n = let nFloat = float_of_int n in let dx = (b -. a) /. (float_of_int n) in let somme = ref (evalHorner p a +. evalHorner p b) in for i = 1 to n-1 do let j = float_of_int i in somme := 2.*.(evalHorner p (dx*.j))+. !somme; done; (b -. a) /. (2. *. nFloat) *. !somme;; (*******************************************************************) let facteurLagrange x = [(1., 1); (-.x, 0)];; (*******************************************************************) let rec lagrange_f_i pts x_i = match pts with [] -> polynomeUnite | (x, y)::t -> if x = x_i then lagrange_f_i t x_i else multPolynomes (facteurLagrange x) (lagrange_f_i t x_i);; (*******************************************************************) (* p_i = y_i * f_i(x)/f_i(x_i) *) (*******************************************************************) let lagrange_p_i pts (x_i, y_i) = let f_i = lagrange_f_i pts x_i in multConst f_i (y_i /. (evalHorner f_i x_i));; (*******************************************************************) (* p(x) = somme des p_i(x) *) (*******************************************************************) let interpolationLagrange pts = let rec add_p_i l = (match l with [] -> polynomeNul | h::t -> addPolynomes (lagrange_p_i pts h) (add_p_i t) ) in add_p_i pts;; (*******************************************************************) let interpolationVandermonde pts = let n = List.length pts in let x = make_vect n 0. in let y = make_matrix n 1 in let rec initTab l m = match l with | [] -> () | (xi, yi)::tail -> ( x.(m) <- xi; y.(m).(0) <- yi; initTab tail (m+1) ) in initTab pts 0; let a = make_matrix n n in fill_matrix a (function (i, j) -> if (j = 0) then 1. else a.(i).(j-1) *. x.(i)); let x = solve_system a y in let rec coeffs k = if (k >= 0) then (x.(k).(0), k)::(coeffs (k-1)) else [] in coeffs (n-1);; (*******************************************************************) (* u_(n+1) = u_n - f(x)/f'(x) *) (*******************************************************************) let rec newton p x = function 0 -> x | n -> newton p (x -. (evalHorner p x)/.(evalHorner (derivePolynome p) x)) (n - 1);; (*******************************************************************) let degre = function | [] -> -1 | (c, e)::_ -> if (e <> 0) then e else if (estNul c) then -1 else 0;; (*******************************************************************) let rec divisionPolynomes a b = if (degre a < degre b) then (polynomeNul, a) else ( match a, b with (ca, ea)::_, (cb, eb)::_ -> let q = [(ca/.cb, ea-eb)] in let r = subPolynomes a (multPolynomes b q) in let (qrec, rrec) = divisionPolynomes r b in ((addPolynomes q qrec), rrec) );; (*******************************************************************) let rec pgcd a b = if (degre b <= 0) then a else let (_, r) = divisionPolynomes a b in pgcd b r;; let r = multPolynomes (facteurLagrange (-.1.)) (facteurLagrange 4.);; strOfPolynome r;; let p = multPolynomes (facteurLagrange 5.) (multPolynomes r (facteurLagrange 2.));; let q = multPolynomes (facteurLagrange 6.) (multPolynomes r (facteurLagrange 1.));; strOfPolynome (pgcd p q);; (*******************************************************************) let rec euclideE a b lastV = if (degre b <= 0) then (a, polynomeNul, lastV) else let (q, r) = divisionPolynomes a b in let (p, u, v)= euclideE b r lastV in (d, v, subPolynomes u (multPolynomes v q));; strOfPolynome (euclideE p q polynomeNul);; (*******************************************************************) let sturmNext n npUn = let q, r = divisionPolynomes n npUn in multConst r (-1.);; let sturmSuite p = let rec sturmSuiteAux u0 u1 = let u2 = sturmNext u0 u1 in u2::( if (degre u2 <= 0) then [] else sturmSuiteAux u1 u2 ) in let dp = multConst (derivePolynome p) (-1.) in p::dp::(sturmSuiteAux p dp);; (*******************************************************************) let rec nbSignChange s x = match s with [] -> 0 | h::[] -> 0 | h1::h2::t -> ( let h1x = evalHorner h1 x in let h2x = evalHorner h2 x in if (h1x < 0. && h2x > 0.) or (h1x > 0. && h2x < 0.) then 1 else 0 ) + (nbSignChange (h2::t) x);; (*******************************************************************) let rec map f = function [] -> [] | h::t -> (f h)::(map f t);; (*******************************************************************) let rec trouveRacines s inf sup eps = let sinf = nbSignChange s inf in let ssup = nbSignChange s sup in if (sinf = ssup) then [] else let mid = (inf +. sup) /. 2. in if (absFloat (sup -. inf) < eps) then [mid] else (trouveRacines s inf mid eps)@(trouveRacines s mid sup eps);; (*******************************************************************) let majoreRacines p = let maxFloat a b = if (a < b) then b else a in match p with (a_n, _)::_ -> let rec maxCoeff l = ( match l with | (h, _)::[] -> absFloat h | (h, _)::t -> maxFloat (absFloat h) (maxCoeff t) | _ -> failwith "oups" ) in (maxCoeff p)/.(absFloat a_n) +. 1. | _ -> failwith "oups";; (*******************************************************************) let sturm p eps = let s = sturmSuite p in let a = majoreRacines p in trouveRacines s (-.a) a eps;; (*******************************************************************) let rec monstronome n = match n with 0 -> polynomeUnite | n-> multPolynomes (expPolynome (facteurLagrange (float_of_int n)) (n mod 4 + 1)) (monstronome (n-1));; let rec gronome = function 0 -> [(-.1., 0)] | n-> (float_of_int ((if n mod 2 = 0 then (-n) else n) / 2) , n)::(gronome (n-1));; for i=1 to 8 do let p = monstronome i in print_float (integrePolynome p 0. 5.); print_string "\n"; print_float (rectangles p 0. 5. 100000); print_string "\n"; print_float (trapezes p 0. 5. 100000); print_string "\n"; print_string (strOfPolynome p); print_string "\n"; print_float (majoreRacines p); print_string "\n"; map (function e -> print_float e; print_string ";") (sturm p 1.e-10); print_string "\n"; done;;

Previous Up