立方根の計算 in Haskell DouKaku?

xは0以上1000未満の実数です。 y * y * y = xになるような実数y(立方根)を小数点以下12桁以上の正確さで求める関数cube_rootを作って下さい。

ただし、このお題の趣旨は実数区間での探索なので、立方根関数があっても使ってはいけません。指数関数と対数関数も禁止します。

Pythonで表現した入出力の例：

>>> cube_root(10.0)
2.1544346900318834
>>> _ ** 3
9.9999999999999947
>>> cube_root(100.0)
4.6415888336127793
>>> _ ** 3
100.00000000000003

そのままHaskellで

cbrt :: Double -> Double -> Double
cbrt e x = iter 1.0 x x
  where iter old new x
           | good old new = new
           | otherwise    = iter new (improve new x) x
        improve g x  = (x / g^2 + 2 * g) / 3
        good old new = abs (1 - old / new) < e

cubicRoot :: Double -> Double
cubicRoot = cbrt 1.0e-13

{-
*Main> cubicRoot 1.0
1.0
*Main> it ^ 3
1.0
*Main> cubicRoot 10.0
2.154434690031884
*Main> it ^ 3 
10.000000000000002
*Main> cubicRoot 100  
4.641588833612778
*Main> it ^ 3
99.99999999999997
*Main> cubicRoot 1000
10.0
*Main> it ^ 3
1000.0
*Main> cubicRoot 1.0e-30
1.0e-10
*Main> it ^ 3
1.0e-30
-}

0.0だけは別判定しないとだめですね。

starling f g x = f x (g x)
cubicRoot 0 = 0
cubicRoot x = snd . head . filter ((e >) . abs . subtract 1 . uncurry (/)) 
            . starling zip tail . flip iterate 1 . improve $ x
  where improve x g = (x / g^2 + 2 * g) / 3
         e = 1.0e-12

自動微分ライブラリを使ってみました。 HackageDBにあります。

このライブラリのderiv関数を使うと導関数を「数値的に」ではなく「記号的に」求める事ができます。つまりf'(x) = {f(x + h) - f(x)} / hという近似を利用するのではなく、 (x^3)' = 3x^2 の様に直接導関数を生成します。

円周率と自然対数も求めてみました。

import Data.Number.Dif

newton f init eps =
  let f' = deriv f
      next = init - (unDif f init) / (f' init)
  in if abs (next - init) < eps
    then init
    else newton f next eps


-- aの立方根は x^3 - a = 0の解
cbrt a = newton (\x -> x**3 - a) 1 1e-13

main = do
  putStrLn $ "cbrt 10   = " ++ show (cbrt 10)
  putStrLn $ "cbrt 100  = " ++ show (cbrt 100)
  putStrLn $ "cbrt 1000 = " ++ show (cbrt 1000)
  putStrLn $ "pi        = " ++ show (newton sin 3 1e-13)
  putStrLn $ "e         = " ++ show (newton (\x -> log x - 1) 2 1e-13)

{-
*Main> :main
cbrt 10   = 2.154434690031884
cbrt 100  = 4.641588833612779
cbrt 1000 = 10.000000000000004
pi        = 3.141592653589793
e         = 2.7182818284590455
-}

DouKaku?

立方根の計算

Posted feedbacks - Haskell

Index

Feed

Other

Link

Pathtraq