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;;