packages feed

AERN-Real (empty) → 0.9.0

raw patch · 23 files changed

+4187/−0 lines, 23 filesdep +basedep +binarydep +containerssetup-changed

Dependencies added: base, binary, containers

Files

+ AERN-Real.cabal view
@@ -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+
+ LICENCE view
@@ -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.
+ Setup.lhs view
@@ -0,0 +1,3 @@+#!/usr/bin/env runhaskell+> import Distribution.Simple+> main = defaultMain
+ src/Data/Number/ER.hs view
@@ -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
+ src/Data/Number/ER/BasicTypes.hs view
@@ -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+
+ src/Data/Number/ER/ExtendedInteger.hs view
@@ -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"+
+ src/Data/Number/ER/Misc.hs view
@@ -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)+
+ src/Data/Number/ER/PlusMinus.hs view
@@ -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
+ src/Data/Number/ER/Real.hs view
@@ -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+
+ src/Data/Number/ER/Real/Approx.hs view
@@ -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+        +        +
+ src/Data/Number/ER/Real/Approx/Elementary.hs view
@@ -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+    +    +    +
+ src/Data/Number/ER/Real/Approx/Interval.hs view
@@ -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
+ src/Data/Number/ER/Real/Approx/Sequence.hs view
@@ -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)+       +
+ src/Data/Number/ER/Real/Arithmetic/Elementary.hs view
@@ -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 <http://en.wikipedia.org/wiki/Computing_π>)+    +    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)) +    +
+ src/Data/Number/ER/Real/Arithmetic/Integration.hs view
@@ -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+-}        +        +
+ src/Data/Number/ER/Real/Arithmetic/Newton.hs view
@@ -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
+ src/Data/Number/ER/Real/Arithmetic/Taylor.hs view
@@ -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	+
+ src/Data/Number/ER/Real/Base.hs view
@@ -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
+ src/Data/Number/ER/Real/Base/CombinedMachineAP.hs view
@@ -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+
+ src/Data/Number/ER/Real/Base/Float.hs view
@@ -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+
+ src/Data/Number/ER/Real/Base/MachineDouble.hs view
@@ -0,0 +1,87 @@+{-# INCLUDE <fenv.h> #-}+{-# 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+    +
+ src/Data/Number/ER/Real/Base/Rational.hs view
@@ -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++++
+ src/Data/Number/ER/Real/DefaultRepr.hs view
@@ -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+