newton.hs

-- 与えられた座標列のn階差分商を求める関数
diffmethod :: [(Double, Double)] -> Double
diffmethod [x0] = snd x0
diffmethod [x0,x1] = (snd x1 - snd x0) / (fst x1 - fst x0)
diffmethod lst = (diffmethod (tail lst) - diffmethod (init lst)) 
               / (fst (last lst) - fst (head lst))

-- 第一引数に与えられた座標列から第二引数のy座標を補完し返す関数
newtonPoly :: [(Double, Double)] -> Double -> Double
newtonPoly [] _  = 0
newtonPoly lst x = newtonPoly (init lst) x 
                 + product (map (\a -> x - fst a) (init lst)) 
                 * diffmethod lst

-- 元の関数
originf x = sin x

-- 0.5刻みで0から2πまでの座標列を生成
sample = map (\x -> (x, originf x)) [0.0,0.5..2*pi]

main = do
  -- 0.1刻みで0から2πまでのそれぞれの点においての補間関数と実際の関数の値を表示
  let p = map (\x -> concat [show x, ", ", show $ newtonPoly sample x,", ", show $ originf x ]) [0.0,0.1..2*pi]
  mapM_ putStrLn p