cdar-mBound (empty) → 0.1.0.0
raw patch · 15 files changed
+3291/−0 lines, 15 filesdep +QuickCheckdep +basedep +cdar-mBoundsetup-changed
Dependencies added: QuickCheck, base, cdar-mBound, containers, criterion, deepseq, integer-gmp, parallel, smallcheck, tasty, tasty-hunit, tasty-quickcheck, tasty-smallcheck
Files
- LICENSE +30/−0
- README.md +66/−0
- Setup.hs +2/−0
- bench/bench.hs +205/−0
- cdar-mBound.cabal +158/−0
- examples/ExamplesRump.hs +36/−0
- examples/Logistic.hs +204/−0
- examples/Main.hs +11/−0
- src/Data/CDAR.hs +23/−0
- src/Data/CDAR/Approx.hs +1882/−0
- src/Data/CDAR/Classes.hs +16/−0
- src/Data/CDAR/Dyadic.hs +264/−0
- src/Data/CDAR/Extended.hs +53/−0
- src/Data/CDAR/IntegerLog.hs +43/−0
- test/Test.hs +298/−0
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2014, Jens Blanck++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of Jens Blanck nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ README.md view
@@ -0,0 +1,66 @@+# CDAR (cdar-mBound)++Implementation of computable real numbers, which is also often referred to as Exact Real Arithmetic.++Please see [section at the bottom](#variation-mbound-vs-master) how this branch (mBound) differs from the master branch.++## Computable Real Numbers++Computable real numbers form a countable subset of the real numbers. It contains all real numbers that can be described by a finite program. This includes all numbers that are commonly used in mathematics, including π and e. It is also closed under all field operations and all common transcendental functions such as exponential and trigonometric functions.++### Interval Arithmetic and Computable Real numbers++This implementation used interval arithmetic internally but that doesn't mean that this is _Interval Arithmetic_. Computable real numbers are _exact_, so we are working under the assumption that arbitrary good approximations of input values can be found. This is obviously different from interval arithmetic where one of the main starting points is that input values may not be exact.++### This implementation++Is based on Centred Dyadic Approximations as described in [Blanck 2006](http://cs.swan.ac.uk/~csjens/pdf/centred.pdf).++It is also heavily inspired by the [iRRAM](http://irram.uni-trier.de/) implementation by Norbert Müller. In particular, it uses nested intervals rather than Cauchy sequences with a computable modulus function.++This is implementation should have comparable efficiency to implementations using Cauchy sequences for shallow expressions. However, for deeply nested expressions, such as iterated functions, it should be significantly faster.++Each computable real can be viewed as an infinite list of rapidly shrinking intervals. Thus, laziness of Haskell is used to keep the computations finite.++### Other Related stuff on Hackage++* [**MPFR**](https://hackage.haskell.org/package/hmpfr-0.4.3/docs/Data-Number-MPFR.html) Arbitrary precision floating point numbers with specified rounding modes. While arbitrary precision can be used it is still fixed for a computation so this is still floating point arithmetic, but with larger precision.++* [**AERN2**](https://hackage.haskell.org/package/aern2-real) Computable real numbers and continuous functions using Cauchy sequences. ++ * AERN2 uses on the CDAR type `Approx` as safely-rounded multi-precision floating-point numbers and builds interval arithmetic on top of them.++* [**ireal**](http://hackage.haskell.org/package/ireal) Computable real numbers using Cauchy sequences with a fixed modulus.++* [**constrible**](http://hackage.haskell.org/package/constructible) From the description this appears to be the real closure of the field of rational numbers. This allows for decidable equality, but excludes transcendental functions.++* [**exact-real**](http://hackage.haskell.org/package/exact-real) Computable real numbers using Cauchy sequences with a fixed modulus.++* **ERA** (Can't find a link at the moment) Computable real number using Cauchy sequences with computable modulus.++## Comparison with ireal and AERN2++(https://github.com/michalkonecny/haskell-reals-comparison)++## Installation++Should build under `stack`.++## Motivation++Although the terminology Exact Real Arithmetic promises the ability to compute arbitrarily the result with huge precision, this is not the real strength of this implementation. It is certainly possible to compute 10000 digits of π, and it doesn't take very long. But, even more useful is the fact that if you ask for 10 digits from your resulting computable real number, then you will get 10 correct digits even if the computation has to compute intermediate values with much higher precision.++Main usage cases:+* To provide _exact_ results for important computations.+* To provide numerical computations without having to understand the problems of floating point computations. See [Goldberg 1991](http://dl.acm.org/citation.cfm?id=103163).+* To provide numerical analysis without having to do error analysis.++A word of warning though: Some operations are, by necessity, partial. In particular, comparisons are partial, and so is 1/x near 0.++<!-- ## Examples -->++## Variation mBound vs master++In the mBounds branch, each `Approx` has an additional integer component, called `mBound`, a bound for the number of bits in the integer component `m`. Thus there is a limit on the bit size of the results of arithmetic operations.+In the original variation, the size of `m` could grow arbitrarily for large numbers. +`mBound` plays a similar role for `Approx` numbers as mantissa size does for floating-point numbers.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ bench/bench.hs view
@@ -0,0 +1,205 @@+import Control.Concurrent+import Control.Monad ()+import Criterion.Main++import Data.Bits+import Data.Ratio+import Data.CDAR++logMap :: (Fractional a) => Rational -> a -> a+logMap c x = (fromRational c)*x*(1-x)++logMap2 :: Approx -> Approx+logMap2 Bottom = Bottom+logMap2 (Approx mb m e s) =+ let one = bit (-s)+ t = abs (bit (-s-1) - m)+ in boundErrorTerm $+ if t >= e+ then approxMB mb (m*(one-m)-e^2) (e * abs (one-2*m)) (2*s)+ else approxMB mb (bit (-2*s - 2)) ((t + e)^2) (2*s)++orbit :: Fractional a => Rational -> [a]+orbit c = iterate (logMap c) (fromRational (1%8))++orbit2 :: Rational -> [CR]+orbit2 c = iterate (((fromRational c)*) . CR . fmap logMap2 . unCR) (fromRational (1%8))++suite :: [Benchmark]+suite = [+ bgroup "Logistic 4" [+ bench "Double" $ nf (!! 10000) (orbit 4 :: [Double])+ , bench "CR" $ nf (require 10 . (!! 10000)) (orbit 4 :: [CR])+ , bench "CR (2)" $ nf (require 10 . (!! 10000)) (orbit2 4 :: [CR])+ ],+ bgroup "Logistic 3.5" [+ bench "Double" $ nf (!! 10000) (orbit (7%2) :: [Double])+ , bench "CR" $ nf (require 10 . (!! 10000)) (orbit (7%2) :: [CR])+ , bench "CR (2)" $ nf (require 10 . (!! 10000)) (orbit2 (7%2) :: [CR])+ ],+ bgroup "Logistic 3" [+ bench "Double" $ nf (!! 10000) (orbit 3 :: [Double])+ , bench "CR" $ nf (require 10 . (!! 10000)) (orbit 3 :: [CR])+ , bench "CR (2)" $ nf (require 10 . (!! 10000)) (orbit2 3 :: [CR])+ ]+ ]++newSuite :: [Benchmark]+newSuite =+ [ bgroup "expappr"+ [ bench "double" $ nf exp (1 :: Double)+ , bench "B40" $ nf (expBinarySplittingA 40) 1+ , bench "T40" $ nf (expTaylorA 40) 1+ , bench "T'40" $ nf expTaylorA' (setMB 40 1)+ , bench "CR40" $ nf (require 40 . expCR) 1+ , bench "B400" $ nf (expBinarySplittingA 400) 1+ , bench "T400" $ nf (expTaylorA 400) 1+ , bench "T'400" $ nf expTaylorA' (setMB 400 1)+ , bench "CR400" $ nf (require 400 . expCR) 1+ , bench "B4000" $ nf (expBinarySplittingA 4000) 1+ , bench "T4000" $ nf (expTaylorA 4000) 1+ , bench "T'4000" $ nf (expTaylorA') (setMB 4000 1)+ , bench "CR4000" $ nf (require 4000 . expCR) 1+ ]+ , bgroup "logappr"+ [ bench "doubleLog" $ nf log (1.5 :: Double)+ , bench "40B" $ nf (logBinarySplittingA 40) (Approx 10 3 0 (-1))+ , bench "40T" $ nf (logTaylorA 40) (Approx 10 3 0 (-1))+ -- , bench "40agm" $ nf (logAgmA (-40)) (Approx 10 3 0 (-1))+ , bench "400B" $ nf (logBinarySplittingA 400) (Approx 10 3 0 (-1))+ , bench "400T" $ nf (logTaylorA 400) (Approx 10 3 0 (-1))+ -- , bench "400agm" $ nf (logAgmA (-400)) (Approx 10 3 0 (-1))+ , bench "4000B" $ nf (logBinarySplittingA 4000) (Approx 10 3 0 (-1))+ , bench "4000T" $ nf (logTaylorA 4000) (Approx 10 3 0 (-1))+ -- , bench "4000agm" $ nf (logAgmA (-4000)) (Approx 10 3 0 (-1))+ ]+ , bgroup "exp"+ [ bench "double" $ nf exp (1 :: Double)+ , bench "40" $ nf (require 40 . exp) 1+ , bench "400" $ nf (require 400 . exp) 1+ ]+ , bgroup "log"+ [ bench "double" $ nf log (2 :: Double)+ , bench "40" $ nf (require 40 . log) 2+ , bench "400" $ nf (require 400 . log) 2+ ]+ , bgroup "sin"+ [ bench "double" $ nf sin (1 :: Double)+ , bench "40" $ nf (require 40 . sin) 1+ , bench "40T" $ nf sinTaylorA (setMB 40 1)+ , bench "40CR" $ nf (require 40 . sinCR) 1+ , bench "400" $ nf (require 400 . sin) 1+ , bench "400T" $ nf sinTaylorA (setMB 400 1)+ , bench "400CR" $ nf (require 400 . sinCR) 1+-- , bench "4000" $ nf (require 4000 . sin) 1+ , bench "4000T" $ nf sinTaylorA (setMB 4000 1)+ , bench "4000CR" $ nf (require 4000 . sinCR) 1+ ]+ , bgroup "cos"+ [ bench "double" $ nf cos (1 :: Double)+ , bench "40" $ nf (require 40 . cos) 1+ , bench "400" $ nf (require 400 . cos) 1+ ]+ , bgroup "atan"+ [ bench "double" $ nf atan (1 :: Double)+ , bench "40" $ nf (require 40 . atan) 1+ , bench "400" $ nf (require 400 . atan) 1+ ]+ , bgroup "pi"+ [ bench "double" $ nf (\_ -> pi :: Double) (1 :: Double)+ , bench "40" $ nf (\_ -> require 40 $ pi) 1+ , bench "400" $ nf (\_ -> require 400 $ pi) 1+ ]+ , env setupEnvPi $ \ ~(pi40,pi400,pi4000) ->+ bgroup "ccc2015"+ [ bench "pi D" $ nf (\_ -> pi :: Double) (1 :: Double)+ , bench "pi 40" $ nf (\_ -> require 40 $ pi) 1+ , bench "pi 400" $ nf (\_ -> require 400 $ pi) 1+ , bench "pi 4000" $ nf (\_ -> require 4000 $ pi) 1+ , bench "+ D" $ nf (\x -> x+x) (pi :: Double)+ , bench "+ 40" $ nf (\x -> x+x) pi40+ , bench "+ 400" $ nf (\x -> x+x) pi400+ , bench "+ 4000" $ nf (\x -> x+x) pi4000+ , bench "* D" $ nf (\x -> x*x) (pi :: Double)+ , bench "* 40" $ nf (\x -> x*x) pi40+ , bench "* 400" $ nf (\x -> x*x) pi400+ , bench "* 4000" $ nf (\x -> x*x) pi4000+ , bench "* D" $ nf (1/) (pi :: Double)+ , bench "* 40" $ nf (1/) pi40+ , bench "* 400" $ nf (1/) pi400+ , bench "* 4000" $ nf (1/) pi4000+ , bench "exp D" $ nf exp (1 :: Double)+ , bench "exp 40" $ nf (require 40 . exp) 1+ , bench "exp 400" $ nf (require 400 . exp) 1+ , bench "exp 4000" $ nf (require 4000 . exp) 1+ , bench "log D" $ nf log (2 :: Double)+ , bench "log 40" $ nf (require 40 . log) 2+ , bench "log 400" $ nf (require 400 . log) 2+ , bench "log 4000" $ nf (require 4000 . log) 2+ , bench "sin D" $ nf sin (1 :: Double)+ , bench "sin 40" $ nf (require 40 . sinCR) 1+ , bench "sin 400" $ nf (require 400 . sinCR) 1+ , bench "sin 4000" $ nf (require 4000 . sinCR) 1+ , bench "cos D" $ nf cos (1 :: Double)+ , bench "cos 40" $ nf (require 40 . cosCR) 1+ , bench "cos 400" $ nf (require 400 . cosCR) 1+ , bench "cos 4000" $ nf (require 4000 . cosCR) 1+ , bench "atan D" $ nf atan (1 :: Double)+ , bench "atan 40" $ nf (require 40 . atanCR) 1+ , bench "atan 400" $ nf (require 400 . atanCR) 1+ , bench "atan 4000" $ nf (require 4000 . atanCR) 1+ ]+ , env setupEnv $ \ ~(pi1,pi2) ->+ bgroup "elementary Approx"+ [ bench "+ double" $ nf (\x -> x+x) (pi :: Double)+ , bench "* double" $ nf (\x -> x*x) (pi :: Double)+ , bench "rec double" $ nf (1/) (pi :: Double)+ , bench "sqrt double" $ nf (sqrt) (pi :: Double)+ , bench "+ 50" $ nf (\x -> x+x) pi1+ , bench "* 50" $ nf (\x -> x*x) pi1+ , bench "rec 50" $ nf recipA (setMB 50 pi1)+ , bench "sqrt 50" $ nf sqrtA (setMB 50 pi1)+ , bench "log 50" $ nf (logA) (setMB 50 pi1)+ , bench "+ 1000" $ nf (\x -> x+x) pi2+ , bench "* 1000" $ nf (\x -> x*x) pi2+ , bench "rec 1000" $ nf recipA (setMB 1000 pi2)+ , bench "sqrt 1000" $ nf sqrtA (setMB 1000 pi2)+ , bench "sqrtRec 1000" $ nf (sqrtRecA 1000) pi2+ , bench "log 1000" $ nf logA (setMB 1000 pi2)+ ]+ ]++setupEnv :: IO (Approx, Approx)+setupEnv = return . (\a -> (limitAndBound 50 a, a)) . limitAndBound 1000 . require 1000 $ pi+setupEnvPi :: IO (Approx, Approx, Approx)+setupEnvPi = return . (\a -> (limitAndBound 40 a, limitAndBound 400 a, limitAndBound 4000 a)) . require 4000 $ pi++threadSuite :: MVar Approx -> MVar Approx -> [Benchmark]+threadSuite u v =+ [ bgroup "thread"+ [ bench "communicate" $ nfIO (do putMVar u (Approx 30 145324626 123 (-30)); a <- takeMVar v; return a)+ ]+ ]++main :: IO ()+main = defaultMain $ newSuite -- ++ suite++{- Are threads making criterion confused, times seem ok, but reported as unreliable.++threadSuite :: MVar Approx -> MVar Approx -> [Benchmark]+threadSuite u v =+ [ bgroup "thread"+ [ bench "communicate" $ nfIO (do putMVar u pi2; a <- takeMVar v; return a)+ ]+ ]++main :: IO ()+main = do+ u <- newEmptyMVar+ v <- newEmptyMVar+ forkIO $ server u v+ defaultMain $ threadSuite u v+ where server u v = forever $ do+ a <- takeMVar u+ putMVar v (a + 1)+-}
+ cdar-mBound.cabal view
@@ -0,0 +1,158 @@+-- Initial cdar.cabal generated by cabal init. For further documentation, +-- see http://haskell.org/cabal/users-guide/++name: cdar-mBound+version: 0.1.0.0+synopsis: Exact real arithmetic using Centred Dyadic Approximations+description: Please see https://github.com/michalkonecny/cdar/tree/mBound#readme+license: BSD3+license-file: LICENSE+author: Jens Blanck, Michal Konecny+maintainer: mikkonecny@gmail.com+-- copyright: +category: Data+build-type: Simple+extra-source-files: README.md+cabal-version: >=1.10++source-repository head+ type: git+ location: https://github.com/michalkonecny/cdar+ branch: mBound++library+ exposed-modules:+ Data.CDAR+ , Data.CDAR.Classes+ , Data.CDAR.Extended+-- , Data.CDAR.POrd+ , Data.CDAR.Approx+ , Data.CDAR.Dyadic+-- , Data.CDAR.Interval+ , Data.CDAR.IntegerLog+ -- other-modules: + other-extensions: BangPatterns, GADTs, TypeSynonymInstances, FlexibleInstances, DeriveFunctor, MagicHash+ build-depends:+ base >= 4.7 && < 4.15,+ containers,deepseq,integer-gmp+ -- base >=4.7 && <4.10+ -- , containers >=0.5 && <0.6+ -- , deepseq >=1.3 && <1.5+ -- , integer-gmp >=0.5 && <1.1+ -- , QuickCheck >=2.6 && <3+ -- , parallel >=3.2 && <3.3+ hs-source-dirs: src+ default-language: Haskell2010+ ghc-options: -Wall -fno-warn-type-defaults -O2++-- test-suite test+-- type: exitcode-stdio-1.0+-- main-is: QT.hs+-- build-depends: base, QuickCheck, cdar+-- hs-source-dirs: test+-- default-language: Haskell2010++test-suite test+ type: exitcode-stdio-1.0+ main-is: Test.hs+ build-depends:+ base, containers, deepseq, integer-gmp, QuickCheck, smallcheck, tasty-hunit, tasty-quickcheck, tasty-smallcheck, tasty, cdar-mBound+ -- base >=4.7 && <4.10+ -- , containers >=0.5 && <0.6+ -- , deepseq >=1.3 && <1.5+ -- , integer-gmp >=0.5 && <1.1+ -- , QuickCheck >=2.6 && <3+ -- , smallcheck >= 1.1 && < 1.2+ -- , tasty-hunit >= 0.8 && < 0.10+ -- , tasty-quickcheck >= 0.8 && < 0.9+ -- , tasty-smallcheck >= 0.8 && < 0.9+ -- , tasty >= 0.8 && < 0.12+ -- , cdar >= 0.1+ -- , parallel >=3.2 && <3.3+ hs-source-dirs: test+ default-language: Haskell2010++-- Benchmark bench-transcend+-- type: exitcode-stdio-1.0+-- main-is: bench-transcend.hs+-- build-depends:+-- base >=4.7 && <4.9+-- , containers >=0.5 && <0.6+-- , deepseq >=1.3 && <1.5+-- , integer-gmp >=0.5 && <1.1+-- , cdar >= 0.1+-- , criterion >= 1.1 && <1.2+-- , parallel >=3.2 && <3.3+-- hs-source-dirs: bench+-- default-language: Haskell2010++Benchmark bench+ type: exitcode-stdio-1.0+ main-is: bench.hs+ build-depends:+ base, containers, deepseq, integer-gmp, cdar-mBound, criterion, parallel+ -- base >=4.7 && <4.10+ -- , containers >=0.5 && <0.6+ -- , deepseq >=1.3 && <1.5+ -- , integer-gmp >=0.5 && <1.1+ -- , cdar >= 0.1+ -- , criterion >= 1.1 && <1.2+ -- , parallel >=3.2 && <3.3+ hs-source-dirs: bench+ default-language: Haskell2010+ ghc-options: -O2 -rtsopts++-- Benchmark bench-thread+-- type: exitcode-stdio-1.0+-- main-is: bench-thread.hs+-- build-depends:+-- base >=4.7 && <4.10+-- , containers >=0.5 && <0.6+-- , deepseq >=1.3 && <1.5+-- , integer-gmp >=0.5 && <1.1+-- , cdar >= 0.1+-- , criterion >= 1.1 && <1.2+-- , parallel >=3.2 && <3.3+-- hs-source-dirs: bench+-- default-language: Haskell2010+-- ghc-options: -threaded -O2 -rtsopts++-- Benchmark BinSplitPar+-- type: exitcode-stdio-1.0+-- main-is: BinSplitPar.hs+-- build-depends:+-- base >=4.7 && <4.10+-- , containers >=0.5 && <0.6+-- , deepseq >=1.3 && <1.5+-- , integer-gmp >=0.5 && <1.1+-- , cdar >= 0.1+-- , criterion >= 1.1 && <1.2+-- , parallel >=3.2 && <3.3+-- , monad-par +-- hs-source-dirs: bench+-- default-language: Haskell2010+-- ghc-options: -threaded -O2 -rtsopts++Executable examples+ hs-source-dirs: examples+ main-is: Main.hs+ other-modules: ExamplesRump+ build-depends:+ base, containers, deepseq, integer-gmp, cdar-mBound+ -- base >= 4 && < 5+ -- , cdar >= 0.1+ -- , containers >= 0.5 && < 0.6+ -- , deepseq >=1.3 && <1.5+ -- , integer-gmp >=0.5 && <1.1+ -- , cdar >= 0.1+ -- , criterion >= 1.1 && <1.2+ -- , parallel >=3.2 && <3.3+ default-language: Haskell2010++Executable logistic+ hs-source-dirs: examples+ main-is: Logistic.hs+ build-depends: base, cdar-mBound+ default-language: Haskell2010++--SinBorwein, Time,
+ examples/ExamplesRump.hs view
@@ -0,0 +1,36 @@+{-# LANGUAGE NoMonomorphismRestriction #-}+{-|+Two examples given to me by S Rump that challenge the precision of floating+point computations.++The first example gives the wrong sign with Double arithmetic, but even the+first approximation of a CReal is correct.++The second example gives the correct value with Double, but we can see that it+is difficult as it takes some extra precision to produce the right result.++The usage of showCRealN is to show the underlying computations; it is not+intended way to use CReal.+-}+module ExamplesRump (examplesRump) where++import Data.CDAR++a :: Num t => t+a = 77617+b :: Num t => t+b = 33096+p :: Fractional a => a+p = 206987/2048+q :: Fractional a => a+q = 119504/2048+r :: Fractional a => a+r = p^^3*(p^^16+6561*q^^16-17496*p^^2*q^^14+20412*p^^4*q^^12-13608*p^^6*q^^10+5670*p^^8*q^^8-1512*p^^10*q^^6+252*p^^12*q^^4-24*p^^14*q^^2) - q++examplesRump :: IO ()+examplesRump = do+ print $ 21*b*b-2*a*a+55*b*b*b*b-10*a*a*b*b+a/(2*b)+ putStrLn . showCRN 3 $ 21*b*b-2*a*a+55*b*b*b*b-10*a*a*b*b+a/(2*b)+ print r+ putStrLn . showCRN 3 $ r+
+ examples/Logistic.hs view
@@ -0,0 +1,204 @@+{-# LANGUAGE NoMonomorphismRestriction #-}+module Main where++import Data.CDAR+import Data.Bits+import Data.Ratio+import System.Environment (getArgs)+import Text.Printf++main :: IO ()+main = do+ args <- getArgs+ let n = read $ args !! 0 :: Int+ let c = if length args > 1 then+ let (i,f) = break (== '.') $ args !! 1+ f' = tail f+ l = length f'+ in (read (i++f') :: CR) / (10^l)+ else scale 1961 (-9)+ printf "Constant is %s.\n" $ showA . require 10 $ c + printf "\nComputing %d iterations of the map.\nIt may take some seconds...\n" n+ printf "\nAfter %d iterations the value is %s\n" n $ showA . limitSize 34 . require 34 $ (orbitC2 c)!!n++logMap :: (Fractional a) => a -> a+logMap x = 4*x*(1-x)++logMap2 :: Approx -> Approx+logMap2 Bottom = Bottom+logMap2 (Approx mb m e s) =+ let one = bit (-s)+ t = abs (bit (-s-1) - m)+ in boundErrorTerm $+ if t >= e+ then approxMB mb (m*(one-m)-e^2) (e * abs (one-2*m)) (2*s)+ else approxMB mb (bit (-2*s - 2)) ((t + e)^2) (2*s)++orbit :: Fractional a => [a]+orbit = iterate logMap (fromRational (1%8))++orbit1 :: [CR]+orbit1 = iterate (polynomial [0,4,4]) (fromRational (1%8))++orbit2 :: [CR]+orbit2 = iterate (CR . fmap ((4*) . logMap2) . unCR) (fromRational (1%8))++logMapC :: Fractional a => a -> a -> a+logMapC c x = c*x*(1-x)++orbitC :: CR -> [CR]+orbitC c = iterate (logMapC c) (1/8)++orbitC1 :: CR -> [CR]+orbitC1 c = iterate ((c*) . polynomial [0,1,-1]) (1/8)++orbitC2 :: CR -> [CR]+orbitC2 c = iterate ((c*) . CR . fmap logMap2 . unCR) (1/8)++{-++-- Using orbitC2:++jens@cspcjb:~/cdar> time stack exec logistic 10 3.82 +Constant is 3.82000000000000000000000~.++Computing 10 iterations of the map.+It may take some seconds...++After 10 iterations the value is 0.1650561779~++real 0m0.114s+user 0m0.094s+sys 0m0.022s+jens@cspcjb:~/cdar> time stack exec logistic 100 3.82+Constant is 3.82000000000000000000000~.++Computing 100 iterations of the map.+It may take some seconds...++After 100 iterations the value is 0.4309357855~++real 0m0.113s+user 0m0.087s+sys 0m0.033s+jens@cspcjb:~/cdar> time stack exec logistic 1000 3.82+Constant is 3.82000000000000000000000~.++Computing 1000 iterations of the map.+It may take some seconds...++After 1000 iterations the value is 0.8687524058~++real 0m0.121s+user 0m0.093s+sys 0m0.030s+jens@cspcjb:~/cdar> time stack exec logistic 10000 3.82+Constant is 3.82000000000000000000000~.++Computing 10000 iterations of the map.+It may take some seconds...++After 10000 iterations the value is 0.9216483032~++real 0m0.342s+user 0m0.307s+sys 0m0.037s+jens@cspcjb:~/cdar> time stack exec logistic 50000 3.82+Constant is 3.82000000000000000000000~.++Computing 50000 iterations of the map.+It may take some seconds...++After 50000 iterations the value is 0.2527618912~++real 0m6.753s+user 0m6.679s+sys 0m0.075s+jens@cspcjb:~/cdar> time stack exec logistic 100000 3.82+Constant is 3.82000000000000000000000~.++Computing 100000 iterations of the map.+It may take some seconds...++After 100000 iterations the value is 0.6195746474~++real 0m41.659s+user 0m41.501s+sys 0m0.159s++---++-- Using orbitC:++jens@cspcjb:~/cdar> time stack exec logistic 10 3.82 +Constant is 3.82000000000000000000000~.++Computing 10 iterations of the map.+It may take some seconds...++After 10 iterations the value is 0.1650561779~++real 0m0.110s+user 0m0.093s+sys 0m0.023s+jens@cspcjb:~/cdar> time stack exec logistic 100 3.82+Constant is 3.82000000000000000000000~.++Computing 100 iterations of the map.+It may take some seconds...++After 100 iterations the value is 0.4309357855~++real 0m0.113s+user 0m0.089s+sys 0m0.031s+jens@cspcjb:~/cdar> time stack exec logistic 1000 3.82+Constant is 3.82000000000000000000000~.++Computing 1000 iterations of the map.+It may take some seconds...++After 1000 iterations the value is 0.8687524058~++real 0m0.128s+user 0m0.081s+sys 0m0.049s+jens@cspcjb:~/cdar> time stack exec logistic 10000 3.82+Constant is 3.82000000000000000000000~.++Computing 10000 iterations of the map.+It may take some seconds...++After 10000 iterations the value is 0.9216483032~++real 0m1.044s+user 0m1.009s+sys 0m0.036s+jens@cspcjb:~/cdar> time stack exec logistic 50000 3.82+Constant is 3.82000000000000000000000~.++Computing 50000 iterations of the map.+It may take some seconds...++After 50000 iterations the value is 0.2527618912~++real 0m36.025s+user 0m35.848s+sys 0m0.176s+jens@cspcjb:~/cdar> time stack exec logistic 100000 3.82+Constant is 3.82000000000000000000000~.++Computing 100000 iterations of the map.+It may take some seconds...++After 100000 iterations the value is 0.6195746474~++real 3m47.789s+user 3m46.463s+sys 0m1.324s++---++model name : Intel(R) Core(TM) i7-4790 CPU @ 3.60GHz++-}
+ examples/Main.hs view
@@ -0,0 +1,11 @@+module Main where++import Data.CDAR+import ExamplesRump+--import Criterion.Main+--import Criterion.Config++main :: IO ()+main = do+ putStrLn . showA . require 50 . sin $ pi / 2+ examplesRump
+ src/Data/CDAR.hs view
@@ -0,0 +1,23 @@+{-| Convenience module importing and reexporting the 'CReal' data type of+computable real numbers and the underlying data type including the centred+dyadic approximations 'Approx'.++Modules imported are:++- "Data.CDAR.Approx"+- "Data.CDAR.Classes"+- "Data.CDAR.Dyadic"+- "Data.CDAR.Extended"+- "Data.CDAR.IntegerLog"+-}+module Data.CDAR (module Approx+ ,module Classes+ ,module Dyadic+ ,module Extended+ ,module IntegerLog) where++import Data.CDAR.Approx as Approx+import Data.CDAR.Classes as Classes+import Data.CDAR.Dyadic as Dyadic+import Data.CDAR.Extended as Extended+import Data.CDAR.IntegerLog as IntegerLog
+ src/Data/CDAR/Approx.hs view
@@ -0,0 +1,1882 @@+{-# LANGUAGE BangPatterns,GADTs,TypeSynonymInstances,FlexibleInstances #-}+{-|+= Computable Real Arithmetic+This module provides the data type 'CR' that implements the real closed field of computable real numbers.++== Centred Dyadic Approximations+The computable reals are realised as lists of rapidly shrinking intervals. The intervals used here are centred dyadic intervals, implemented here as the data type 'Approx'.++For more information on the theoretical aspects see <http://cs.swan.ac.uk/~csjens/pdf/centred.pdf>.+-}+module Data.CDAR.Approx (Approx(..)+ ,CR(..)+-- ,errorBits+-- ,errorBound+-- ,defaultPrecision+ ,Precision+ ,showA+ ,showInBaseA+ ,mBound+ ,mapMB+ ,setMB+ ,enforceMB+ ,approxAutoMB+ ,approxMB+ ,approxMB2+ ,endToApprox+ ,lowerBound+ ,upperBound+ ,lowerA+ ,upperA+ ,centre+ ,centreA+ ,radius+ ,diameter+ ,exact+ ,approximatedBy+ ,better+ ,fromDyadic+ ,fromDyadicMB+ ,toApprox+ ,toApproxMB+ ,recipA+ ,divAInteger+ ,modA+ ,divModA+ ,toDoubleA+ ,precision+ ,significance+ ,boundErrorTerm+ ,boundErrorTermMB+ ,limitSize+ ,checkPrecisionLeft+ ,limitAndBound+ ,limitAndBoundMB+ ,unionA+ ,intersectionA+ ,consistentA+ ,poly+ ,pow+ ,powers+ ,sqrtHeronA+ ,sqrtA+ ,sqrtRecA+ ,findStartingValues+ ,sqrtD+ ,shiftD+ ,sqrA+ ,log2Factorials+ ,taylorA+ ,expA+ ,expBinarySplittingA+ ,expTaylorA+ ,expTaylorA'+ ,logA+ ,logInternal+ ,logBinarySplittingA+ ,logTaylorA+ ,sinTaylorA+ ,sinTaylorRed1A+ ,sinTaylorRed2A+ ,sinA+ ,cosA+ ,atanA+ ,sinBinarySplittingA+ ,cosBinarySplittingA+ ,atanBinarySplittingA+ ,atanTaylorA+ ,piRaw+ ,piA+ ,piMachinA+ ,piBorweinA+ ,piAgmA+ ,log2A+ ,lnSuperSizeUnknownPi+ ,logAgmA+ ,agmLn+ ,showCRN+ ,showCR+ ,ok+ ,require+ ,toDouble+ ,fromDouble+ ,fromDoubleAsExactValue+ ,polynomial+ ,taylorCR+ ,atanCR+ ,piCRMachin+ ,piMachinCR+ ,piBorweinCR+ ,piBinSplitCR+ ,ln2+ ,sinCRTaylor+ ,sinCR+ ,cosCR+ ,sqrtCR+ ,expCR+ ) where++import Control.Applicative (ZipList (..))+import Control.DeepSeq+import Control.Exception+import Data.Bits+import Data.Char (isDigit)+import Data.CDAR.Classes+import Data.CDAR.Dyadic+import Data.CDAR.Extended+import Data.CDAR.IntegerLog+import Data.Char (intToDigit)+import Data.List (findIndex, intersperse, transpose, unfoldr, zipWith4)+import Data.Ratio++-- import Debug.Trace++-- |A type synonym. Used to denote number of bits after binary point.+type Precision = Int++{-|+= Centred Dyadic Approximations+There are two constructors for approximations:++- 'Approx' is encodes some finite interval with dyadic endpoints. A real+ number is /approximated/ by the approximation is it belongs to the interval.+- 'Bottom' is the trivial approximation that approximates all real numbers.++The four fields of an @Approx m e s@ should be thought of as:++[@mb@] the midpoint bound, ie maximum bits available for the midpoint+[@m@] the midpoint+[@e@] the error term+[@s@] the exponent++Thus, a value @Approx p m e s@ is to be interpreted as the interval+[(m-e)*2^s, (m+e)*2^s] where |m| <= 2^p.++== Centred intervals+We have opted for a centred representation of the intervals. It is also+possible to represent the endpoints as 'Dyadic' numbers. The rationale for a+centred repersentation is that we often normalise an approximation @Approx p m e+s@ so that @e@ is limited in size. This allows many operations to only work on+one large number @m@.++== Potential for overflow+Since the last field (the exponent) is only an 'Int' it may overflow. This is+an optimisation that was adopted since it seems unlikely that overflow in a 64+bit Int exponent would occur. In a 32 bit system, this is potentially an+issue.++The 'Integer' data type is unbonded, but is, of course, bounded by the+available memory available in the computer. No attempt has been made to check+for exhausted memory.++== Approximations as a Domain++Ordered by reverse inclusion the dyadic intervals encoded by the 'Approx'+approximations (including 'Bottom') constitute the compact elements of a Scott+domain /D/. (It is a substructure of the (algebraic) interval domain.)+We will identify our approximations with the compact elements of /D/.++Increasing sequences in /D/ have suprema. A sequence /converges/ if the length+of the approximations tend to zero. The supremum of a converging sequence is a+singleton set containing a real number. Let ρ be the map taking a converging+sequence to the unique real number in the supremum. The computations on+(computable) real numbers is via this representation map ρ.++There is no check that the sequences we have are in fact increasing, but we+are assuming that all sequences are pairwise consistent. We can thus create an+increasing sequence by considering the sequence of finite suprema. For+correctness, we have to ensure that all operations done on consistent+sequences result in consistent sequences. If non-consistent sequences are+somehow input we can make no guarantees at all about the computed value.++Note, that we cannot ensure that converging sequences are mapped to converging+sequences because of properties of computable real arithmetic. In particular,+at any discuntinuity, it is impossible to compute a converging sequence.+-}+data Approx = Approx Int Integer Integer Int+ | Bottom+ deriving (Read,Show)++instance NFData Approx where+ rnf Bottom = ()+ rnf (Approx b m e s) = rnf b `seq` rnf m `seq` rnf e `seq` rnf s++instance Scalable Approx where+ scale Bottom _ = Bottom+ scale (Approx b m e s) n = Approx b m e (s+n)++instance Scalable CR where+ scale (CR x) n = CR $ flip scale n <$> x++{-|+=The Computable Real data type++Computable reals are realised as infinite sequences of centred dyadic+representations.++All approximations in such a sequence should be pairwise consistent, i.e.,+have a non-empty intersection. However, there is no check that this is+actually the case.++If the diameter of the approximations tend to zero we say that the sequences+converges to the unique real number in the intersection of all intervals.+Given the domain /D/ of approximations described in 'Approx', we have a+representation (a retraction) ρ from the converging sequences in /D/ to ℝ.+Some operations on computable reals are partial, notably equality and+rnf m `seq` there is no guarantee that a+rnf m `seq` +rnf m `seq` +rnf m `seq` there is a bound on how much effort is+rnf m `seq` mation. For involved computations it is+rnf m `seq` proximations are trivial, i.e.,+rnf m `seq` ually converge, it will generate proper+rnf m `seq` initial trivial approximations.+rnf m `seq` +The amount of added effort in each iteration is rather substantial so the+expected precision of approximations increase very quickly.++==The actual data type++In fact, the type 'CR' is a newtype of 'ZipList' 'Approx' in the+implementation of infinite sequences of approximations, as that allows for+applicative style. Hopefully, it is not needed to access the internal+representation of 'CR' directly.+-}+newtype CR = CR {unCR :: ZipList Approx}++-- |Number of bits that error term is allowed to take up. A larger size allows+-- for more precise but slightly more costly computations. The value here is+-- suggested by test runs.+errorBits :: Int+errorBits = 10++errorBound :: Integer+errorBound = 2^errorBits++errorBitsMB :: Int+errorBitsMB = 1++errorBoundMB :: Integer+errorBoundMB = 2^errorBitsMB+++-- |The default cutoff for diverging computations. May well be chosen much+-- smaller. 31 corresponds to about 10 decimal places.+defaultPrecision :: Precision+defaultPrecision = 31++{-|++Gives a decimal representation of an approximation. It tries to give as many+decimal digits as possible given the precision of the approximation. The+representation may be wrong by 1 ulp (unit in last place). If the value is not+exact the representation will be followed by @~@.++The representation is not always intuitive:++>>> showA (Approx 1 1 0)+"1.~"++The meaning of the above is that it is 1, but then the added @~@ (which must+be after the decimal point) means that the last position may be off by 1,+i.e., it could be down to 0 or up to 2. And [0,2] is indeed the range encoded+by the above approximation.+-}+showA :: Approx -> String+showA = showInBaseA 10++-- |Similar to 'showA' but can generate representations in other bases (<= 16).+{- am is the absolute value of the significand+ b corresponds to the value 1 with respect to the shift s -- this is used to find the digits in the auxiliary functions+ i is the integral part of am+ f is the fractional part of am+ i' and f' are the integral and fractional parts relevant for near zero approximations+ e' is the error term shifted appropriately when s positive, also set to at least 1+ (otherwise odd bases will yield infinite expansions)+-}+showInBaseA :: Int -> Approx -> String+showInBaseA _ Bottom = "⊥"+showInBaseA base (Approx _ m e s)+ | e == 0 && (even base || s >= 0)+ = sign ++ showExactA base b i f+ | am < e = "±" ++ showNearZeroA base b i' f'+ | otherwise = sign ++ showInexactA base b i f e'+ where b = bit (max 0 (-s))+ am = abs m+ i = shift am s+ e' = max 1 $ shift e (max 0 s)+ f = am .&. (b-1)+ i' = shift (am+e) s+ f' = (am+e) .&. (b-1)+ sign = if m < 0 then "-" else ""++showExactA :: Int -> Integer -> Integer -> Integer -> String+showExactA base b i f =+ let g i' = let (q,r) = quotRem i' (fromIntegral base)+ in if i' == 0 then Nothing+ else Just (intToDigit (fromIntegral r), q)+ ip = reverse (unfoldr g i)+ h f' = let (q,r) = quotRem ((fromIntegral base)*f') b+ in if f' == 0 then Nothing+ else Just (intToDigit (fromIntegral q), r)+ fp = unfoldr h f+ in (if null ip then "0" else ip)+ ++ (if null fp then "" else ".")+ ++ fp++showNearZeroA :: Int -> Integer -> Integer -> Integer -> String+showNearZeroA base b i f =+ let s = showExactA base b i f+ t = takeWhile (flip elem "0.~") s+ u = takeWhile (/= '.') s+ in if null t+ then replicate (length u) '~'+ else t ++ "~"++showInexactA :: Int -> Integer -> Integer -> Integer -> Integer -> String+showInexactA base b i f e =+ let g (0,b',f') = if b' < f'+e+ then Just ('1', (0, (fromIntegral base)*b', f'))+ else Nothing+ g (n,b',f') = let (q,r) = quotRem n (fromIntegral base)+ z = (q, (fromIntegral base)*b', r*b'+f')+ in if e+f' <= b'+ then Just (intToDigit (fromIntegral r), z)+ else if e <= min f' b'+ then Just (intToDigit ((fromIntegral r + 1) `rem` (fromIntegral base)), z)+ else Just ('~', z)+ intRev = unfoldr g (i,b,f)+ noFrac = case intRev of+ [] -> False+ (x:_) -> x == '~'+ int = if null intRev then "0" else reverse intRev+ h (f',err) = let (q,r) = quotRem ((fromIntegral base)*f') b+ err' = (fromIntegral base)*err+ z = (r, err')+ in if err' + r <= b+ then Just (intToDigit (fromIntegral q), z)+ else if err' <= min r b+ then Just (intToDigit ((fromIntegral q + 1) `rem` (fromIntegral base)), z)+ else Nothing+ frac = unfoldr h (f,e)+ in int ++ if noFrac+ then ""+ else "." ++ frac ++ "~"++{-|+ Give the "bound on midpoint bit-size" component of an 'Approx'.++ The midpoint coponent should always be bounded by this as follows:+ @ abs m <= 2^mb@.+-}+mBound :: Approx -> Int+mBound (Approx mb _ _ _) = mb+mBound Bottom = error "mBound Bottom"++approxAutoMB :: Integer -> Integer -> Int -> Approx+approxAutoMB m e s = Approx mb m e s+ where+ ame = (abs m) + e+ mb | ame <= 4 = 2+ | otherwise = 1 + integerLog2 (ame - 1)+++mapMB :: (Int -> Int) -> Approx -> Approx+mapMB f (Approx mb m e s) = approxMB (f mb) m e s+mapMB _f Bottom = Bottom++setMB :: Int -> Approx -> Approx+setMB mb = mapMB (const mb)+++approxMB :: Int -> Integer -> Integer -> Int -> Approx+approxMB mb m e s = + enforceMB $ Approx mb m e s++approxMB2 :: Int -> Int -> Integer -> Integer -> Int -> Approx+approxMB2 mb1 mb2 m e s = + enforceMB $ Approx (mb1 `max` mb2) m e s++enforceMB :: Approx -> Approx+enforceMB Bottom = Bottom+enforceMB a@(Approx mb m e s)+ | m_size <= mb = a+ | abs m <= 1 = a+ | otherwise = Approx mb m' e'' (s + d)+ where+ m_size = 1+integerLog2 (abs m - 1) -- |m| <= 2^m_size+ d = m_size - mb+ m' = unsafeShiftR m d -- we have: m' * 2^d <= m+ e' = 1 + (unsafeShiftR (e-1) d) -- we have: 0 <= e <= e' * 2^d+ e''| m == unsafeShiftL m' d = e' -- no loss of information+ | otherwise = 1 + e'+++-- |Construct a centred approximation from the end-points.+endToApprox :: Int -> Extended Dyadic -> Extended Dyadic -> Approx+endToApprox mb (Finite l) (Finite u)+ | u < l = Bottom -- Might be better with a signalling error.+ | otherwise =+ let a@(m:^s) = scale (l + u) (-1)+ (n:^t) = u-a+ r = min s t+ m' = unsafeShiftL m (s-r)+ n' = unsafeShiftL n (t-r)+ in (approxMB mb m' n' r)+endToApprox _ _ _ = Bottom++-- Interval operations+-- |Gives the lower bound of an approximation as an 'Extended' 'Dyadic' number.+lowerBound :: Approx -> Extended Dyadic+lowerBound (Approx _ m e s) = Finite ((m-e):^s)+lowerBound Bottom = NegInf++-- |Gives the upper bound of an approximation as an 'Extended' 'Dyadic' number.+upperBound :: Approx -> Extended Dyadic+upperBound (Approx _ m e s) = Finite ((m+e):^s)+upperBound Bottom = PosInf++-- |Gives the lower bound of an 'Approx' as an exact 'Approx'.+lowerA :: Approx -> Approx+lowerA Bottom = Bottom+lowerA (Approx mb m e s) = Approx mb (m-e) 0 s++-- |Gives the upper bound of an 'Approx' as an exact 'Approx'.+upperA :: Approx -> Approx+upperA Bottom = Bottom+upperA (Approx mb m e s) = Approx mb (m+e) 0 s++-- |Gives the mid-point of an approximation as a 'Maybe' 'Dyadic' number.+centre :: Approx -> Maybe Dyadic+centre (Approx _ m _ s) = Just (m:^s)+centre _ = Nothing++-- |Gives the centre of an 'Approx' as an exact 'Approx'.+centreA :: Approx -> Approx+centreA Bottom = Bottom+centreA (Approx mb m _e s) = Approx mb m 0 s++-- |Gives the radius of an approximation as a 'Dyadic' number. Currently a+-- partial function. Should be made to return an 'Extended' 'Dyadic'.+radius :: Approx -> Extended Dyadic+radius (Approx _ _ e s) = Finite (e:^s)+radius _ = PosInf++-- |Gives the lower bound of an approximation as an 'Extended' 'Dyadic' number.+diameter :: Approx -> Extended Dyadic+diameter (Approx _ _ e s) = Finite $ 2 * (e:^s)+diameter _ = PosInf++-- |Returns 'True' if the approximation is exact, i.e., it's diameter is 0.+exact :: Approx -> Bool+exact (Approx _ _ 0 _) = True+exact _ = False++-- |Checks if a number is approximated by an approximation, i.e., if it+-- belongs to the interval encoded by the approximation.+approximatedBy :: Real a => a -> Approx -> Bool+_ `approximatedBy` Bottom = True+r `approximatedBy` d =+ let r' = toRational r+ in toRational (lowerBound d) <= r' && r' <= toRational (upperBound d)++-- |A partial order on approximations. The first approximation is better than+-- the second if it is a sub-interval of the second.+better :: Approx -> Approx -> Bool+d `better` e = lowerBound d >= lowerBound e &&+ upperBound d <= upperBound e++-- |Turns a 'Dyadic' number into an exact approximation.+fromDyadic :: Dyadic -> Approx+fromDyadic (m:^s) = approxAutoMB m 0 s++-- |Turns a 'Dyadic' number into an exact approximation.+fromDyadicMB :: Int -> Dyadic -> Approx+fromDyadicMB mb (m:^s) = approxMB mb m 0 s++-- |Two approximations are equal if they encode the same interval.+instance Eq Approx where+ (Approx _ m e s) == (Approx _ n f t)+ | s >= t = let k = s-t+ in unsafeShiftL m k == n && unsafeShiftL e k == f+ | s < t = let k = t-s+ in m == unsafeShiftL n k && e == unsafeShiftL f k+ Bottom == Bottom = True+ _ == _ = False++-- |Not a sensible instance. Just used to allow to allow enumerating integers+-- using \'..\' notation.+instance Enum Approx where+ toEnum n = Approx 64 (fromIntegral n) 0 0+ fromEnum (Approx _ m _ s) = fromIntegral $ shift m s+ fromEnum Bottom = 0++instance Num Approx where+ (Approx mb1 m e s) + (Approx mb2 n f t)+ | s >= t = let k = s-t+ in approxMB2 mb1 mb2 (unsafeShiftL m k + n) (unsafeShiftL e k + f) t+ | s < t = let k = t-s+ in approxMB2 mb1 mb2 (m + unsafeShiftL n k) (e + unsafeShiftL f k) s+ _ + _ = Bottom+ (Approx mb1 m e s) * (Approx mb2 n f t)+ | am >= e && an >= f && a > 0 = approxMB2 mb1 mb2 (a+d) (ab+ac) u+ | am >= e && an >= f && a < 0 = approxMB2 mb1 mb2 (a-d) (ab+ac) u+ | am < e && n >= f = approxMB2 mb1 mb2 (a+b) (ac+d) u+ | am < e && -n >= f = approxMB2 mb1 mb2 (a-b) (ac+d) u+ | m >= e && an < f = approxMB2 mb1 mb2 (a+c) (ab+d) u+ | -m >= e && an < f = approxMB2 mb1 mb2 (a-c) (ab+d) u+ | a == 0 = approxMB2 mb1 mb2 (0) (ab+ac+d) u+ | am < e && an < f && a > 0 && ab > ac = approxMB2 mb1 mb2 (a+ac) (ab+d) u+ | am < e && an < f && a > 0 && ab <= ac = approxMB2 mb1 mb2 (a+ab) (ac+d) u+ | am < e && an < f && a < 0 && ab > ac = approxMB2 mb1 mb2 (a-ac) (ab+d) u+ | am < e && an < f && a < 0 && ab <= ac = approxMB2 mb1 mb2 (a-ab) (ac+d) u+ where am = (abs m)+ an = (abs n)+ a = m*n+ b = m*f+ c = n*e+ d = e*f+ ab = (abs b)+ ac = (abs c)+ u = s+t+ _ * _ = Bottom+ negate (Approx mb m e s) = Approx mb (-m) e s+ negate Bottom = Bottom+ abs (Approx mb m e s)+ | m' < e = let e' = m'+e+ in Approx mb e' e' (s-1)+ | otherwise = Approx mb m' e s+ where m' = abs m+ abs Bottom = Bottom+ signum (Approx _ m e _)+ | e == 0 = Approx 64 (signum m) 0 0+ | abs m < e = Approx 64 0 1 0+ | abs m == e = Approx 64 (signum m) 1 (-1)+ | otherwise = Approx 64 (signum m) 0 0+ signum Bottom = Approx 64 0 1 0+ fromInteger i = mapMB (max 64) $ approxAutoMB i 0 0++-- |Convert a rational number into an approximation of that number with+-- 'Precision' bits correct after the binary point.+toApprox :: Precision -> Rational -> Approx+toApprox t r = approxAutoMB m e (-t - 1)+ where+ m = (2 * r_scaled_rounded)+ r_scaled_rounded = round r_scaled+ r_scaled = r*2^^t+ e | r_scaled == fromIntegral r_scaled_rounded = 0+ | otherwise = 1++-- |Convert a rational number into an approximation of that number with+-- 'mBound' significant bits correct.+toApproxMB :: Int -> Rational -> Approx+toApproxMB mb r = + (Approx mb (numerator r) 0 0) / (Approx mb (denominator r) 0 0)++-- |Not a proper Fractional type as Approx are intervals.+instance Fractional Approx where+ fromRational = toApprox defaultPrecision+ recip = recipA+ a1 / a2 = a1 * (recipA $ setMB mb a2)+ where+ mb = mBound a1 `max` mBound a2++-- |Compute the reciprocal of an approximation. The number of bits after the+-- binary point is bounded by the 'midpoint bound' if the input is exact.+-- Otherwise, a good approximation with essentially the same significance as+-- the input will be computed.+recipA :: Approx -> Approx+recipA Bottom = Bottom+recipA (Approx mb m e s)+ | e == 0 && m /= 0+ = let s' = integerLog2 (abs m)+ in if abs m == bit s'+ then+ Approx mb (signum m) 0 (-s-s')+ else+ Approx mb+ (round (bit (mb+s') % m))+ 1+ (-mb-s-s')+ | (abs m) > e = let d = m*m-e*e+ d2 = unsafeShiftR d 1+ s' = integerLog2 d + 2 * errorBits+ in boundErrorTerm $ approxMB mb+ ((unsafeShiftL m s' + d2) `div` d)+ ((unsafeShiftL e s' + d2+1) `div` d + 1)+ (-s-s')+ -- (abs m) > e = let d = m*m-e*e+ -- s' = 2 * (integerLog2 m + errorBits)+ -- in boundErrorTerm $ Approx+ -- (round (unsafeShiftL m s'%(d)))+ -- (ceiling (1%2 + unsafeShiftL e s'%(d)))+ -- (-s-s')+ | otherwise = Bottom++-- |Divide an approximation by an integer.+divAInteger :: Approx -> Integer -> Approx+divAInteger Bottom _ = Bottom+divAInteger (Approx mb m e s) n =+ let t = integerLog2 n+ in approxMB mb+ (round (unsafeShiftL m t % n))+ (ceiling (unsafeShiftL e t % n))+ s++-- |Compute the modulus of two approximations.+modA :: Approx -> Approx -> Approx+modA (Approx mb1 m e s) (Approx mb2 n f t) =+ let r = min s t+ (d,m') = divMod (unsafeShiftL m (s-r)) (unsafeShiftL n (t-r))+ e' = scale e (s-r) + abs d * scale f (t-r)+ in approxMB2 mb1 mb2 m' e' r+modA _ _ = Bottom++-- |Compute the integer quotient (although returned as an 'Approx' since it+-- may be necessary to return 'Bottom' if the integer quotient can't be+-- determined) and the modulus as an approximation of two approximations.+divModA :: Approx -> Approx -> (Approx, Approx)+divModA (Approx mb1 m e s) (Approx mb2 n f t) =+ let r = min s t+ (d,m') = divMod (unsafeShiftL m (s-r)) (unsafeShiftL n (t-r))+ e' = e + abs d * f+ in (fromIntegral d, approxMB2 mb1 mb2 m' e' r)+divModA _ _ = (Bottom, Bottom)++-- |Not a proper Ord type as Approx are intervals.+instance Ord Approx where+ compare (Approx _ m e s) (Approx _ n f t)+ | e == 0 && f == 0 = compare (m:^s) (n:^t)+ | abs ((m:^s)-(n:^t)) > (e:^s)+(f:^t) = compare (m:^s) (n:^t)+ | otherwise = error "compare: comparisons are partial on Approx"+ compare _ _ = error "compare: comparisons are partial on Approx"++-- |The 'toRational' function is partial since there is no good rational+-- number to return for the trivial approximation 'Bottom'.+--+-- Note that converting to a rational number will give only a single rational+-- point. Do not expect to be able to recover the interval from this value.+instance Real Approx where+ toRational (Approx _ m e s) = approxRational+ (toRational (m:^s))+ (toRational (e:^s))+ toRational _ = undefined++-- |Convert the centre of an approximation into a 'Maybe' 'Double'.+toDoubleA :: Approx -> Maybe Double+toDoubleA = fmap (fromRational . toRational) . centre+++-- |Computes the precision of an approximation. This is roughly the number of+-- correct bits after the binary point.+precision :: Approx -> Extended Precision+precision (Approx _ _ 0 _) = PosInf+precision (Approx _ _ e s) = Finite $ - s - (integerLog2 e) - 1+precision Bottom = NegInf++-- |Computes the significance of an approximation. This is roughly the number+-- of significant bits.+significance :: Approx -> Extended Int+significance (Approx _ _ 0 _) = PosInf+significance (Approx _ 0 _ _) = NegInf+significance (Approx _ m 1 _) = Finite $ integerLog2 (abs m) - 1+significance (Approx _ m e _) =+ Finite $ (integerLog2 (abs m)) - (integerLog2 (e-1)) - 1+significance Bottom = NegInf++{-|+This function bounds the error term of an 'Approx'.++It is always the case that @x `'better'` boundErrorTerm x@.++Consider an approximation @Approx mb m e s@. If @e@ has /k/ bits then that+essentially expresses that the last /k/ bits of @m@ are unknown or garbage. By+scaling both @m@ and @e@ so that @e@ has a small number of bits we save on+memory space and computational effort to compute operations. On the other+hand, if we remove too many bits in this way, the shift in the mid-point of the+interval becomes noticable and it may adversely affect convergence speed of+computations. The number of bits allowed for @e@ after the operation is+determined by the constant 'errorBits'.++== Domain interpretation and verification++For this implementation to be correct it is required that this function is+below the identity on the domain /D/ of 'Approx'. For efficiency it is+desirable to be as close to the identity as possible.++This function will map a converging sequence to a converging sequence.+-}+boundErrorTerm :: Approx -> Approx+boundErrorTerm Bottom = Bottom+boundErrorTerm a@(Approx mb m e s)+ | e < errorBound = a+ | otherwise =+ let k = integerLog2 e + 1 - errorBits+ t = testBit m (k-1)+ m' = unsafeShiftR m k+ -- may overflow and use errorBits+1+ e' = unsafeShiftR (e + bit (k-1)) k + 1+ in if t+ then Approx mb (m'+1) e' (s+k)+ else Approx mb m' e' (s+k)++boundErrorTermMB :: Approx -> Approx+boundErrorTermMB Bottom = Bottom+boundErrorTermMB a@(Approx _ m e s)+ | e < errorBoundMB = a+ | otherwise =+ let k = integerLog2 e + 1 - errorBitsMB+ t = testBit m (k-1)+ m' = unsafeShiftR m k+ -- may overflow and use errorBits+1+ e' = unsafeShiftR (e + bit (k-1)) k + 1+ in if t+ then approxAutoMB (m'+1) e' (s+k)+ else approxAutoMB m' e' (s+k)++{-|+Limits the size of an approximation by restricting how much precision an+approximation can have.++It is always the case that @x `'better'` limitSize x@.++This is accomplished by restricting the exponent of the approximation from+below. In other words, we limit the precision possible.++It is conceivable to limit the significance of an approximation rather than+the precision. This would be an interesting research topic.++== Domain interpretation and verification++For this implementation to be correct it is required that this function is+below the identity on the domain /D/ of 'Approx'. For efficiency it is+desirable to be as close to the identity as possible.++This function will NOT map a converging sequence to a converging sequence for+a fixed precision argument. However, if the function is applied with+increasing precision for a converging sequence, then this will give a+converging sequence.+-}+limitSize :: Precision -> Approx -> Approx+limitSize _ Bottom = Bottom+limitSize l a@(Approx mb m e s)+ | k > 0 = Approx (mb-k)+ ((if testBit m (k-1) then (+1) else id) (unsafeShiftR m k))+ (1 + (unsafeShiftR (e + bit (k-1)) k))+ (-l)+ | otherwise = a+ where k = (-s)-l++-- |Throws an exception if the precision of an approximation is not larger+-- than the deafult minimum.+checkPrecisionLeft :: Approx -> Approx+checkPrecisionLeft a+ | precision a > pure defaultPrecision = a+ | otherwise = throw $ LossOfPrecision++-- |Bounds the error term and limits the precision of an approximation.+--+-- It is always the case that @x `'better'` limitAndBound x@.+limitAndBound :: Precision -> Approx -> Approx+limitAndBound limit =+ limitSize limit . boundErrorTerm++limitAndBoundMB :: Precision -> Approx -> Approx+limitAndBoundMB limit =+ limitSize limit . boundErrorTermMB++-- | Find the hull of two approximations.+unionA :: Approx -> Approx -> Approx+unionA Bottom _ = Bottom+unionA _ Bottom = Bottom+unionA a@(Approx mb1 _ _ _) b@(Approx mb2 _ _ _) =+ endToApprox (mb1 `max` mb2) (lowerBound a `min` lowerBound b) (upperBound a `max` upperBound b)++-- | Find the intersection of two approximations.+intersectionA :: Approx -> Approx -> Approx+intersectionA Bottom a = a+intersectionA a Bottom = a+intersectionA a@(Approx mb1 _ _ _) b@(Approx mb2 _ _ _) =+ if l <= u+ then endToApprox (mb1 `max` mb2) l u+ else error "Trying to take intersection of two non-overlapping Approx."+ where l = (lowerBound a `max` lowerBound b)+ u = (upperBound a `min` upperBound b)++-- | Determine if two approximations overlap.+consistentA :: Approx -> Approx -> Bool+consistentA Bottom _ = True+consistentA _ Bottom = True+consistentA a b = (lowerBound a `max` lowerBound b) <= (upperBound a `min` upperBound b)++-- |Given a list of polynom coefficients and a value this evaluates the+-- polynomial at that value.+--+-- This works correctly only for exact coefficients.+--+-- Should give a tighter bound on the result since we reduce the dependency+-- problem.+poly :: [Approx] -> Approx -> Approx+poly [] _ = 0+poly _ Bottom = Bottom+poly as x@(Approx mb _ _ _) =+ let --poly' :: [Dyadic] -> Dyadic -> Dyadic+ poly' as' x' = sum . zipWith (*) as' $ pow x'+ ms = map ((maybe (error "Can't compute poly with Bottom coefficients") id) . centre) as+ (Just c) = centre x+ (m':^s) = poly' ms c+ ds = zipWith (*) (tail as) (map fromIntegral ([1,2..] :: [Int]))+ (Finite b) = upperBound . abs $ poly' ds x+ (Finite (e':^_)) = fmap (b*) $ radius x+ -- exponent above will be same as s+ in approxMB mb m' e' s++-- |Gives a list of powers of a number, i.e., [1,x,x^2,...].+pow :: (Num a) => a -> [a]+pow x = iterate (* x) 1++-- |Computes lists of binomial coefficients. [[1], [1,1], [1,2,1], [1,3,3,1], ...]+binomialCoefficients :: (Num a) => [[a]]+binomialCoefficients =+ let f ss = 1 : zipWith (+) ss (tail ss) ++ [1]+ in iterate f [1]++-- |Computes powers of approximations. Should give tighter intervals than the+-- general 'pow' since take the dependency problem into account. However, so+-- far benchmarking seems to indicate that the cost is too high, but this may+-- depend on the application.+powers :: Approx -> [Approx]+powers (Approx mb m e s) =+ let ms = pow m+ es = pow e+ f = reverse . zipWith (*) ms . reverse . zipWith (*) es+ sumAlt [] = (0,0)+ sumAlt (x:[]) = (x,0)+ sumAlt (x:y:xs) = let (a,b) = sumAlt xs in (a+x,b+y)+ g s' (m', e') = approxMB mb m' e' s'+ in zipWith g (iterate (+s) 0) $ map (sumAlt . f) binomialCoefficients+powers _ = repeat Bottom++{-|+Old implementation of sqrt using Heron's method. Using the current method+below proved to be more than twice as fast for small arguments (~50 bits) and+many times faster for larger arguments.+-}+sqrtHeronA :: Precision -> Approx -> Approx+sqrtHeronA _ Bottom = Bottom+sqrtHeronA k a@(Approx mb m e s)+ | -m > e = error "Attempting sqrt of Approx containing only negative numbers."+ | m < e = Bottom+ | e == 0 = let (n:^t) = shiftD (-k) $ sqrtD (-k-2) (m:^s)+ in mapMB (max mb) $ approxAutoMB n 1 t+ | m == e = let (n:^t) = sqrtD (s `quot` 2 -errorBits) ((m+e):^s)+ n' = (n+2) `quot` 2+ in approxMB mb n' n' t+ | otherwise = let (Finite p) = significance a+ s' = s `quot` 2 - p - errorBits+ l@(n:^t) = sqrtD s' ((m-e):^s)+ (n':^t') = sqrtD' s' ((m+e):^s) l+ in endToApprox mb (Finite ((n-1):^t)) (Finite ((n'+1):^t'))++{-|+Compute the square root of an approximation.++This and many other operations on approximations is just a reimplementation of+interval arithmetic, with an extra argument limiting the effort put into the+computation. This is done via the precision argument.++The resulting approximation should approximate the image of every point in the+input approximation.+-}+sqrtA :: Approx -> Approx+sqrtA Bottom = Bottom+sqrtA x@(Approx _ 0 0 _) = x+sqrtA x@(Approx mb _ _ _) + | upperA x < 1 = + sqrtRecA k (recipA $ setMB k x)+ | otherwise =+ -- limitAndBoundMB mb $ + x * sqrtRecA k x+ where+ k = 2*mb + 2++{-|+This uses Newton's method for computing the reciprocal of the square root.+-}+sqrtRecA :: Precision -> Approx -> Approx+sqrtRecA _ Bottom = Bottom+sqrtRecA k a@(Approx mb m e s)+ | -m > e = error "Attempting sqrtRec of Approx containing only negative numbers."+ | m < e = Bottom+ | e == 0 = let (n:^t) = shiftD (-k) $ sqrtRecD (-k-2) (m:^s)+ in mapMB (max mb) $ approxAutoMB n 1 t+ | m == e = let (n:^t) = sqrtRecD (s `quot` 2 -errorBits) ((m+e):^s)+ n' = (n+2) `quot` 2+ in approxMB mb n' n' t+ | otherwise = let (Finite p) = significance a+ s' = s `quot` 2 - p - errorBits+ (n:^t) = sqrtRecD s' ((m-e):^s) -- upper bound of result+ -- We have tried to use sqrtRecD' with the above value as+ -- a first approximation to the result. However, the low+ -- endpoint may be too far away as a starting value to+ -- ensure convergence to the right root. It's possible+ -- that if we swap the order we would be fine. But as it+ -- is, this computes a new first approximation.+ (n':^t') = sqrtRecD s' ((m+e):^s) -- lower bound of result+ in endToApprox mb (Finite ((n'-1):^t')) (Finite ((n+1):^t))+++{-|+The starting values for newton iterations can be found using the auxiliary function findStartingValues below.++For example, to generate the starting values for sqrtRecD above using three leading bits for both odd and even exponents the following was used:++> findStartingValues ((1/) . sqrt) [1,1.125..2]+[Approx 4172150648 1 (-32),Approx 3945434766 1 (-32),Approx 3752147976 1 (-32),Approx 3584793264 1 (-32),Approx 3438037830 1 (-32),Approx 3307969824 1 (-32),Approx 3191645366 1 (-32),Approx 3086800564 1 (-32)]+> mapM_ (putStrLn . showInBaseA 2 . limitSize 6) it+0.111110~+0.111011~+0.111000~+0.110101~+0.110011~+0.110001~+0.110000~+0.101110~+> findStartingValues ((1/) . sqrt) [2,2.25..4]+[Approx 2950156016 1 (-32),Approx 2789843678 1 (-32),Approx 2653169278 1 (-32),Approx 2534831626 1 (-32),Approx 2431059864 1 (-32),Approx 2339087894 1 (-32),Approx 2256834080 1 (-32),Approx 2182697612 1 (-32)]+> mapM_ (putStrLn . showInBaseA 2 . limitSize 6) it+0.101100~+0.101010~+0.101000~+0.100110~+0.100100~+0.100011~+0.100010~+0.100001~+-}+findStartingValues :: (Double -> Double) -> [Double] -> [Approx]+findStartingValues f = map (fromRational . toRational . (/2)) . (\l -> zipWith (+) l (tail l)) . map f++-- |Square an approximation. Gives the exact image interval, as opposed to+-- multiplicating a number with itself which will give a slightly larger+-- interval due to the dependency problem.+sqrA :: Approx -> Approx+sqrA Bottom = Bottom+sqrA (Approx mb m e s)+ | am > e = approxMB mb (m^(2 :: Int) + e^(2 :: Int)) (2*am*e) (2*s)+ | otherwise = let m' = (am + e)^(2 :: Int) in approxMB mb m' m' (2*s-1)+ where am = abs m+-- Binary splitting++{-|+Binary splitting summation of linearly convergent series as described in+/'Fast multiprecision evaluation of series of rational numbers'/ by B Haible+and T Papanikolaou, ANTS-III Proceedings of the Third International Symposium+on Algorithmic Number Theory Pages 338-350, 1998.++The main idea is to balance the computations so that more operations are+performed with values of similar size. Using the underlying fast+multiplication algorithms this will give performance benefits.++The algorithm parallelises well. However, a final division is needed at the+end to compute /T\/BQ/ which amount to a substantial portion of the+computation time.+-}+abpq :: Num a => [Integer] -> [Integer] -> [a] -> [a] -> Int -> Int -> (a, a, Integer, a)+abpq as bs ps qs n1 n2+ | n == 1 = (ps !! n1, qs !! n1, bs !! n1, fromIntegral (as !! n1) * ps !! n1)+ | n < 6 = let as' = take n $ drop n1 as+ bs' = take n $ drop n1 bs+ ps' = take n $ drop n1 ps+ qs' = take n $ drop n1 qs+ pbs = product bs'+ bs'' = map (pbs `div`) bs'+ ps'' = scanl1 (*) ps'+ qs'' = scanr1 (*) (tail qs' ++ [1])+ in (ps'' !! (n-1), product qs', pbs+ , sum $ zipWith4 (\a b p q -> fromIntegral a * fromIntegral b * p * q)+ as' bs'' ps'' qs'')+ | n > 1 =+ let (pl, ql, bl, tl) = abpq as bs ps qs n1 m+ (pr, qr, br, tr) = abpq as bs ps qs m n2+ in (pl * pr, ql * qr, bl * br, fromIntegral br * qr * tl + fromIntegral bl * pl * tr)+ | otherwise = error "Non-expected case in binary splitting"+ where+ n = n2 - n1+ m = (n1 + n2 + 1) `div` 2++ones :: Num a => [a]+ones = repeat 1++{-|+Computes the list [lg 0!, lg 1!, lg 2!, ...].+-}+-- To be changed to Stirling formula if that is faster+log2Factorials :: [Int]+log2Factorials = map integerLog2 . scanl1 (*) $ 1:[1..]++-- Straighforward Taylor summation++{-|+Computes a sum of the form ∑ aₙ/qₙ where aₙ are approximations and qₙ are+integers. Terms are added as long as they are larger than the current+precision bound. The sum is adjusted for the tail of the series. For this to+be correct we need the the terms to converge geometrically to 0 by a factor of+at least 2.+-}+taylor :: Precision -> [Approx] -> [Integer] -> Approx+taylor res as qs =+ let res' = res + errorBits+ f a q = limitAndBound res' $ a * recipA (setMB (mBound a) $ fromIntegral q)+ mb = zipWith f as qs+ (cs,(d:_)) = span nonZeroCentredA mb -- This span and the sum on the next line do probably not fuse!+ in fudge (sum cs) d++-- | A list of factorial values.+fac :: Num a => [a]+fac = map fromInteger $ 1 : scanl1 (*) [1..]++-- | A list of the factorial values of odd numbers.+oddFac :: Num a => [a]+oddFac = let f (_:x:xs) = x:f xs+ f _ = error "Impossible"+ in f fac++{-+evenFac :: Num a => [a]+evenFac = let f (x:_:xs) = x:f xs+ f _ = error "Impossible"+ in f fac+-}++-- | Checks if the centre of an approximation is not 0.+nonZeroCentredA :: Approx -> Bool+nonZeroCentredA Bottom = False+nonZeroCentredA (Approx _ 0 _ _) = False+nonZeroCentredA _ = True++-- This version is faster especially far smaller precision.++{-|+Computes the sum of the form ∑ aₙxⁿ where aₙ and x are approximations.++Terms are added as long as they are larger than the current precision bound.+The sum is adjusted for the tail of the series. For this to be correct we need+the the terms to converge geometrically to 0 by a factor of at least 2.+-}+taylorA :: Precision -> [Approx] -> Approx -> Approx+taylorA res as x =+ fudge sm d+ where+ (sm, d) = sumAndNext . takeWhile (nonZeroCentredA . fst) . addNext . map (limitAndBound res) $ zipWith (*) as (pow x)+ sumAndNext = aux 0+ where+ aux a [(b,dd)] = (a+b,dd)+ aux a ((b,_):rest) = aux (a+b) rest+ aux _ _ = undefined+ addNext (x1:x2:xs) = (x1,x2):(addNext (x2:xs))+ addNext _ = error "taylorA: end of initite sequence" ++{- Exponential computed by standard Taylor expansion after range reduction.+-}++{-|+The exponential of an approximation. There are three implementation using+Taylor expansion here. This is just choosing between them.++More thorough benchmarking would be desirable.++Is faster for small approximations < ~2000 bits.+-}+expA :: Approx -> Approx+expA = expTaylorA'++-- | Exponential by binary splitting summation of Taylor series.+expBinarySplittingA :: Precision -> Approx -> Approx+expBinarySplittingA _ Bottom = Bottom+expBinarySplittingA res a@(Approx mb m e s) =+ let s' = s + integerLog2 m+ -- r' chosen so that a' below is smaller than 1/2+ r' = floor . sqrt . fromIntegral . max 5 $ res+ r = s' + r'+ -- a' is a scaled by 2^k so that 2^(-r') <= a' < 2^(-r'+1)+ a' = Approx mb m e (s-r)+ mb' = mb+res+ -- (Finite c) = min (significance a) (Finite res)+ (Just n) = findIndex (>= res+r) $ zipWith (+) log2Factorials [0,r'..]+ (p, q, b, t) = abpq ones+ ones+ (map (setMB mb') $ 1:repeat a')+ (map (setMB mb') $ 1:[1..])+ 0+ n+ nextTerm = a * p * recipA (fromIntegral n * q)+ ss = iterate (boundErrorTerm . sqrA) $ fudge (t * recipA (fromIntegral b*q)) nextTerm+ in ss !! r++-- | Exponential by summation of Taylor series.+expTaylorA :: Precision -> Approx -> Approx+expTaylorA _ Bottom = Bottom+expTaylorA res (Approx mb m e s) =+ let s' = s + integerLog2 m+ -- r' chosen so that a' below is smaller than 1/2+ r' = floor . sqrt . fromIntegral . max 5 $ res+ r = max 0 $ s' + r'+ -- a' is a scaled by 2^k so that 2^(-r') <= a' < 2^(-r'+1)+ a' = (Approx (mb + res) m e (s-r))+ t = taylor+ (res + r)+ (iterate (a'*) 1)+ (scanl1 (*) $ 1:[1..])+ in (!! r) . iterate (boundErrorTermMB . sqrA) $ t+ +-- | Exponential by summation of Taylor series.+expTaylorA' :: Approx -> Approx+expTaylorA' Bottom = Bottom+expTaylorA' a + | upperA a < 0 = recipA $ aux (-a)+ | otherwise = aux a+ where+ aux Bottom = Bottom+ aux (Approx mb 0 0 _) = Approx mb 1 0 0+ aux (Approx mb m 0 s) =+ let s' = s + (integerLog2 m)+ -- r' chosen so that a' below is smaller than 1/2+ r' = floor . sqrt . fromIntegral . max 5 $ mb+ r = max 0 $ s' + r'+ mb'_ = mb + r + (integerLog2 m) + 1+ mb' = (120*mb'_) `div` 100+ -- a' is a scaled by 2^k so that 2^(-r') <= a' < 2^(-r'+1)+ a' = (Approx mb' m 0 (s-r))+ t = boundErrorTermMB $ taylorA mb' (map (recipA . setMB mb') fac) a'+ in (!! r) . iterate (boundErrorTermMB . sqrA) $ t+ aux a2 = aux (lowerA a2) `unionA` aux (upperA a2)+ +{- Logarithms computed by ln x = 2*atanh ((x-1)/(x+1)) after range reduction.+-}++{-|++Computing the logarithm of an approximation. This chooses the fastest implementation.++More thorough benchmarking is desirable.++Binary splitting is faster than Taylor. AGM should be used over ~1000 bits.+-}+logA :: Approx -> Approx+-- This implementation asks for the dyadic approximation of the endpoints, we+-- should instead use that, after the first range reduction, the derivative is+-- less than 3/2 on the interval, so it easy to just compute one expensive+-- computation. We could even make use of the fact that the derivative on the+-- interval x is bounded by 1/x to get a tighter bound on the error.+logA Bottom = Bottom+logA x@(Approx _ m e _)+ | m > e && upperA x < 1 = -(logInternal (recipA x))+ | m > e = logInternal x+-- let (n :^ t) = logD (negate p) $ (m-e) :^ s+-- (n' :^ t') = logD (negate p) $ (m+e) :^ s+-- in endToApprox (Finite ((n-1):^t)) (Finite ((n'+1):^t'))+ | otherwise = Bottom++logInternal :: Approx -> Approx+logInternal Bottom = error "LogInternal: impossible"+logInternal (Approx mb m e s) =+ let t' = (negate mb) - 10 - max 0 (integerLog2 m + s) -- (5 + size of argument) guard digits+ r = s + integerLog2 (3*m) - 1+ x = scale (m :^ s) (-r) -- 2/3 <= x' <= 4/3+ y = divD' t' (x - 1) (x + 1) -- so |y| <= 1/5+ (n :^ s') = flip scale 1 $ atanhD t' y+ (e' :^ s'') = divD' t' (e:^(s-r)) x -- Estimate error term.+ res = approxMB mb n (scale (e' + 1) (s'' - s')) s'+ in boundErrorTerm $ res + fromIntegral r * log2A t'++-- | Logarithm by binary splitting summation of Taylor series.+logBinarySplittingA :: Precision -> Approx -> Approx+logBinarySplittingA _ Bottom = Bottom+logBinarySplittingA res a@(Approx mb m e s) =+ if m <= e then Bottom -- only defined for strictly positive arguments+ else+ let r = s + integerLog2 (3*m) - 1+ a' = Approx (mb+res) m e (s-r) -- a' is a scaled by a power of 2 so that 2/3 <= |a'| <= 4/3+ u = a' - 1+ v = a' + 1+ u2 = sqrA u+ v2 = sqrA v+ Finite res' = min (significance a) (Finite res)+ n = ceiling . (/2) $ fromIntegral (-res')/(log 0.2/log 2) - 1+ (_, q, b, t) = abpq (repeat 2)+ [1,3..]+ (u:repeat u2)+ (v:repeat v2)+ 0+ n+ nextTerm = recipA (setMB (mb+res) 5) ^^ (2*n+1)+ in boundErrorTerm $ fudge (t * recipA (fromIntegral b*q) + fromIntegral r * log2A (-res)) nextTerm++-- | Logarithm by summation of Taylor series.+logTaylorA :: Precision -> Approx -> Approx+logTaylorA _ Bottom = Bottom+logTaylorA res (Approx mb m e s) =+ if m <= e then Bottom -- only defined for strictly positive arguments+ else+ let res' = res + errorBits+ r = s + integerLog2 (3*m) - 1+ a' = approxMB mb m e (s-r) -- a' is a scaled by a power of 2 so that 2/3 <= a' <= 4/3+ u = a' - 1+ v = a' + 1+ x = u * recipA v -- so |u/v| <= 1/5+ x2 = boundErrorTerm $ sqrA x+ t = taylor+ res'+ (iterate (x2*) x)+ [1,3..]+ in boundErrorTerm $ 2 * t + fromIntegral r * log2A (-res')++-- Sine computed by Taylor expansion after 2 stage range reduction.++-- | Computes sine by summation of Taylor series after two levels of range reductions.+sinTaylorA :: Approx -> Approx+sinTaylorA Bottom = Approx 64 0 1 0 -- [-1,1]+sinTaylorA a@(Approx mb _ e _) + | e == 0 = sinTaylorRed2A aRed+ | otherwise = sL `unionA` sR -- aRed is in the interval [-π/2,π/2] where sine is monotone+ where+ (aRed, (maRedL, maRedR)) = sinTaylorRed1A a+ sL =+ case maRedL of+ Nothing -> Approx mb (-1) 0 0 -- aRed probably contains -pi/2+ Just aRedL -> sinTaylorRed2A aRedL+ sR =+ case maRedR of+ Nothing -> Approx mb 1 0 0 -- aRed probably contains +pi/2+ Just aRedR -> sinTaylorRed2A aRedR++-- | First level of range reduction for sine. Brings it into the interval [-π/2,π/2].+sinTaylorRed1A :: Approx -> (Approx, (Maybe Approx, Maybe Approx))+sinTaylorRed1A Bottom = (Bottom, (Nothing, Nothing))+sinTaylorRed1A a@(Approx mb _ _ _) = + let _pi = piA (mb+10)+ _halfPi = scale _pi (-1)+ x = setMB mb . (subtract _halfPi) . abs . (_pi -) . abs . (subtract _halfPi) . modA a $ 2*_pi+ xL = lowerA x+ xR = upperA x+ _halfPiL = lowerA _halfPi+ in (x, + (if (- _halfPiL) <= xL then Just xL else Nothing, -- guarantee -π/2 <= xL+ if xR <= _halfPiL then Just xR else Nothing)) -- guarantee xR <= π/2+ +-- | Second level of range reduction for sine.+sinTaylorRed2A :: Approx -> Approx+sinTaylorRed2A Bottom = Approx 64 0 1 0 -- [-1,1]+sinTaylorRed2A a@(Approx mb m _ s) = + let k = max 0 (integerLog2 m + s + (floor . sqrt $ fromIntegral mb))+ a' = a * recipA ((setMB mb 3)^k)+ a2 = negate $ sqrA a'+ t = taylorA mb (map (recipA . setMB mb) oddFac) a2+ step x = boundErrorTerm $ x * (3 - 4 * sqrA x)+ in (!! k) . iterate step . boundErrorTerm $ t * a'++-- | Computes the sine of an approximation. Chooses the best implementation.+sinA :: Approx -> Approx+sinA = sinTaylorA++-- | Computes the cosine of an approximation. Chooses the best implementation.+cosA :: Approx -> Approx+cosA Bottom = Approx 64 0 1 0 -- [-1,1]+cosA x@(Approx mb _ _ _) = sinA ((Approx 1 1 0 (-1)) * piA (mb+2) - x)++-- | Computes the sine of an approximation by binary splitting summation of Taylor series.+--+-- Begins by reducing the interval to [0,π/4].+sinBinarySplittingA :: Precision -> Approx -> Approx+sinBinarySplittingA _ Bottom = Bottom+sinBinarySplittingA res a =+ let _pi = piBorweinA res+ (Approx mb' m' e' s') = 4 * a * recipA _pi+ (k,m1) = m' `divMod` bit (-s')+ a2 = _pi * fromDyadicMB mb' (1:^(-2)) * (Approx mb' m1 e' s')+ in case k `mod` 8 of+ 0 -> sinInRangeA res a2+ 1 -> cosInRangeA res (_pi * fromDyadicMB mb' (1:^(-2)) - a2)+ 2 -> cosInRangeA res a2+ 3 -> sinInRangeA res (_pi * fromDyadicMB mb' (1:^(-2)) - a2)+ 4 -> - sinInRangeA res a2+ 5 -> - cosInRangeA res (_pi * fromDyadicMB mb' (1:^(-2)) - a2)+ 6 -> - cosInRangeA res a2+ 7 -> - sinInRangeA res (_pi * fromDyadicMB mb' (1:^(-2)) - a2)+ _ -> error "Impossible"++-- | Computes the cosine of an approximation by binary splitting summation of Taylor series.+--+-- Begins by reducing the interval to [0,π/4].+cosBinarySplittingA :: Precision -> Approx -> Approx+cosBinarySplittingA _ Bottom = Bottom+cosBinarySplittingA res a =+ let _pi = piBorweinA res+ (Approx mb' m' e' s') = 4 * a * recipA _pi+ (k,m1) = m' `divMod` bit (-s')+ a2 = _pi * fromDyadicMB mb' (1:^(-2)) * (Approx mb' m1 e' s')+ in case k `mod` 8 of+ 0 -> cosInRangeA res a2+ 1 -> sinInRangeA res (_pi * fromDyadicMB mb' (1:^(-2)) - a2)+ 2 -> - sinInRangeA res a2+ 3 -> - cosInRangeA res (_pi * fromDyadicMB mb' (1:^(-2)) - a2)+ 4 -> - cosInRangeA res a2+ 5 -> - sinInRangeA res (_pi * fromDyadicMB mb' (1:^(-2)) - a2)+ 6 -> sinInRangeA res a2+ 7 -> cosInRangeA res (_pi * fromDyadicMB mb' (1:^(-2)) - a2)+ _ -> error "Impossible"+++-- | Computes the arc tangent of an approximation. Chooses the best implementation.+atanA :: Precision -> Approx -> Approx+atanA = atanTaylorA++-- | Computes the arc tangent of an approximation by binary splitting summation of Taylor series.+atanBinarySplittingA :: Precision -> Approx -> Approx+atanBinarySplittingA _ Bottom = Bottom+atanBinarySplittingA res a =+ let rr x = x * recipA (1 + sqrtA (1 + sqrA x))+ a' = rr . rr . rr $ a -- range reduction so that |a'| < 1/4+ a2 = - sqrA a'+ res' = case (significance a) of+ (Finite _r) -> min res _r+ _ -> res+-- Finite res' = min (significance a) (Finite res)+ n = (res' + 1) `div` 2+ (_, q, b, t) = abpq ones+ [1,3..]+ (a':repeat a2)+ (repeat 1)+ 0+ n+ nextTerm = Approx (mBound a + res) 1 0 (-2*n)+ in boundErrorTerm . (8*) $ fudge (t * recipA (fromIntegral b*q)) nextTerm++-- + Bottom+-- + Deal with sign -- Because of next line, not worthwhile+-- + if lowerbound(abs a) > 2 then pi/2 - atan (1/a) -- Don't want to do this, what if 0 \in a?+-- + else+-- - r = min res (significance a)+-- - k = floor (sqrt r) / 2 `min` 2 (guarantee |x| < 1/2)+-- - x = rr^k(a)+-- - taylor (r + k + 5) (-x^2) [1,3..]+-- - (x*)+-- - same error as x+-- - (2^k *)++atanTaylorA :: Precision -> Approx -> Approx+atanTaylorA _ Bottom = Bottom+atanTaylorA res a@(Approx mb _ _ _) =+ let (Finite r) = min (pure res) (significance a)+ k = min (floor (sqrt (fromIntegral r)) `div` 2) 2+ res' = (mb `max` res) + k + 5+ rr _x = _x * recipA (1 + sqrtA (1 + sqrA _x))+ x = boundErrorTerm $ iterate rr (setMB res' a) !! k+ x2 = negate (sqrA x)+ t = boundErrorTerm $ x * taylorA res' (map (recipA . setMB res') [1,3..]) x2+ in scale t k++-- > let x = fromDouble (-0.2939788524332769)+-- > require 10 $ x+-- Approx (-5295852201093248) 1 (-54)+-- > require 10 . tan $ atan x+-- Approx (-10845905307838971904) 907 (-65)+-- > scale (-5295852201093248) 11+-- -10845905307838971904+--+-- problemet är att 1 måste skalas till 2^11, men blev bara 907+--+-- Men problemet verkar vara i tan, inte i atan.+++-- | Computes the arc tangent of an approximation by summation of Taylor series.+-- atanTaylorA :: Precision -> Approx -> Approx+-- atanTaylorA _ Bottom = Bottom+-- atanTaylorA res a =+-- let rr x = x * recipA res (1 + sqrtA res (1 + sqrA x))+-- a' = rr . rr . rr $ a -- range reduction so that |a'| < 1/4+-- a2 = - sqrA a'+-- res' = case (significance a) of+-- (Finite _r) -> min res _r+-- _ -> res+-- -- Finite res' = min (significance a) (Finite res)+-- t = taylorA res' (map (recipA res') [1,3..]) a2+-- in boundErrorTerm . (8*) $ t++{-+swapSinCos :: Precision -> Approx -> Approx+swapSinCos res a = sqrtA res $ 1 - sqrA a+-}++-- Computes sine if second argument is in the range [0,pi/4]+sinInRangeA :: Precision -> Approx -> Approx+sinInRangeA _ Bottom = Bottom+sinInRangeA res a =+ let n = res `div` 2 -- need to improve this estimate (is valid from res>=80)+ (_, q, b, t) = abpq ones+ ones+ (a:repeat (- sqrA a))+ (1:[2*i*(2*i+1) | i <- [1..]] :: [Approx])+ 0+ n+ nextTerm = fromDyadicMB (mBound a + res) (1:^(-res))+ in boundErrorTerm $ fudge (t * recipA (fromIntegral b*q)) nextTerm++-- Computes cosine if second argument is in the range [0,pi/4]+cosInRangeA :: Precision -> Approx -> Approx+cosInRangeA _ Bottom = Bottom+cosInRangeA res a =+ let n = res `div` 2 -- need to improve this estimate (is valid from res>=80)+ (_, q, b, t) = abpq ones+ ones+ (1:repeat (- sqrA a))+ (1:[2*i*(2*i-1) | i <- [1..]] :: [Approx])+ 0+ n+ nextTerm = fromDyadicMB (mBound a + res) (1:^(-res))+ in boundErrorTerm $ fudge (t * recipA (fromIntegral b*q)) nextTerm++{-|+Computes a sequence of approximations of π using binary splitting summation of+Ramanujan's series. See Haible and Papanikolaou 1998.+-}+piRaw :: [Approx]+piRaw = unfoldr f (1, (1, 1, 1, 13591409))+ where as = [13591409,13591409+545140134..]+ bs = ones+ ps = (1:[-(6*n-5)*(2*n-1)*(6*n-1) | n <- [1,2..]])+ qs = (1:[n^3*640320^2*26680 | n <- [1,2..]])+ f (i, (pl, ql, bl, tl)) = + let i2 = i*2+ (pr, qr, br, tr) = abpq as bs ps qs i i2+ n = 21+47*(i-1)+ x = fromIntegral tl * recipA (setMB n $ fromIntegral (bl*ql))+ x1 = fudge x $ fromDyadicMB n (1:^(-n))+ x2 = boundErrorTermMB $ sqrtA (setMB n 1823176476672000) * recipA x1+ in Just ( x2+ , (i2, (pl * pr, ql * qr, bl * br, fromIntegral br * qr * tl + fromIntegral bl * pl * tr))+ )++-- | Computes an 'Approx' of π of the given precision.+piA :: Precision -> Approx+piA res = limitAndBound res . head $ dropWhile ((< pure res) . precision) piRaw++{-|+Second argument is noise to be added to first argument. Used to allow for the+error term when truncating a series.+-}+fudge :: Approx -> Approx -> Approx+fudge a (Approx _ 0 0 _) = a+fudge (Approx mb m 0 s) (Approx mb' m' e' s') =+ approxMB2 mb mb' (m `shift` (s - s')) (abs m' + e' + 1) s'+fudge (Approx mb m e s) (Approx mb' m' e' s') =+ let m'' = 1 + (abs m' + e') `shift` (s' - s + 1)+ in approxMB2 mb mb' m (e+m'') s+fudge _ _ = Bottom++--++-- | Compute π using Machin's formula. Lifted from computation on dyadic numbers.+piMachinA :: Precision -> Approx+piMachinA t = let (m:^s) = piMachinD (-t) in approxAutoMB m 1 s++-- | Compute π using AGM as described in Borwein and Borwein's book 'Pi and+-- the AGM'. Lifted from computation on dyadic numbers.+piBorweinA :: Precision -> Approx+piBorweinA t = let (m:^s) = piBorweinD (-t) in approxAutoMB m 1 s+++-- | Compute π using AGM as described in Borwein and Borwein's book 'Pi and+-- the AGM'.+piAgmA :: Precision -> Approx -> Approx+piAgmA t x_@(Approx mb_ _ _ _) = + let -- t' = t - 10+ mb = mb_ + t+ a = setMB mb 1+ x = setMB mb x_+ b = boundErrorTerm $ (2*x*recipA (x^2-1))^2+ ss = agmA t a b+ c = boundErrorTerm . (1-) . (*recipA (1-b^2)) . agm2 . agm1 $ ss+ d = sqrtA (1+b)+ b2 = b^2+ b3 = b2*b+ b4 = b2^2+ l = boundErrorTerm $ + (((Approx mb 1 0 (-1))*b-(Approx mb 3 0 (-4))*b2+(Approx mb 9 0 (-5))*b3)*c*d-1/(1+b)+(2+b2)/d) + / ((2+(Approx mb 1 0 (-1))*b2+(Approx mb 9 0 (-5))*b4)*c+b2)+ u = boundErrorTerm $ + ((Approx mb 1 0 (-1))*b*c*d-1/(1+b)+(2+b2+(Approx mb 3 0 (-3))*b3+(Approx mb 9 0 (-3))*b4)/d) + / ((2+(Approx mb 1 0 (-1))*b2)*c+b2+(Approx mb 9 0 (-3))*b4)+ r = boundErrorTerm $ unionA l u+ e = boundErrorTerm $ unionA+ ((2+(Approx mb 1 0 (-1))*b2)*r-(Approx mb 1 0 (-1))*b*d)+ ((2+(Approx mb 1 0 (-1))*b2+(Approx mb 9 0 (-5))*b4)*r+ -((Approx mb 1 0 (-1))*b-(Approx mb 3 0 (-4))*b2+(Approx mb 9 0 (-5))*b3)*d)+ _pi = boundErrorTerm $ unionA (2*(snd (last ss))*e) (2*(fst (last ss))*e)+ in _pi+piAgmA _ _ = Bottom++-- | Compute approximations of ln 2. Lifted from computation on dyadic numbers.+log2A :: Precision -> Approx+log2A t = let (m:^s) = ln2D t in approxAutoMB m 1 s+++-- AGM++-- | Compute logarithms using AGM as described in Borwein and Borwein's book 'Pi and+-- the AGM'. An approximation of pi is produced as a by-product.+lnSuperSizeUnknownPi :: Precision -> Approx -> (Approx,Approx)+lnSuperSizeUnknownPi t x_@(Approx mb_ _ _ _) =+ let --t' = t - 10+ mb = mb_ + t+ a = setMB mb 1+ x = setMB mb x_+ b = boundErrorTerm $ (2*x*recipA (x^2-1))^2+ ss = agmA t a b+ (an,bn) = last ss+ c = boundErrorTerm . (1-) . (*recipA (1-b^2)) . agm2 . agm1 $ ss+ d = sqrtA (1+b)+ b2 = b^2+ b3 = b2*b+ b4 = b2^2+ l = boundErrorTerm $+ (((Approx mb 1 0 (-1))*b-(Approx mb 3 0 (-4))*b2+(Approx mb 9 0 (-5))*b3)*c*d-1/(1+b)+(2+b2)/d)+ / ((2+(Approx mb 1 0 (-1))*b2+(Approx mb 9 0 (-5))*b4)*c+b2)+ u = boundErrorTerm $+ ((Approx mb 1 0 (-1))*b*c*d-1/(1+b)+(2+b2+(Approx mb 3 0 (-3))*b3+(Approx mb 9 0 (-3))*b4)/d)+ / ((2+(Approx mb 1 0 (-1))*b2)*c+b2+(Approx mb 9 0 (-3))*b4)+ r = boundErrorTerm $ unionA l u+ e = boundErrorTerm $ unionA+ ((2+(Approx mb 1 0 (-1))*b2)*r-(Approx mb 1 0 (-1))*b*d)+ ((2+(Approx mb 1 0 (-1))*b2+(Approx mb 9 0 (-5))*b4)*r+ -((Approx mb 1 0 (-1))*b-(Approx mb 3 0 (-4))*b2+(Approx mb 9 0 (-5))*b3)*d)+ _pi = boundErrorTerm $ unionA (2*bn*e) (2*an*e)+ in (r,_pi) --[a,b,c,d,b2,b3,b4,l,u,r,e,pi]+lnSuperSizeUnknownPi _ Bottom = (Bottom, Bottom)++-- | Compute logarithms using AGM as described in Borwein and Borwein's book 'Pi and+-- the AGM'. An approximation of pi is needed as an extra argument.+lnSuperSizeKnownPi :: Precision -> Approx -> Approx -> Approx+lnSuperSizeKnownPi _t _pi Bottom = Bottom+lnSuperSizeKnownPi t _pi x_@(Approx mb_ _ _ _) =+ let --t' = t - 10+ mb = mb_ + t+ a = setMB mb 1+ x = setMB mb x_+ b = boundErrorTerm $ (2*x*recipA (x^2-1))^2+ b2 = b^2+ b3 = b2*b+ b4 = b2^2+ b1sqrt = sqrtA (1+b)+ step (_a,_b) = (boundErrorTerm $ Approx mb 1 0 (-1) * (_a+_b)+ ,boundErrorTerm $ sqrtA (_a*_b))+ close (_a,_b) = approximatedBy 0 $ _a-_b+ ((an,bn):_) = dropWhile (not . close) $ iterate step (a,b)+ i = boundErrorTerm $ unionA (_pi*recipA (2*an)) (_pi*recipA (2*bn))+ l = (i + ((Approx mb 1 0 (-1))*b-(Approx mb 3 0 (-4))*b2+(Approx mb 9 0 (-5))*b3)*b1sqrt)+ / (2 + (Approx mb 1 0 (-1))*b2 + (Approx mb 9 0 (-5))*b4)+ u = (i + (Approx mb 1 0 (-1))*b*b1sqrt) / (2 + (Approx mb 1 0 (-1))*b2)+ in boundErrorTerm $ unionA l u++lnLarge :: Precision -> Approx -> Approx+lnLarge _t Bottom = Bottom+lnLarge t x_@(Approx mb_ _ _ _) =+ let (Finite k) = min (significance x) (Finite (-t))+ x = setMB (mb_ + t) x_+ _pi = piBorweinA t+ iL2 = integerLog2+ fI = fromIntegral+ n = max 0 . (1+) . (+(iL2 (fI k)-2)) . negate . iL2 . fI . iL2 . truncate $ toRational x+ (Approx mb2 m e s) = lnSuperSizeKnownPi t _pi $ x^(2^n)+ in Approx mb2 m e (s-n)++lnSmall :: Precision -> Approx -> Approx+lnSmall _ Bottom = Bottom+lnSmall t x_@(Approx mb_ m _ s) =+ let (Finite t') = min (significance x) (Finite (-t))+ x = setMB (mb_ + t) x_+ _pi = piBorweinA t'+ iL2 = integerLog2+ -- fI = fromIntegral+ k = (-t) `div` 4 - iL2 m - s+ logx2k = lnSuperSizeKnownPi (-t') _pi $ x * 2^k+ log2k = lnSuperSizeKnownPi (-t') _pi $ 2^k+ in logx2k - log2k++-- | Compute logarithms using AGM as described in Borwein and Borwein's book 'Pi and+-- the AGM'. +-- TODO: adapt to mBound+logAgmA :: Precision -> Approx -> Approx+logAgmA t x+ | significance x < pure 5 = Bottom+ | 0 `approximatedBy` x = Bottom+ | signum x == (-1) = error "Trying to take logarithm of purely negative Approx."+ | lowerBound x > pure 2 = lnLarge t x+ | upperBound x < pure 3 = lnSmall t x+ | otherwise = error "Logic fault in logAgmA."+++agmA :: Precision -> Approx -> Approx -> [(Approx,Approx)]+agmA t a b = let --t' = t - 5+ step (_a,_b) = (boundErrorTerm $ Approx t 1 0 (-1) * (a+b), boundErrorTerm $ sqrtA (_a*_b))+ close (_a,_b) = approximatedBy 0 $ _a-_b+ in (\(as, bs) -> as ++ take 1 bs) . break close $ iterate step (a,b)++sqDiff :: Approx -> Approx -> Approx+sqDiff a b = boundErrorTerm $ a^2 - b^2++agm1 :: [(Approx, Approx)] -> [Approx]+agm1 = zipWith (*) [Approx 4 1 0 i | i <- [-1,0..]] . map (uncurry sqDiff)++agm2 :: [Approx] -> Approx+agm2 xs = sum (init xs) + unionA 0 (2 * last xs)++-- | Compute logarithms using AGM as described in Borwein and Borwein's book 'Pi and+-- the AGM'.+agmLn :: Precision -> Approx -> Approx+agmLn t x_@(Approx mb_ _ _ _) = + let --t' = t - 10+ mb = mb_ + t+ a = setMB mb 1+ x = setMB mb x_+ b = boundErrorTerm $ (2*x*recipA (x^2-1))^2+ ss = agmA t a b+ -- (an,bn) = last ss+ c = boundErrorTerm . (1-) . (*recipA (1-b^2)) . agm2 . agm1 $ ss+ d = sqrtA (1+b)+ b2 = b^2+ b3 = b2*b+ b4 = b2^2+-- l = boundErrorTerm $ (((Approx 1 0 (-1))*b-(Approx 3 0 (-4))*b2+(Approx 9 0 (-5))*b3)*c*d-recipA t' (1+b)+(2+b2)*recipA t' d) * recipA t' ((2+(Approx 1 0 (-1))*b2+(Approx 9 0 (-5))*b4)*c+b2)+-- u = boundErrorTerm $ ((Approx 1 0 (-1))*b*c*d-recipA t' (1+b)+(2+b2+(Approx 3 0 (-3))*b3+(Approx 9 0 (-3))*b4)*recipA t' d) *recipA t' ((2+(Approx 1 0 (-1))*b2)*c+b2+(Approx 9 0 (-3))*b4)+ l = boundErrorTerm $ (((Approx mb 1 0 (-1))*b-(Approx mb 3 0 (-4))*b2+(Approx mb 9 0 (-5))*b3)*c*d-1/(1+b)+(2+b2)/d) / ((2+(Approx mb 1 0 (-1))*b2+(Approx mb 9 0 (-5))*b4)*c+b2)+ u = boundErrorTerm $ ((Approx mb 1 0 (-1))*b*c*d-1/(1+b)+(2+b2+(Approx mb 3 0 (-3))*b3+(Approx mb 9 0 (-3))*b4)/d) / ((2+(Approx mb 1 0 (-1))*b2)*c+b2+(Approx mb 9 0 (-3))*b4)+ r = boundErrorTerm $ unionA l u+ e = boundErrorTerm $ unionA+ ((2+(Approx mb 1 0 (-1))*b2)*r-(Approx mb 1 0 (-1))*b*d)+ ((2+(Approx mb 1 0 (-1))*b2+(Approx mb 9 0 (-5))*b4)*r-((Approx mb 1 0 (-1))*b-(Approx mb 3 0 (-4))*b2+(Approx mb 9 0 (-5))*b3)*d)+ _pi = boundErrorTerm $ unionA (2*(snd (last ss))*e) (2*(fst (last ss))*e)+ in r --[a,b,c,d,b2,b3,b4,l,u,r,e,_pi]+agmLn _t _ = Bottom+ +-- The CR implementation++type Resources = Int++startLimit :: Int+startLimit = 80++bumpLimit :: Int -> Int+bumpLimit n = n * 3 `div` 2++resources :: ZipList Resources+resources = ZipList $ iterate bumpLimit startLimit++-- Should not use show as it would be impossible to write a corresponding read instance.+-- instance Show CR where+-- show = show . require 40++op2withResource ::+ (Approx -> Approx -> Approx) -> + (Approx -> Approx) -> + (Approx -> Approx -> Precision -> Approx)+op2withResource op2 post a b l =+ post $ limitAndBound l (op2 a (enforceMB $ mapMB (max l) b))++op1withResource ::+ (Approx -> Approx) -> + (Approx -> Approx) -> + (Approx -> Precision -> Approx)+op1withResource op1 post a l =+ post $ (op1 (enforceMB $ mapMB (max l) a))++instance Num CR where+ (CR x) + (CR y) = CR $ op2withResource (+) (ok 10) <$> x <*> y <*> resources+ (CR x) * (CR y) = CR $ op2withResource (*) (ok 10) <$> x <*> y <*> resources+ negate (CR x) = CR $ negate <$> x+ abs (CR x) = CR $ abs <$> x+ signum (CR x) = CR $ signum <$> x+ fromInteger n = + CR $ (\ a l -> ok 10 $ enforceMB $ setMB l a) <$> pure (fromInteger n) <*> resources++instance Fractional CR where+ recip (CR x) = CR $ op1withResource recipA id <$> x <*> resources+ fromRational x = CR $ toApprox <$> resources <*> pure x++instance Eq CR where+ (==) = error "CR does not have a total equality."++instance Ord CR where+ compare = error "CR does not have a total ordering."++instance Real CR where+ toRational = toRational . require 40++-- | Shows the internal representation of a 'CR'. The first /n/+-- approximations are shown on separate lines.+showCRN :: Int -> CR -> String+showCRN n (CR x) = concat . intersperse "\n" . map showA . take n . getZipList $ x++-- | Shows a 'CR' with the desired precision.+showCR :: Int -> CR -> String+showCR p = showA . require p++-- There is no show instance of 'CR' since the representation would be infinite. We can therefore not satisfy (read . show) = id.++-- | Reads a floating point representation of a real number and interprets+-- that as a CR. Does not currently allow for the same format output by+-- 'showCR'.+instance Read CR where+ readsPrec _ input =+ let (intPart, rest) = span isDigit input+ in if null rest || head rest /= '.'+ then [(CR $ pure (fromInteger (read intPart :: Integer)), rest)]+ else let (fracPart, rest') = span isDigit (tail rest)+ in [((CR $ pure (fromInteger (read (intPart ++ fracPart) :: Integer))) / 10^(length fracPart), rest')]++-- | Check that an approximation has at least /d/ bits of precision. This is+-- used to bail out of computations where the size of approximation grow+-- quickly.+ok :: Int -> Approx -> Approx+ok d a = if precision a > fromIntegral d then a else Bottom++-- | Given a 'CR' this functions finds an approximation of that number to+-- within the precision required.+require :: Int -> CR -> Approx+require d (CR x) = head . dropWhile (== Bottom) . getZipList $ ok d <$> x++-- | Gives a 'Double' approximation of a 'CR' number.+toDouble :: CR -> Maybe Double+toDouble = toDoubleA . require (54+errorBits)++fromDouble :: Double -> CR+fromDouble x =+ let (m, s) = decodeFloat x+ -- When the mantissa of a floating point is interpreted as a whole number+ -- instead of as a fraction in the IEEE 754 encoding the exponent 972+ -- corresponds to 1024, which is what IEEE 754 use to encode infinity and+ -- NaN.+ in if (m == 972) then CR $ pure Bottom+ else CR $ pure (Approx 64 m 1 s)++fromDoubleAsExactValue :: Double -> CR+fromDoubleAsExactValue x =+ let (m, s) = decodeFloat x+ -- When the mantissa of a floating point is interpreted as a whole number+ -- instead of as a fraction in the IEEE 754 encoding the exponent 972+ -- corresponds to 1024, which is what IEEE 754 use to encode infinity and+ -- NaN.+ in if (m == 972) then CR $ pure Bottom+ else CR $ pure (Approx 64 m 0 s)++transposeZipList :: [ZipList a] -> ZipList [a]+transposeZipList = ZipList . transpose . map getZipList++-- | Evaluate a polynomial, given as a list of its coefficients, at the given point.+polynomial :: [CR] -> CR -> CR+polynomial as (CR x) = + CR $ (\as' x' l -> ok 10 . limitAndBound l $ poly as' x') <$> transposeZipList (map unCR as) <*> x <*> resources++-- | Computes the sum of a Taylor series, given as a list of its coefficients,+-- at the given point.+taylorCR :: [CR] -> CR -> CR+taylorCR as (CR x) =+ CR $ (\as' x' l -> sum . takeWhile nonZeroCentredA . map (limitAndBound l) $ zipWith (*) as' (pow x'))+ <$> transposeZipList (map unCR as) <*> x <*> resources++epsilon :: CR+epsilon = CR $ Approx 10 0 1 . negate <$> resources++-- | The square root function. Lifted from square root application on 'Approx'+-- approximations.+sqrtCR :: CR -> CR+sqrtCR (CR x) = CR $ op1withResource sqrtA id <$> x <*> resources++alternateSign :: Num a => [a] -> [a]+alternateSign = zipWith (*) (cycle [1,-1])++atanCR :: CR -> CR+atanCR x =+ let rr y = y / (1 + sqrt (1 + x^2))+ x' = rr . rr . rr $ x -- range reduction so that |a'| < 1/4+ x2 = - x'^2+ (CR t) = epsilon + x * taylorCR (map ((1/) . fromIntegral) [1,3..]) x2+ in CR $ boundErrorTerm . (8*) <$> t+-- let x2 = x^2+-- in epsilon + x * taylor (map (1/) . alternateSign . map fromInteger $ [1,3..]) x2++-- | π computed using Machin's formula. Computed directly on 'CR'.+piCRMachin :: CR+piCRMachin = 4*(4*atanCR (1/5)-atanCR (1/239))++-- | π computed using Machin's formula. Computed on 'Approx' approximations.+piMachinCR :: CR+piMachinCR = CR $ piMachinA . negate <$> resources++-- | π computed using Borwein's formula. Computed on 'Approx' approximations.+piBorweinCR :: CR+piBorweinCR = CR $ piBorweinA . negate <$> resources++-- | π computed using binary splitting. Computed on 'Approx' approximations.+piBinSplitCR :: CR+piBinSplitCR = CR $ limitAndBound <$> resources <*> (require <$> resources <*> ZipList (repeat (CR $ ZipList piRaw)))++-- | The constant ln 2.+ln2 :: CR+ln2 = CR $ log2A . negate <$> resources++-- | The exponential computed using Taylor's series. Computed directly on+-- 'CR'. Will have poor behaviour on larger inputs as no range reduction is+-- performed.+expCR :: CR -> CR+expCR = (+ epsilon) . taylorCR (map (1/) $ fac)++halfPi :: CR+halfPi = scale pi (-1)++sinRangeReduction :: CR -> CR+sinRangeReduction (CR x) = (subtract halfPi) . abs . (pi -) . abs . (subtract halfPi) . CR $ modA <$> x <*> unCR (2 * pi)++sinRangeReduction2 :: CR -> CR+sinRangeReduction2 (CR x) =+ let k = (\a -> case a of + (Approx _ m _ s) -> max 0 $ 8 * (integerLog2 m + s + 3) `div` 5+ Bottom -> 0) <$> x+ (CR y) = sinCRTaylor ((CR x) / (CR $ fromIntegral . (3^) <$> k))+ step z = z*(3-4*z^2)+ in CR $ (\y' k' l -> limitAndBound l $ foldr ($) y' (replicate k' step)) <$> y <*> k <*> resources++-- | The sine function computed using Taylor's series. Computed directly on+-- 'CR'. Will have poor behaviour on larger inputs as no range reduction is+-- performed.+sinCRTaylor :: CR -> CR+sinCRTaylor x = let x2 = x^2+ in epsilon + x * taylorCR (map (1/) $ alternateSign oddFac) x2++-- | The sine function computed using Taylor's series. Computed directly on+-- 'CR'.+sinCR :: CR -> CR+sinCR = sinRangeReduction2 . sinRangeReduction++-- | The cosine function computed using Taylor's series. Computed directly on+-- 'CR'.+cosCR :: CR -> CR+cosCR = sinCR . (halfPi -)++instance Floating CR where+ sqrt (CR x) = CR $ op1withResource sqrtA id <$> x <*> resources+ pi = piBinSplitCR+ exp (CR x) = CR $ op1withResource expA id <$> x <*> resources+ log (CR x) = CR $ op1withResource logA id <$> x <*> resources+ sin (CR x) = CR $ op1withResource sinA id <$> x <*> resources+ cos (CR x) = CR $ op1withResource cosA id <$> x <*> resources+ asin x = 2 * (atan (x / (1 + (sqrt (1 - x^2)))))+ acos x = halfPi - asin x+ atan (CR x) = CR $ atanA <$> resources <*> x+ sinh x = ((exp x) - (exp $ negate x)) / 2+ cosh x = ((exp x) + (exp $ negate x)) / 2+ tanh x = let t = exp (2*x) in (t-1)/(t+1)+ asinh x = log (x + sqrt (x^2 + 1))+ acosh x = CR $ op1withResource logA id <$> unCR (x + sqrt (x^2 - 1)) <*> resources+ atanh x = (CR $ op1withResource logA id <$> unCR ((1+x) / (1-x)) <*> resources) / 2
+ src/Data/CDAR/Classes.hs view
@@ -0,0 +1,16 @@+{-|+Defines the 'Scalable' type class.+-}+module Data.CDAR.Classes where++import Data.Bits++-- | 'Scalable' allows scaling numerical data types by powers of 2.+class Scalable a where+ scale :: a -> Int -> a++-- | The 'Integer' instance.+instance Scalable Integer where+ scale x n+ | n >= 0 = unsafeShiftL x n+ | otherwise = unsafeShiftR (x + bit (-n-1)) (-n)
+ src/Data/CDAR/Dyadic.hs view
@@ -0,0 +1,264 @@+{- |The Dyadic module provides dyadic numbers. Syntax for dyadic number+ 'a * 2 ^ s' is 'a :^ s'. The exponent 's' is an 'Int', but the 'a' is an+ arbitrary 'Integer'.+-}+module Data.CDAR.Dyadic (Dyadic(..),shiftD,sqrtD,sqrtD',sqrtRecD,sqrtRecD',initSqrtRecD,initSqrtRecDoubleD,piMachinD,piBorweinD,ln2D,logD,atanhD,divD,divD') where++import Data.Ratio+import Data.Bits+import Data.CDAR.Classes+import Data.CDAR.IntegerLog++-- Dyadic numbers++infix 8 :^++-- |The Dyadic data type.+data Dyadic = Integer :^ Int deriving (Read,Show)++instance Eq Dyadic where+ a :^ s == b :^ t | s <= t = a == unsafeShiftL b (t-s)+ | otherwise = unsafeShiftL a (s-t) == b++instance Ord Dyadic where+ compare (a :^ s) (b :^ t) | s <= t = compare a (unsafeShiftL b (t-s))+ | otherwise = compare (unsafeShiftL a (s-t)) b++instance Scalable Dyadic where+ scale (a :^ s) n = a :^ (s+n)++instance Num Dyadic where+ a :^ s + b :^ t + | s <= t = (a + scale b (t-s)) :^ s+ | otherwise = (scale a (s-t) + b) :^ t+ a :^ s * b :^ t = (a * b) :^ (s+t)+ negate (a :^ s) = (negate a) :^ s+ abs (a :^ s) = (abs a) :^ s+ signum (a :^ _) = (signum a):^ 0+ fromInteger i = i :^ 0++instance Real Dyadic where+ toRational (a :^ s) = (toRational a)*2^^s++-- | Shift a dyadic number to a given base and round in case of right shifts.+shiftD :: Int -> Dyadic -> Dyadic+shiftD t (m:^s) =+ if t <= s+ then unsafeShiftL m (s-t) :^ t+ else unsafeShiftR (m + bit (t-s-1)) (t-s) :^ t++-- Square root++divD :: Dyadic -> Dyadic -> Dyadic+divD a (n:^t) = let (m:^_) = shiftD (2*t) a+ in round (m % n) :^ t++{- |Computes the square root of a Dyadic number to the specified base. The+ Newton-Raphson method may overestimates the square root, but the+ overestimate is bounded by 1 ulp. For example, sqrtD 0 2 will give 2,+ whereas the closest integer to the square root is 1. Need double precision+ to guarantee correct rounding, which was not considered worthwhile.++ This is actually Heron's method and is no longer used in Approx as it is+ faster to use sqrtRecD.+-}+sqrtD :: Int -> Dyadic -> Dyadic+sqrtD t x = sqrtD' t x $ initSqrtD x+ where+ -- Initial approximation of square root to about 3 bits+ initSqrtD :: Dyadic -> Dyadic+ initSqrtD (0:^_) = (0:^0)+ initSqrtD (m:^s) = let i = integerLog2 m+ n = shift m (2-i)+ s' = (s+i) `div` 2 - 3+ in if odd (s+i)+ then (n+8):^s'+ else (n+4):^s'++-- |Square root with initial approximation as third argument.+sqrtD' :: Int -> Dyadic -> Dyadic -> Dyadic+sqrtD' t x@(m:^_) y+ | m == 0 = 0:^0+ | m > 0 = converge . iterate (newtonStep x t) $ shiftD t y+ | otherwise = error "Attempting sqrt of negative dyadic number."+ where+ -- One step of Newton iteration to find square root.+ -- The base of a need to be the same as t', and the result will have base t'.+ newtonStep :: Dyadic -> Int -> Dyadic -> Dyadic+ newtonStep d t' a = shiftD t' $ (1:^(-1)) * (a + divD d a)++-- |Reciprocal of square root using Newton's iteration.+sqrtRecD :: Int -> Dyadic -> Dyadic+sqrtRecD t a = sqrtRecD' t a $ initSqrtRecDoubleD a++-- |Gives initial values for the iteration. Based on the three most+-- significant bits of the argument. Should consider to use machine floating+-- point to give initial value.+initSqrtRecD :: Dyadic -> Dyadic+initSqrtRecD (0 :^ _) = error "Divide by zero in initSqrtRecD"+initSqrtRecD (m :^ s) =+ let i = integerLog2 m+ n = shift m (3-i)+ s' = (-i-s-1) `div` 2 - 5+ in if even (s+i)+ then ([62,59,56,53,51,49,48,46]!!(fromIntegral n-8)) :^ s'+ else ([44,42,40,38,36,35,34,33]!!(fromIntegral n-8)) :^ s'++-- Using Double to find a good first approximation to SqrtRec. Extracts 52+-- bits from the dyadic number and encodes a float in the range [1/2,2).+-- Decodes the result of the Double computation as a dyadic taking care of the+-- exponent explicitly (as the range of Double exponents is too small).+initSqrtRecDoubleD :: Dyadic -> Dyadic+initSqrtRecDoubleD (0 :^ _) = error "Divide by zero in initSqrtRecDoubleD"+initSqrtRecDoubleD (m :^ s) =+ let i = integerLog2 m+ n = shift m (52-i)+ s' = (-i-s) `div` 2 - 52+ (m',_) = decodeFloat . (1/) . sqrt $ encodeFloat n (if even (s+i) then (-52) else (-51))+ t = if m' /= bit 52 && even (s+i) then s'-1 else s'+ in m' :^ t++-- |Reciprocal of square root using Newton's iteration with inital value as third argument.+sqrtRecD' :: Int -> Dyadic -> Dyadic -> Dyadic+sqrtRecD' t a x0 =+ let step x = shiftD t $ x + shiftD t (x * (1 - shiftD t (x * x * a)) * (1 :^ (-1)))+ xs = iterate step x0+ in converge xs++-- isqrt :: Integral t => t -> t+-- isqrt 0 = 0+-- isqrt 1 = 1+-- isqrt n = head $ dropWhile (\x -> x*x > n) $ iterate (\x -> (x + n `div` x) `div` 2) (n `div` 2)++-- Assuming a seqence of dyadic numbers with the same exponent converges+-- quadratically. Then this returns the first element where it is known that+-- the sequence will become essentially constant (differing by at most 1). If+-- the numbers contain sufficiently many bits we need only wait until the+-- difference has less than half as many bits (with a small margin).+converge :: [Dyadic] -> Dyadic+converge ((0:^s):_) = (0:^s)+converge ((m:^_):ds@(d@(n:^_):_)) =+ let a = integerLog2 m+ b = integerLog2 (abs (m - n))+ in if a <= 6+ then if abs (m-n) <= 1 then d else converge ds+ else if 2 * b <= a - 3 then d else converge ds+converge _ = error "List terminating in converge."++divD' :: Int -> Dyadic -> Dyadic -> Dyadic+divD' p a b = let (m:^_) = shiftD (2*p) a+ (n:^_) = shiftD p b+ in round (m%n) :^ p++-- |Compute dyadic values close to pi by Machin's formula.+piMachinD :: Int -> Dyadic+piMachinD t = let t' = t-10-integerLog2 (fromIntegral (abs t))+ a = map ((:^t') . round . (bit (-t') %)) $ iterate ((-25)*) (5 :: Integer)+ b = map ((:^t') . round . (bit (-t') %)) $ iterate ((-57121)*) (239 :: Integer)+ c = map ((:^t') . round . (bit (-t') %)) ([1,3..] :: [Integer])+ d = takeWhile (/= 0) . map (shiftD t') $ zipWith (*) a c+ e = takeWhile (/= 0) . map (shiftD t') $ zipWith (*) b c+ in shiftD t $ 4 * (4 * sum d - sum e)++-- |Compute dyadic values close to pi by Borwein's formula.+piBorweinD :: Int -> Dyadic+piBorweinD t = let t' = t-10-integerLog2 (fromIntegral (abs t))+ s = sqrtD t' 2+ d = head . dropWhile (\(_,y,_) -> y /= 0) $ iterate (f t') (6-4*s,s-1,1 :: Int)+ in shiftD t . (\(a,_,_) -> divD' t' 1 a) $ d+ where+ f :: Int -> (Dyadic, Dyadic, Int) -> (Dyadic, Dyadic, Int)+ f l (a,y,k) = let u = sqrtD l . sqrtD l $ 1-y^(4::Int)+ y' = shiftD l $ divD (1-u) (1+u)+ a' = shiftD l $ a*(1+y')^(4::Int)-(bit (2*k+1):^0)*y'*(1+y'+y'*y')+ in (a',y',k+1)++-- |Compute dyadic values close to ln 2.+ln2D :: Int -> Dyadic+ln2D t = let t' = t - 10 - 2 * integerLog2 (fromIntegral (abs t))+ a = map ((:^t') . round . (bit (-t') %)) $ iterate (3*) (3 :: Integer)+ b = map ((:^t') . round . (bit (-t') %)) $ iterate (4*) (4 :: Integer)+ c = map ((:^t') . round . (bit (-t') %)) ([1..] :: [Integer])+ d = zipWith (+) a b+ e = takeWhile (/= 0) . map (shiftD t') $ zipWith (*) d c+ in shiftD t $ sum e++logD :: Int -> Dyadic -> Dyadic+logD t x@(a :^ s) =+ if a <= 0 then error "logD: Non-positive argument"+ else let t' = t - 10 - max 0 (integerLog2 a + s) -- (5 + size of argument) guard digits+ r = s + integerLog2 (3*a) - 1+ x' = scale x (-r) -- 2/3 <= x' <= 4/3+ y = divD' t' (x' - 1) (x' + 1) -- so |y| <= 1/5+ z = fromIntegral r * ln2D t'+ in shiftD t . (+z) . flip scale 1 $ atanhD t' y++atanhD :: Int -> Dyadic -> Dyadic+atanhD t x@(a :^ s) =+ if integerLog2 (abs a) + s >= 0 then error "atanhD: Argument outside domain, (-1,1)"+ else let + -- number of guard digits is 5+k where k depends on the precision [how do we know if this is enough?]+ t' = t - 5 - integerLog2 (abs $ fromIntegral t) `div` 2+ g _x _y = shiftD t' (_x * _y)+ x2 = g x x+ b = iterate (g x2) 1+ c = (map (divD' t' 1 . fromIntegral) [1,3..])+ in shiftD t . g x . sum . takeWhile (/= 0) $ zipWith g b c++{-+agmD :: Int -> Dyadic -> Dyadic -> Dyadic+agmD t a b = let t' = t - 5+ agmStep (c,d) = ((1:^(-1)) * (c+d), sqrtD t' (c*d))+ close (c,d) = abs (c-d) < 1:^t+ in fst . last . takeWhile (not . close) $ iterate agmStep (shiftD t' a, shiftD t' b)++theta2D :: Int -> Dyadic -> Dyadic+theta2D t q = let q2 = shiftD t $ q*q+ q4 = sqrtD t . sqrtD t $ q+ a = iterate (shiftD t . (q2*)) 1+ b = scanl1 (\c d -> shiftD t (c*d)) a+ in (2*) . shiftD t . (q4*) . sum $ takeWhile (/= 0) b++theta3D :: Int -> Dyadic -> Dyadic+theta3D t q = let q2 = shiftD t $ q*q+ -- q4 = sqrtD t . sqrtD t $ q+ a = iterate (shiftD t . (q2*)) q+ b = scanl1 (\c d -> shiftD t (c*d)) a+ in (1+) . (2*) . shiftD t . sum $ takeWhile (/= 0) b++lnBigD :: Int -> Dyadic -> Dyadic+lnBigD t x = let t' = t-10-integerLog2 (fromIntegral (abs t))+ t2 = theta2D t' x+ t3 = theta3D t' x+ a = agmD t' (t2*t2) (t3*t3)+ in shiftD t $ divD' t' (piBorweinD t') a+-}+{-+agm a b = let step (a,b) = (0.5 * (a+b), sqrt (a*b))+ close (a,b) = abs (a-b) < 1e-6+ in fst . last . takeWhile (not . close) $ iterate step (a, b)++theta2 q = let q2 = q^2+ q4 = sqrt . sqrt $ q+ a = iterate (q2*) 1+ b = scanl1 (*) a+ in (2*) . (q4*) . sum $ takeWhile (> 1e-6) b++theta3 q = let q2 = q^2+ q4 = sqrt . sqrt $ q+ a = iterate (q2*) q+ b = scanl1 (*) a+ in (1+) . (2*) . sum $ takeWhile (> 1e-6) b+-}++{-+-- checking that converging dyadic numbers only differ in last bit+check ((1:^_):(1:^_):_) = False+check ((1:^_):xs) = check xs+check ((-1:^_):(-1:^_):_) = False+check ((-1:^_):xs) = check xs+check ((0:^_):xs) = check xs+check ((_:^_):_) = False+check [] = True+check [_] = True+-}
+ src/Data/CDAR/Extended.hs view
@@ -0,0 +1,53 @@+{- |The Extended module allows real-valued numeric data types to be extended by+ positive and negative infinity.+-}+{-# LANGUAGE DeriveFunctor #-}+module Data.CDAR.Extended (Extended(..)) where++import Control.Monad++-- |Extended numbers are either finite numbers or one of the two infinities.+data Extended a = PosInf | NegInf | Finite a deriving (Eq,Read,Show,Functor)++instance Applicative Extended where+ pure = Finite+ (Finite f) <*> (Finite x) = Finite $ f x+ (Finite _) <*> PosInf = PosInf+ (Finite _) <*> NegInf = NegInf+ PosInf <*> _ = PosInf+ NegInf <*> _ = NegInf++instance Monad Extended where+ return = pure+ (Finite x) >>= f = f x+ PosInf >>= _ = PosInf+ NegInf >>= _ = NegInf++instance Ord a => Ord (Extended a) where+ compare PosInf PosInf = EQ+ compare NegInf NegInf = EQ+ compare _ PosInf = LT+ compare NegInf _ = LT+ compare PosInf _ = GT+ compare _ NegInf = GT+ compare (Finite a) (Finite b) = compare a b++instance (Ord a, Num a) => Num (Extended a) where+ -- PosInf + NegInf should be undefined, but here it is the first argument+ (+) = liftM2 (+)+ -- 0 * ???Inf should be undefined, but here it is PosInf+ a * PosInf = if a < 0 then NegInf else PosInf+ PosInf * a = if a < 0 then NegInf else PosInf+ a * NegInf = if a < 0 then PosInf else NegInf+ NegInf * a = if a < 0 then PosInf else NegInf+ a * b = (*) <$> a <*> b+ negate = fmap negate+ abs = fmap abs+ signum PosInf = Finite 1+ signum NegInf = Finite (-1)+ signum a = signum <$> a+ fromInteger i = Finite $ fromInteger i++instance Real a => Real (Extended a) where+ toRational (Finite x) = toRational x+ toRational _ = undefined
+ src/Data/CDAR/IntegerLog.hs view
@@ -0,0 +1,43 @@+{-# LANGUAGE MagicHash #-}++{- |The IntegerLog module provides a fast integer logarithm in base 2. Current implementation relies in the integer-gmp package.+-}+module Data.CDAR.IntegerLog (integerLog2) where++import GHC.Integer.Logarithms -- only reason for needing integer-gmp package+import GHC.Exts++{-# INLINE integerLog2 #-}+{- |The 'integerLog2' function gives the floor of the logarithm of a positive integer. For non-positive arguments the returned result is meaningless. Possible overflow of returned value is ignored.++It could be argued that NegInf should be returned for 0 if the return type had been 'Extended' 'Int'.+-}+integerLog2 :: Integer -> Int+integerLog2 n = I# (integerLog2# n)++{- Should be able to rely on some provided integer logarithms from ghc 7.2+ onwards. A Haskell2010 version previously used follows. ++integerLog2 :: Integer -> Int+integerLog2 i =+ if i < 2 then 0+ else integerLog2' i 1 2++integerLog2' :: Integer -> Int -> Int -> Int+integerLog2' i l h =+ if i < bit h then+ integerLog2'' i l h+ else+ integerLog2' i h (2*h)++integerLog2'' :: Integer -> Int -> Int -> Int+integerLog2'' i l h =+ if (h-l) <= 1 then+ l+ else+ if i < bit m then+ integerLog2'' i l m+ else+ integerLog2'' i m h+ where m = (l+h) `div` 2+-}
+ test/Test.hs view
@@ -0,0 +1,298 @@+{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses, StandaloneDeriving, DeriveGeneric #-}+import Test.Tasty+import Test.Tasty.SmallCheck as SC+import Test.Tasty.QuickCheck as QC+import Test.Tasty.HUnit+import Test.SmallCheck.Series as SC+import GHC.Generics++import Control.Applicative (ZipList (..))+import Control.Monad (liftM, liftM3, liftM4)+import Data.Functor++import Data.List+import Data.Ord+import Data.Fixed (mod')+import Data.CDAR+++-- Main test++main = defaultMain tests++tests :: TestTree+tests = testGroup "Tests" [dyadic, localOption (mkTimeout 500000) approx, creal, unitTests]+ -- [dyadic+ -- ,approx+ -- ,creal+ -- ,unitTests]++-- General auxiliary properties.++idempotent :: Eq a => (a -> a) -> a -> Bool+idempotent f = \x -> f(x) == f(f(x))++autoinverse :: Eq a => (a -> a) -> a -> Bool+autoinverse f = \x -> x == f(f(x))++dyadic :: TestTree+dyadic = testGroup "Dyadic" [propertiesD]++deriving instance Generic Dyadic+instance (Monad m) => Serial m Dyadic++instance Arbitrary Dyadic where+ arbitrary = do a <- arbitrary+ s <- choose (-1000,100)+ return (a :^ s)++propertiesD :: TestTree+propertiesD = testGroup "Properties of Dyadic" [scPropD, qcPropD]++scPropD = testGroup "(checked by SmallCheck)"+ [ SC.testProperty "compare" $+ \x y -> compare (x :: Dyadic) (y :: Dyadic) == (compare (toRational x) (toRational y))+ , SC.testProperty "multiplication" $+ \a b -> toRational (a*b :: Dyadic) == (toRational a)*(toRational b)+ , SC.testProperty "addition" $+ \a b -> toRational (a+b :: Dyadic) == (toRational a)+(toRational b)+ , SC.testProperty "negate" $+ \a -> (a :: Dyadic) + negate a == 0+ , SC.testProperty "negation exists" $+ \a b -> (a :: Dyadic) + b == 0 SC.==> a == negate b+ , SC.testProperty "negation autoinverse" $ autoinverse (negate :: Dyadic -> Dyadic)+ , SC.testProperty "abs idempotent" $ idempotent (abs :: Dyadic -> Dyadic)+ , SC.testProperty "signum idempotent" $ idempotent (signum :: Dyadic -> Dyadic)+ , SC.testProperty "signum in {-1,0,1}" $ \x -> signum (x :: Dyadic) `elem` [-1,0,1]+ , SC.testProperty "negate abs signum" $+ \a -> (a :: Dyadic) + (negate (signum a) * (abs a)) == 0+ , SC.testProperty "read show" $+ \a -> (a :: Dyadic) == read (show a)+ , SC.testProperty "shift left identity" $ \d@(a :^ s) (SC.NonNegative n) -> shiftD (s - n) d == d+ , SC.testProperty "sqrt" $ \d@(a :^ s) -> sqrtD (min 0 (2*s)) (d*d) == abs d+ ]++qcPropD = testGroup "(checked by QuickCheck)"+ [ QC.testProperty "compare" $ + \x y -> compare (x :: Dyadic) (y :: Dyadic) == (compare (toRational x) (toRational y))+ , QC.testProperty "multiplication" $+ \a b -> toRational (a*b :: Dyadic) == (toRational a)*(toRational b)+ , QC.testProperty "addition" $+ \a b -> toRational (a+b :: Dyadic) == (toRational a)+(toRational b)+ , QC.testProperty "negate" $+ \a -> (a :: Dyadic) + negate a == 0+ , QC.testProperty "negate abs signum" $+ \a -> (a :: Dyadic) + (negate (signum a) * (abs a)) == 0:^0+ , QC.testProperty "read show" $+ \a -> (a :: Dyadic) == read (show a)+ , QC.testProperty "negation autoinverse" $ autoinverse (negate :: Dyadic -> Dyadic)+ , QC.testProperty "abs idempotent" $ idempotent (abs :: Dyadic -> Dyadic)+ , QC.testProperty "signum idempotent" $ idempotent (signum :: Dyadic -> Dyadic)+ , QC.testProperty "signum in {-1,0,1}" $ \x -> signum (x :: Dyadic) `elem` [-1,0,1]+ , QC.testProperty "shift left identity" $ \d@(a :^ s) (QC.NonNegative n) -> shiftD (s - n) d == d+ , QC.testProperty "sqrt" $ \d@(a :^ s) -> sqrtD (min 0 (2*s)) (d*d) == abs d+ ]++{-+Multiplication and addition of Approx is not tested. Those two operations should return the exact Approx of the result.++There are no test making sure that the loss of information is acceptable for reciprocal, exp, log, ... .+-}+approx :: TestTree+approx = testGroup "Approx" [scPropA, qcPropA]++instance Monad m => Serial m Approx where+ series = cons0 Bottom \/ (quadToApprox <$> series)+ where+ quadToApprox (SC.Positive mb, m, SC.NonNegative e, s) = approxMB mb m e s++instance Arbitrary Approx where+ arbitrary = frequency [(10,+ liftM3 approxAutoMB+ arbitrary + (return 0)+ (choose (-100,20))),+ (25,+ liftM4 approxMB+ (choose (10,100000))+ arbitrary + (elements [1,2,3,254,255,256,257,510,511,512,513,767,768,1020,1021,1022,1023,1024,1025])+ (choose (-100,20))),+ (60,+ liftM4 approxMB+ (choose (10,100000))+ arbitrary + (choose (0,10))+ (choose (-100,20))),+ (5, return Bottom)]++scPropA :: TestTree+scPropA = testGroup "(checked by smallcheck)"+ [ SC.testProperty "read . show = id" $ \a -> (a :: Approx) == read (show a)+ , SC.testProperty "diameter" $ \a -> diameter a == upperBound a - lowerBound a+ , SC.testProperty "fromDyadic exact and centre" $ \d -> let a = fromDyadic d in exact a && Just d == centre a+ , SC.testProperty "abs negate sign" $ \a -> approximatedBy 0 $ (a ::Approx) + (negate (signum a) * (abs a))+ , SC.testProperty "a+(-a) contains 0" $ \a -> 0 `approximatedBy` ((a :: Approx) + negate a)+ , SC.testProperty "toApprox r contains r" $ \r -> approximatedBy r $ toApprox 20 (r :: Rational)+ , SC.testProperty "1/(1/a) contains a" $ \a -> not (0 `approximatedBy` a) SC.==> + a `better` (1/(1/a))+ , SC.testProperty "a*(1/a) contains 1" $ \a -> not (0 `approximatedBy` a) SC.==> 1 `approximatedBy` (a*(1/a))+ , SC.testProperty "boundErrorTerm" $+ \a -> case a of+ Approx _ _ _ _ -> let b@(Approx _ _ e _) = boundErrorTerm a+ in (a `better` b) && (e < 1024) -- 1024 is errorBound when errorBits is 10+ Bottom -> better a (boundErrorTerm a)+ , SC.testProperty "limitSize" $+ \a -> case a of+ Approx _ _ _ _ -> let b@(Approx _ _ _ s) = limitSize 2 a+ in (a `better` b) && (s >= -2)+ Bottom -> better a (limitSize 2 a)+ , SC.testProperty "sqrt" $ \a -> let b = abs a in better b $ (sqrtA b)^2+ , SC.testProperty "recipA" $ \a -> let b = abs a in 0 `approximatedBy` b || (lowerBound (recipA b) * upperBound b <= 1+ && upperBound (recipA b) * lowerBound b >= 1)+ ]++qcPropA :: TestTree+qcPropA = testGroup "(checked by quickcheck)"+ [ QC.testProperty "read . show = id" $ \a -> (a :: Approx) == read (show a)+ , QC.testProperty "diameter" $ \a -> diameter a == upperBound a - lowerBound a+ , QC.testProperty "fromDyadic exact and centre" $ \d -> let a = fromDyadic d in exact a && Just d == centre a+ , QC.testProperty "abs negate sign" $ \a -> approximatedBy 0 $ (a ::Approx) + (negate (signum a) * (abs a))+ , QC.testProperty "addition precision" $ \a b -> + collect (precision (boundErrorTerm (a+b)) - (precision a `min` precision b))+ True+{-+ , QC.testProperty "multiplication precision" $ \a b -> + collect (precision (boundErrorTerm (a*b)) - (precision a `min` precision b))+ True+-}+ , QC.testProperty "a+(-a) contains 0" $ \a -> 0 `approximatedBy` ((a :: Approx) + negate a)+ , QC.testProperty "toApprox r contains r" $ \r -> approximatedBy r $ toApprox 20 (r :: Rational)+ , QC.testProperty "1/(1/a) contains a" $ \a -> not (0 `approximatedBy` a) QC.==> + a `better` (1/(1/a))+ , QC.testProperty "a*(1/a) contains 1" $ \a -> not (0 `approximatedBy` a) QC.==> 1 `approximatedBy` (a*(1/a))+ , QC.testProperty "boundErrorTerm" $+ \a -> case a of + Approx _ _ _ _ -> let b@(Approx _ _ e _) = boundErrorTerm a+ in (a `better` b) && (e < 1024) -- 1024 is errorBound when errorBits is 10+ Bottom -> better a (boundErrorTerm a)+ , QC.testProperty "limitSize" $+ \a -> case a of + Approx _ _ _ _ -> let b@(Approx _ _ _ s) = limitSize 2 a+ in (a `better` b) && (s >= -2)+ Bottom -> better a (limitSize 2 a)+ , QC.testProperty "values" $ \a -> let types = a :: Approx in collect a True+ , QC.testProperty "sqr . sqrt" $ \a -> let b = abs a in b `better` sqrA (sqrtA b)+ , QC.testProperty "sqrt . sqr" $ \a -> abs a `better` sqrtA (sqrA a)+ , QC.testProperty "recipA" $ \a -> let b = abs a in 0 `approximatedBy` b || (lowerBound (recipA b) * upperBound b <= 1+ && upperBound (recipA b) * lowerBound b >= 1)+ ]++creal :: TestTree+creal = testGroup "CReal" [propertiesCR]++genCRApprox :: Gen Approx -> Gen CR+genCRApprox = liftM (CR . pure)++genCR :: Gen CR+genCR = do+ x <- choose (0,1) :: Gen Double+ let (m,_) = decodeFloat x+ s <- choose (-125,-42) :: Gen Int+ s' <- choose (-56,-51) :: Gen Int+ frequency [(2, return . CR . pure $ approxAutoMB m 0 s)+ ,(2, return . CR . pure $ approxAutoMB m 0 s')+ ,(2, return . CR . pure $ approxAutoMB (-m) 0 s)+ ,(2, return . CR . pure $ approxAutoMB (-m) 0 s')+ ,(1, return . CR . pure $ approxAutoMB m 1 s)+ ,(1, return . CR . pure $ approxAutoMB m 1 s')+ ,(1, return . CR . pure $ approxAutoMB (-m) 1 s)+ ,(1, return . CR . pure $ approxAutoMB (-m) 1 s')+ ]++genCROpen :: CR -> CR -> Gen CR+genCROpen a b = do+ x <- choose (0,1) :: Gen Double+ return $ a + (b-a) * fromDoubleAsExactValue x++genCRClosed :: CR -> CR -> Gen CR+genCRClosed a b = do+ x <- choose (0,1) :: Gen Double+ frequency [(96, return $ a + (b-a) * fromDoubleAsExactValue x)+ ,(2, return a)+ ,(2, return b)+ ]++instance Arbitrary CR where+ arbitrary = genCR++instance Show CR where+ show (CR x) = showA . head . getZipList $ x++checkCRN :: Int -> CR -> CR -> Bool+checkCRN n (CR x) (CR y) = and $ zipWith better (take n $ getZipList x) (take n $ getZipList y)++checkFun :: (CR -> CR) -> CR -> Bool+checkFun f x = let a = require 1 x+ l = lowerA a+ u = upperA a+ fx = require 1 (f x)+ in better (require 200 . f . CR . pure $ l) fx &&+ better (require 200 . f . CR . pure $ u) fx++propertiesCR :: TestTree+propertiesCR = testGroup "Properties of CReal" [scPropCR, qcPropCR]++scPropCR = testGroup "(checked by smallCheck)"+ []++{-+Some of the test should hold for larger domains than tested over. Consider,+for example, atanh . tanh: with an input of -1000, tanh will give essentially+-1 with some error bound. Given that atanh is computed by a logarithm applied+to (x+1)/(x-1), which is essentially 0, we may exhaust memory before a good+enough result is computed.+-}+qcPropCR = testGroup "(checked by quickCheck)"+ [ QC.testProperty "values" $ \x -> let types = x :: CR in collect (show x) True+ , QC.testProperty "Pythagorean identity" $ QC.forAll (genCRClosed (-pi) pi) $ \x -> checkCRN 5 1 ((sin x)^2 + (cos x)^2)+ , QC.testProperty "^2 . sqrt" $ QC.forAll (genCRClosed 0 100) $ \x -> checkCRN 5 x ((sqrt x)^2)+ , QC.testProperty "sqrt . ^2" $ \x -> checkCRN 5 (abs x) (sqrt (x^2))+ , QC.testProperty "log . exp" $ \x -> checkCRN 5 x (log (exp x))+ , QC.testProperty "exp . log" $ QC.forAll (genCROpen 0 100) $ \x -> checkCRN 5 x (exp (log x))+ , QC.testProperty "asin . sin" $ QC.forAll (genCROpen (-pi/2) (pi/2)) $ \x -> checkCRN 5 x (asin (sin x))+ , QC.testProperty "sin . asin" $ QC.forAll (genCROpen (-1) 1) $ \x -> checkCRN 5 x (sin (asin x))+ , QC.testProperty "acos . cos" $ QC.forAll (genCROpen 0 pi) $ \x -> checkCRN 5 x (acos (cos x))+ , QC.testProperty "cos . acos" $ QC.forAll (genCROpen (-1) 1) $ \x -> checkCRN 5 x (cos (acos x))+ , QC.testProperty "atan . tan" $ QC.forAll (genCROpen (-pi/2) (pi/2)) $ \x -> checkCRN 5 x (atan (tan x))+ , QC.testProperty "tan . atan" $ \x -> checkCRN 1 x (tan (atan x))+ , QC.testProperty "asinh . sinh" $ \x -> checkCRN 5 x (asinh (sinh x))+ , QC.testProperty "sinh . asinh" $ \x -> checkCRN 5 x (sinh (asinh x))+ , QC.testProperty "acosh . cosh" $ QC.forAll (genCRClosed 0 10) $ \x -> checkCRN 5 x (acosh (cosh x))+ , QC.testProperty "cosh . acosh" $ QC.forAll (genCRClosed 1 10) $ \x -> checkCRN 5 x (cosh (acosh x))+ , QC.testProperty "atanh . tanh" $ QC.forAll (genCRClosed (-100) 100) $ \x -> checkCRN 5 x (atanh (tanh x))+ , QC.testProperty "tanh . atanh" $ QC.forAll (genCROpen (-1) 1) $ \x -> checkCRN 5 x (tanh (atanh x))+ , QC.testProperty "recip . recip" $ \x -> checkCRN 5 x (1/(1/x))+ , QC.testProperty "sqrt" $ QC.forAll (genCRClosed 0 100) $ checkFun sqrt+ , QC.testProperty "sin" $ checkFun sin+ , QC.testProperty "cos" $ checkFun cos+ , QC.testProperty "tan" $ checkFun tan+ , QC.testProperty "exp" $ QC.forAll (genCRClosed (-100) 100) $ checkFun exp+ , QC.testProperty "log" $ QC.forAll (genCROpen 0 100) $ checkFun log+ , QC.testProperty "asin" $ QC.forAll (genCRClosed (-1) 1) $ checkFun asin+ , QC.testProperty "acos" $ QC.forAll (genCRClosed (-1) 1) $ checkFun asin+ , QC.testProperty "atan" $ checkFun atan+ , QC.testProperty "sinh" $ QC.forAll (genCRClosed (-100) 100) $ checkFun sinh+ , QC.testProperty "asinh" $ checkFun asinh+ , QC.testProperty "cosh" $ QC.forAll (genCRClosed 0 10) $ checkFun cosh+ , QC.testProperty "acosh" $ QC.forAll (genCRClosed 1 10) $ checkFun acosh+ , QC.testProperty "tanh" $ QC.forAll (genCRClosed (-10) 10) $ checkFun tanh+ , QC.testProperty "atanh" $ QC.forAll (genCROpen (-1) 1) $ checkFun atanh+ ]++unitTests = testGroup "Unit tests"+ [ testCase "showA 1" $ showA (Approx 3 7 2 2) @?= "3~"+ , testCase "showA 2" $ showA (Approx 3 7 2 3) @?= "~~"+ , testCase "showA (Approx 99 2 0)" $ showA (Approx 8 99 2 0) @?= "10~"+ , testCase "sin π/2 = 1 (1000 bits)" $ approximatedBy 1 (require 1000 . sin $ pi/2) @?= True+ ]