Dichotomie

-- Signature la plus générale
--dicho  :: (Fractional t, Num u, Ord t, Ord u) =>
--          (t -> u) -> t -> t -> t -> t
 
-- Signature dans le cas d'une fonction sur 64 bits mais alors mpfr not possible
-- dicho :: (Double -> Double) -> Double -> Double -> Double -> Double
dicho      f          a    b    p =
    -- si l'intervalle est de largeur <= à la précision
    if (abs (b - a) <= p) then
        -- je renvoie m (le milieu de [a,b])
        m
    else
        -- si a et m sont de part et d'autre de la solution
        if (f a) * (f m) < 0 then
            -- je recherche la solution entre a et m
            dicho f a m p
        else
            -- sinon je la recherche entre m et b
            dicho f m b p
    -- m étant le milieu de [a,b]
    where m = 0.5 * (a + b)
 
 
-- Problème ! si p est de l'ordre de l'eps machine, le test échoue
-- (Rule of thumb : éviter les while (ou équivalent))
-- On essaie de prévoir le nombre nécessaire d'itérations pour passer à un for
-- On cherche le nb d'étapes à partir duquel (b0 - a0)/2^n devient <= p
-- On trouve n >= ????
 
          
dicho2 f a b p =
    let n = ceiling ((log ((b - a) / p)) / (log 2)) in
    let aux aa bb k =
            let m = 0.5 * (aa + bb) in
            if k == n then
                m
            else
                if (f aa) * (f m) < 0 then
                    aux aa m (k + 1)
                else
                    aux m bb (k + 1)
    in aux a b 0