Add special data type for passing parameters for root-finding functions

Currently only `ridders' make use of that data type. Also add data type
for specifying error tolerance

Numeric/RootFinding.hs

@@ -10,10 +10,15 @@
 --
 -- Haskell functions for finding the roots of real functions of real arguments.
 module Numeric.RootFinding
-    (
+    ( -- * Data types
       Root(..)
     , fromRoot
+    , Tolerance(..)
+    , withinTolerance
+    -- * Ridders algorithm
+    , RiddersParam(..)
     , ridders
+    -- * Newton-Raphson algorithm
     , newtonRaphson
     -- * References
     -- $references
@@ -22,12 +27,18 @@ module Numeric.RootFinding
 import Control.Applicative              (Alternative(..), Applicative(..))
 import Control.Monad                    (MonadPlus(..), ap)
 import Data.Data                        (Data, Typeable)
+import Data.Default.Class
 #if __GLASGOW_HASKELL__ > 704
 import GHC.Generics                     (Generic)
 #endif
-import Numeric.MathFunctions.Comparison (within)
+import Numeric.MathFunctions.Comparison (within,eqRelErr)
+import Numeric.MathFunctions.Constants  (m_epsilon)
 
 
+----------------------------------------------------------------
+-- Data types
+----------------------------------------------------------------
+
 -- | The result of searching for a root of a mathematical function.
 data Root a = NotBracketed
             -- ^ The function does not have opposite signs when
@@ -81,6 +92,54 @@ fromRoot _ (Root a) = a
 fromRoot a _        = a
 
 
+-- | Error tolerance for finding root. It describes when root finding
+--   algorithm should stop trying to improve approximation.
+data Tolerance
+  = RelTol !Double
+    -- ^ Relative error tolerance. Given @RelTol ε@ two values are
+    --   considered approximately equal if
+    --   \[ |a - b| / |\operatorname{max}(a,b)} < \vareps \]
+  | AbsTol !Double
+    -- ^ Absolute error tolerance. Given @AbsTol δ@ two values are
+    --   considered approximately equal if \[ |a - b| < \delta \].
+    --   Note that @AbsTol 0@ could be used to require to find
+    --   approximation within machine precision.
+  deriving (Eq, Read, Show, Typeable, Data
+#if __GLASGOW_HASKELL__ > 704
+           , Generic
+#endif
+           )
+
+withinTolerance :: Tolerance -> Double -> Double -> Bool
+withinTolerance (RelTol eps) a b = eqRelErr eps a b
+-- NOTE: `<=` is needed to allow 0 absolute tolerance which is used to
+--       describe precision of 1ulp
+withinTolerance (AbsTol tol) a b = abs (a - b) <= tol
+
+----------------------------------------------------------------
+-- Ridders algorithm
+----------------------------------------------------------------
+
+-- | Parameters for 'ridders' root finding
+data RiddersParam = RiddersParam
+  { riddersMaxIter :: !Int
+    -- ^ Maximum number of iterations.
+  , riddersTol     :: !Tolerance
+    -- ^ Error tolerance for root approximation.
+  }
+  deriving (Eq, Read, Show, Typeable, Data
+#if __GLASGOW_HASKELL__ > 704
+           , Generic
+#endif
+           )
+
+instance Default RiddersParam where
+  def = RiddersParam
+        { riddersMaxIter = 100
+        , riddersTol     = RelTol (4 * m_epsilon)
+        }
+
+
 -- | Use the method of Ridders[Ridders1979] to compute a root of a
 --   function. It doesn't require derivative and provide quadratic
 --   convergence (number of significant digits grows quadratically
@@ -90,7 +149,7 @@ fromRoot a _        = a
 --   and upper bounds of the search (i.e. the root must be
 --   bracketed). If there's more that one root in the bracket
 --   iteration will converge to some root in the bracket.
-ridders :: Double
+ridders :: RiddersParam
         -- ^ Absolute error tolerance. Iterations will be stopped when
         --   difference between root and estimate is less than
         --   tolerance or when precision couldn't be improved further
@@ -100,7 +159,7 @@ ridders :: Double
         -> (Double -> Double)
         -- ^ Function to find the roots of.
         -> Root Double
-ridders tol (lo,hi) f
+ridders p (lo,hi) f
     | flo == 0    = Root lo
     | fhi == 0    = Root hi
     -- root is not bracketed
@@ -114,14 +173,14 @@ ridders tol (lo,hi) f
     --
     go !a !fa !b !fb !i
         -- Root is bracketed within 1 ulp. No improvement could be made
-        | within 1 a b       = Root a
+        | within 1 a b = Root a
         -- Root is found. Check that f(m) == 0 is nessesary to ensure
         -- that root is never passed to 'go'
-        | fm == 0            = Root m
-        | fn == 0            = Root n
-        | d < tol            = Root n
+        | fm == 0      = Root m
+        | fn == 0      = Root n
+        | withinTolerance (riddersTol p) a b = Root n
         -- Too many iterations performed. Fail
-        | i >= (100 :: Int)  = SearchFailed
+        | i >= riddersMaxIter p              = SearchFailed
         -- Ridder's approximation coincide with one of old bounds or
         -- went out of (a,b) range due to numerical problems. Revert
         -- to bisection
@@ -133,7 +192,6 @@ ridders tol (lo,hi) f
         | fn*fa < 0          = go a fa n fn (i+1)
         | otherwise          = go n fn b fb (i+1)
       where
-        d    = abs (b - a)
         dm   = (b - a) * 0.5
         -- Mean point
         !m   = a + dm
@@ -143,6 +201,11 @@ ridders tol (lo,hi) f
         !fn  = f n
 
 
+
+----------------------------------------------------------------
+-- Newton-Raphson algorithm
+----------------------------------------------------------------
+
 -- | Solve equation using Newton-Raphson iterations.
 --
 -- This method require both initial guess and bounds for root. If

math-functions.cabal

@@ -39,11 +39,12 @@ library
     DeriveGeneric
 
   ghc-options:          -Wall -O2
-  build-depends:        base >=4.5 && <5
+  build-depends:        base                >= 4.5 && < 5
                       , deepseq
-                      , vector >= 0.7
+                      , data-default-class  >= 0.1.2.0
+                      , vector              >= 0.7
                       , primitive
-                      , vector-th-unbox
+                      , vector-th-unbox     >= 0.2.1.6
   if flag(system-expm1) || !os(windows)
     cpp-options: -DUSE_SYSTEM_EXPM1
   exposed-modules:
@@ -83,19 +84,19 @@ test-suite tests
     Tests.SpecFunctions
     Tests.SpecFunctions.Tables
     Tests.Sum
-  build-depends:
-    math-functions,
-    base >=4.5 && <5,
-    deepseq,
-    primitive,
-    vector >= 0.7,
-    vector-th-unbox,
-    erf,
-    HUnit      >= 1.2,
-    QuickCheck >= 2.7,
-    test-framework,
-    test-framework-hunit,
-    test-framework-quickcheck2
+  build-depends:    base >=4.5 && <5
+                  , math-functions
+                  , data-default-class  >= 0.1.2.0
+                  , deepseq
+                  , primitive
+                  , vector >= 0.7
+                  , vector-th-unbox
+                  , erf
+                  , HUnit      >= 1.2
+                  , QuickCheck >= 2.7
+                  , test-framework
+                  , test-framework-hunit
+                  , test-framework-quickcheck2
 
 source-repository head
   type:     git

tests/Tests/RootFinding.hs

@@ -1,6 +1,7 @@
 -- |
 module Tests.RootFinding ( tests ) where
 
+import Data.Default.Class
 import Test.Framework
 import Test.Framework.Providers.HUnit
 
@@ -11,9 +12,11 @@ import Tests.Helpers
 tests :: Test
 tests = testGroup "Root finding"
   [ testGroup "Ridders"
-    [ testAssertion "sin x - 0.525 [exact]"  $ testRiddersSin0_525 0
-    , testAssertion "sin x - 0.525 [1e-12]"  $ testRiddersSin0_525 1e-12
-    , testAssertion "sin x - 0.525 [1e-6]"   $ testRiddersSin0_525 1e-6
+    [ testAssertion "sin x - 0.525 [exact]"     $ testRiddersSin0_525 (AbsTol 0)
+    , testAssertion "sin x - 0.525 [abs 1e-12]" $ testRiddersSin0_525 (AbsTol 1e-12)
+    , testAssertion "sin x - 0.525 [abs 1e-6]"  $ testRiddersSin0_525 (AbsTol 1e-6)
+    , testAssertion "sin x - 0.525 [rel 1e-12]" $ testRiddersSin0_525 (RelTol 1e-12)
+    , testAssertion "sin x - 0.525 [rel 1e-6]"  $ testRiddersSin0_525 (RelTol 1e-6)
     ]
   , testGroup "Newton-Raphson"
     [ testAssertion "sin x - 0.525 [1e-12]"  $ testNewtonSin0_525 1e-12
@@ -23,10 +26,12 @@ tests = testGroup "Root finding"
   where
     -- Exact root for equation: sin x - 0.525 = 0
     exactRoot = 0.5527151130967832
+    --
     testRiddersSin0_525 tol
-      = abs (r - exactRoot) <= tol
+      = withinTolerance tol r exactRoot
       where
-        Root r = ridders tol (0, pi/2) (\x -> sin x - 0.525)
+        Root r = ridders def{riddersTol = tol} (0, pi/2) (\x -> sin x - 0.525)
+    --
     testNewtonSin0_525 tol
       = abs (r - exactRoot) <= tol * r
       where