diff --git a/exact-real.cabal b/exact-real.cabal --- a/exact-real.cabal +++ b/exact-real.cabal @@ -1,5 +1,5 @@ name: exact-real -version: 0.10.0 +version: 0.11.0 synopsis: Exact real arithmetic description: A type to represent exact real number using a fast binary Cauchy sequence @@ -17,6 +17,7 @@ readme.md stack.yaml cabal-version: >=1.10 +tested-with: GHC >= 7.10 && < 8 source-repository head type: git @@ -25,6 +26,7 @@ library exposed-modules: Data.CReal + Data.CReal.Converge Data.CReal.Internal build-depends: base >= 4.8 && < 4.9, diff --git a/readme.md b/readme.md --- a/readme.md +++ b/readme.md @@ -2,6 +2,10 @@ ========== [![Hackage](https://budueba.com/hackage/exact-real)](https://hackage.haskell.org/package/exact-real) +[![Stackage Nightly](http://stackage.org/package/exact-real/badge/nightly)](http://stackage.org/nightly/package/exact-real) +[![Build Status](https://travis-ci.org/expipiplus1/exact-real.svg)](https://travis-ci.org/expipiplus1/exact-real) +[![Build status](https://ci.appveyor.com/api/projects/status/sviohnqtyo9y7goa?svg=true)](https://ci.appveyor.com/project/expipiplus1/exact-real) +[![Coverage Status](https://coveralls.io/repos/expipiplus1/exact-real/badge.svg?branch=codecov&service=github)](https://coveralls.io/github/expipiplus1/exact-real?branch=codecov) Exact real arithmetic implemented by fast binary Cauchy sequences. diff --git a/src/Data/CReal/Converge.hs b/src/Data/CReal/Converge.hs new file mode 100644 --- /dev/null +++ b/src/Data/CReal/Converge.hs @@ -0,0 +1,139 @@ +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeFamilies #-} + +-- | The Converge type class. +module Data.CReal.Converge + ( Converge(..) + ) where + +import Control.Arrow ((&&&)) +import Data.Coerce (coerce) +import Data.CReal.Internal (CReal(..), atPrecision, crMemoize) +import Data.Function (on) +import Data.Proxy (Proxy) +import GHC.TypeLits (someNatVal, SomeNat(..)) + +-- $setup +-- >>> :set -XFlexibleContexts +-- >>> import Data.CReal.Internal + +-- | If a type is an instance of Converge then it represents a stream of values +-- which are increasingly accurate approximations of a desired value +class Converge a where + -- | The type of the value the stream converges to. + type Element a + + -- | 'converge' is a function that returns the value the stream is converging + -- to. If given a stream which doens't converge to a single value then + -- 'converge' will not terminate. + -- + -- If the stream is empty then it should return nothing. + -- + -- >>> let initialGuess = 1 :: Double + -- >>> let improve x = (x + 121 / x) / 2 + -- >>> converge (iterate improve initialGuess) + -- Just 11.0 + -- + -- >>> converge [] :: Maybe [Int] + -- Nothing + converge :: a -> Maybe (Element a) + + -- 'convergeErr' is a function that returns the value the stream is + -- converging to. It also takes a function err which returns a value which + -- varies monotonically with the error of the value in the stream. This can + -- be used to ensure that when 'convergeErr' terminates when given a + -- non-converging stream or a stream which enters a cycle close to the + -- solution. See the documentation for the CReal instance for a caveat with + -- that implementation. + -- + -- It's often the case that streams generated with approximation + -- functions such as Newton's method will generate worse approximations for + -- some number of steps until they find the "zone of convergence". For these + -- cases it's necessary to drop some values of the stream before handing it + -- to convergeErr. + -- + -- For example trying to find the root of the following funciton @f@ with a + -- poor choice of starting point. Although this doesn't find the root, it + -- doesn't fail to terminate. + -- >>> let f x = x ^ 3 - 2 * x + 2 + -- >>> let f' x = 3 * x ^ 2 - 2 + -- >>> let initialGuess = 0.1 + -- >>> let improve x = x - f x / f' x + -- >>> let err x = abs (f x) + -- >>> convergeErr err (iterate improve initialGuess) + -- Just 0.1 + convergeErr :: Ord (Element a) => (Element a -> Element a) -> a -> Maybe (Element a) + +-- | Every list of equatable values is an instance of 'Converge'. 'converge' +-- returns the first element which is equal to the succeeding element in the +-- list. If the list ends before the sequence converges the last value is +-- returned. +instance {-# OVERLAPPABLE #-} Eq a => Converge [a] where + type Element [a] = a + + converge = lastMay . takeWhilePairwise (/=) + {-# INLINE converge #-} + + convergeErr err xs = fmap snd . lastMay . takeWhilePairwise ((>) `on` fst) $ es + where es = (err &&& id) <$> xs + {-# INLINE convergeErr #-} + +-- | The overlapping instance for @'CReal' n@ has a slightly different +-- behavior. The instance for 'Eq' will cause 'converge' to return a value when +-- the list converges to within 2^-n (due to the 'Eq' instance for @'CReal' n@) +-- despite the precision the value is requested at by the surrounding +-- computation. This instance will return a value approximated to the correct +-- precision. +-- +-- It's important to note when the error function reaches zero this function +-- behaves like 'converge' as it's not possible to determine the precision at +-- which the error function should be evaluated at. +-- +-- Find where log x = π using Newton's method +-- >>> let initialGuess = 1 +-- >>> let improve x = x - x * (log x - pi) +-- >>> let Just y = converge (iterate improve initialGuess) +-- >>> showAtPrecision 10 y +-- "23.1406" +-- >>> showAtPrecision 50 y +-- "23.1406926327792686" +instance {-# OVERLAPPING #-} Converge [CReal n] where + type Element [CReal n] = CReal n + + converge [] = Nothing + converge xs = + Just $ crMemoize (\p -> + case someNatVal (toInteger p) of + Nothing -> error "Data.CReal.Converge p should be non negative" + Just (SomeNat (_ :: Proxy p')) -> + let modifyPrecision = coerce :: [CReal n] -> [CReal p'] + in (last . takeWhilePairwise (/=) . modifyPrecision $ xs) `atPrecision` p) + {-# INLINE converge #-} + + convergeErr _ [] = Nothing + convergeErr err xs = + Just $ crMemoize (\p -> + case someNatVal (toInteger p) of + Nothing -> error "Data.CReal.Converge p should be non negative" + Just (SomeNat (_ :: Proxy p')) -> + let modifyPrecision = coerce :: [CReal n] -> [CReal p'] + modifyFunPrecision = coerce :: (CReal n -> CReal n) -> CReal p' -> CReal p' + es = (modifyFunPrecision err &&& id) <$> modifyPrecision xs + continue (e1, x1) (e2, x2) = if e1 == 0 then x1 /= x2 else e1 > e2 + in (snd . last . takeWhilePairwise continue $ es) `atPrecision` p) + {-# INLINE convergeErr #-} + +takeWhilePairwise :: (a -> a -> Bool) -> [a] -> [a] +takeWhilePairwise p (x1:x2:xs) = if x1 `p` x2 + then x1 : takeWhilePairwise p (x2:xs) + else [x1] +takeWhilePairwise _ xs = xs + +lastMay :: [a] -> Maybe a +lastMay [] = Nothing +lastMay xs = Just (last xs) + + diff --git a/test/Test.hs b/test/Test.hs --- a/test/Test.hs +++ b/test/Test.hs @@ -1,14 +1,19 @@ {-# LANGUAGE TemplateHaskell #-} {-# LANGUAGE DataKinds #-} +{-# LANGUAGE FlexibleContexts #-} {-# OPTIONS_GHC -fno-warn-orphans #-} module Main (main) where +import Data.List (inits) +import Data.Maybe (fromJust) import Data.Ratio ((%)) import Test.Tasty (testGroup, TestTree) import Test.Tasty.QuickCheck (Positive(..), testProperty, (===), Property, (==>), (.&&.), testProperty) import Test.Tasty.TH (defaultMainGenerator) +import Test.Tasty.HUnit (Assertion, (@=?), testCase) +import Data.CReal.Converge import Data.CReal.Internal import Data.CReal.Extra () @@ -66,6 +71,61 @@ prop_showNumDigits (Positive places) x = let s = rationalToDecimal places x in length (dropWhile (/= '.') s) === places + 1 + +-- +-- Testing Data.CReal.Converge +-- + +case_convergeErrEmptyCReal :: Assertion +case_convergeErrEmptyCReal = convergeErr undefined [] @=? (Nothing :: Maybe (CReal 0)) + +case_convergeErrEmptyUnit :: Assertion +case_convergeErrEmptyUnit = convergeErr undefined [] @=? (Nothing :: Maybe ()) + +case_convergeEmptyCReal :: Assertion +case_convergeEmptyCReal = converge [] @=? (Nothing :: Maybe (CReal 0)) + +case_convergeEmptyUnit :: Assertion +case_convergeEmptyUnit = converge [] @=? (Nothing :: Maybe ()) + +prop_convergeCollatzInteger :: Positive Integer -> Property +prop_convergeCollatzInteger (Positive x) = converge (iterate collatz x) === Just 1 + where collatz :: Integer -> Integer + collatz c | c == 1 = 1 + | even c = c `div` 2 + | otherwise = c * 3 + 1 + + +case_convergePointNineRecurringCReal :: Assertion +case_convergePointNineRecurringCReal = (Just 1 :: Maybe (CReal Precision)) @=? + converge (read <$> pointNineRecurring) + where pointNineRecurring = ("0.9" ++) <$> inits (repeat '9') + +prop_convergeErrSqrtCReal :: Positive (CReal Precision) -> Property +prop_convergeErrSqrtCReal (Positive x) = sqrt' (x ^ (2::Int)) === x + where sqrt' x' = let initialGuess = x' + improve y = (y + x' / y) / 2 + err y = abs (x' - y * y) + in fromJust $ convergeErr err (tail $ iterate improve initialGuess) + +-- Test that the behavior when error is too small is correct +prop_convergeErrSmallSqrtCReal :: Positive (CReal Precision) -> Property +prop_convergeErrSmallSqrtCReal (Positive x) = sqrt' (x ^ (2::Int)) === x + where sqrt' x' = let initialGuess = x' + improve y = (y + x' / y) / 2 + err y = abs (x' - y * y) / 128 + in fromJust $ convergeErr err (tail $ iterate improve initialGuess) + +prop_convergeErrSqrtInteger :: Positive Integer -> Property +prop_convergeErrSqrtInteger (Positive x) = sqrt' (x ^ (2::Int)) === x + where sqrt' x' = let initialGuess = x' + improve y = (y + x' `quot` y) `quot` 2 + err y = abs (x' - y * y) + in fromJust $ convergeErr err (tail $ iterate improve initialGuess) + +-- +-- +-- {-# ANN test_boundedFunctions "HLint: ignore Use camelCase" #-} test_boundedFunctions :: [TestTree]