(* $Id: polynomials.ml,v 1.27 2008/01/30 20:05:00 averell Exp $ *) open Format module PolyOnField (R : Signatures.RING) = struct open Array type coeffRing = R.t type t = coeffRing array type poly_t = t let unk = "X" let zero = [| |] let one = [| R.one |] let x = [| R.zero; R.one|] let monome coeff deg = if R.is_zero coeff then zero else let res = make (deg+1) R.zero; in res.(deg) <- coeff; res let normal polx = let i = ref (length polx - 1) in while (!i >=0) && (R.( =/ ) polx.(!i) R.zero) do i := !i - 1 done; if !i = (length polx - 1) then polx else sub polx 0 (!i + 1) let degree polx = (* if (length polx) <> length (normal polx) then ( printf "--%i-%i--!." (length polx) (length (normal polx)); failwith "Non normalized polynomials") else *) (length polx) - 1 let degree_intern polx = (length polx) - 1 let coeff polx n = if n > degree polx then R.zero else polx.(n) let lcoeff polx = let dpolx = degree polx in if dpolx = -1 then R.zero else polx.(dpolx) let multXn polx n = if degree polx = -1 then zero else append (make n R.zero) polx let multX polx = multXn polx 1 let of_R r = normal [|r|] let of_int i = normal (of_R (R.of_int i)) let of_Rlist rl = normal (of_list rl) let to_Rlist pol = to_list pol let multCoeff polx r = (map (R.( */ ) r) polx) let rec ( +/ ) polx poly = if degree_intern polx < degree_intern poly then poly +/ polx else let res = copy polx in for i=0 to degree_intern poly do res.(i) <- R.(+/) res.(i) poly.(i) done; normal res let negate = map R.negate let ( -/ ) polx poly = polx +/ (negate poly) let is_zero pol = (length pol) = 0 let ( =/ ) a b = is_zero (a -/ b) (* ( = ) : t -> t -> bool *) (* not quite correct but fast*) (* multiplication using naive algorithm *) let naiveMult polx poly = if degree_intern polx = -1 || degree_intern poly = -1 then zero else let res = make ((degree_intern polx) + (degree_intern poly) + 1) R.zero in for ix=0 to (degree_intern polx) do for iy=0 to (degree_intern poly) do res.(ix+iy) <- R.( +/ ) res.(ix+iy) ( R.( */ ) polx.(ix) poly.(iy)) done; done; normal res let cut ar lCut = let lar = length ar in if lar <= lCut then (ar, zero) else ((sub ar 0 lCut), (sub ar lCut (lar - lCut))) (* multiplication using Karatsuba algorithm *) let rec ( */ ) polx poly = let karaThreshold = 16 in (* Naive / Karatsuba methode *) let dx = degree_intern polx and dy = degree_intern poly in if dx <= karaThreshold || dy <= karaThreshold then naiveMult polx poly else let chooseDeg = (max dx dy) in let lCut = (chooseDeg+1) / 2 in let (x0, x1) = cut polx lCut and (y0, y1) = cut poly lCut in let x0y0 = x0 */ y0 and x1y1 = x1 */ y1 and xxyy = (x0 +/ x1) */ (y0 +/ y1) in ( (x0y0 +/ (multXn x1y1 (2 * lCut))) +/ (multXn (xxyy -/ (x0y0 +/ x1y1)) lCut) ) let repSquaringInv one mult inv x a = let rec powacc accu x a = if a = 0 then accu else if a > 0 then if a mod 2 = 0 then powacc accu (mult x x) (a/2) else powacc (mult accu x) (mult x x) (a/2) else powacc one (inv x) (-a) in powacc one x a let repSquaring one mult = repSquaringInv one mult (function x -> failwith "Non invertible") let pow = repSquaring one ( */ ) let diff polx = let res = make (degree polx) R.zero in let rec diff_aux i iR = if i >= (length res) then res else (res.(i) <- R.( */ ) iR polx.(i+1); diff_aux (i+1) (R.( +/ ) iR R.one)) in normal (diff_aux 0 R.one) let print polx = if degree polx = -1 then R.print R.zero else let first = ref true in printf "@["; if not (R.(=/) polx.(0) R.zero) then begin first := false; printf "(@["; R.print polx.(0); printf "@])" end; for ix=1 to degree polx do if not (R.(=/) polx.(ix) R.zero) then begin if !first then (first := false; printf "@[";) else printf "@ + @["; if not (R.(=/) polx.(ix) R.one) then begin printf "(@["; R.print polx.(ix); printf "@])"; end; if ix > 1 then Format.printf "%s^%i@]" unk ix else Format.printf "%s@]" unk end; done; printf "@]" let eval pol x = let rec eval_rec accu i = if i = -1 then accu else eval_rec (R.( +/ ) pol.(i) (R.( */ ) x accu)) (i-1) in eval_rec R.zero (degree pol) let subs pol y = let rec eval_rec accu i = if i = -1 then accu else eval_rec ((of_R pol.(i)) +/ (x */ accu)) (i-1) in normal (eval_rec zero (degree pol)) end