Skip to content

Instantly share code, notes, and snippets.

@frostburn
Last active April 9, 2018 05:02
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save frostburn/5123472 to your computer and use it in GitHub Desktop.
Save frostburn/5123472 to your computer and use it in GitHub Desktop.
Cube root in Haskell
--{-# LANGUAGE TemplateHaskell #-}
--import Test.QuickCheck
--import Test.QuickCheck.All
--prop_cbrt :: Double -> Bool
--prop_cbrt x = abs (x - cbrt x ^ 3) < 1e-12
--main = ($quickCheckAll)
module Cbrt (cbrt) where
--Newton's iteration.
newton :: Fractional a => (a -> a) -> (a -> a) -> a -> [a]
newton f f' x0 = iterate (\x -> x - f x / f' x) x0
--Assuming that a sequence of is converging to a fixed value,
--take the first value that fails to improve in precision.
converge :: [Double] -> Double
converge xs = snd . head $ dropWhile ((< 0) . fst) $ zip deltaDeltas xs where
deltaDeltas = zipWith (-) (tail deltas) deltas
deltas = map abs $ zipWith (-) (tail xs) xs
--Alternative implementation
converge' :: [Double] -> Double
converge' (x0:x1:x2:xs) = if (change >= lastChange) then x2 else converge' (x1:x2:xs) where
change = abs (x2 - x1)
lastChange = abs (x1 - x0)
--Given y find x such that x^3 == y using Newton's method.
cbrt' :: Double -> Double
cbrt' y = converge $ newton (\x -> x*x*x - y) (\x -> 3 * x*x) y
--Find the cube root of x.
cbrt :: Double -> Double
cbrt 0 = 0
cbrt x = if (abs x > 0.01)
then cbrt' x
else (1/) $ cbrt' (1 / x)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment