module UseDif where

--import NR
import Dif

-- dy/dx = y
-- y(0) = 1
el :: Dif Double
el = mkDif 1 el

-- dy/dx = y^2 + 1
-- y(0) = 0
yl :: Dif Double
yl = mkDif 0 (yl^2 + 1)

-- The Lambert function W(z)
-- W(z)*exp(W(z)) = z
-- z = W(z)*exp(W(z))
--
-- dz/dW = 1*exp(W) + W*exp(W) = 
--       = exp(W) * (1 + W)
-- dW/dz = exp(-W) / (1 + W)
-- W(0) = 0
-- So the derivatives at 0 are
wl :: (Floating a) => Dif a
wl = mkDif 0 (exp(-wl) / (1 + wl))

relative :: (Num a, Ord a) => a -> [a] -> a
relative eps (a:b:rest) | abs(a-b) < eps * abs b = b
                        | otherwise              = relative eps (b:rest)
relative _ _ = error "relative"

fac :: (Num a) => Integer -> a
fac 0 = 1
fac n = fromInteger n * fac(n-1)

sums :: (Num a) => [a] -> [a]
sums = scanl (+) 0

getAllDerivatives :: (Num a) => Dif a -> [a]
getAllDerivatives = map val . iterate df

taylors :: (Fractional a) => a -> Dif a -> a -> [a]
taylors a d x = sums $ zipWith (\ n fn -> fn / fac n * (x-a)^n) [0..] (getAllDerivatives d)

--maclaurins :: (Fractional a) => Dif a -> a -> [a]
--maclaurins = taylors 0

lambert :: (Floating a, Ord a) => a -> a -> a
lambert eps = relative eps . taylors 0 wl

--------

-- newtonRaphson eps f z gives a solution to
-- f x = 0, using Newton-Raphson starting at z,
-- and using relative approx eps
newtonRaphson :: (Fractional a, Ord a) => a -> (Dif a -> Dif a) -> a -> a
newtonRaphson eps f z = relative eps $ iterate (val . next . dVar) z
  where next x = x - f x / df (f x)