diff --git a/AERN-Real.cabal b/AERN-Real.cabal new file mode 100644 --- /dev/null +++ b/AERN-Real.cabal @@ -0,0 +1,81 @@ +Name: AERN-Real +Version: 0.9.0 +Cabal-Version: >= 1.2 +Build-Type: Simple +License: BSD3 +License-File: LICENCE +Author: Michal Konecny +Copyright: (c) 2007-2008 Michal Konecny, Amin Farjudian, Jan Duracz +Maintainer: Michal Konecny +Stability: experimental +Category: Data, Math +Synopsis: datatypes and abstractions for approximating exact real numbers +Tested-with: GHC ==6.8.2 +Description: + Datatypes and abstractions for approximating exact real numbers + and a basic arithmetic over such approximations. The design is + inspired to some degree by Mueller's iRRAM and Lambov's RealLib + (both are C++ libraries for exact real arithmetic). + . + Abstractions are provided via 4 type classes: + . + * ERRealBase: abstracts floating point numbers + . + * ERApprox: abstracts neighbourhoods of real numbers + . + * ERIntApprox: abstracts neighbourhoods of real numbers that are known to be intervals + . + * ERApproxElementary: abstracts real number approximations that support elementary operations + . + For ERRealBase we give several implementations. The default is + an arbitrary precision floating point type that uses Double + for lower precisions and an Integer-based simulation for higher + precisions. Rational numbers can be used as one of the alternatives. + Augustsson's Data.Number.BigFloat can be easily wrapped as an instance + of ERRealBase except that it uses a different method to control precision. + . + ERIntApprox is implemented via outwards-rounded arbitrary precision interval arithmetic. + Any instance of ERRealBase can be used for the endpoints of the intervals. + . + ERApproxElementary is implemented generically for any implementation + of ERIntApprox. This way some of the most common elementary operations are provided, + notably: sqrt, exp, log, sin, cos, atan. These operations converge + to an arbitrary precision and also work well over larger intervals without + excessive wrapping. + . + There is also some support for generic Taylor series, interval Newton method + and simple numerical integration. + +Flag containers-in-base + +Library + hs-source-dirs: src + if flag(containers-in-base) + Build-Depends: + base < 3, binary >= 0.4 + else + Build-Depends: + base >= 3, containers, binary >= 0.4 + Exposed-modules: + Data.Number.ER, + Data.Number.ER.Real, + Data.Number.ER.Real.DefaultRepr, + Data.Number.ER.Real.Base.MachineDouble, + Data.Number.ER.Real.Base.CombinedMachineAP, + Data.Number.ER.Real.Base.Rational, + Data.Number.ER.Real.Base.Float, + Data.Number.ER.Real.Base, + Data.Number.ER.Real.Arithmetic.Elementary, + Data.Number.ER.Real.Arithmetic.Integration, + Data.Number.ER.Real.Arithmetic.Taylor, + Data.Number.ER.Real.Arithmetic.Newton, + Data.Number.ER.Real.Approx.Sequence, + Data.Number.ER.Real.Approx.Elementary, + Data.Number.ER.Real.Approx.Interval, + Data.Number.ER.Real.Approx, + Data.Number.ER.PlusMinus, + Data.Number.ER.BasicTypes, + Data.Number.ER.Misc, + Data.Number.ER.ExtendedInteger + Extensions: DeriveDataTypeable, ForeignFunctionInterface, ScopedTypeVariables + diff --git a/LICENCE b/LICENCE new file mode 100644 --- /dev/null +++ b/LICENCE @@ -0,0 +1,30 @@ +Copyright (c) 2007-2008 Michal Konecny, Amin Farjudian, Jan Duracz + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions +are met: + +1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + +2. 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. + +3. Neither the name of the author nor the names of his contributors + may be used to endorse or promote products derived from this software + without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE 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 AUTHORS 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. diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,3 @@ +#!/usr/bin/env runhaskell +> import Distribution.Simple +> main = defaultMain diff --git a/src/Data/Number/ER.hs b/src/Data/Number/ER.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER.hs @@ -0,0 +1,25 @@ +{-| + Module : Data.Number.ER + Description : top level of the exactreals framework + Copyright : (c) Michal Konecny + License : LGPL + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : non-portable (requires fenv.h) + + This namespace is the root for the AERN family of packages. + AERN stands for Approximated Exact Real Numbers. + All AERN packages build on the package AERN-Real. + + Module "Data.Number.ER.Real" contains an overview + of the AERN-Real package. + +-} +module Data.Number.ER +( + module Data.Number.ER.Real +) +where + +import Data.Number.ER.Real diff --git a/src/Data/Number/ER/BasicTypes.hs b/src/Data/Number/ER/BasicTypes.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/BasicTypes.hs @@ -0,0 +1,92 @@ +{-| + Module : Data.Number.ER.BasicTypes + Description : generic types for exact real number processing + Copyright : (c) Michal Konecny + License : LGPL + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : portable + + generic types for exact real number processing +-} +module Data.Number.ER.BasicTypes +where + +import qualified Data.Number.ER.ExtendedInteger as EI + +import qualified Data.Map as Map + +{-| + Precision represents an upper bound on the measure of + an approximation viewed as a set; + not to be confused with the precision of + an 'Data.Number.ER.Real.Base.Float.ERFloat' and similar. + + In an approximation comprising a number of + instances of 'Data.Number.ER.Real.Base.ERRealBase', + we will refer to the bit-precision of these base components + as the 'Granularity' of the approximation. +-} +type Precision = EI.ExtendedInteger + +{-| + The bit size of the floating point numbers (or similar) + used internally in real number and function approximations. +-} +type Granularity = Int + +prec2gran :: Precision -> Granularity +prec2gran = fromInteger . toInteger + +{-| + This type synonym should be used for funciton parameter(s) + that guide the convergence of the function's result to + a perfect (exact) result. + + The name should remind us + that there is no universally valid relationship between + this integer the quality (precision) of the result. + The only condition usually assumed is that in the limit + when the effort index rises to infinity, the result + should be exact. +-} +type EffortIndex = Integer + +effIx2gran :: EffortIndex -> Granularity +effIx2gran = fromInteger . toInteger + +effIx2prec :: EffortIndex -> Precision +effIx2prec = fromInteger . toInteger + +effIx2int :: EffortIndex -> Int +effIx2int = fromInteger . toInteger + +int2effIx :: Int -> EffortIndex +int2effIx = fromInteger . toInteger + +prec2effIx :: Precision -> EffortIndex +prec2effIx = fromInteger . toInteger + +gran2effIx :: Granularity -> EffortIndex +gran2effIx = fromInteger . toInteger + +{-| + A variable identifier for axes in function domains, polynomials etc. +-} +type VarID = Int +defaultVar :: VarID +defaultVar = 0 + +{-| + A many-dimensional point or interval. +-} +type Box ira = Map.Map VarID ira + +{-| using 'defaultVar' -} +unaryDom :: ira -> Box ira +unaryDom r = Map.singleton defaultVar r + +noinfoDom :: Box ira +noinfoDom = Map.empty + diff --git a/src/Data/Number/ER/ExtendedInteger.hs b/src/Data/Number/ER/ExtendedInteger.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/ExtendedInteger.hs @@ -0,0 +1,125 @@ +{-| + Module : Data.Number.ER.ExtendedInteger + Description : integer with infinities + Copyright : (c) Michal Konecny + License : LGPL + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : portable + + An arbitrary sized integer type with additional +infinity and -infinity. + + To be imported qualified, usually with prefix EI. +-} +module Data.Number.ER.ExtendedInteger +( + ExtendedInteger(..), + isInfinite, binaryLog, take +) +where + +import Prelude hiding (isInfinite, take) +import qualified Prelude + +data ExtendedInteger + = MinusInfinity | Finite Integer | PlusInfinity + deriving (Eq) + +isInfinite :: ExtendedInteger -> Bool +isInfinite MinusInfinity = True +isInfinite PlusInfinity = True +isInfinite _ = False + +{-| + the smallest integer i for which 2^i <= abs n +-} +binaryLog :: ExtendedInteger -> ExtendedInteger +binaryLog PlusInfinity = PlusInfinity +binaryLog MinusInfinity = PlusInfinity +binaryLog (Finite n) + | n < 0 = binaryLog (Finite (- n)) + | n == 0 = MinusInfinity + | otherwise = -- (n > 0) + -- how to do this fast? + intBinaryLog n + +intBinaryLog n + | n > 1 = 1 + (intBinaryLog (n `div` 2)) + | n == 1 = 0 + +instance Show ExtendedInteger where + show MinusInfinity = "-InfInt" + show PlusInfinity = "+InfInt" + show (Finite i) = show i + +take :: ExtendedInteger -> [a] -> [a] +take MinusInfinity _ = error "takeEI called with MinusInfinity" +take PlusInfinity list = list +take (Finite n) list = Prelude.take (fromInteger n) list + +instance Ord ExtendedInteger where + compare MinusInfinity MinusInfinity = EQ + compare MinusInfinity _ = LT + compare _ MinusInfinity = GT + compare PlusInfinity PlusInfinity = EQ + compare PlusInfinity _ = GT + compare _ PlusInfinity = LT + compare (Finite i1) (Finite i2) = + compare i1 i2 + +instance Num ExtendedInteger where + fromInteger i = Finite i + {- abs -} + abs MinusInfinity = PlusInfinity + abs PlusInfinity = PlusInfinity + abs (Finite i) = Finite $ abs i + {- signum -} + signum ei + | ei < 0 = -1 + | ei > 0 = 1 + | otherwise = 0 + {- negate -} + negate (Finite i) = Finite (-i) + negate MinusInfinity = PlusInfinity + negate PlusInfinity = MinusInfinity + {- addition -} + PlusInfinity + MinusInfinity = + error "cannot add PlusInfinity and MinusInfinity" + MinusInfinity + PlusInfinity = + error "cannot add PlusInfinity and MinusInfinity" + PlusInfinity + ei = PlusInfinity + ei + PlusInfinity = PlusInfinity + MinusInfinity + ei = MinusInfinity + ei + MinusInfinity = MinusInfinity + (Finite i1) + (Finite i2) = Finite $ i1 + i2 + {- multiplication -} + ei1 * ei2 | ei1 > ei2 = ei2 * ei1 + MinusInfinity * ei + | ei < 0 = PlusInfinity + | ei > 0 = MinusInfinity + | otherwise = error "cannot multiply MinusInfinity and 0" + ei * PlusInfinity + | ei < 0 = MinusInfinity + | ei > 0 = PlusInfinity + | otherwise = error "cannot multiply PlusInfinity and 0" + (Finite i1) * (Finite i2) = Finite $ i1 * i2 + +instance Enum ExtendedInteger where + toEnum i = Finite $ toInteger i + fromEnum (Finite i) = fromInteger i + fromEnum _ = error "infinite integers cannot be enumerated" + +instance Real ExtendedInteger where + toRational (Finite i) = toRational i + toRational _ = error "infinite integers cannot be converted to rational" + +instance Integral ExtendedInteger where + quotRem (Finite i) (Finite m) = + (Finite a, Finite b) + where + (a,b) = quotRem i m + quotRem _ _ = error "cannot make a quotient involving an infinite integer" + toInteger (Finite i) = i + toInteger _ = error "infinite integers cannot be converted to Integer" + diff --git a/src/Data/Number/ER/Misc.hs b/src/Data/Number/ER/Misc.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/Misc.hs @@ -0,0 +1,268 @@ +{-| + Module : Data.Number.ER.Misc + Description : general purpose extras + Copyright : (c) Michal Konecny + License : LGPL + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : portable + + Miscelaneous utilities (eg related to Ordering, pairs, booleans, strings) +-} +module Data.Number.ER.Misc where + +import List +import System.IO.Unsafe + +unsafePrint msg val = + unsafePerformIO $ + do + putStrLn $ "unsafe: " ++ msg + return val + +{-| + Compose as when defining the lexicographical ordering. +-} +compareCompose :: Ordering -> Ordering -> Ordering +compareCompose EQ o = o +compareCompose o _ = o + +{-| + Compose as when defining the lexicographical ordering. +-} +compareComposeMany :: [Ordering] -> Ordering +compareComposeMany [] = EQ +compareComposeMany (EQ:os) = compareComposeMany os +compareComposeMany (o:_) = o + +{-| + The lexicographical ordering. +-} +compareLex :: (Ord a) => [a] -> [a] -> Ordering +compareLex [] _ = LT +compareLex _ [] = GT +compareLex (x:xs) (y:ys) + | x == y = compareLex xs ys + | otherwise = compare x y + +mapFst :: (a1 -> a2) -> (a1,b) -> (a2,b) +mapFst f (a,b) = (f a,b) +mapSnd :: (b1 -> b2) -> (a,b1) -> (a,b2) +mapSnd f (a,b) = (a,f b) +mapPair :: (a1 -> a2, b1 -> b2) -> (a1,b1) -> (a2,b2) +mapPair (f1, f2) (a,b) = (f1 a, f2 b) +mapPairHomog :: (a1 -> a2) -> (a1,a1) -> (a2,a2) +mapPairHomog f = mapPair (f,f) + +unpair :: [(a,a)] -> [a] +unpair = (\(l1,l2) -> l1 ++ l2) . unzip + +bool2maybe :: Bool -> Maybe () +bool2maybe True = Just () +bool2maybe False = Nothing + +dropLast :: Int -> [a] -> [a] +dropLast n list = reverse $ drop n (reverse list) + +{-| + eg + +> concatWith "," ["a","b"] = "a,b" + +-} +concatWith :: + String {-^ a connective -} -> + [String] -> + String +concatWith sep [] = "" +concatWith sep [str] = str +concatWith sep (str : strs) = str ++ sep ++ (concatWith sep strs) + +{-| + eg + +> replicateSeveral [(2,"a"),(1,"b")] = "aab" + +-} +replicateSeveral :: [(Int,a)] -> [a] +replicateSeveral [] = [] +replicateSeveral ((n,e):rest) = + replicate n e ++ (replicateSeveral rest) + +{-| + eg + +> countDuplicates "aaba" = [(2,"a"),(1,"b"),(1,"a")] + +-} +countDuplicates :: + Eq a => + [a] -> + [(Int,a)] +countDuplicates list = + map (\ g -> (length g, head g)) $ group list + +{-| + eg + +> allCombinations +> [ +> (1,['a']), +> (2,['b','c']), +> (3,['d','e','f']) +> ] = +> [ +> [(1,'a'),(2,'b'),(3,'d')], +> [(1,'a'),(2,'b'),(3,'e')], +> [(1,'a'),(2,'b'),(3,'f')], +> [(1,'a'),(2,'c'),(3,'d')], +> [(1,'a'),(2,'c'),(3,'e')], +> [(1,'a'),(2,'c'),(3,'f')] +> ] +-} +allCombinations :: + [(k,[v])] -> [[(k,v)]] +allCombinations [] = [[]] +allCombinations ((k, vals) : rest) = + concat $ map (\ v -> map ((k,v):) restCombinations) vals + where + restCombinations = + allCombinations rest + +allPairsCombinations :: + [(k,(v,v))] -> [[(k,v)]] +allPairsCombinations [] = [[]] +allPairsCombinations ((k, (v1,v2)) : rest) = + (map ((k, v1) :) restCombinations) + ++ + (map ((k, v2) :) restCombinations) + where + restCombinations = + allPairsCombinations rest + + +{-| + eg + +> allPairsCombinationsEvenOdd +> [ +> (1,('a0','a1'), +> (2,('b0','b1'), +> (3,('c0','c1') +> ] = +> ([ +> [(1,'a0'),(2,'b0'),(3,'c0')], +> [(1,'a0'),(2,'b1'),(3,'c1')], +> [(1,'a1'),(2,'b1'),(3,'c0')], +> [(1,'a1'),(2,'b0'),(3,'c1')] +> ] +> ,[ +> [(1,'a0'),(2,'b0'),(3,'c1')], +> [(1,'a0'),(2,'b1'),(3,'c0')], +> [(1,'a1'),(2,'b0'),(3,'c0')], +> [(1,'a1'),(2,'b1'),(3,'c1')] +> ] +> ) +-} +allPairsCombinationsEvenOdd :: + [(k,(v,v))] {-^ the first value is even, the second odd -} -> + ([[(k,v)]], [[(k,v)]]) +allPairsCombinationsEvenOdd [] = ([[]], []) +allPairsCombinationsEvenOdd ((k, (evenVal,oddVal)) : rest) = + ( + (map ((k, evenVal) :) restCombinationsEven) + ++ + (map ((k, oddVal) :) restCombinationsOdd) + , + (map ((k, evenVal) :) restCombinationsOdd) + ++ + (map ((k, oddVal) :) restCombinationsEven) + ) + where + (restCombinationsEven, restCombinationsOdd) = + allPairsCombinationsEvenOdd rest + + + +{- numeric -} + +intLog :: + (Num n1, Num n2, Ord n1) => + n1 {-^ base -} -> + n1 {-^ x -} -> + n2 +intLog b n + | n > 0 = p2 + where + (p2, pe2) = findSlow (p1, pe1) (p1 + 1, pe1 * b) + (p1, pe1) = findFast (1, b) (2, b*b) + findFast (p, pe) (pp, ppe) + | ppe < n = findFast (pp, ppe) (2 * pp, ppe * ppe) + | otherwise = (p, pe) + findSlow (p, pe) (pp, ppe) + | ppe < n = findSlow (pp, ppe) (pp + 1, ppe * b) + | otherwise = (pp, ppe) + +{-| + Directionally rounded versions of @+,*,sum,prod@. +-} +plusUp, plusDown, timesUp, timesDown :: + (Num t) => + t -> t -> t +sumUp, sumDown, productDown, productUp :: + (Num t) => + [t] -> t +plusUp = (+) +plusDown c1 c2 = - ((- c1) - c2) +sumUp = foldl plusUp 0 +sumDown = foldl plusDown 0 +timesUp = (*) +timesDown c1 c2 = - ((- c1) * c2) +productUp = foldl timesUp 1 +productDown = foldl timesDown 1 + +{- parsing -} +readMaybe :: (Read a) => String -> Maybe a +readMaybe s = + case reads s of + [] -> Nothing + (val,_) : _ -> Just val + + +{- sequences -} +listUpdate :: Int -> a -> [a] -> [a] +listUpdate i newx (x:xs) + | i == 0 = newx : xs + | i > 0 = x : (listUpdate (i - 1) newx xs) + + +listHasMatch :: (a -> Bool) -> [a] -> Bool +listHasMatch f s = + foldl (\b a -> b && (f a)) False s + +--{-| types encoding natural numbers -} +--class TypeNumber n +-- where +-- getTNData :: n +-- getTNNumber :: n -> Int +-- +--data TN_0 = TN_0 +--tn_0 = TN_0 +--data TN_SUCC tn_prev = TN_SUCC tn_prev +-- +--type TN_ONE = TN_SUCC TN_0 +--tn_1 = TN_SUCC TN_0 +-- +--instance (TypeNumber TN_0) +-- where +-- getTNData = TN_0 +-- getTNNumber _ = 0 +-- +--instance +-- (TypeNumber tn_prev) => +-- (TypeNumber (TN_SUCC tn_prev)) +-- where +-- getTNData = TN_SUCC getTNData +-- getTNNumber (TN_SUCC p) = 1 + (getTNNumber p) + diff --git a/src/Data/Number/ER/PlusMinus.hs b/src/Data/Number/ER/PlusMinus.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/PlusMinus.hs @@ -0,0 +1,47 @@ +{-# LANGUAGE DeriveDataTypeable #-} +{-| + Module : Data.Number.ER.PlusMinus + Description : mini sign datatype + Copyright : (c) Michal Konecny + License : LGPL + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : portable + + A mini enumeration to represent the sign of different numbers and approximations. +-} +module Data.Number.ER.PlusMinus where + +import Data.Typeable +import Data.Generics.Basics +import Data.Binary +--import BinaryDerive + +data PlusMinus = Minus | Plus + deriving (Eq, Ord, Typeable, Data) + +instance Show PlusMinus where + show Plus = "+" + show Minus = "-" + +{- the following has been generated by BinaryDerive -} +instance Binary PlusMinus where + put Minus = putWord8 0 + put Plus = putWord8 1 + get = do + tag_ <- getWord8 + case tag_ of + 0 -> return Minus + 1 -> return Plus + _ -> fail "no parse" +{- the above has been generated by BinaryDerive -} + +signNeg Plus = Minus +signNeg Minus = Plus + +signMult Plus s = s +signMult Minus s = signNeg s + +signToNum Plus = 1 +signToNum Minus = -1 diff --git a/src/Data/Number/ER/Real.hs b/src/Data/Number/ER/Real.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/Real.hs @@ -0,0 +1,70 @@ +{-| + Module : Data.Number.ER.Real + Description : overview of AERN-Real + Copyright : (c) Michal Konecny + License : LGPL + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : non-portable (requires fenv.h) + + Datatypes and abstractions for approximating exact real numbers + and a basic arithmetic over such approximations. The design is + inspired to some degree by Mueller's iRRAM and Lambov's RealLib + (both are C++ libraries for exact real arithmetic). + + Abstractions are provided via 4 type classes: + + * 'B.ERRealBase': abstracts floating point numbers + + * 'RA.ERApprox': abstracts neighbourhoods of real numbers + + * 'RA.ERIntApprox': abstracts neighbourhoods of real numbers that are known to be intervals + + * 'RAEL.ERApproxElementary': abstracts real number approximations that support elementary operations + + For ERRealBase we give several implementations. The default is + an arbitrary precision floating point type that uses Double + for lower precisions and an Integer-based simulation for higher + precisions. Rational numbers can be used as one of the alternatives. + Augustsson's Data.Number.BigFloat can be easily wrapped as an instance + of ERRealBase except that it uses a different method to control precision. + + ERIntApprox is implemented via outwards-rounded arbitrary precision interval arithmetic. + Any instance of ERRealBase can be used for the endpoints of the intervals. + + ERApproxElementary is implemented generically for any implementation + of ERIntApprox. This way some of the most common elementary operations are provided, + notably: sqrt, exp, log, sin, cos, atan. These operations converge + to an arbitrary precision and also work well over larger intervals without + excessive wrapping. + + There is also some support for generic Taylor series, interval Newton method + and simple numerical integration. + +-} +module Data.Number.ER.Real +( + B.ERRealBase, + RA.ERApprox, + RA.ERIntApprox, + RAEL.ERApproxElementary, + module Data.Number.ER.Real.DefaultRepr, + module Data.Number.ER.Real.Approx.Sequence, + module Data.Number.ER.Real.Arithmetic.Taylor, + module Data.Number.ER.Real.Arithmetic.Newton, + module Data.Number.ER.Real.Arithmetic.Integration, + module Data.Number.ER.BasicTypes +) +where + +import Data.Number.ER.Real.DefaultRepr +import Data.Number.ER.BasicTypes +import qualified Data.Number.ER.Real.Base as B +import qualified Data.Number.ER.Real.Approx as RA +import qualified Data.Number.ER.Real.Approx.Elementary as RAEL +import Data.Number.ER.Real.Approx.Sequence +import Data.Number.ER.Real.Arithmetic.Taylor +import Data.Number.ER.Real.Arithmetic.Newton +import Data.Number.ER.Real.Arithmetic.Integration + diff --git a/src/Data/Number/ER/Real/Approx.hs b/src/Data/Number/ER/Real/Approx.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/Real/Approx.hs @@ -0,0 +1,303 @@ +{-| + Module : Data.Number.ER.Real.Approx + Description : classes abstracting exact reals + Copyright : (c) Michal Konecny + License : LGPL + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : portable + + Definitions of classes that describe what is + required from arbitrary precision approximations + of exact real numbers. + + We introduce two levels of abstraction for these + approximations: + + * 'ERApprox' = + a *set* of approximated numbers whose size is + measured using some fixed measure + + * 'ERIntApprox' = + an *interval* of real numbers with finitely + representable endpoints + + To be imported qualified, usually with the synonym RA. +-} +module Data.Number.ER.Real.Approx +( + ERApprox(..), + ERIntApprox(..), + bounds2ira, + effIx2ra, + splitIRA, +-- checkShrinking, +-- eqSingletons, + exactMiddle, + maxExtensionR2R +) +where + +import Data.Number.ER.BasicTypes +import qualified Data.Number.ER.ExtendedInteger as EI + +import Data.Typeable + +{-| + A type whose elements represent sets that can be used + to approximate a single extended real number with arbitrary precision. +-} +class (Fractional ra, Ord ra) => ERApprox ra where + getPrecision :: ra -> Precision + {-^ + Precision is a measure of the set size. + + The default interpretation: + + * If the diameter of the set is d, then the precision + should be near floor(- log_2 d). + -} + getGranularity :: ra -> Granularity + -- ^ the lower the granularity the bigger the rounding errors + setGranularity :: Granularity -> ra -> ra + -- ^ increase or safely decrease granularity + setMinGranularity :: Granularity -> ra -> ra + -- ^ ensure granularity is not below the first arg + isEmpty :: ra -> Bool + -- ^ true if this represents a computational error + isBottom :: ra -> Bool + -- ^ true if this holds no information + isExact :: ra -> Bool + -- ^ true if this is a singleton + isDisjoint :: ra -> ra -> Bool + isDisjoint a b = isEmpty $ a /\ b + isBounded :: ra -> Bool + -- ^ true if the approximation excludes infinity + bottomApprox :: ra + -- ^ the bottom element - any number + emptyApprox :: ra + -- ^ the top element - error + refines :: ra -> ra -> Bool + -- ^ first arg is a subset of the second arg + (/\) :: ra -> ra -> ra + -- ^ join; combining two approximations of the same number + intersectMeasureImprovement :: + EffortIndex -> ra -> ra -> (ra, ra) + {-^ + Like intersection but the second component: + + * measures improvement of the intersection relative to the first of the two approximations + + * is a positive number: 1 means no improvement, 2 means doubled precision, etc. + -} + equalReals :: ra -> ra -> Maybe Bool + -- ^ nothing if overlapping and not singletons + compareReals :: ra -> ra -> Maybe Ordering + -- ^ nothing if overlapping and not singletons + leqReals :: ra -> ra -> Maybe Bool + -- ^ nothing if overlapping on interior or by a wrong endpoint + equalApprox :: ra -> ra -> Bool + -- ^ syntactic comparison + compareApprox :: ra -> ra -> Ordering + -- ^ syntactic linear ordering + double2ra :: Double -> ra + showApprox :: + Int {-^ number of relevant decimals to show -} -> + Bool {-^ should show granularity -} -> + Bool {-^ should show internal representation details -} -> + ra {-^ the approximation to show -} -> + String + +{-| + Assuming the arguments are singletons, equality is decidable. +-} +eqSingletons :: (ERApprox ra) => ra -> ra -> Bool +eqSingletons s1 s2 = + case equalReals s1 s2 of + Just b -> b + _ -> False + +{-| + Assuming the arguments are singletons, @<=@ is decidable. +-} +leqSingletons :: (ERApprox ra) => ra -> ra -> Bool +leqSingletons s1 s2 = + case compareReals s1 s2 of + Just EQ -> True + Just LT -> True + _ -> False + +{-| + Assuming the arguments are singletons, @<@ is decidable. +-} +ltSingletons :: (ERApprox ra) => ra -> ra -> Bool +ltSingletons s1 s2 = + case compareReals s1 s2 of + Just LT -> True + _ -> False + +{-| + For a finite sequence of real approximations, determine + whether it is a shrinking sequence. +-} +checkShrinking :: + (ERApprox ra) => + [ra] -> + Maybe (ra, ra) +checkShrinking [] = Nothing +checkShrinking [_] = Nothing +checkShrinking (a : b : rest) + | b `refines` a = checkShrinking (b : rest) + | otherwise = Just (a,b) + + +{-| + A type whose elements represent sets that can be used + to approximate a recursive set of closed extended real number intervals + with arbitrary precision. +-} +--class (ERApprox sra) => SetOfRealsApprox sra where +-- (\/) :: sra -> sra -> sra -- ^ union; either approximation could be correct + +{-| + A type whose elements represent real *intervals* that can be used + to approximate a single extended real number with arbitrary precision. + + Sometimes, these types can be used to approximate + a closed extended real number interval with arbitrary precision. + Nevetheless, this is not guaranteed. +-} +class (ERApprox ira) => ERIntApprox ira + where + doubleBounds :: ira -> (Double, Double) + floatBounds :: ira -> (Float, Float) + integerBounds :: ira -> (EI.ExtendedInteger, EI.ExtendedInteger) + bisectDomain :: + Maybe ira {-^ point to split at -} -> + ira {-^ interval to split -} -> + (ira, ira) -- ^ left and right, overlapping on a singleton + defaultBisectPt :: ira -> ira + -- | returns thin approximations of endpoints, in natural order + bounds :: ira -> (ira, ira) + {-| + meet, usually constructing interval from approximations of its endpoints + + This does not need to be the meet of the real intervals + but it has to be a maximal element in the set of all + ira elements that are below the two parameters. + -} + (\/) :: ira -> ira -> ira + +bounds2ira :: + (ERIntApprox ira) => + ira -> ira -> ira +bounds2ira = (\/) + +{- old stuff that will probably never be resurrected: + +-- It is intended that ra and ira are the same type. +-- We distinguish them so that we can conveniently +-- switch between two levels of abstraction when +-- working with values of this one type. +-- +-- Given some ra or ira, the other type is determined uniquely. + +-- -- | coercion to more concrete view (allows a more intentional computation) +-- ra2ira :: ra -> ira +-- -- | coercion to more abstract view (guarantees certain extensionality and convergence properties) +-- ira2ra :: ira -> ra + +-- -- | coercion +-- ira2sra :: ira -> sra +-- sraCover :: sra -> ira +-- sraAllIntervals :: sra -> [ira] -- ^ disjoint, in natural order +-} + +-- +--bounds2ira :: +-- (ERIntApprox ira) => +-- ra -> +-- ra -> +-- ira +--bounds2ira leftRA rightRA = +-- (ra2ira leftRA) \/ (ra2ira rightRA) + +effIx2ra :: + (ERApprox ra) => + EffortIndex -> ra +effIx2ra = fromInteger . toInteger + +{-| + Split an interval to a sequence of intervals whose union is the + original interval using a given sequence of cut points. + The cut points are expected to be in increasing order and contained + in the given interval. Violations of this rule are tolerated. +-} +splitIRA :: + (ERIntApprox ira) => + ira {-^ an interval to be split -} -> + [ira] {-^ approximations of the cut points in increasing order -} -> + [ira] +splitIRA interval splitPoints = + doSplit [] end pointsRev + where + (start, end) = bounds interval + pointsRev = reverse $ start : splitPoints + doSplit previousSegments nextRight [] = previousSegments + doSplit previousSegments nextRight (nextLeft : otherPoints) = + doSplit (nextLeft \/ nextRight : previousSegments) nextLeft otherPoints + +{-| + * Return the endpoints of the interval as well as the exact midpoint. + + * To be able to do this, there may be a need to increase granularity. + + * All three singleton intervals are set to the same new granularity. +-} +exactMiddle :: + (ERIntApprox ira) => + ira -> + (ira,ira,ira,Granularity) +exactMiddle dom = + case isExact domM of + True -> + (domL, domM, domR, gran) + False -> + (domLhg, domMhg, domRhg, higherGran) + where + (domL, domR) = bounds dom + gran = max (getGranularity domL) (getGranularity domR) + domM = (domL + domR) / 2 + higherGran = gran + 1 + domLhg = setMinGranularity higherGran domL + domRhg = setMinGranularity higherGran domR + domMhg = (domLhg + domRhg) / 2 + + +{-| + This produces a function that computes the maximal extension of the + given function. A maximal extension function has the property: + f(I) = { f(x) | x in I }. Here we get this property only for the + limit function for ix tending to infinity. +-} +maxExtensionR2R :: + (ERIntApprox ira) => + (EffortIndex -> ira -> [ira]) + {-^ returns a safe approximation of all extrema within the interval -} -> + (EffortIndex -> ira -> ira) + {-^ a function behaving well on sequences that intersect to a point -} -> + (EffortIndex -> ira -> ira) + {- ^ a function behaving well on sequences that intersect to a non-empty interval -} +maxExtensionR2R getExtremes f ix x + | getPrecision x < effIx2prec ix = + (f ix xL) \/ (f ix xR) \/ + (foldl (\/) emptyApprox $ getExtremes ix x) + -- x is thin enough (?), don't bother evaluating by endpoints and extrema: + | otherwise = + f ix x + where + (xL, xR) = bounds x + + + diff --git a/src/Data/Number/ER/Real/Approx/Elementary.hs b/src/Data/Number/ER/Real/Approx/Elementary.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/Real/Approx/Elementary.hs @@ -0,0 +1,60 @@ +{-| + Module : Data.Number.ER.Real.Approx.Elementary + Description : abstraction of exact reals capable of elementary operations + Copyright : (c) Michal Konecny + License : LGPL + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : portable + + To be imported qualified, usually with the synonym RAEL. +-} +module Data.Number.ER.Real.Approx.Elementary +( + ERApproxElementary(..) +) +where + +import Prelude hiding (exp, log, sin, cos) + +import qualified Data.Number.ER.Real.Approx as RA +import Data.Number.ER.BasicTypes + +import Data.Number.ER.Real.Arithmetic.Elementary + +{-| + A class defining various common real number operations + in a approximation-aware fashion, ie introducing effort indices. + + All operations here have default implementations based on + "Data.Number.ER.Real.Arithmetic.Elementary". +-} +class (RA.ERIntApprox ra) => (ERApproxElementary ra) + where + abs :: EffortIndex -> ra -> ra + abs ix = Prelude.abs + min :: EffortIndex -> ra -> ra -> ra + min ix = Prelude.min + max :: EffortIndex -> ra -> ra -> ra + max ix = Prelude.max + exp :: EffortIndex -> ra -> ra + exp = erExp_IR + log :: EffortIndex -> ra -> ra + log = erLog_IR + (**) :: EffortIndex -> ra -> ra -> ra + (**) ix b e = exp ix $ e * (log ix b) + pi :: EffortIndex -> ra + pi = erPi_R + sin :: EffortIndex -> ra -> ra + sin = erSine_IR + cos :: EffortIndex -> ra -> ra + cos = erCosine_IR + tan :: EffortIndex -> ra -> ra + tan ix r = (sin ix r) / (cos ix r) + atan :: EffortIndex -> ra -> ra + atan = erATan_IR + + + + diff --git a/src/Data/Number/ER/Real/Approx/Interval.hs b/src/Data/Number/ER/Real/Approx/Interval.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/Real/Approx/Interval.hs @@ -0,0 +1,531 @@ +{-# LANGUAGE DeriveDataTypeable #-} +{-| + Module : Data.Number.ER.Real.Approx.Interval + Description : safe interval arithmetic + Copyright : (c) Michal Konecny + License : LGPL + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : portable + + This module defines an arbitrary precision interval type and + most of its interval arithmetic operations. +-} +module Data.Number.ER.Real.Approx.Interval +( + ERInterval(..), + normaliseERInterval +) +where + +import qualified Data.Number.ER.Real.Approx as RA +import qualified Data.Number.ER.Real.Approx.Elementary as RAEL +import qualified Data.Number.ER.Real.Base as B +import qualified Data.Number.ER.ExtendedInteger as EI + +import Data.Number.ER.BasicTypes +import Data.Number.ER.Misc + +import Data.Ratio + +import Data.Typeable +import Data.Generics.Basics +import Data.Binary +--import BinaryDerive + +{-| + Type for arbitrary precision interval arithmetic. +-} +data ERInterval base = + ERIntervalEmpty -- ^ usually represents computation error (top element in the interval domain) + | ERIntervalAny -- ^ represents no knowledge of result (bottom element in the interval domain) + | ERInterval + { + erintv_left :: base, + erintv_right :: base + } + deriving (Typeable, Data) + +{- the following has been generated by BinaryDerive -} +instance (Binary a) => Binary (ERInterval a) where + put ERIntervalEmpty = putWord8 0 + put ERIntervalAny = putWord8 1 + put (ERInterval a b) = putWord8 2 >> put a >> put b + get = do + tag_ <- getWord8 + case tag_ of + 0 -> return ERIntervalEmpty + 1 -> return ERIntervalAny + 2 -> get >>= \a -> get >>= \b -> return (ERInterval a b) + _ -> fail "no parse" +{- the above has been generated by BinaryDerive -} + + +{-| + convert to a normal form, ie: + + * no NaNs as endpoints + + * @l <= r@ + + * no (-Infty, +Infty) +-} +normaliseERInterval :: + (B.ERRealBase b) => + ERInterval b -> ERInterval b +normaliseERInterval (ERInterval minusInfty plusInfty) + | B.isPlusInfinity plusInfty && B.isPlusInfinity (- minusInfty) = + ERIntervalAny +normaliseERInterval (ERInterval nan1 nan2) + | B.isERNaN nan1 && B.isERNaN nan2 = + ERIntervalAny +normaliseERInterval (ERInterval nan r) + | B.isERNaN nan = + ERInterval (- B.plusInfinity) r +normaliseERInterval (ERInterval l nan) + | B.isERNaN nan = + ERInterval l (B.plusInfinity) +normaliseERInterval (ERInterval l r) + | l > r = ERIntervalEmpty +normaliseERInterval i = i + +{-| + erintvPrecision returns an approximation of the number of bits required + to represent the mantissa of a normalised size of the interval: + + + > - log_2 ((r - l) / (1 + abs(r) + abs(l))) + + Notice that this is +Infty for singleton and empty intervals + and -Infty for the whole real line. +-} +erintvPrecision :: + (B.ERRealBase b) => + ERInterval b -> EI.ExtendedInteger +erintvPrecision (ERInterval l r) = + - (B.getApproxBinaryLog $ (r - l)/(1 + abs(r) + abs(l))) +erintvPrecision ERIntervalEmpty = EI.PlusInfinity +erintvPrecision ERIntervalAny = EI.MinusInfinity + +erintvGranularity :: + (B.ERRealBase b) => + ERInterval b -> Int +erintvGranularity ERIntervalAny = 0 +erintvGranularity ERIntervalEmpty = 0 +erintvGranularity (ERInterval l r) = + min (B.getGranularity l) (B.getGranularity r) + +{- syntactic comparisons -} + +{-| + a syntactic equality test +-} +erintvEqualApprox :: + (B.ERRealBase b) => + ERInterval b -> ERInterval b -> Bool +erintvEqualApprox (ERInterval l1 r1) (ERInterval l2 r2) = + l1 == l2 && r1 == r2 +erintvEqualApprox ERIntervalEmpty ERIntervalEmpty = True +erintvEqualApprox ERIntervalAny ERIntervalAny = True +erintvEqualApprox _ _ = False + +{-| + a syntactic linear order +-} +erintvCompareApprox :: + (B.ERRealBase b) => + ERInterval b -> ERInterval b -> Ordering +erintvCompareApprox ERIntervalEmpty ERIntervalEmpty = EQ +erintvCompareApprox ERIntervalEmpty _ = LT +erintvCompareApprox _ ERIntervalEmpty = GT +erintvCompareApprox ERIntervalAny ERIntervalAny = EQ +erintvCompareApprox ERIntervalAny _ = LT +erintvCompareApprox _ ERIntervalAny = GT +erintvCompareApprox (ERInterval l1 r1) (ERInterval l2 r2) = + case compare l1 l2 of + EQ -> compare r1 r2 + res -> res + +{- semantic comparisons -} + +{-| + Compare for equality two intervals interpreted as approximations for + 2 single real numbers. When equality or inequality cannot + be established, return Nothing (ie bottom). +-} +erintvEqualReals :: + (B.ERRealBase b) => + ERInterval b -> + ERInterval b -> + Maybe Bool +erintvEqualReals ERIntervalEmpty _ = Nothing +erintvEqualReals _ ERIntervalEmpty = Nothing +erintvEqualReals ERIntervalAny _ = Nothing +erintvEqualReals _ ERIntervalAny = Nothing +erintvEqualReals (ERInterval l1 r1) (ERInterval l2 r2) + | l1 == r1 && l2 == r2 && l1 == l2 = Just True + | r1 < l2 || l1 > r2 = Just False + | otherwise = Nothing + +{-| + Compare in natural order two intervals interpreted as approximations for + 2 single real numbers. When equality or inequality cannot + be established, return Nothing (ie bottom). +-} +erintvCompareReals :: + (B.ERRealBase b) => + ERInterval b -> + ERInterval b -> + Maybe Ordering +erintvCompareReals ERIntervalEmpty _ = Nothing +erintvCompareReals _ ERIntervalEmpty = Nothing +erintvCompareReals ERIntervalAny _ = Nothing +erintvCompareReals _ ERIntervalAny = Nothing +erintvCompareReals i1@(ERInterval l1 r1) i2@(ERInterval l2 r2) + | r1 < l2 = Just LT + | l1 > r2 = Just GT + | l1 == r1 && l2 == r2 && l1 == l2 = Just EQ + | otherwise = Nothing + +{-| + Compare in natural order two intervals interpreted as approximations for + 2 single real numbers. When relaxed equality cannot + be established nor disproved, return Nothing (ie bottom). +-} +erintvLeqReals :: + (B.ERRealBase b) => + ERInterval b -> + ERInterval b -> + Maybe Bool +erintvLeqReals ERIntervalEmpty _ = Nothing +erintvLeqReals _ ERIntervalEmpty = Nothing +erintvLeqReals ERIntervalAny _ = Nothing +erintvLeqReals _ ERIntervalAny = Nothing +erintvLeqReals i1@(ERInterval l1 r1) i2@(ERInterval l2 r2) + | r1 <= l2 = Just True + | l1 > r2 = Just False + | otherwise = Nothing + + +{-| + + Default splitting: + + > [-Infty,+Infty] |-> [-Infty,0] [0,+Infty] + + > [-Infty,x] |-> [-Infty,2*x-1] [2*x-1, x] (x <= 0) + + > [-Infty,x] |-> [-Infty,0] [0, x] (x > 0) + + > [x,+Infty] |-> [x,2*x+1] [2*x+1,+Infty] (x => 0) + + > [x,+Infty] |-> [x,0] [0,+Infty] (x < 0) + + > [x,y] |-> [x, (x+y)/2] [(x+y)/2, y] + + > empty |-> empty empty +-} +erintvDefaultBisectPt :: + (B.ERRealBase b) => + Granularity -> + (ERInterval b) -> + (ERInterval b) +erintvDefaultBisectPt gran ERIntervalAny = 0 +erintvDefaultBisectPt gran ERIntervalEmpty = ERIntervalEmpty +erintvDefaultBisectPt gran (ERInterval l r) = + ERInterval m m + where + m + | B.isPlusInfinity r = + if l < 0 + then 0 + else 2 * (B.setMinGranularity gran l) + 1 + | B.isPlusInfinity (-l) = + if r > 0 + then 0 + else 2 * (B.setMinGranularity gran r) - 1 + | otherwise = + ((B.setMinGranularity gran l) + r)/2 + + +erintvBisect :: + (B.ERRealBase b, RealFrac b) => + Granularity -> + (Maybe (ERInterval b)) -> + (ERInterval b) -> + (ERInterval b, ERInterval b) +erintvBisect gran maybePt i = + (l RA.\/ m, m RA.\/ r) + where + (l,r) = RA.bounds i + m = + case maybePt of + Just m -> m + Nothing -> erintvDefaultBisectPt gran i + +instance (B.ERRealBase b) => Eq (ERInterval b) where + i1 == i2 = + case erintvEqualReals i1 i2 of + Nothing -> + error $ + "ERInterval: Eq: comparing overlapping intervals:\n" ++ + show i1 ++ "\n" ++ + show i2 + Just b -> b + +instance (B.ERRealBase b) => Ord (ERInterval b) where + compare i1 i2 = + case erintvCompareReals i1 i2 of + Nothing -> + error $ + "ERInterval: Ord: comparing overlapping intervals:\n" ++ + show i1 ++ "\n" ++ + show i2 + Just r -> r + {- max: + (Default implementation is wrong in this case: + eg compare is not defined for overlapping intervals.) + -} + max i1@(ERInterval l1 r1) i2@(ERInterval l2 r2) = + normaliseERInterval $ ERInterval (max l1 l2) (max r1 r2) + max ERIntervalEmpty _ = ERIntervalEmpty + max _ ERIntervalEmpty = ERIntervalEmpty + max ERIntervalAny ERIntervalAny = ERIntervalAny + max ERIntervalAny (ERInterval l r) = ERInterval l B.plusInfinity + max (ERInterval l r) ERIntervalAny = ERInterval l B.plusInfinity + {- min: -} + min i1@(ERInterval l1 r1) i2@(ERInterval l2 r2) = + normaliseERInterval $ ERInterval (min l1 l2) (min r1 r2) + min ERIntervalEmpty _ = ERIntervalEmpty + min _ ERIntervalEmpty = ERIntervalEmpty + min ERIntervalAny ERIntervalAny = ERIntervalAny + min ERIntervalAny (ERInterval l r) = ERInterval (- B.plusInfinity) r + min (ERInterval l r) ERIntervalAny = ERInterval (- B.plusInfinity) r + +instance (B.ERRealBase b) => Show (ERInterval b) + where + show = erintvShow 6 True False + +erintvShow numDigits showGran showComponents interval = + showERI interval + where + showERI ERIntervalEmpty = "[NONE]" + showERI ERIntervalAny = "[ANY]" + showERI (ERInterval l r) + | l == r = "<" ++ showBase l ++ ">" + | otherwise = + "[" ++ showBase l ++ "," ++ showBase r ++ "]" + showBase = B.showDiGrCmp numDigits showGran showComponents + +instance (B.ERRealBase b) => Num (ERInterval b) where + fromInteger n = + normaliseERInterval $ ERInterval (fromInteger n) (fromInteger n) + {- abs -} + abs (ERInterval l r) + | l < 0 && r > 0 = ERInterval 0 (max (-l) r) + | r <= 0 = ERInterval (-r) (-l) + | otherwise = ERInterval l r + abs ERIntervalAny = ERInterval 0 B.plusInfinity + abs ERIntervalEmpty = ERIntervalEmpty + {- signum -} + signum i@(ERInterval l r) + | l < 0 && r > 0 = ERInterval (-1) 1 -- need many-valuedness via sequences of intervals + | r < 0 = ERInterval (-1) (-1) + | l > 0 = ERInterval 1 1 + | l == 0 && r == 0 = i + | l == 0 = ERInterval 0 1 + | r == 0 = ERInterval (-1) 0 + signum ERIntervalAny = ERInterval (-1) 1 + signum ERIntervalEmpty = ERIntervalEmpty + {- negate -} + negate (ERInterval l r) = (ERInterval (-r) (-l)) + negate ERIntervalEmpty = ERIntervalEmpty + negate ERIntervalAny = ERIntervalAny + {- addition -} + (ERInterval l1 r1) + (ERInterval l2 r2) = + normaliseERInterval $ + ERInterval + (-((-l1) + (-l2))) -- reverse the rounding mode + (r1 + r2) + ERIntervalAny + i2 = ERIntervalAny + i1 + ERIntervalAny = ERIntervalAny + ERIntervalEmpty + i2 = ERIntervalEmpty + i1 + ERIntervalEmpty = ERIntervalEmpty + {- multiplication -} + (ERInterval l1 r1) * (ERInterval l2 r2) + | haveNan = ERIntervalAny + | otherwise = + normaliseERInterval $ + ERInterval minProd maxProd + where + haveNan = or $ map B.isERNaN (prodsL ++ prodsR) + minProd = foldl1 min prodsL + maxProd = foldl1 max prodsR + prodsL = [-((-l1) * l2), -((-l1) * r2), -((-r1) * l2), -((-r1) * r2)] + prodsR = [l1 * l2, l1 * r2, r1 * l2, r1 * r2] + ERIntervalAny * i2 = ERIntervalAny + i1 * ERIntervalAny = ERIntervalAny + ERIntervalEmpty * i2 = ERIntervalEmpty + i1 * ERIntervalEmpty = ERIntervalEmpty + +instance (B.ERRealBase b) => Fractional (ERInterval b) where + fromRational rat = + (fromInteger $ numerator rat) + / (fromInteger $ denominator rat) + {- division -} + (ERInterval l1 r1) / (ERInterval l2 r2) + | l2 < 0 && r2 > 0 = ERIntervalAny + | haveNan = +-- unsafePrint "ERInterval: /: haveNan" $ + ERIntervalAny + | l2 == 0 && r2 > 0 && 1/l2 < 0 = -- minus 0 + (ERInterval l1 r1) / (ERInterval (-l2) r2) -- correct it to +0 + | r2 == 0 && l2 < 0 && 1/r2 > 0 = -- plus 0 + (ERInterval l1 r1) / (ERInterval l2 (-r2)) -- correct it to -0 + | otherwise = + normaliseERInterval $ + ERInterval minDiv maxDiv + where + haveNan = or $ map B.isERNaN (divsL ++ divsR) + minDiv = foldl1 min divsL + maxDiv = foldl1 max divsR + divsL = [-(l1 / (-l2)), -(l1 / (-r2)), -(r1 / (-l2)), -(r1 / (-r2))] + divsR = [l1 / l2, l1 / r2, r1 / l2, r1 / r2] + ERIntervalAny / i2 = ERIntervalAny + i1 / ERIntervalAny = ERIntervalAny + ERIntervalEmpty / i2 = ERIntervalEmpty + i1 / ERIntervalEmpty = ERIntervalEmpty + +instance (B.ERRealBase b, RealFrac b) => RA.ERApprox (ERInterval b) where + getPrecision i = erintvPrecision i + getGranularity i = erintvGranularity i + {- setMinGranularity -} + setMinGranularity gr (ERInterval l r) = + normaliseERInterval $ + (ERInterval (- (B.setMinGranularity gr (-l))) (B.setMinGranularity gr r)) + setMinGranularity _ i = i + {- setGranularity -} + setGranularity gr (ERInterval l r) = + normaliseERInterval $ + (ERInterval (- (B.setGranularity gr (-l))) (B.setGranularity gr r)) + setGranularity _ i = i + {- isDisjoint -} + isDisjoint i1 i2 = RA.isEmpty $ i1 RA./\ i2 + {- bottomApprox -} + bottomApprox = ERIntervalAny + {- emptyApprox -} + emptyApprox = ERIntervalEmpty + {- isEmpty -} + isEmpty ERIntervalEmpty = True + isEmpty _ = False + {- isBottom -} + isBottom ERIntervalAny = True + isBottom (ERInterval l r) = + B.isPlusInfinity r && B.isPlusInfinity (-l) + isBottom _ = False + {- isExact -} + isExact ERIntervalEmpty = False + isExact ERIntervalAny = False + isExact (ERInterval l r) = l == r + {- isBounded -} + isBounded ERIntervalEmpty = True + isBounded ERIntervalAny = False + isBounded (ERInterval l r) = + (- B.plusInfinity) < l && r < B.plusInfinity + {- intersection -} + ERIntervalEmpty /\ i = ERIntervalEmpty + i /\ ERIntervalEmpty = ERIntervalEmpty + ERIntervalAny /\ i = i + i /\ ERIntervalAny = i + (ERInterval l1 r1) /\ (ERInterval l2 r2) = + normaliseERInterval $ + ERInterval (max l1 l2) (min r1 r2) + {- intersectMeasureImprovement -} + intersectMeasureImprovement _ ERIntervalEmpty i = (ERIntervalEmpty, 1) + intersectMeasureImprovement _ i ERIntervalEmpty = (ERIntervalEmpty, 1) + intersectMeasureImprovement _ ERIntervalAny i = (i, 1) + intersectMeasureImprovement _ i ERIntervalAny = (i, 1) + intersectMeasureImprovement ix i1 i2 = + (isec, impr) + where + isec = i1 RA./\ i2 + impr + | 0 `RA.refines` isecWidth && 0 `RA.refines` i1Width = 1 -- 0 -> 0 is no improvement + | otherwise = i1Width / isecWidth + i1Width = i1H - i1L + isecWidth = isecH - isecL + (isecL, isecH) = RA.bounds $ RA.setMinGranularity gran isec + (i1L, i1H) = RA.bounds $ RA.setMinGranularity gran i1 + gran = effIx2gran ix + {- refines -} + refines _ ERIntervalAny = True + refines ERIntervalEmpty _ = True + refines ERIntervalAny _ = False + refines _ ERIntervalEmpty = False + refines (ERInterval l1 r1) (ERInterval l2 r2) = + l2 <= l1 && r1 <= r2 + {- semantic comparisons -} + equalReals = erintvEqualReals + compareReals = erintvCompareReals + leqReals = erintvLeqReals + {- non-semantic comparisons -} + equalApprox = erintvEqualApprox + compareApprox = erintvCompareApprox + {- conversion from Double -} + double2ra d = + ERInterval b b + where + b = B.fromDouble d + {- formatting -} + showApprox = erintvShow + +instance (B.ERRealBase b, RealFrac b) => RA.ERIntApprox (ERInterval b) + where + doubleBounds ERIntervalAny = (- infinity, infinity) + where + infinity = 1/0 + doubleBounds ERIntervalEmpty = + error "SuiteERInterval: iraDoubleBounds: empty interval" + doubleBounds (ERInterval l r) = + (B.toDouble l, B.toDouble r) + floatBounds ERIntervalAny = (- infinity, infinity) + where + infinity = 1/0 + floatBounds ERIntervalEmpty = + error "SuiteERInterval: iraFloatBounds: empty interval" + floatBounds (ERInterval l r) = + (B.toFloat l, B.toFloat r) + integerBounds ERIntervalAny = (- infinity, infinity) + where + infinity = EI.PlusInfinity + integerBounds ERIntervalEmpty = + error "SuiteERInterval: iraIntegerBounds: empty interval" + integerBounds (ERInterval l r) = + (- (mkEI (- l)), mkEI r) + where + mkEI f + | B.isPlusInfinity f = EI.PlusInfinity + | B.isPlusInfinity (-f) = EI.MinusInfinity + | otherwise = ceiling f + defaultBisectPt dom = erintvDefaultBisectPt (RA.getGranularity dom + 1) dom + bisectDomain maybePt dom = + erintvBisect (RA.getGranularity dom + 1) maybePt dom + {- \/ -} + ERIntervalEmpty \/ i = i + i \/ ERIntervalEmpty = i + ERIntervalAny \/ _ = ERIntervalAny + _ \/ ERIntervalAny = ERIntervalAny + (ERInterval l1 r1) \/ (ERInterval l2 r2) = + normaliseERInterval $ + ERInterval (min l1 l2) (max r1 r2) + {- RA.bounds -} + bounds ERIntervalAny = + (ERInterval (-B.plusInfinity) (-B.plusInfinity), + ERInterval B.plusInfinity B.plusInfinity) + bounds ERIntervalEmpty = (ERIntervalEmpty, ERIntervalEmpty) + bounds (ERInterval l r) = + (ERInterval l l, ERInterval r r) + +instance (B.ERRealBase b, RealFrac b) => RAEL.ERApproxElementary (ERInterval b) +-- all operations here have appropriate default implementations diff --git a/src/Data/Number/ER/Real/Approx/Sequence.hs b/src/Data/Number/ER/Real/Approx/Sequence.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/Real/Approx/Sequence.hs @@ -0,0 +1,220 @@ +{-| + Module : Data.Number.ER.Real.Approx.Sequence + Description : exact reals via convergent sequences + Copyright : (c) Michal Konecny + License : LGPL + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : portable + + Types and methods related to explicit + convergent sequences of real number approximations. +-} +module Data.Number.ER.Real.Approx.Sequence +( + ConvergRealSeq, + makeFastConvergRealSeq, + convertFuncRA2Seq, + convertBinFuncRA2Seq, + convergRealSeqElem, + showConvergRealSeq, + showConvergRealSeqAuto +) +where + +import qualified Data.Number.ER.Real.Approx as RA +import Data.Number.ER.BasicTypes + +import Data.Maybe +import Data.Ratio + +{-| + A converging sequence of real number approximations. + + * Every finite subsequence has a non-empty intersection. + + * The limit should be a singleton. +-} +data ConvergRealSeq ra = + ConvergRealSeq (EffortIndex -> ra) + +convergRealSeqElem :: (ConvergRealSeq ra) -> EffortIndex -> ra +convergRealSeqElem (ConvergRealSeq sq) ix = sq ix + +{-| + Using this operator, a unary funtion working over + approximations can be converted to one that works + over exact numbers represented through a sequence + of approximations. +-} +convertFuncRA2Seq :: + (EffortIndex -> ra -> ra) -> + (ConvergRealSeq ra) -> + (ConvergRealSeq ra) +convertFuncRA2Seq f (ConvergRealSeq argSeq) = + ConvergRealSeq resultSeq + where + resultSeq ix = + f ix (argSeq ix) + +{-| + The same as above, where f is binary +-} +convertBinFuncRA2Seq :: + (EffortIndex -> ra -> ra -> ra) -> + (ConvergRealSeq ra) -> + (ConvergRealSeq ra) -> + (ConvergRealSeq ra) + +convertBinFuncRA2Seq f (ConvergRealSeq arg1) (ConvergRealSeq arg2) = + ConvergRealSeq resultSeq + where + resultSeq ix = + f ix (arg1 ix) (arg2 ix) + +{-| + Turn an arbitrary convergent sequence into one with + a guaranteed convergence rate - the precision (as defined + by 'RA.ERApprox.RA.getPrecision') of x_ix is at least ix. +-} +makeFastConvergRealSeq :: + (RA.ERApprox ra) => + (ConvergRealSeq ra) -> + (ConvergRealSeq ra) +makeFastConvergRealSeq (ConvergRealSeq argSeq) = + ConvergRealSeq fastSeq + where + fastSeq ix = + head $ catMaybes $ map (precisionOK . argSeq) indexSeries + where + indexSeries = + -- take 5 $ -- upper bound on iteration - for testing + binGeomSeries (max 1 ix) + precisionOK ra + | RA.getPrecision ra >= (effIx2prec ix) = Just ra + | otherwise = Nothing + +{-| + binGeomSeries n is the geometric series + [ n, 2n, 4n, 8n, ...] +-} +binGeomSeries + :: (Num a) + => a + -> [a] +binGeomSeries n = + n : (binGeomSeries (2 * n)) + +instance (RA.ERApprox ra) => Show (ConvergRealSeq ra) + where + show = showConvergRealSeq 6 True True 10 -- cheating here, should throw an error + + +{-| + Show function for ConvergRealSeq's with full arguments. +-} +showConvergRealSeq + :: (RA.ERApprox ra) + => Int + -> Bool + -> Bool + -> Precision + -> (ConvergRealSeq ra) + -> String + +showConvergRealSeq numDigits showGran showComponents prec r = + RA.showApprox numDigits showGran showComponents $ + convergRealSeqElem (makeFastConvergRealSeq r) (prec2effIx prec) + + +{-| + Show function for ConvergRealSeq's with all parameters fixed + except for number of digits +-} +showConvergRealSeqAuto + :: (RA.ERApprox ra) + => Int + -> (ConvergRealSeq ra) + -> String +showConvergRealSeqAuto numDigits argSeq = + showConvergRealSeq numDigits True True prec argSeq + where + prec = effIx2prec $ ceiling $ (fromInteger $ toInteger numDigits) * 3.3219280948873626 + + + +instance + (RA.ERApprox ra) + => Eq (ConvergRealSeq ra) + where + r1 == r2 = + iterateRA_A raEq 2 [r1, r2] + where + raEq _ ([a1,a2]) = RA.equalReals a1 a2 + +instance + (RA.ERApprox ra) + => Ord (ConvergRealSeq ra) + where + compare r1 r2 = + iterateRA_A eraComp 2 [r1, r2] + where + eraComp _ ([a1,a2]) = RA.compareReals a1 a2 + +pointwiseConvergRealSeq1 f (ConvergRealSeq sq) = + ConvergRealSeq (f . sq) +pointwiseConvergRealSeq2 f (ConvergRealSeq sq1) (ConvergRealSeq sq2) = + ConvergRealSeq (\ix -> f (sq1 ix) (sq2 ix)) + +instance + (RA.ERApprox ra) + => Num (ConvergRealSeq ra) + where + fromInteger n = ConvergRealSeq sq + where + sq ix = + RA.setMinGranularity (effIx2gran ix) $ fromInteger n + abs = pointwiseConvergRealSeq1 $ abs + signum = pointwiseConvergRealSeq1 $ signum + negate = pointwiseConvergRealSeq1 $ negate + (+) = pointwiseConvergRealSeq2 $ (+) + (*) = pointwiseConvergRealSeq2 $ (*) + +instance + (RA.ERApprox ra) + => Fractional (ConvergRealSeq ra) + where + fromRational q = ConvergRealSeq sq + where + sq ix = + (RA.setMinGranularity (effIx2gran ix) num) / denom + num = fromInteger $ numerator q + denom = fromInteger $ denominator q + recip = pointwiseConvergRealSeq1 $ recip + +{-| + Take a converging sequence of partial functions F_i that operate on + real approximations and turn it into a function F that operates on converging sequences. + F looks for some members of the real approximation sequences + and an i so that F_i is defined for the chosen approximations + and returns its result. +-} +iterateRA_A + :: (EffortIndex -> [ra] -> Maybe a) + -- ^ a sequence of partial functions based on approximations + -> EffortIndex -- ^ a starting index to use when searching sequences + -> ([ConvergRealSeq ra] -> a) + -- ^ a total function based on sequences + +iterateRA_A fn_RA startIx args = + head $ catMaybes $ map ((uncurry fn_RA) . args_Prec) indexSeries + where + indexSeries = +-- take 5 $ -- upper bound on iteration - for testing + binGeomSeries $ max 1 startIx + -- [(max 1 startIx)..] + args_Prec currentIndex = + (currentIndex, map (\ arg -> convergRealSeqElem arg currentIndex) args) + + diff --git a/src/Data/Number/ER/Real/Arithmetic/Elementary.hs b/src/Data/Number/ER/Real/Arithmetic/Elementary.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/Real/Arithmetic/Elementary.hs @@ -0,0 +1,602 @@ +{-| + Module : Data.Number.ER.Real.Arithmetic.Elementary + Description : some elementary functions + Copyright : (c) Michal Konecny, Amin Farjudian, Jan Duracz + License : LGPL + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : portable + + Some important elementary functions for real approximations + and their maximal extensions for interval approximations. +-} +module Data.Number.ER.Real.Arithmetic.Elementary +( + -- * specialised exponentiation + erSqr_R, + erSqr_IR, + erPow_R, + erPow_IR, + erSqrt_R, + erSqrt_IR, + erRoot_R, + erRoot_IR, + -- * exponentiation and logarithm + erExp_R, + erExp_IR, + erLog_R, + erLog_IR, + -- * trigonometrics + erSine_R, + erSine_IR, + erCosine_R, + erCosine_IR, + erATan_R, + erATan_IR, + erPi_R +) +where + +import qualified Data.Number.ER.Real.Approx as RA +import Data.Number.ER.BasicTypes + +import Data.Number.ER.Real.Arithmetic.Taylor +-- import Data.Number.ER.Real.Arithmetic.Newton + +import Data.Number.ER.Misc + +{- + sqr +-} + +erSqr_IR :: + (RA.ERIntApprox ira) => + EffortIndex -> + ira -> ira +erSqr_IR = erSqr_R + +erSqr_R :: + (RA.ERIntApprox ira) => + EffortIndex -> + ira -> ira +erSqr_R ix a + | RA.isEmpty a = + RA.emptyApprox + | otherwise = + max 0 $ a' * a' + where + a' = RA.setMinGranularity gran a + gran = effIx2gran ix + +{- + integer exponentiation x ^ p +-} + +erPow_IR :: + (RA.ERIntApprox ira) => + EffortIndex -> + Integer -> + ira -> ira +erPow_IR = erPow_R + +erPow_R :: + (RA.ERIntApprox ira) => + EffortIndex -> + Integer -> + ira -> ira +erPow_R ix p a + | RA.isEmpty a = + RA.emptyApprox + | p < 0 = + 1 / erPow_R ix (-p) a + | p == 0 = + 1 + | even p = + erPow_R ix (div p 2) (erSqr_R ix a') + | otherwise = + a' * (erPow_R ix (div (p - 1) 2) (erSqr_R ix a')) + where + a' = RA.setMinGranularity gran a + gran = effIx2gran ix + +{- + sqrt +-} + +erSqrt_R :: + (RA.ERIntApprox ira) => + EffortIndex -> ira -> ira +erSqrt_R = erSqrtNewton_R + +erSqrt_IR :: + (RA.ERIntApprox ira) => + EffortIndex -> ira -> ira +erSqrt_IR = + RA.maxExtensionR2R + sqrtExtrema + (\ ix x -> erSqrt_R ix x) + +sqrtExtrema ix x + | 0 `RA.refines` x = [0] + | otherwise = [] + + + +erSqrtContFr_R :: + (RA.ERIntApprox ira) => + EffortIndex -> ira -> ira +erSqrtContFr_R ix a + | aR == 0 = 0 + | aL == 1/0 = 1/0 + | aR < 0 = RA.emptyApprox + | otherwise = + contFrIter (ix + 3) $ + RA.setMinGranularity gran $ max 0 (0 RA.\/ a) + where + gran = effIx2gran ix + (aL, aR) = RA.bounds a + aM1 = a - 1 + + contFrIter i x_i + | i == 0 = + x_i + | otherwise = + 1 + (aM1 / (x_iPlus1 + 1)) + where + x_iPlus1 = contFrIter (i - 1) x_i + +erSqrtNewton_R :: + (RA.ERIntApprox ira) => + EffortIndex -> ira -> ira +erSqrtNewton_R ix a + | RA.isEmpty a = RA.emptyApprox + | aR == 0 = 0 + | aL == 1/0 = 1/0 + | aR < 0 = RA.emptyApprox + | otherwise = + x_i RA.\/ (a/x_i) + where + gran = effIx2gran ix + (aL, aR) = RA.bounds a + aM1 = a - 1 + + x_i = + newtonIter ((ix `div` 100) + 5) $ + RA.setMinGranularity gran $ max 0 aR + newtonIter i x_i + | i == 0 = x_i + | otherwise = + snd $ RA.bounds $ + (x_iMinus1 + a / (x_iMinus1)) / 2 + where + x_iMinus1 = newtonIter (i - 1) x_i + +{- + pth root x ^ (1/p) +-} + +erRoot_R :: + (RA.ERIntApprox ira) => + EffortIndex -> Integer -> ira -> ira +erRoot_R = erRootNewton_R + +erRoot_IR :: + (RA.ERIntApprox ira) => + EffortIndex -> Integer -> ira -> ira +erRoot_IR ix p = + RA.maxExtensionR2R + (rootExtrema p) + (\ ix x -> erRoot_R ix p x) $ + ix + +rootExtrema p ix x + | 0 `RA.refines` x && even p = [0] + | otherwise = [] + +erRootNewton_R :: + (RA.ERIntApprox ira) => + EffortIndex -> Integer -> ira -> ira +erRootNewton_R ix p a + | RA.isEmpty a = RA.emptyApprox + | aR == 0 = 0 + | aL == 1/0 = 1/0 + | aR < 0 && even p = RA.emptyApprox + | aR < 0 = - erRootNewton_R ix p (-a) + | p > 0 = + x_i RA.\/ (a/x_i_pow_p_minus_1) + | otherwise = + 1 / (erRootNewton_R ix (-p) a) -- TODO: check extremes + where + gran = effIx2gran ix + (aL, aR) = RA.bounds a + aM1 = a - 1 + pIRA = fromInteger p + pIRA_minus_1 = pIRA - 1 + + (x_i, x_i_pow_p_minus_1) = + newtonIter (ix + 5) $ + RA.setMinGranularity gran $ max 0 aR + + newtonIter i x_0 + | i == 0 = + (x_0, x_0_pow_p_minus_1) + | otherwise = + (x_i, x_i_pow_p_minus_1) + + where + (x_iMinus1, x_iMinus1_pow_p_minus_1) = + newtonIter (i - 1) x_0 + x_i = + snd $ RA.bounds $ + (pIRA_minus_1 * x_iMinus1 + a / x_iMinus1_pow_p_minus_1) / pIRA + x_i_pow_p_minus_1 = + erPow_R ix (p - 1) x_i + x_0_pow_p_minus_1 = + erPow_R ix (p - 1) x_0 + +{- + e^x and log +-} + +erExp_R :: + (RA.ERIntApprox ira) => + EffortIndex -> ira -> ira + +erExp_R = erExp_Tay_Opt_R + +{- + exp as derived from Taylor series is already a maximal extension + because it does not suffer from the wrapping effect - all + functions used are monotone - all Taylor coeffs are positive +-} +erExp_IR :: + (RA.ERIntApprox ira) => + EffortIndex -> ira -> ira + +erExp_IR ix x + | 0 `RA.refines` x || (-1/0) `RA.refines` x= + RA.maxExtensionR2R + (\ ix x -> []) + (\ ix x -> erExp_R ix x) + ix x + | otherwise = + erExp_R ix x + + +{- Log using Newton -} + +erLog_R :: + (RA.ERIntApprox ira) => + EffortIndex -> ira -> ira + +erLog_R = + logDivSeries_R +-- erLog_IR -- intervals are more efficient for log than singletons + +erLog_IR :: + (RA.ERIntApprox ira) => + EffortIndex -> ira -> ira + +erLog_IR = + RA.maxExtensionR2R + logExtrema + (\ ix x -> logDivSeries_R ix x) + +logExtrema ix x + | 0 `RA.refines` x = [-1/0] + | otherwise = [] + +{-| log using a fast converging series, designed to be used with singletons -} +logDivSeries_R :: + (RA.ERIntApprox ira) => EffortIndex -> ira -> ira +logDivSeries_R ix x + | RA.isEmpty posx = RA.emptyApprox + | posx `RA.refines` 0 = -1/0 + | posx `RA.refines` (1/0) = 1/0 + | otherwise = + case RA.compareReals posx 1 of + Just LT -> + - (logDivSeries_R ix (recip posx)) + _ -> + nearLogx + 2 * t * (series ix (RA.setMinGranularity gran 1)) + where + gran = effIx2gran ix + posx = (RA.setMinGranularity gran x) RA./\ (0 RA.\/ (1/0)) + nearLogx = + 0.69314718055994530941 * (fromInteger $ intLog 2 $ xCeiling) + remNearLogx = + posx / (erExp_R ix nearLogx) -- should be very close to 1 + xCeiling = + snd $ RA.integerBounds posx + t = + ((remNearLogx - 1) / (remNearLogx + 1)) -- the range of this expression is [-1,1] + RA./\ ((-1) RA.\/ 1) -- correction of wrapping + tsqare = abs $ t * t -- the range is [0,1] + series termsCount currentDenominator + | termsCount > 0 = + (recip currentDenominator) + tsqare * (series (termsCount - 1) (currentDenominator + 2)) + | otherwise = + (recip currentDenominator) + * (1 RA.\/ (recip $ 1 - tsqare)) -- [1,1/(1-t^2)] is a valid error bound + +--{- log using Newton -} +-- +--logNewton_RA +-- :: (RA.ERIntApprox ira) +-- => EffortIndex +-- -> ra -- must not be below 1 +-- -> ra +-- +--logNewton_RA i x = +-- case compareReals posx 1 of +-- Just LT -> +-- - (logNewton_RA i (recip posx)) +-- _ -> +-- erNewton_FullArgs +-- ( \ i y -> (erExp_RA i y) - posx, erExp_RA) +-- (RA.setMinGranularity gran nearLogx) +-- (RA.setMinGranularity gran 1) +-- (fromInteger $ toInteger i) +-- i +-- where +-- gran = effIx2gran i +-- posx = +-- RA.setMinGranularity gran x /\ (ira2ra $ 0 RA.\/ (1/0)) +-- nearLogx = +-- 0.69314718055994530941 * (fromInteger $ intLog 2 $ xCeiling) +-- xCeiling +-- | RA.isEmpty posx = 1 -- choice of constant irrelevant +-- | otherwise = +-- snd $ RA.iraIntegerBounds $ ra2ira posx + + +{- + sin(x) and cos(x) +-} + +erSine_R :: + (RA.ERIntApprox ira) => + EffortIndex -> ira -> ira + +erSine_R ix x + | (1/0) `RA.refines` x || (-1/0) `RA.refines` x = + (-1) RA.\/ 1 + | otherwise = + erSine_Tay_Opt_R ix x + +erCosine_R :: + (RA.ERIntApprox ira) => + EffortIndex -> ira -> ira + +erCosine_R ix x + | (1/0) `RA.refines` x || (-1/0) `RA.refines` x = + (-1) RA.\/ 1 + | otherwise = + erCosine_Tay_Opt_R ix x + + +{- Sine using generic Taylor (see Taylor for an optimised version) -} + +erSine_Tay_R :: + (RA.ERIntApprox ira) => + EffortIndex -> ira -> ira + +erSine_Tay_R ix x + | (1/0) `RA.refines` x || (-1/0) `RA.refines` x = + (-1) RA.\/ 1 + | otherwise = + erTaylor_R ix sine_coefSeq sine_error 0 x + +sine_coefSeq :: + (RA.ERIntApprox ira) => + Int -> ira + +sine_coefSeq n + | n `mod` 4 == 0 = 0 + | n `mod` 4 == 1 = 1 + | n `mod` 4 == 2 = 0 + | n `mod` 4 == 3 = -1 + +sine_error n = (-1) RA.\/ 1 + +{- maximal extensions -} + +erSine_IR :: + (RA.ERIntApprox ira) => + EffortIndex -> ira -> ira + +erSine_IR = + RA.maxExtensionR2R sineExtremes erSine_R + +erCosine_IR :: + (RA.ERIntApprox ira) => + EffortIndex -> ira -> ira + +erCosine_IR = + RA.maxExtensionR2R cosineExtremes erCosine_R + +sineExtremes ix x + | RA.isBounded x = + alternatingExtremes 1 (-1) ix scaledX + | otherwise = [-1,1] + where + scaledX = (x / (erPi_R ix)) - 0.5 + +cosineExtremes ix x + | RA.isBounded x = + alternatingExtremes 1 (-1) ix scaledX + | otherwise = [-1,1] + where + scaledX = (x / (erPi_R ix)) + +alternatingExtremes extr0 extr1 ix scaledX + | extremesCount >= 2 = [extr0,extr1] + | extremesCount == 1 && even minExtremeN = [extr0] + | extremesCount == 1 = [extr1] + | otherwise = [] + where + extremesCount = 1 + maxExtremeN - minExtremeN + (xFloor, xCeiling) = RA.integerBounds scaledX + minExtremeN = + case RA.compareReals (fromInteger $ toInteger xFloor) scaledX of + Just LT -> (xFloor + 1) + _ -> xFloor + maxExtremeN = + case RA.compareReals scaledX (fromInteger $ toInteger xCeiling) of + Just LT -> xCeiling - 1 + _ -> xCeiling + + +{- + tan(x), atan(x) and pi +-} + +erATan_R :: + (RA.ERIntApprox ira) => + EffortIndex -> ira -> ira + +erATan_R = atanEuler_R + +erATan_IR :: + (RA.ERIntApprox ira) => + EffortIndex -> ira -> ira + +erATan_IR = + RA.maxExtensionR2R atanExtremes erATan_R + +atanExtremes ix x = [] + +{- atan using Euler's series: + x / (1 + x^2) * (1 + t*2*1/(2*1 + 1)*(1 + t*2*2/(2*2 + 1)*(1 + ... (1 + t*2*n/(2*n+1)*(1 + ...))))) + where + t = x^2/(1 + x^2) + + where the tail (1 + t*2*n/(2*n+1)*(1 + ...)) is inside the interval: + [1 + (x^2*2n/(2n + 1)), 1 + x^2] +-} + +atanEuler_R :: + (RA.ERIntApprox ira) => + EffortIndex -> ira -> ira + +atanEuler_R ix x + | RA.isEmpty x = RA.emptyApprox + | otherwise = + (x / xSquarePlus1) * (series ix (RA.setMinGranularity gran 2)) + where + gran = effIx2gran ix + series termsCount coeffBase + | termsCount > 0 = + 1 + xSquareOverXSquarePlus1 * coeff * (series (termsCount - 1) (coeffBase + 2)) + | otherwise = + 1 + xSquare * (coeff RA.\/ 1) + where + coeff = coeffBase / (coeffBase + 1) + xSquare = abs $ x * x + xSquarePlus1 = xSquare + 1 + xSquareOverXSquarePlus1 = xSquare / xSquarePlus1 + +--{- atan using Newton -} +-- +--atanNewton_RA :: +-- (RA.ERIntApprox ira) => +-- EffortIndex -> ra -> ra +-- +--atanNewton_RA i x = +-- erNewton_FullArgs +-- ( \ i y -> (erTan_RA i y) - x, erTanDeriv_RA) +-- (RA.setMinGranularity (effIx2gran i) (x)) +-- (RA.setMinGranularity (effIx2gran i) 1) +-- (fromInteger $ toInteger i) +-- i + +{- tan -} + +erTan_R :: + (RA.ERIntApprox ira) => + EffortIndex -> ira -> ira + +erTan_R ix x = + (erSine_R ix x) / (erCosine_R ix x) + +erTanDeriv_R ix x = + recip $ abs $ cosx * cosx + where + cosx = erCosine_R ix x + + +{- pi -} + +{-| + pi using Bellard's formula + + Convergence properties: + + * shrinking sequence + + * rate at least 2^(-i). + +-} +erPi_R :: + (RA.ERIntApprox ira) => + EffortIndex -> ira +erPi_R = piBellard_R + +-- | pi using atan +piAtan_R :: + (RA.ERIntApprox ira) => + EffortIndex -> ira +piAtan_R ix = + (*) 4 $ atanEuler_R ix 1 + +{-| + pi using Bellard's formula + (see ) + + Convergence properties: + + * shrinking sequence + + * rate at least 2^(-i). + +-} +piBellard_R :: + (RA.ERIntApprox ira) => + EffortIndex -> ira +piBellard_R ix = + r1over64 * (sum $ reverse $ bellardTerms 0 (10 + (ix `div` 10)) (1,z,z)) + {- + sum from the smallest to the largest + (got this trick from Martin Escardo who said he got it from Andrej Bauer) + + the rounding error dominates the truncation error to such + a degree that the truncation error can be safely left out + + each bellard term contributes 10 binary digits that the following terms + do not influence + -} + where + gran = max 0 (effIx2gran ix) + 10 + r1over64 = (RA.setMinGranularity gran 1) / 64 + r1over1024 = (RA.setMinGranularity gran 1) / 1024 + z = RA.setMinGranularity gran 0 + bellardTerms n nMax (mult, r4n, r10n) + | n >= nMax = [] + | otherwise = + termN : rest + where + rest = + bellardTerms (n + 1) nMax (- mult * r1over1024, r4n + 4, r10n + 10) + termN = + mult * bellardSum + bellardSum = + -- sum from the smallest to the largest + (recip $ r10n + 9) + - (recip $ r4n + 3) + - 4 * ((recip $ r10n + 7) + (recip $ r10n + 5)) + - (64 / (r10n + 3)) + - (32 / (r4n + 1)) + + (256 / (r10n + 1)) + + diff --git a/src/Data/Number/ER/Real/Arithmetic/Integration.hs b/src/Data/Number/ER/Real/Arithmetic/Integration.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/Real/Arithmetic/Integration.hs @@ -0,0 +1,141 @@ +{-| + Module : Data.Number.ER.Real.Arithmetic.Integration + Description : simple integration methods + Copyright : (c) Michal Konecny + License : LGPL + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : portable + + Simple integration methods for Haskell functions operating + on real number approximations. +-} +module Data.Number.ER.Real.Arithmetic.Integration +( + integrateCont, +-- integrateDiff, + integrateCont_R, + integrateContAdapt_R +) +where + +import qualified Data.Number.ER.Real.Approx as RA +import Data.Number.ER.BasicTypes +import Data.Number.ER.Real.Approx.Sequence +import Data.Number.ER.Real.Arithmetic.Elementary + +testIntegr1 :: + (RA.ERIntApprox ira) => + (ConvergRealSeq ira) +testIntegr1 = integrateCont erExp_IR 0 1 + +integrateCont :: + (RA.ERIntApprox ira) => + (EffortIndex -> ira -> ira) -> + (ConvergRealSeq ira) -> (ConvergRealSeq ira) -> (ConvergRealSeq ira) + +integrateCont f = convertBinFuncRA2Seq $ integrateContAdapt_R f + +integrateDiff :: + (RA.ERIntApprox ira) => + (EffortIndex -> ira -> ira, EffortIndex -> ira -> ira) -> + (ConvergRealSeq ra) -> (ConvergRealSeq ra) -> (ConvergRealSeq ra) + +integrateDiff f = convertBinFuncRA2Seq $ integrateDiffAdapt_RA f + + +{-| + naive integration, using a partition of 2 * prec equally sized intervals +-} +integrateCont_R :: + (RA.ERIntApprox ira) => + (EffortIndex -> ira -> ira) -> + EffortIndex -> (ira) -> (ira) -> (ira) +integrateCont_R f ix a b = + sum $ map rectArea rectangles + where + rectArea (width, height) = width * height + rectangles = + zip (repeat width) $ map (f ix) covering + width = (b - a) / (fromInteger rectCount) + rectCount = 2 * (fromInteger $ toInteger gran) + gran = effIx2gran ix + covering = getCoveringIntervals division + getCoveringIntervals ( pt1 : pt2 : rest ) = + ((pt1) RA.\/ (pt2)) : (getCoveringIntervals $ pt2 : rest) + getCoveringIntervals _ = [] + division = map getEndPoint $ [0..rectCount] + getEndPoint n = + a + ((fromInteger n) * width) + +{-| + integration using divide and conquer adaptive partitioning +-} +integrateContAdapt_R :: + (RA.ERIntApprox ira) => + (EffortIndex -> ira -> ira) -> + EffortIndex -> (ira) -> (ira) -> (ira) +integrateContAdapt_R f ix a b = + sum rectangleAreas + where + rectangleAreas = + getRs ix a b + getRs subix l r + | RA.getPrecision area >= prec = [area] + | otherwise = + (getRs nsubix l m) ++ (getRs nsubix m r) + where + prec = foldl1 min [effIx2prec ix, RA.getPrecision l, RA.getPrecision r] + area = (r - l) * (f subix (l RA.\/ r)) + nsubix = subix + 1 + m = (l + r)/2 + + +{-| + integration using divide and conquer adaptive partitioning + making use of the derivative of the integrated function +-} +integrateDiffAdapt_RA :: + (RA.ERIntApprox ira) => + (EffortIndex -> ira -> ira, EffortIndex -> ira -> ira) -> + EffortIndex -> (ra) -> (ra) -> (ra) +integrateDiffAdapt_RA f prec a b = + error "TODO" + +{- + sum rectangleAreas + where + rectangleAreas = + getRs prec a b + getRs p l r + | getPrecision area >= prec = [area] + | otherwise = + (getRs np l m) ++ (getRs np m r) + where + np = p + 1 + m = (l + r)/2 +-- area = areaDiff + area = areaRect /\ areaDiff + -- merge the information given by the rectangle method + -- with the information given by the derivative method + areaRect = w * fVal -- same as in integrateContAdapt_R + (fVal, fDeriv) = applyRdiffR f p (l \/ r) + w = r - l + areaDiff + | isExact fDeriv = w * (fl + fr) / 2 -- derivative is constant and perfectly known + | otherwise = areaLow \/ areaHigh + fl = fst $ applyRdiffR f (2 * p) l + fr = fst $ applyRdiffR f (2 * p) r + -- interestingly, we have to request fl, fr with higher precision than + -- we requested fDeriv so that the derivative would be of any use + -- with these values - replace (2 * p) by p and it will not converge! + -- area computed by a scary formula: + areaLow = t + w * (fl * dHigh - fr * dLow) / dDiff + areaHigh = - t - w * (fl * dLow - fr * dHigh) / dDiff -- swap dHigh and dLow + t = (w^2 * dLow * dHigh + (fr - fl)^2)/(2 * dDiff) + dDiff = dHigh - dLow + (dLow, dHigh) = getBounds fDeriv +-} + + diff --git a/src/Data/Number/ER/Real/Arithmetic/Newton.hs b/src/Data/Number/ER/Real/Arithmetic/Newton.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/Real/Arithmetic/Newton.hs @@ -0,0 +1,201 @@ +{-| + + Module : Data.Number.ER.Real.Arithmetic.Taylor + Description : interval Newton method + Copyright : (c) Amin Farjudian, Michal Konecny + License : LGPL + + Maintainer : mik@konecny.aow.cz + Stability : alpha + Portability : portable + + Interval Newton's method for root finding. + + To be used for obtaining functions out of their inverse(s) over various + intervals. +-} +module Data.Number.ER.Real.Arithmetic.Newton +( + erNewton_FullArgs, + erNewton_mdfd_FullArgs +) +where + +import qualified Data.Number.ER.Real.Approx as RA +import Data.Number.ER.BasicTypes +import Data.Number.ER.Real.Arithmetic.Taylor + +erNewton_FullArgs + :: (RA.ERIntApprox ira) + => (EffortIndex -> ira -> ira, EffortIndex -> ira -> ira) -- ^ a function and its derivative + -> ira -- ^ a starting point + -> ira -- ^ a lower bound of the absolute value of the derivative over the working interval + -> Int -- ^ number of iterations + -> EffortIndex -- ^ the initial index to use for argument function and its derivative + -> ira -- ^ the result + +erNewton_FullArgs (f ,df) startPt minDrv iterCnt i = + erNewton_FullArgs_aux startPt startOtherPt iterCnt + where + erNewton_FullArgs_aux newtonPt otherPt iterCnt + | (iterCnt <= 0 || RA.getPrecision result >= prec) = + result + | otherwise = + erNewton_FullArgs_aux newNewtonPt newOtherPt (iterCnt - 1) + where + result = + newtonPt RA.\/ otherPt + prec = effIx2prec i + newNewtonPt = + midPoint $ RA.bounds $ + (newtonPt - ( (f i newtonPt) / (( df i newtonPt)))) + -- /\ (ira2ra ((ra2ira minDrv) \/ 100000000))))) + newOtherPt = otherEndPoint newNewtonPt + startOtherPt = otherEndPoint startPt + otherEndPoint a = a - ((f i a) / minDrv) -- /\ (0 \/ 10000000) + + +{-| + This auxiliary function returns the average of two ra's +-} +midPoint + :: (RA.ERIntApprox ira) + => (ira ,ira) + -> ira +midPoint (x1, x2) = (x1 + x2) / 2 + + +{-| Modified Newton Method + Notes: + + 1. It has a cubic convergence speed, as opposed to the original Newton's + square convergence speed. + + 2. It does not deal with multiple roots. + + 3. Per iteration, it makes two queries on the derivative, so it best + suits the cases where computation of the derivative is at most as + expensive as the function itself. +-} +erNewton_mdfd_FullArgs + :: (RA.ERIntApprox ira) + => (EffortIndex -> ira -> ira, EffortIndex -> ira -> ira) -- ^ a function and its derivative + -> ira -- ^ a starting point + -> ira -- ^ The minimum of absolute value of derivative over the working interval + -> Int -- ^ number of iterations + -> EffortIndex -- ^ It triggers the initial index to be called by the argument function and its derivative. + -> ira -- ^ the result + +erNewton_mdfd_FullArgs (f ,df) startPt minDrv iterCnt i = + erNewton_FullArgs_aux startPt startOtherPt iterCnt + where + erNewton_FullArgs_aux newtonPt otherPt iterCnt + | iterCnt <= 0 = newtonPt RA.\/ otherPt + | otherwise = erNewton_FullArgs_aux newNewtonPt newOtherPt (iterCnt - 1) + where + valueAtNewOtherPt = f i newOtherPt + derivAtNewtonPt = df i newOtherPt + unblurredDeriv = midPoint $ RA.bounds $ derivAtNewtonPt + intermediaryPt = midPoint $ RA.bounds $ newtonPt - valueAtNewOtherPt / (2 * derivAtNewtonPt) + derivAtIntermediaryPt = df i intermediaryPt + newNewtonPt = + midPoint $ RA.bounds $ + (newtonPt - ( valueAtNewOtherPt / derivAtIntermediaryPt)) + newOtherPt = otherEndPoint newNewtonPt + startOtherPt = otherEndPoint startPt + otherEndPoint a = a - ((f i a) / minDrv) + +erNewton_mdfd + :: (RA.ERIntApprox ira) + => (EffortIndex -> ira -> ira, EffortIndex -> ira -> ira) -- ^ a function and its derivative + -> ira -- ^ a starting point + -> ira -- ^ The minimum of absolute value of derivative over the working interval + -> EffortIndex -- ^ It triggers the initial index to be called by the argument function and its derivative. + -> ira -- ^ the result + +erNewton_mdfd (f ,df) startPt minDrv i = + erNewton_mdfd_FullArgs (f, df) startPt minDrv (fromInteger $ toInteger $ i) i + +--apNewton_mdfd +-- :: (RA.ERIntApprox ira) +-- => (EffortIndex -> ra -> ra, EffortIndex -> ra -> ra) -- ^ a function and its derivative +-- -> ra -- ^ a starting point +-- -> ra -- ^ The minimum of absolute value of derivative over the working interval +-- -> EffortIndex -- ^ It triggers the initial index to be called by the argument function and its derivative. Moreover, the number of iterations are predefined by this argument. +-- -> ra -- ^ the result +-- +--apNewton_mdfd (f, df) startPt minDrv i = +-- erNewton_mdfd_FullArgs +-- + +--id_RA +-- :: (RA.ERIntApprox ira) +-- => EffortIndex -> ira -> ira +-- +--id_RA i x = x +-- +--const_one_RA +-- :: (RA.ERIntApprox ira) +-- => EffortIndex -> ira -> ira +-- +--const_one_RA i x = (setMinGranularity (effIx2gran i) 1) +-- +-- +--test_erNewton_FullArgs_01_RA +-- :: (RA.ERIntApprox ira) +-- => EffortIndex -> ira -> ira +-- +--test_erNewton_FullArgs_01_RA i x = erNewton_FullArgs_01 ( id_RA, const_one_RA) x 10 i +-- +--test_erNewton_FullArgs_01 +-- :: (RA.ERIntApprox ira) +-- => (ConvergRealSeq ira) -> (ConvergRealSeq ira) +-- +--test_erNewton_FullArgs_01 = convertFuncRA2Seq test_erNewton_FullArgs_01_RA +-- +--exp_Ra_minus_r_RA +-- :: (RA.ERIntApprox ira) +-- => EffortIndex -> ira -> ira -> ira +-- +--exp_Ra_minus_r_RA i r x = (erExp_RA i x) - r +-- +--exp_Ra_minus_r +-- :: (RA.ERIntApprox ira) +-- => (ConvergRealSeq ira) -> (ConvergRealSeq ira) -> (ConvergRealSeq ira) +-- +--exp_Ra_minus_r = convertBinFuncRA2Seq exp_Ra_minus_r_RA +-- +--logNewton_RA_02 +-- :: (RA.ERIntApprox ira) +-- => EffortIndex -> ira -> ira +-- +--logNewton_RA_02 i x = +-- erNewton_FullArgs_02 +-- ( \ i y -> (erExp_RA i y) - x, erExp_RA) +-- (setMinGranularity (effIx2gran i) (2)) +-- (setMinGranularity (effIx2gran i) 1) +-- i +-- +--logNewton_02 +-- :: (RA.ERIntApprox ira) +-- => (ConvergRealSeq ira) -> (ConvergRealSeq ira) +-- +--logNewton_02 = convertFuncRA2Seq logNewton_RA_02 + + +--logNewton_mdfd_RA +-- :: (RA.ERIntApprox ira) +-- => EffortIndex -> ira -> ira + +--logNewton_mdfd_RA i x = +-- erNewton_mdfd_FullArgs +-- ( \ i y -> (erExp_RA i y) - x, erExp_RA) +-- (setMinGranularity (effIx2gran i) (2)) +-- (setMinGranularity (effIx2gran i) 1) +-- i + +--logNewton_mdfd +-- :: (RA.ERIntApprox ira) +-- => (ConvergRealSeq ira) -> (ConvergRealSeq ira) +-- +--logNewton_mdfd = convertFuncRA2Seq logNewton_mdfd_RA diff --git a/src/Data/Number/ER/Real/Arithmetic/Taylor.hs b/src/Data/Number/ER/Real/Arithmetic/Taylor.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/Real/Arithmetic/Taylor.hs @@ -0,0 +1,177 @@ +{-| + Module : Data.Number.ER.Real.Arithmetic.Taylor + Description : implementation of Taylor series + Copyright : (c) Amin Farjudian, Michal Konecny + License : LGPL + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : portable + + Taylor series related functions. +-} +module Data.Number.ER.Real.Arithmetic.Taylor where + +import qualified Data.Number.ER.Real.Approx as RA +import qualified Data.Number.ER.ExtendedInteger as EI +import Data.Number.ER.BasicTypes + + +erTaylor_R + :: (RA.ERIntApprox ira) + => EffortIndex + -> (Int -> ira) -- ^ coefficients of the Taylor series + -> (Int -> ira) -- ^ function to estimate the n'th derivative between a and x + -> ira -- ^ centre of the Taylor Expansion + -> ira + -> ira +erTaylor_R ix coefSeq derivBounds a x = + erTaylor_R_FullArgs coefSeq derivBounds n a gran x + where + n = fromInteger ix + gran = fromInteger $ toInteger $ ix + +erTaylor_R_FullArgs + :: (RA.ERIntApprox ira) + => (Int -> ira) -- ^ coefficients of the Taylor series + -> (Int -> ira) -- ^ function to estimate the n'th derivative between a and x + -> Int -- ^ use this many elements of the series (+ account for error appropriately) + -> ira -- ^ centre of the Taylor Expansion + -> Granularity -- ^ make all constants have this granularity, thus influencing rounding errors + -> ira + -> ira +erTaylor_R_FullArgs coefSeq derivBounds n a gran x = + rec_apTaylor (RA.setMinGranularity gran 0) 0 + where + rec_apTaylor i j + | n > j = (coefSeq(j)) + + ((x - a)/(i+1)) * (rec_apTaylor (i+1) (j+1)) + | n == j = derivBounds n + | otherwise = + error "Data.Number.ER.Real.Arithmetic.Taylor.hs: erTaylor_RA_FullArgs: The index n cannot be negative" + +{-| + A Taylor series for exponentiation. +-} +erExp_Tay_Opt_R + :: (RA.ERIntApprox ira) + => EffortIndex + -> ira + -> ira +erExp_Tay_Opt_R ix x + | RA.isEmpty x = RA.emptyApprox + | x `RA.refines` (-1/0) = 0 -- -infty is not handled well by the Taylor formula + | otherwise = 1 + (te ix x (RA.setMinGranularity gran 1)) + where + gran = effIx2gran ix + te steps x i + | steps > 0 = + (x/i) * (1 + (te (steps - 1) x (i + 1))) + | steps == 0 = + errorBound + where + errorBound = + (x/i) * ithDerivBound + ithDerivBound + | xCeiling == EI.MinusInfinity = -- certainly -infty: + 0 + | xCeiling < 0 = -- certainly negative: + pow26xFloor RA.\/ 1 + | xFloor > 0 = -- certainly positive: + 1 RA.\/ pow28xCeiling + | otherwise = -- could contain 0: + pow26xFloor RA.\/ pow28xCeiling + where + (xFloor, xCeiling) = RA.integerBounds x + pow26xFloor + | xFloor == EI.MinusInfinity = + 0 + | otherwise = + ((RA.setMinGranularity gran 26)/10) ^^ xFloor + -- lower estimate of e^x + pow28xCeiling + | xCeiling == EI.PlusInfinity = + (1/0) + | otherwise = + ((RA.setMinGranularity gran 28)/10) ^^ xCeiling + -- upper estimate of e^x + +{- + The sine and cosine are implemented in almost exactly the same way +-} + +{-| + A Taylor series for sine. +-} +erSine_Tay_Opt_R + :: (RA.ERIntApprox ira) + => EffortIndex + -> ira + -> ira +erSine_Tay_Opt_R ix x = taylor_seg ix x (RA.setMinGranularity gran 1) + where + gran = effIx2gran ix + taylor_seg i x n -- 'i' for iterator + | i > 0 = x - ((x*x)/((n+1)*(n+2))) * (taylor_seg (i-2) x (n+2)) + | otherwise = errorRegion + where + errorRegion = (-1) RA.\/ (1) + +{-| + A Taylor series for cosine. +-} +erCosine_Tay_Opt_R + :: (RA.ERIntApprox ira) + => EffortIndex + -> ira + -> ira +erCosine_Tay_Opt_R ix x = taylor_seg ix x (RA.setMinGranularity gran 1) + where + gran = effIx2gran ix + taylor_seg i x n -- 'i' for iterator + | i > 0 = 1 - ((x*x)/(n*(n+1))) * (taylor_seg (i-2) x (n+2)) + | otherwise = errorRegion + where + errorRegion = (-1) RA.\/ (1) + + + +{-| Natural Logarithm: The following is a code for obtaining natural + logarithm using taylor series. Note that it only works for + x in [ 1, 2]. For other values, a scaling by factors of e^q is + best to do, i.e. if x is not in [1,2], then find some ratioal number + q where exp(q) * x is in [1,2]. Then you have: + log ( exp(q) * m) = q + log(m) +-} + +{-| Coefficients of the taylor series around a=1 -} +--logTayCoefs +-- :: (RA.ERIntApprox ira) +-- => Int -- up to how many terms of the Taylor series is desired +-- -> Int +-- -> ra +-- +--logTayCoefs n i +---- | i < 0 = error "ERTaylor.logTayCoefs: Negative n for the n-th term of Taylor series for logarithm" +-- | i == 0 = 0 +-- | i == n = fromInteger $ toInteger $ amFact(n-1) +-- | otherwise = fromInteger $ toInteger $ ((-1)^(i-1) * amFact(i-1)) +-- where +-- amFact (m) = product [1..m] + +--logTay_RA +-- :: (RA.ERIntApprox ira) +-- => EffortIndex +-- -> ra +-- -> ra +-- +--logTay_RA i = erTaylor_RA_FullArgs (logTayCoefs $fromInteger $ toInteger $ i) +-- (100000) (setMinGranularity (effIx2gran i) 1) (effIx2gran i) +-- +-- +--logTay +-- :: (RA.ERIntApprox ira) +-- => (ConvergRealSeq ra) +-- -> (ConvergRealSeq ra) +--logTay = convertFuncRA2Seq logTay_RA + diff --git a/src/Data/Number/ER/Real/Base.hs b/src/Data/Number/ER/Real/Base.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/Real/Base.hs @@ -0,0 +1,58 @@ +{-| + Module : Data.Number.ER.Real.Base + Description : class abstracting floats + Copyright : (c) Michal Konecny + License : LGPL + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : portable + + Abstraction over various fixed and floating point types as well + as rational numbers. + + This module should be included qualified as is often given the local + synonym B. +-} +module Data.Number.ER.Real.Base +( + module Data.Number.ER.BasicTypes, + ERRealBase(..) +) +where + +import Data.Number.ER.BasicTypes +import qualified Data.Number.ER.ExtendedInteger as EI + +import Data.Typeable + +{-| + This class is an abstraction of a subset of real numbers + with upwards rounded operations. +-} +class (Fractional rb, Ord rb) => ERRealBase rb + where + defaultGranularity :: rb -> Granularity + getApproxBinaryLog :: rb -> EI.ExtendedInteger + getGranularity :: rb -> Granularity + setMinGranularity :: Granularity -> rb -> rb + setGranularity :: Granularity -> rb -> rb + {-| + if @a@ is rounded to @ao@ then @|a-ao| <= getBaseMaxRounding ao@ + -} + getMaxRounding :: rb -> rb + isERNaN :: rb -> Bool + erNaN :: rb + isPlusInfinity :: rb -> Bool + plusInfinity :: rb + minusInfinity :: rb + minusInfinity = - plusInfinity + fromDouble :: Double -> rb + toDouble :: rb -> Double + fromFloat :: Float -> rb + toFloat :: rb -> Float + showDiGrCmp :: + Int {- ^ number of decimal digits to show -} -> + Bool {-^ whether to show granularity -} -> + Bool {-^ whether to show internal structure -} -> + rb -> String diff --git a/src/Data/Number/ER/Real/Base/CombinedMachineAP.hs b/src/Data/Number/ER/Real/Base/CombinedMachineAP.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/Real/Base/CombinedMachineAP.hs @@ -0,0 +1,238 @@ +{-# LANGUAGE DeriveDataTypeable #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-| + Module : Data.Number.ER.Real.Base.CombinedMachineAP + Description : auto-switching hardware-software floats + Copyright : (c) Michal Konecny + License : LGPL + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : non-portable (requires fenv.h) + + Arbitrary precision floating point numbers that use + machine double up to its precision. When a higher + granularity is required, it automatically switches + to another floating point type. +-} +module Data.Number.ER.Real.Base.CombinedMachineAP +( + ERMachineAP, + doubleDigits +) +where + +import qualified Data.Number.ER.Real.Base as B +import qualified Data.Number.ER.ExtendedInteger as EI +import Data.Number.ER.Real.Base.MachineDouble +import Data.Number.ER.Real.Base.Float +import Data.Number.ER.BasicTypes +import Data.Number.ER.Misc + +import Data.Typeable +import Data.Generics.Basics +import Data.Binary +--import BinaryDerive + +import Data.Ratio + +data ERMachineAP b = + ERMachineAPMachineDouble + { + machapfltDoubleGranularity :: Granularity + {-^ this has to be between 1 and 'doubleDigits' -} + , + machapfltDouble :: Double + } + | + ERMachineAPB + { + machapfltB :: b + } + deriving (Typeable, Data) + +doubleDigits = floatDigits (0 :: Double) + +{- the following has been generated by BinaryDerive -} +instance (Binary b) => (Binary (ERMachineAP b)) where + put (ERMachineAPMachineDouble a b) = putWord8 0 >> put a >> put b + put (ERMachineAPB a) = putWord8 1 >> put a + get = do + tag_ <- getWord8 + case tag_ of + 0 -> get >>= \a -> get >>= \b -> return (ERMachineAPMachineDouble a b) + 1 -> get >>= \a -> return (ERMachineAPB a) + _ -> fail "no parse" +{- the above has been generated by BinaryDerive -} + +lift1ERMachineAP :: + (Double -> Double) -> + (b -> b) -> + (ERMachineAP b -> ERMachineAP b) +lift1ERMachineAP fD fB (ERMachineAPMachineDouble g d) = + ERMachineAPMachineDouble g (fD d) +lift1ERMachineAP fD fB (ERMachineAPB b) = + ERMachineAPB (fB b) + +op1ERMachineAP :: + (Double -> a) -> + (b -> a) -> + (ERMachineAP b -> a) +op1ERMachineAP fD fB (ERMachineAPMachineDouble g d) = + fD d +op1ERMachineAP fD fB (ERMachineAPB b) = + fB b + +lift2ERMachineAP :: + (B.ERRealBase b) => + (Double -> Double -> Double) -> + (b -> b -> b) -> + (ERMachineAP b -> ERMachineAP b -> ERMachineAP b) +lift2ERMachineAP fD fB + (ERMachineAPMachineDouble g1 d1) (ERMachineAPMachineDouble g2 d2) = + ERMachineAPMachineDouble (max g1 g2) (fD d1 d2) +lift2ERMachineAP fD fB + (ERMachineAPMachineDouble g1 d1) (ERMachineAPB b2) = + ERMachineAPB $ fB (B.fromDouble d1) b2 +lift2ERMachineAP fD fB + (ERMachineAPB b1) (ERMachineAPMachineDouble g2 d2) = + ERMachineAPB $ fB b1 (B.fromDouble d2) +lift2ERMachineAP fD fB + (ERMachineAPB b1) (ERMachineAPB b2) = + ERMachineAPB $ fB b1 b2 + +op2ERMachineAP :: + (B.ERRealBase b) => + (Double -> Double -> a) -> + (b -> b -> a) -> + (ERMachineAP b -> ERMachineAP b -> a) +op2ERMachineAP fD fB + (ERMachineAPMachineDouble g1 d1) (ERMachineAPMachineDouble g2 d2) = + fD d1 d2 +op2ERMachineAP fD fB + (ERMachineAPMachineDouble g1 d1) (ERMachineAPB b2) = + fB (B.fromDouble d1) b2 +op2ERMachineAP fD fB + (ERMachineAPB b1) (ERMachineAPMachineDouble g2 d2) = + fB b1 (B.fromDouble d2) +op2ERMachineAP fD fB + (ERMachineAPB b1) (ERMachineAPB b2) = + fB b1 b2 + +instance (B.ERRealBase b) => Show (ERMachineAP b) + where + show = showERMachineAP 6 True True + +showERMachineAP numDigits showGran showComponents = + showEMA + where + maybeGran gr + | showGran = "{g=" ++ show gr ++ "}" + | otherwise = "" + maybeComps + | showComponents = "{Double}" + | otherwise = "" + showEMA (ERMachineAPMachineDouble gr d) = + show d ++ (maybeGran gr) ++ maybeComps + showEMA (ERMachineAPB b) = + B.showDiGrCmp numDigits showGran showComponents b + + +instance (B.ERRealBase b) => Eq (ERMachineAP b) + where + (==) = op2ERMachineAP (==) (==) + +instance (B.ERRealBase b) => Ord (ERMachineAP b) + where + compare = op2ERMachineAP compare compare + + + +instance (B.ERRealBase b) => Num (ERMachineAP b) + where + fromInteger n + | gran < doubleDigits = + ERMachineAPMachineDouble gran $ fromInteger n + | otherwise = + ERMachineAPB b + where + gran = B.getGranularity b + b = fromInteger n + abs = lift1ERMachineAP abs abs + signum = lift1ERMachineAP signum signum + negate = lift1ERMachineAP negate negate + (+) = lift2ERMachineAP (+) (+) + (*) = lift2ERMachineAP (*) (*) + +instance (B.ERRealBase b) => Fractional (ERMachineAP b) + where + fromRational rat = + (fromInteger $ numerator rat) + / (fromInteger $ denominator rat) + recip = lift1ERMachineAP recip recip + (/) = lift2ERMachineAP (/) (/) + +instance (B.ERRealBase b, Real b) => Real (ERMachineAP b) + where + toRational = op1ERMachineAP toRational toRational + +instance (B.ERRealBase b, RealFrac b) => RealFrac (ERMachineAP b) + where + properFraction (ERMachineAPMachineDouble g d) = + (a, ERMachineAPMachineDouble g remainder) + where + (a,remainder) = properFraction d + properFraction (ERMachineAPB b) = + (a, ERMachineAPB remainder) + where + (a,remainder) = properFraction b + + +instance (B.ERRealBase b) => B.ERRealBase (ERMachineAP b) + where + defaultGranularity _ = (B.defaultGranularity (0 :: b)) + getApproxBinaryLog = + op1ERMachineAP doubleBinaryLog B.getApproxBinaryLog + where + doubleBinaryLog d + | d < 0 = + error $ "ERMachineAP: getApproxBinaryLog: negative argument " ++ show d + | d == 0 = EI.MinusInfinity + | d >= 1 = + fromInteger $ intLog 2 $ ceiling d + | d < 1 = + negate $ fromInteger $ intLog 2 $ ceiling $ recip d + getGranularity (ERMachineAPB b) = B.getGranularity b + getGranularity (ERMachineAPMachineDouble gr _) = gr + setMinGranularity gran (ERMachineAPMachineDouble g d) + | gran > doubleDigits = + ERMachineAPB $ B.setMinGranularity gran $ B.fromDouble d + | otherwise = + ERMachineAPMachineDouble gran d + setMinGranularity gran (ERMachineAPB b) = ERMachineAPB $ B.setMinGranularity gran b + setGranularity gran (ERMachineAPMachineDouble g d) + | gran > doubleDigits = + ERMachineAPB $ B.setGranularity gran $ B.fromDouble d + | otherwise = + ERMachineAPMachineDouble gran d + setGranularity gran (ERMachineAPB b) + | gran <= doubleDigits = + ERMachineAPMachineDouble gran $ B.toDouble b + | otherwise = + ERMachineAPB $ B.setGranularity gran b + getMaxRounding _ = + error "ERMachineAP: getMaxRounding not implemented yet" + isERNaN = op1ERMachineAP isNaN B.isERNaN + erNaN = B.fromDouble (0/0) + isPlusInfinity = + op1ERMachineAP (== 1/0) B.isPlusInfinity + plusInfinity = B.fromDouble $ 1/0 + fromDouble d = + ERMachineAPMachineDouble (B.defaultGranularity (0 :: b)) d + toDouble = op1ERMachineAP id B.toDouble + fromFloat f = + ERMachineAPMachineDouble (B.defaultGranularity (0 :: b)) $ + fromRational $ toRational f + toFloat = op1ERMachineAP (fromRational . toRational) B.toFloat + showDiGrCmp = showERMachineAP + diff --git a/src/Data/Number/ER/Real/Base/Float.hs b/src/Data/Number/ER/Real/Base/Float.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/Real/Base/Float.hs @@ -0,0 +1,508 @@ +{-# LANGUAGE DeriveDataTypeable #-} +{-| + Module : Data.Number.ER.Real.Base + Description : arbitrary precision floating point numbers + Copyright : (c) Michal Konecny + License : LGPL + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : portable + + This module defines an arbitrary precision floating point type and + its operations. It should be viewed more abstractly as an instance + of 'B.ERRealBase' when used as interval endpoints. +-} +module Data.Number.ER.Real.Base.Float +( + ERFloat +) +where + +import qualified Data.Number.ER.ExtendedInteger as EI +import Data.Number.ER.PlusMinus +import Data.Number.ER.Misc +import Data.Number.ER.BasicTypes +import qualified Data.Number.ER.Real.Base as B + +import Data.Ratio + +import Data.Typeable +import Data.Generics.Basics +import Data.Binary +-- import BinaryDerive + +--debugMsg = unsafePrint +debugMsg msg = id + + +{-| +A floating point number with a given but arbitrary precision represented by its 'Granularity'. + + * base: 2. + + * granularity specifies the bit-size of both the significand and the exponent + + * special values: NaN, signed Infinity and signed Zero + + * no denormalised numbers + + * operations unify the granularity of their operands to the maximum 'Granularity' + + * Rounding is always towards +Infinity. + For field operations, the rounded result is as close as possible to the exact result. +-} +data ERFloat = + ERFloatNaN -- any number / bottom + { + apfltGran :: Granularity -- >0 + } + | ERFloatInfty + { + apfltGran :: Granularity, -- >0 + apfltSign :: PlusMinus + } + | ERFloatZero + { + apfltGran :: Granularity, -- >0 + apfltSign :: PlusMinus + } + | ERFloat + { + -- represents: + -- sign * (1 + (mant/2^gran)) * (2 ^ exp) + apfltGran :: Granularity, -- >0 granularity + apfltSign :: PlusMinus, + apfltMant :: Integer, -- 0 .. (2^gran - 1) + apfltExp :: Integer -- -2^gran..2^gran + } + deriving (Typeable, Data) + +zero = ERFloatZero 10 Plus + +{- the following has been generated by BinaryDerive -} +instance Binary ERFloat where + put (ERFloatNaN a) = putWord8 0 >> put a + put (ERFloatInfty a b) = putWord8 1 >> put a >> put b + put (ERFloatZero a b) = putWord8 2 >> put a >> put b + put (ERFloat a b c d) = putWord8 3 >> put a >> put b >> put c >> put d + get = do + tag_ <- getWord8 + case tag_ of + 0 -> get >>= \a -> return (ERFloatNaN a) + 1 -> get >>= \a -> get >>= \b -> return (ERFloatInfty a b) + 2 -> get >>= \a -> get >>= \b -> return (ERFloatZero a b) + 3 -> get >>= \a -> get >>= \b -> get >>= \c -> get >>= \d -> return (ERFloat a b c d) + _ -> fail "no parse" +{- the above has been generated by BinaryDerive -} + + +{-| normalisation + + * ensures that the components are within their regions + + * possibly turning the number into a zero or infinity +-} +normaliseERFloat :: ERFloat -> ERFloat +normaliseERFloat flt@(ERFloat gr s m e) + | m < 0 = + normaliseERFloat $ + ERFloat gr s (2*m + grmax) (e - 1) + | m >= grmax = + normaliseERFloat $ + ERFloat gr s ((m - grmax + (rndCorr s)) `div` 2) (e + 1) + | e > grmax = + case s of + Plus -> ERFloatInfty gr Plus + Minus -> minERFloat gr -- round upwards! + | e < -grmax = + case s of + Plus -> ulpERFloat gr -- round upwards! + Minus -> ERFloatZero gr Minus + | otherwise = flt + where + grmax = 2^gr +normaliseERFloat flt = flt + +ulpERFloat gr = + ERFloat gr Plus 0 (-2^gr) + +minERFloat gr = + ERFloat gr Minus (grmax - 1) grmax + where + grmax = 2^gr + +maxERFloat gr = + ERFloat gr Plus (grmax - 1) grmax + where + grmax = 2^gr + +rndCorr Plus = 1 +rndCorr Minus = 0 + +increaseERFloatExp e flt@(ERFloat gr s m eOld) = + ERFloat gr s mNew e + where + mNew = + -grmax + ((m + grmax + (rndCorr s) * (ediff - 1)) `div` ediff) + -- .^^^^^^^^^^^^^^^^^^^^^^^^^ round upwards + grmax = 2^gr + ediff = 2^(e - eOld) -- assuming e >= eOld +increaseERFloatExp _ flt = flt + +decreaseERFloatExp e flt@(ERFloat gr s m eOld) = + ERFloat gr s mNew e + where + mNew = + -grmax + ediff * (m + grmax) + grmax = 2^gr + ediff = 2^(eOld - e) -- assuming e <= eOld +decreaseERFloatExp _ flt = flt + + +apFloatExponent :: ERFloat -> EI.ExtendedInteger + +apFloatExponent (ERFloatInfty _ _) = EI.PlusInfinity +apFloatExponent (ERFloatZero _ _) = EI.MinusInfinity +apFloatExponent (ERFloatNaN _) = EI.PlusInfinity -- includes infinity +apFloatExponent flt = EI.Finite $ apfltExp flt + + +setERFloatGranularity :: + Granularity -> ERFloat -> ERFloat +setERFloatGranularity gr flt@(ERFloat oldGr s m e) + | gr > 0 = + normaliseERFloat $ ERFloat gr s newM e + | otherwise = + flt + where + newM = + (m * (2^gr) + + ((rndCorr s)*(2^oldGr - 1))) -- round upwards! + `div` (2^oldGr) +setERFloatGranularity gr f = f { apfltGran = gr } + +setERFloatMinGranularity :: + Granularity -> ERFloat -> ERFloat +setERFloatMinGranularity gr flt + | gr > oldGr = + setERFloatGranularity gr flt + | otherwise = flt + where + oldGr = apfltGran flt + +apfltGranularity = apfltGran + +{-^ see the documentation of 'ERRealBase.getBaseMaxRounding' -} +apfltGetMaxRounding :: + ERFloat -> ERFloat +apfltGetMaxRounding f@(ERFloatNaN _) = f +apfltGetMaxRounding f@(ERFloatInfty _ _) = f +apfltGetMaxRounding (ERFloatZero gr _) = + ERFloat gr Plus 0 (-(2^gr)) +apfltGetMaxRounding (ERFloat gr s m e) = + ERFloat gr Plus 0 (max (e - (toInteger gr)) (-(2^gr))) + +instance Show ERFloat where + show = showERFloat 6 True False + + +showERFloat numDigits showGran showComponents flt = + showERF flt + where + maybeGran gr + | showGran = "{g=" ++ show gr ++ "}" + | otherwise = "" + showERF (ERFloatNaN gr) = "NaN" ++ (maybeGran gr) + showERF (ERFloatZero gr pm) = show pm ++ "0" ++ (maybeGran gr) + showERF (ERFloatInfty gr pm) = show pm ++ "oo" ++ (maybeGran gr) + showERF f@(ERFloat gr s m e) = + decimal ++ (maybeGran gr) ++ maybeComps + where + maybeComps + | showComponents = "{val="++ show (s,m,e) ++ "}" + | otherwise = "" + decimal = + show s + ++ show digit1 ++ "." ++ (concat $ map show $ take numDigits digits) + ++ "E" ++ show dexp + dexp = dexpBound - zerosCount + digit1 : digits = + drop zerosCount preDigits + dexpBound -- upper bound of dexp: f/10^dexpBound < 1 + | e > 0 = intLog 10 (2^e) + | e <= 0 = 2 - (intLog 10 (2^(-e))) + (zerosCount, preDigits) = + getDigits 0 $ (abs $ setERFloatGranularity numBinDigits f) / (ten ^^ dexpBound) + ten = setERFloatGranularity numBinDigits 10 + numBinDigits = 4 * numDigits + getDigits prevZeros ff + | digit == 0 = (zerosCount, digit : digits) + | otherwise = (prevZeros, digit : digits) + where + digit :: Integer + digit = truncate ff + (zerosCount, digits) = + getDigits zerosNow ((ff - (fromInteger digit)) * ten) + zerosNow + | digit == 0 = prevZeros + 1 + | otherwise = 0 + + +{- + Beware: cannot use List.elem with ERFloat because of + the intensional nature of Eq (eg ERFloatNaN /= ERFloatNaN). +-} +instance Eq ERFloat where + (ERFloatNaN _) == _ = + False + -- error "cannot compare NaN" + _ == (ERFloatNaN _) = + False + -- error "cannot compare NaN" + (ERFloatZero _ _) == (ERFloatZero _ _) = True + (ERFloatInfty _ pm1) == (ERFloatInfty _ pm2) = (pm1 == pm2) + f1@(ERFloat gr1 s1 m1 e1) == f2@(ERFloat gr2 s2 m2 e2) + | gr1 < gr2 = + (setERFloatGranularity gr2 f1) == f2 + | gr1 > gr2 = + f1 == (setERFloatGranularity gr1 f2) + | otherwise = + s1 == s2 && m1 == m2 && e1 == e2 + _ == _ = False + +isERFloatNaN (ERFloatNaN _) = True +isERFloatNaN _ = False + +instance Ord ERFloat where + {- compare NaN -} + compare _ (ERFloatNaN _) = + error "ERFloat: comparing NaN - aborting" + compare (ERFloatNaN _) _ = + error "ERFloat: comparing NaN - aborting" + {- compare infty -} + compare (ERFloatInfty gr1 pm1) (ERFloatInfty gr2 pm2) = + compare pm1 pm2 + compare _ (ERFloatInfty _ Plus) = LT + compare _ (ERFloatInfty _ Minus) = GT + compare (ERFloatInfty _ Plus) _ = GT + compare (ERFloatInfty _ Minus) _ = LT + {- compare zero -} + compare (ERFloatZero gr1 pm1) (ERFloatZero gr2 pm2) = EQ + compare (ERFloatZero _ _) (ERFloat _ Plus _ _) = LT + compare (ERFloatZero _ _) (ERFloat _ Minus _ _) = GT + compare (ERFloat _ Minus _ _) (ERFloatZero _ _) = LT + compare (ERFloat _ Plus _ _) (ERFloatZero _ _) = GT + {- compare regular -} + compare (ERFloat _ Minus _ _) (ERFloat _ Plus _ _) = LT + compare (ERFloat _ Plus _ _) (ERFloat _ Minus _ _) = GT + compare (ERFloat gr1 Plus m1 e1) (ERFloat gr2 _ m2 e2) + | e1 < e2 = LT + | e1 > e2 = GT + | gr1 == gr2 = compare m1 m2 + | otherwise = compare ((2^gr2)*m1) ((2^gr1)*m2) + compare f1@(ERFloat _ Minus _ _) f2@(ERFloat _ _ _ _) = + compare (-f2) (-f1) + +instance Num ERFloat where + fromInteger n + | n == 0 = ERFloatZero (B.defaultGranularity zero) Plus + | n < 0 = + normaliseERFloat $ ERFloat gr Minus m e + | otherwise = + normaliseERFloat $ ERFloat gr Plus m e + where + gr = fromInteger e + e = max (toInteger (B.defaultGranularity zero)) $ (intLog 2 $ abs n) - 1 + m = (abs n) - 2^gr + abs f@(ERFloatNaN _) = f + abs f = f { apfltSign = Plus } + signum f@(ERFloatNaN _) = f + signum (ERFloatZero gr Plus) = setERFloatMinGranularity gr 1 + signum (ERFloatZero gr Minus) = setERFloatMinGranularity gr (-1) + signum (ERFloatInfty gr Plus) = setERFloatMinGranularity gr 1 + signum (ERFloatInfty gr Minus) = setERFloatMinGranularity gr (-1) + signum flt = + case apfltSign flt of { Plus -> 1; Minus -> -1 } + negate (ERFloat gr s m e) = ERFloat gr (signNeg s) m e + negate (ERFloatZero gr s) = ERFloatZero gr (signNeg s) + negate (ERFloatInfty gr s) = ERFloatInfty gr (signNeg s) + negate f@(ERFloatNaN _) = f + {- addition -} + f1 + f2 -- ensure equal granularity: + | gr1 > gr2 = f1 + (setERFloatGranularity gr1 f2) + | gr1 < gr2 = (setERFloatGranularity gr2 f1) + f2 + where + gr1 = apfltGran f1 + gr2 = apfltGran f2 + f@(ERFloatNaN _) + _ = f + _ + f@(ERFloatNaN _) = f + (ERFloatZero _ _) + f = f + f + (ERFloatZero _ _) = f + (ERFloatInfty gr Plus) + (ERFloatInfty _ Minus) = + debugMsg ("ERFloat: infty - infty -> NaN\n") $ + ERFloatNaN gr + (ERFloatInfty gr Minus) + (ERFloatInfty _ Plus) = + debugMsg ("ERFloat: -infty + infty -> NaN\n") $ + ERFloatNaN gr + f@(ERFloatInfty gr s) + _ = f + _ + f@(ERFloatInfty gr s) = f + f1@(ERFloat gr s1 m1 e1) + f2@(ERFloat _ s2 m2 e2) + -- equalise the exponents: + | e1 < e2 = f1 + (decreaseERFloatExp e1 f2) + | e1 > e2 = (decreaseERFloatExp e2 f1) + f2 + -- ensure positive comes before negative: + | s1 == Minus && s2 == Plus = + f2 + f1 + -- opposite signs: + | s1 == Plus && s2 == Minus && m1 == m2 = + ERFloatZero gr Plus + | s1 == Plus && s2 == Minus && m1 > m2 = + normaliseERFloat $ + ERFloat gr s1 (m1 - m2 - 2^gr) e1 + | s1 == Plus && s2 == Minus && m1 < m2 = + normaliseERFloat $ + ERFloat gr s2 (m2 - m1 - 2^gr) e1 + -- equal signs: + | otherwise = + normaliseERFloat $ + ERFloat gr s1 (m1 + m2 + 2^gr) e1 + {- multiplication -} + -- ensure equal granularity: + f1 * f2 + | gr1 > gr2 = f1 * (setERFloatGranularity gr1 f2) + | gr1 < gr2 = (setERFloatGranularity gr2 f1) * f2 + where + gr1 = apfltGran f1 + gr2 = apfltGran f2 + -- NaN: + f@(ERFloatNaN _) * _ = f + _ * f@(ERFloatNaN _) = f + -- Infty + (ERFloatInfty gr _) * (ERFloatZero _ _) = + debugMsg ("ERFloat: infty * 0 -> NaN\n") $ + ERFloatNaN gr + (ERFloatZero gr _) * (ERFloatInfty _ _) = + debugMsg ("ERFloat: 0 * infty -> NaN\n") $ + ERFloatNaN gr + f * (ERFloatInfty gr s) = ERFloatInfty gr $ signMult s (apfltSign f) + (ERFloatInfty gr s) * f = ERFloatInfty gr $ signMult s (apfltSign f) + -- Zero + (ERFloatZero gr s) * f = ERFloatZero gr $ signMult s (apfltSign f) + f * (ERFloatZero gr s) = ERFloatZero gr $ signMult s (apfltSign f) + -- regular + f1@(ERFloat gr s1 m1 e1) * f2@(ERFloat _ s2 m2 e2) = + normaliseERFloat $ + ERFloat gr s mNew (e1 + e2) + where + s = signMult s1 s2 + mNew = + m1 + m2 + + ((m1 * m2 + (rndCorr s) * (2^gr - 1)) + `div` 2^gr) + +instance Fractional ERFloat where + fromRational rat = +-- error "ERFloat: fromRational cannot be implemented reliably: use apfloatFromRational instead" + (fromInteger $ numerator rat) + / (fromInteger $ denominator rat) + f1 / f2 + | gr1 > gr2 = f1 / (setERFloatGranularity gr1 f2) + | gr1 < gr2 = (setERFloatGranularity gr2 f1) / f2 + where + gr1 = apfltGran f1 + gr2 = apfltGran f2 + f@(ERFloatNaN _) / _ = f + f1 / f2 = + case apfltSign f1 of + Plus -> f1 * (recip f2) + Minus -> (- f1) * (recip (- f2)) -- rounding upwards! + recip f@(ERFloatNaN _) = f + recip (ERFloatZero gr s) = ERFloatInfty gr s + recip (ERFloatInfty gr s) = ERFloatZero gr s + recip (ERFloat gr s m e) = + normaliseERFloat $ + ERFloat gr s mNew (-e) + where + mNew = + (- grmax * m + + (rndCorr s) * (grmax + m -1)) -- rounding upwards! + `div` + (grmax + m) + grmax = 2^gr + + +apfloatFromRational :: + Granularity -> Rational -> ERFloat +apfloatFromRational gran rat = + (setERFloatMinGranularity gran (fromInteger $ numerator rat)) + / (fromInteger $ denominator rat) + + + +instance Real ERFloat where + toRational (ERFloat gr s m e) = + case s of + Plus -> r + Minus -> -r + where + r = (eOn2R) * (1 + mR/(grOn2R)) + mR = toRational m + eOn2R = toRational $ 2 ^^ e + grOn2R = toRational $ 2 ^ gr + toRational (ERFloatZero _ _) = 0 + toRational f = + error $ "cannot covert " ++ show f ++ " to a rational" + +instance RealFrac ERFloat where + properFraction (ERFloatNaN _) = + error "no integral part in ERFloatNaN" + properFraction (ERFloatZero _ _) = + (0, 0) + properFraction (ERFloatInfty _ _) = + error "no integral part in ERFloatInfty" + properFraction f@(ERFloat gr s m e) + | e < 0 = (0, f) + | s == Plus = + (n, frac) + | s == Minus = + (-n, frac) + where + n = fromInteger $ 2^e + (m*(2^e) `div` 2^gr) + frac = f - (fromInteger $ toInteger n) + + +instance B.ERRealBase ERFloat + where + defaultGranularity _ = 10 + getApproxBinaryLog = apFloatExponent + getGranularity = apfltGran + setMinGranularity = setERFloatMinGranularity + setGranularity = setERFloatGranularity + getMaxRounding = apfltGetMaxRounding + isERNaN (ERFloatNaN _) = True + isERNaN _ = False + erNaN = ERFloatNaN (B.defaultGranularity zero) + isPlusInfinity (ERFloatInfty _ Plus) = True + isPlusInfinity _ = False + plusInfinity = ERFloatInfty (B.defaultGranularity zero) Plus + fromDouble d + | isNaN d = ERFloatNaN (B.defaultGranularity zero) + | otherwise = (fromRational . toRational) d + toDouble (ERFloatInfty _ s) = signToNum s * (1/0) + toDouble (ERFloatNaN _) = 0/0 + toDouble flt = + (fromInteger $ numerator rat) / (fromInteger $ denominator rat) + where + rat = toRational flt + fromFloat f + | isNaN f = ERFloatNaN (B.defaultGranularity zero) + | otherwise = (fromRational . toRational) f + toFloat (ERFloatInfty _ s) = signToNum s * (1/0) + toFloat (ERFloatNaN _) = 0/0 + toFloat flt = + (fromInteger $ numerator rat) / (fromInteger $ denominator rat) + where + rat = toRational flt + showDiGrCmp = showERFloat + diff --git a/src/Data/Number/ER/Real/Base/MachineDouble.hs b/src/Data/Number/ER/Real/Base/MachineDouble.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/Real/Base/MachineDouble.hs @@ -0,0 +1,87 @@ +{-# INCLUDE #-} +{-# LANGUAGE ForeignFunctionInterface #-} +{-| + Module : Data.Number.ER.Real.Base.MachineDouble + Description : enabling Double's as interval endpoints + Copyright : (c) Michal Konecny + License : LGPL + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : non-portable (requires fenv.h) + + Make 'Double' an instance of 'B.ERRealBase' as much as possible. +-} +module Data.Number.ER.Real.Base.MachineDouble +( + initMachineDouble +) +where + +import qualified Data.Number.ER.Real.Base as B +import qualified Data.Number.ER.ExtendedInteger as EI +import Data.Number.ER.Misc + +import Foreign.C + +{- + The following section is taken from Oleg Kiselyov's email + http://www.haskell.org/pipermail/haskell/2005-October/016574.html +-} + +type FP_RND_T = CInt -- fenv.h + +eFE_TONEAREST = 0 +eFE_DOWNWARD = 0x400 +eFE_UPWARD = 0x800 +eFE_TOWARDZERO = 0xc00 + +foreign import ccall "fenv.h fegetround" fegetround + :: IO FP_RND_T + +foreign import ccall "fenv.h fesetround" fesetround + :: FP_RND_T -> IO FP_RND_T +{- end of Oleg's code -} + +{-| + Set machine floating point unit to the upwards-directed rounding + mode. + + This procedure has to be executed before using 'Double' + as a basis for interval and polynomial arithmetic defined in this package. +-} +initMachineDouble :: IO () +initMachineDouble = + do + currentRndMode <- fegetround + case currentRndMode == eFE_UPWARD of + True -> + putStrLn "initMachineDouble: already rounding upwards" + False -> + do + fesetround eFE_UPWARD + putStrLn "initMachineDouble: switched to upwards rounding" + +instance B.ERRealBase Double + where + defaultGranularity _ = 53 + getApproxBinaryLog f + | f == 0 = + EI.MinusInfinity + | otherwise = + intLog 2 (abs $ ceiling f) + getGranularity _ = 53 + setMinGranularity _ = id + setGranularity _ = id + getMaxRounding _ = 0 + isERNaN f = isNaN f + erNaN = 0/0 + isPlusInfinity f = isInfinite f && f > 0 + plusInfinity = 1/0 + fromDouble = fromRational . toRational + toDouble = fromRational . toRational + fromFloat = fromRational . toRational + toFloat = fromRational . toRational + showDiGrCmp _numDigits _showGran _showComponents f = show f + + diff --git a/src/Data/Number/ER/Real/Base/Rational.hs b/src/Data/Number/ER/Real/Base/Rational.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/Real/Base/Rational.hs @@ -0,0 +1,242 @@ +{-# LANGUAGE DeriveDataTypeable #-} +{-| + Module : Data.Number.ER.Real.Base.Rational + Description : rational numbers with infinities + Copyright : (c) Michal Konecny + License : LGPL + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : portable + + Unlimited size rational numbers extended with signed infinities and NaN. + + These can serve as endpoints of 'Data.Number.ER.Real.Approx.Interval.ERInterval'. + + To be imported qualified, usually with prefix ERAT. +-} +module Data.Number.ER.Real.Base.Rational +( + ExtendedRational(..) +) +where + +import Prelude hiding (isNaN) + +import qualified Data.Number.ER.Real.Base as B +import qualified Data.Number.ER.ExtendedInteger as EI +import Data.Number.ER.PlusMinus +import Data.Number.ER.Misc + +import Data.Ratio +import Data.Typeable +import Data.Generics.Basics + +import Data.Binary + +data ExtendedRational = + NaN + | Infinity PlusMinus + | Finite Rational + deriving (Typeable, Data) + +{- the following has been generated by BinaryDerive -} +instance Binary ExtendedRational where + put NaN = putWord8 0 + put (Infinity a) = putWord8 1 >> put a + put (Finite a) = putWord8 2 >> put a + get = do + tag_ <- getWord8 + case tag_ of + 0 -> return NaN + 1 -> get >>= \a -> return (Infinity a) + 2 -> get >>= \a -> return (Finite a) + _ -> fail "no parse" +{- the above has been generated by BinaryDerive -} + +eratSign :: ExtendedRational -> PlusMinus +eratSign NaN = error "ExtendedRational: eratSign: NaN" +eratSign (Infinity s) = s +eratSign (Finite r) + | r < 0 = Minus + | otherwise = Plus + +liftToERational1 :: + (Rational -> Rational) -> + (ExtendedRational -> ExtendedRational) +liftToERational1 f (Finite r) = + Finite (f r) + +liftToERational2 :: + (Rational -> Rational -> Rational) -> + (ExtendedRational -> ExtendedRational -> ExtendedRational) +liftToERational2 f (Finite r1) (Finite r2) = + Finite (f r1 r2) + + +instance Show ExtendedRational + where + show = showERational 6 True False + +showERational numDigits _showGran showComponents = + showER + where + showER NaN = "NaN" + showER (Infinity pm) = + show pm ++ "oo" + showER (Finite r) | r == 0 = + "0" + showER (Finite r) = + decimal + ++ (if showComponents then components else "") + where + components = "{" ++ show r ++ "}" + decimal = + show pm + ++ show digit1 ++ "." ++ (concat $ map show $ take numDigits digits) + ++ "E" ++ show dexp + pm | r < 0 = Minus + | otherwise = Plus + dexp = dexpBound - zerosCount + digit1 : digits = + drop zerosCount preDigits + dexpBound = -- upper bound of dexp: f/10^dexpBound < 1 + 2 + (intLog 10 num) - (intLog 10 dnm) + num = numerator absr + dnm = denominator absr + absr = abs r + (zerosCount, preDigits) = + getDigits 0 $ absr / (10 ^^ dexpBound) + getDigits prevZeros rr + | digit == 0 = (zerosCount, digit : digits) + | otherwise = (prevZeros, digit : digits) + where + digit :: Integer + digit = truncate rr + (zerosCount, digits) = + getDigits zerosNow ((rr - (fromInteger digit)) * 10) + zerosNow + | digit == 0 = prevZeros + 1 + | otherwise = 0 + +instance Eq ExtendedRational where + NaN == _ = + False + -- error "cannot compare NaN" + _ == NaN = + False + -- error "cannot compare NaN" + (Infinity pm1) == (Infinity pm2) = (pm1 == pm2) + (Finite r1) == (Finite r2) = r1 == r2 + _ == _ = False + +isNaN NaN = True +isNaN _ = False + +instance Ord ExtendedRational where + {- compare NaN -} + compare _ NaN = + error "comparing NaN - aborting" + compare NaN _ = + error "comparing NaN - aborting" + {- compare infty -} + compare (Infinity pm1) (Infinity pm2) = + compare pm1 pm2 + compare _ (Infinity Plus) = LT + compare _ (Infinity Minus) = GT + compare (Infinity Plus) _ = GT + compare (Infinity Minus) _ = LT + {- compare regular -} + compare (Finite r1) (Finite r2) = compare r1 r2 + +instance Num ExtendedRational where + fromInteger n = Finite (fromInteger n) + abs NaN = NaN + abs (Infinity _) = Infinity Plus + abs r = liftToERational1 abs r + signum NaN = NaN + signum (Infinity Plus) = 1 + signum (Infinity Minus) = -1 + signum r = liftToERational1 signum r + negate NaN = NaN + negate (Infinity s) = Infinity (signNeg s) + negate (Finite r) = Finite (negate r) + {- addition -} + -- NaN + NaN + _ = NaN + _ + NaN = NaN + -- Infty + (Infinity Plus) + (Infinity Minus) = NaN + (Infinity Minus) + (Infinity Plus) = NaN + (Infinity s) + _ = Infinity s + _ + (Infinity s) = Infinity s + -- regular + r1 + r2 = liftToERational2 (+) r1 r2 + {- multiplication -} + -- NaN + NaN * _ = NaN + _ * NaN = NaN + -- Infty + (Infinity _) * (Finite r) | r == 0 = NaN + (Finite r) * (Infinity _) | r == 0 = NaN + r * (Infinity s) = Infinity $ signMult s (eratSign r) + (Infinity s) * r = Infinity $ signMult s (eratSign r) + -- regular + r1 * r2 = liftToERational2 (*) r1 r2 + +instance Fractional ExtendedRational where + fromRational rat = Finite rat + recip NaN = NaN + recip (Infinity s) = 0 + recip (Finite r) + | r == 0 = Infinity Plus + | otherwise = (Finite $ recip r) + +instance Real ExtendedRational where + toRational (Finite r) = r + toRational r = error $ "cannot convert " ++ show r ++ " to Rational" + +instance RealFrac ExtendedRational where + properFraction (Finite r) = + (a, Finite b) + where + (a,b) = properFraction r + properFraction r = + error $ "ExtendedRational: RealFrac: no integral part in " ++ show r + +instance B.ERRealBase ExtendedRational + where + defaultGranularity _ = 10 + getApproxBinaryLog (Finite r) + | r == 0 = + EI.MinusInfinity + | otherwise = + (intLog 2 (abs $ numerator $ r)) + - + (intLog 2 (abs $ denominator $ r)) + getApproxBinaryLog (Infinity _) = EI.PlusInfinity + getApproxBinaryLog (NaN) = error "RationalBase: getApproxBinaryLog: NaN" + getGranularity _ = 0 + setMinGranularity _ = id + setGranularity _ = id + getMaxRounding _ = 0 + isERNaN = isNaN + erNaN = NaN + isPlusInfinity (Infinity Plus) = True + isPlusInfinity _ = False + plusInfinity = Infinity Plus + fromDouble = fromRational . toRational + toDouble (Infinity Plus) = 1/0 + toDouble (Infinity Minus) = -1/0 + toDouble (NaN) = 0/0 + toDouble (Finite r) = fromRational r + fromFloat = fromRational . toRational + toFloat (Infinity Plus) = 1/0 + toFloat (Infinity Minus) = -1/0 + toFloat (NaN) = 0/0 + toFloat (Finite r) = fromRational r + showDiGrCmp = showERational + + + + diff --git a/src/Data/Number/ER/Real/DefaultRepr.hs b/src/Data/Number/ER/Real/DefaultRepr.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/Real/DefaultRepr.hs @@ -0,0 +1,78 @@ +{-| + Module : Data.Number.ER.Real.DefaultRepr + Description : concise names for default real representations + Copyright : (c) Michal Konecny + License : LGPL + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : non-portable (requires fenv.h) + + This module supplies default instances for the real number classes + defined in "Data.Number.ER.Real.Approx". + + These classes express loosely coupled abstraction layers. + To preserve the intended loose coupling, + please use these definitions only in functions that do not import or export + any real numbers or real functions. +-} +module Data.Number.ER.Real.DefaultRepr +( + initMachineDouble, + B, BM, BAP, BMAP, BR, + RA, IRA +) +where + +--import + +import Data.Number.ER.Real.Base.Float +import Data.Number.ER.Real.Base.Rational + +import Data.Number.ER.Real.Approx.Interval + +--import Data.Number.ER.Real.Base.BigFloatBase +import Data.Number.ER.Real.Base.MachineDouble +import Data.Number.ER.Real.Base.CombinedMachineAP + +type BAP = ERFloat + +{-| + Limited granularity, but sometimes up to 100x faster + than ERFloat! + + !!! to be safe, one has to run 'initMachineDouble' +-} +type BM = Double + +{-| + Use machine 'Double' while the granularity is up to its significant bit length + and when the granularity grows beyond that, use 'ERFloat'. + + !!! to be safe, one has to run 'initMachineDouble' +-} +type BMAP = ERMachineAP BAP + +--type BBF = BigFloat Prec50 -- seems incomplete on 25/Jun/2008 + +{-| very inefficient -} +type BR = ExtendedRational + +{-| + the default base type +-} +--type B = BAP +--type B = BM +type B = BMAP +--type B = BR + +{-| + the default instance of 'Data.Number.ER.Real.Approx.ERApprox' +-} +type RA = ERInterval B + +{-| + the default instance of 'Data.Number.ER.Real.Approx.ERIntApprox' +-} +type IRA = ERInterval B +