32 lines
1 KiB
Haskell
32 lines
1 KiB
Haskell
thiele :: [Double] -> [Double] -> Double -> Double
|
|
thiele xs ys = f rho1 (tail xs)
|
|
where
|
|
f _ [] _ = 1
|
|
f r@(r0:r1:r2:rs) (x:xs) v = r2 - r0 + (v - x) / f (tail r) xs v
|
|
rho1 = (!! 1) . (++ [0]) <$> rho
|
|
rho = repeat 0 : repeat 0 : ys : rnext (tail rho) xs (tail xs)
|
|
where
|
|
rnext _ _ [] = []
|
|
rnext r@(r0:r1:rs) x xn =
|
|
let z_ = zipWith
|
|
in z_ (+) (tail r0) (z_ (/) (z_ (-) x xn) (z_ (-) r1 (tail r1))) :
|
|
rnext (tail r) x (tail xn)
|
|
|
|
-- Inverted interpolation function of f
|
|
invInterp :: (Double -> Double) -> [Double] -> Double -> Double
|
|
invInterp f xs = thiele (map f xs) xs
|
|
|
|
main :: IO ()
|
|
main =
|
|
mapM_
|
|
print
|
|
[ 3.21 * inv_sin (sin (pi / 3.21))
|
|
, pi / 1.2345 * inv_cos (cos 1.2345)
|
|
, 7 * inv_tan (tan (pi / 7))
|
|
]
|
|
where
|
|
[inv_sin, inv_cos, inv_tan] =
|
|
uncurry ((. div_pi) . invInterp) <$>
|
|
[(sin, (2, 31)), (cos, (2, 100)), (tan, (4, 1000))]
|
|
-- N points taken uniformly from 0 to Pi/d
|
|
div_pi (d, n) = (* (pi / (d * n))) <$> [0 .. n]
|