AERN-Real 0.9.9 → 0.10.0

raw patch · 41 files changed

+2931/−1698 lines, 41 filesdep +QuickCheckdep +directorydep +filepathdep ~binarydep ~hmpfrPVP ok

version bump matches the API change (PVP)

Dependencies added: QuickCheck, directory, filepath

Dependency ranges changed: binary, hmpfr

API changes (from Hackage documentation)

Files

• AERN-Real.cabal +46/−35
• ChangeLog +10/−0
• demos/Demo.hs +149/−0
• demos/Matrix.hs +385/−0
• demos/Pi.hs +43/−0
• src/Data/Number/ER/BasicTypes.hs +3/−3
• src/Data/Number/ER/BasicTypes/DomainBox.hs +192/−0
• src/Data/Number/ER/BasicTypes/DomainBox/IntMap.hs +207/−0
• src/Data/Number/ER/BasicTypes/ExtendedInteger.hs +125/−0
• src/Data/Number/ER/BasicTypes/PlusMinus.hs +47/−0
• src/Data/Number/ER/BasicTypes/Tests/Generate.hs +92/−0
• src/Data/Number/ER/ExtendedInteger.hs +0/−125
• src/Data/Number/ER/Misc/STM.hs +42/−0
• src/Data/Number/ER/Misc/Tests.hs +54/−0
• src/Data/Number/ER/MiscSTM.hs +0/−31
• src/Data/Number/ER/PlusMinus.hs +0/−47
• src/Data/Number/ER/Real/Approx.hs +209/−135
• src/Data/Number/ER/Real/Approx/Elementary.hs +37/−1
• src/Data/Number/ER/Real/Approx/Interval.hs +245/−286
• src/Data/Number/ER/Real/Approx/OI.hs +56/−0
• src/Data/Number/ER/Real/Approx/Sequence.hs +2/−2
• src/Data/Number/ER/Real/Approx/Tests/Generate.hs +177/−0
• src/Data/Number/ER/Real/Approx/Tests/Properties.hs +266/−0
• src/Data/Number/ER/Real/Approx/Tests/Reporting.hs +167/−0
• src/Data/Number/ER/Real/Approx/Tests/Run.hs +100/−0
• src/Data/Number/ER/Real/Arithmetic/Elementary.hs +83/−54
• src/Data/Number/ER/Real/Arithmetic/LinearSolver.hs +4/−4
• src/Data/Number/ER/Real/Arithmetic/Taylor.hs +26/−10
• src/Data/Number/ER/Real/Base.hs +8/−3
• src/Data/Number/ER/Real/Base/CombinedMachineAP.hs +2/−1
• src/Data/Number/ER/Real/Base/Float.hs +4/−3
• src/Data/Number/ER/Real/Base/MPFR.hs +2/−1
• src/Data/Number/ER/Real/Base/MachineDouble.hs +12/−2
• src/Data/Number/ER/Real/Base/Rational.hs +3/−2
• src/Data/Number/ER/Real/Base/Tests/Generate.hs +90/−0
• src/Data/Number/ER/Real/DomainBox.hs +0/−181
• src/Data/Number/ER/Real/DomainBox/IntMap.hs +0/−196
• tests/Demo.hs +0/−149
• tests/Matrix.hs +0/−384
• tests/Pi.hs +0/−43
• tests/RunERIntervalTests.hs +43/−0

@@ -1,5 +1,5 @@ Name:           AERN-Real-Version:        0.9.9+Version:        0.10.0 Cabal-Version:  >= 1.2 Build-Type:     Simple License:        BSD3@@ -7,7 +7,7 @@ Author:         Michal Konecny (Aston University) Copyright:      (c) 2007-2008 Michal Konecny, Amin Farjudian, Jan Duracz  Maintainer:     mik@konecny.aow.cz-Stability:      experimental+Stability:      beta Category:       Data, Math Synopsis:       arbitrary precision interval arithmetic for approximating exact real numbers Tested-with:    GHC ==6.10.1@@ -25,10 +25,14 @@     .     For an architectural overview, see module "Data.Number.ER.Real".     .-    Simple examples of usage can be found in module @Demo.hs@ in folder @tests@.-Extra-source-files:-    tests/Demo.hs tests/Pi.hs tests/Matrix.hs-Data-files:+    Simple examples of usage can be found in folder @demos@.+    .+    There is a built-in test suite and it can be evoked using+    the module in the folder @tests@.++Extra-Source-Files:+    demos/Demo.hs demos/Pi.hs demos/Matrix.hs+    tests/RunERIntervalTests.hs     ChangeLog  Flag use-hmpfr@@ -38,36 +42,43 @@   hs-source-dirs:  src   if flag(use-hmpfr)       Build-Depends:-        base >= 3, base < 4, containers, binary >= 0.4, html >= 1.0, regex-compat >= 0.71, stm, time, hmpfr == 0.1.3+        base >= 3, base < 4, containers, binary, html >= 1.0, regex-compat >= 0.71, stm, time, QuickCheck == 1.2.0.0, filepath, directory, hmpfr == 0.2       cpp-options: -DUSE_MPFR   else       Build-Depends:-        base >= 3, base < 4, containers, binary >= 0.4, html >= 1.0, regex-compat >= 0.71, stm, time+        base >= 3, base < 4, containers, binary, html >= 1.0, regex-compat >= 0.71, stm, time, QuickCheck == 1.2.0.0, filepath, directory   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.MPFR,-    Data.Number.ER.Real.Base,-    Data.Number.ER.Real.Arithmetic.Elementary,-    Data.Number.ER.Real.Arithmetic.Integration,-    Data.Number.ER.Real.Arithmetic.LinearSolver,-    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.Real.DomainBox,-    Data.Number.ER.Real.DomainBox.IntMap,-    Data.Number.ER.ShowHTML,-    Data.Number.ER.PlusMinus,-    Data.Number.ER.BasicTypes,-    Data.Number.ER.Misc,-    Data.Number.ER.MiscSTM,-    Data.Number.ER.ExtendedInteger- +        Data.Number.ER,+        Data.Number.ER.BasicTypes,+        Data.Number.ER.BasicTypes.DomainBox,+        Data.Number.ER.BasicTypes.DomainBox.IntMap,+        Data.Number.ER.BasicTypes.ExtendedInteger,+        Data.Number.ER.BasicTypes.PlusMinus,+        Data.Number.ER.BasicTypes.Tests.Generate,+        Data.Number.ER.Misc,+        Data.Number.ER.Misc.STM,+        Data.Number.ER.Misc.Tests,+        Data.Number.ER.Real,+        Data.Number.ER.Real.Approx,+        Data.Number.ER.Real.Approx.Elementary,+        Data.Number.ER.Real.Approx.Interval,+        Data.Number.ER.Real.Approx.OI,+        Data.Number.ER.Real.Approx.Sequence,+        Data.Number.ER.Real.Approx.Tests.Generate,+        Data.Number.ER.Real.Approx.Tests.Properties,+        Data.Number.ER.Real.Approx.Tests.Reporting,+        Data.Number.ER.Real.Approx.Tests.Run,+        Data.Number.ER.Real.Arithmetic.Elementary,+        Data.Number.ER.Real.Arithmetic.Integration,+        Data.Number.ER.Real.Arithmetic.LinearSolver,+        Data.Number.ER.Real.Arithmetic.Newton,+        Data.Number.ER.Real.Arithmetic.Taylor,+        Data.Number.ER.Real.Base,+        Data.Number.ER.Real.Base.CombinedMachineAP,+        Data.Number.ER.Real.Base.Float,+        Data.Number.ER.Real.Base.MachineDouble,+        Data.Number.ER.Real.Base.MPFR,+        Data.Number.ER.Real.Base.Rational,+        Data.Number.ER.Real.Base.Tests.Generate,+        Data.Number.ER.Real.DefaultRepr,+        Data.Number.ER.ShowHTML

@@ -1,3 +1,13 @@+0.9.10: 28 July 2009+    * switching to beta status+    * new QuickCheck test suite covering most functionality+    * new support for anti-consistent intervals (eg [2,0])+      (also called directed or improper; using Kaucher arithmetic)+    * new support for inner-rounded interval arithemtic+    * fixed errors in some elementary functions for extreme values+    * fixed performance bug in arctan+    * improved hierarchy of auxiliary modules+ 0.9.9: 23 February 2009     * Small changes needed in other AERN packages:         * New operation for domain boxes: get its dimension.

@@ -0,0 +1,149 @@+{-# LANGUAGE CPP #-}+{-| +    Module      :  Main+    Description :  simple examples of using AERN-Real+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mik@konecny.aow.cz+    Stability   :  experimental+    Portability :  portable++    Simple examples of using AERN-Real+-}+module Main where++import qualified Data.Number.ER.Real as AERN+import Data.Number.ER.Real (ConvergRealSeq(..), convertFuncRA2Seq)++#ifdef USE_MPFR+--type B = AERN.BAP -- use pure Haskell floats+type B = AERN.BMAP -- use combination of double and pure Haskell floats+--type B = AERN.BMPFR -- use MPFR floats+#else+--type B = AERN.BAP -- use pure Haskell floats+type B = AERN.BMAP -- use combination of double and pure Haskell floats+#endif+type RA = AERN.RA B+type IRA = AERN.IRA B+type R = ConvergRealSeq IRA++one :: R+one = 1++two :: R+two = 2++piSeq :: R+piSeq = ConvergRealSeq $ AERN.pi++seqExp = convertFuncRA2Seq $ AERN.exp+seqSine = convertFuncRA2Seq $ AERN.sin+seqCosine = convertFuncRA2Seq $ AERN.cos++main = +    do+    AERN.initialiseBaseArithmetic (0 :: RA)+    putStrLn "****************************"+    putStrLn "Testing interval arithmetic:"+    putStrLn "****************************"+    putStrLn "**** Fractions:"+    putStrLn $+        "(default granularity, show internals) 1/3 =\n  " ++ +        AERN.showApprox 30 True True (1/3 :: RA) +    putStrLn $+        "(granularity 50, show internals) 1/3 =\n  " ++ +        AERN.showApprox 30 True True ((AERN.setGranularityOuter 50 1/3) :: RA) +    putStrLn $+        "(granularity 100, show internals) 1/3 =\n  " ++ +        AERN.showApprox 40 True True ((AERN.setGranularityOuter 100 1/3) :: RA) +    putStrLn $+        "(granularity 100, do not show internals) 1/3 =\n  " ++ +        AERN.showApprox 40 True False ((AERN.setGranularityOuter 100 1/3) :: RA) +    putStrLn $+        "(granularity 100, default show) 1/3 =\n  " ++ +        show ((AERN.setGranularityOuter 100 1/3) :: RA) +    putStrLn "**** Exp:"+    putStrLn $ +        "(effort 5, granularity 50) exp 1 =\n  " ++ +        (show $ AERN.exp 5 (AERN.setGranularityOuter 50 (1::RA)))+    putStrLn $ +        "(effort 10, granularity 50) exp 1 =\n  " ++ +        (show $ AERN.exp 10 (AERN.setGranularityOuter 50 (1::RA)))+    putStrLn $+        "(effort 10, granularity 100) exp 1 =\n  " ++ +        (show $ AERN.exp 10 (AERN.setGranularityOuter 100 (1::RA)))+    putStrLn $ +        "(effort 20, granularity 50) exp 1 =\n  " ++ +        (show $ AERN.exp 20 (AERN.setGranularityOuter 50 (1::RA)))+    putStrLn $+        "(effort 20, granularity 100) exp 1 =\n  " ++ +        (show $ AERN.exp 20 (AERN.setGranularityOuter 100 (1::RA)))+    putStrLn "**** Pi:"+    putStrLn $ +        "(effort 10) pi =\n  " ++ +        (show $ (AERN.pi 10 :: RA))+    putStrLn $ +        "(effort 50) pi =\n  " ++ +        (AERN.showApprox 20 True False $ (AERN.pi 50 :: RA))+    putStrLn $ +        "(effort 100) pi =\n  " ++ +        (AERN.showApprox 35 True False $ (AERN.pi 100 :: RA))+    putStrLn $ +        "(effort 200) pi =\n  " ++ +        (AERN.showApprox 65 True False $ (AERN.pi 200 :: RA))+    putStrLn $ +        "(effort 400) pi =\n  " ++ +        (AERN.showApprox 125 True False $ (AERN.pi 400 :: RA))+    putStrLn "**** Sine:"+    putStrLn $+        "(effort 10, granularity 50) sin 1 =\n  " ++ +        (show $ AERN.sin 10 (AERN.setGranularityOuter 50 (1::RA)))+    putStrLn $+        "(effort 10, granularity 100) sin 1 =\n  " ++ +        (show $ AERN.sin 10 (AERN.setGranularityOuter 100 (1::RA)))+    putStrLn "**** Integration:"+    putStrLn $ +        "(effort 10, granularity 50) integrate exp 0 1 =\n  " ++ +        (show $ AERN.integrateContAdapt_R AERN.exp 10 0 (AERN.setGranularityOuter 50 (1::RA)))+    putStrLn $ +        "(effort 20, granularity 50) integrate exp 0 1 =\n  " ++ +        (show $ AERN.integrateContAdapt_R AERN.exp 20 0 (AERN.setGranularityOuter 50 (1::RA)))+--    putStrLn $ +--        "(effort 30, granularity 50) integrate exp 0 1 =\n  " ++ +--        (show $ AERN.integrateContAdapt_R AERN.exp 30 0 (AERN.setGranularityOuter 50 (1::RA)))+    putStrLn "*****************************"+    putStrLn "Testing convergent sequences:"+    putStrLn "*****************************"+--    putStrLn $ "1 =\n  " ++ show one+--    putStrLn $ "1 + 2 =\n  " ++ (show $ one + two)+    putStrLn "**** Fractions:"+    putStrLn $ +        "(precision 20) 1/3 =\n  " ++ +        (AERN.showConvergRealSeqAuto 20 $ one / 3)+    putStrLn $ +        "(precision 20) 100000000001/300000000000 =\n  " ++ +        (AERN.showConvergRealSeqAuto 20 $ (one + 100000000000)/300000000000 )+    putStrLn $ +        "100000000001/300000000000 =? 1/3:\n  " ++ +        (show $ one/3 == 100000000001/300000000000)+--    putStrLn $ "abs -1 = " ++ (show $ abs (- one))+--    putStrLn $ "neg 2 = " ++ (show $ negate two)+--    putStrLn $ "1 + 2 = " ++ (show $ one + 2)+    putStrLn "**** Elementary:"+    putStrLn $ +        "(precision 30) exp 1 =\n  " ++ +        (AERN.showConvergRealSeqAuto 30 $ seqExp one)+    putStrLn $ +        "(precision 500) pi =\n  " ++ +        (AERN.showConvergRealSeqAuto 500 $ piSeq)+    putStrLn $ +        "(precision 30) cosine(1) =\n  " ++ +        (AERN.showConvergRealSeqAuto 30 $ seqCosine one)    +    putStrLn $+        "(precision 30) sine(1) =\n  " ++ +        (AERN.showConvergRealSeqAuto 30 $ seqSine one)+    putStrLn "**** Integration:"+    putStrLn $ -- very slow for precision > 4+        "(precision 3) integrate exp 0 1 =\n  " ++ +        (AERN.showConvergRealSeqAuto 3 $ AERN.integrateCont AERN.exp 0 one)

@@ -0,0 +1,385 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE DeriveDataTypeable #-}+module Main++where++import qualified Data.Number.ER.Real as AERN+import Data.Number.ER.BasicTypes+import Data.Number.ER.Misc++import Data.Maybe+import qualified Data.List as List+import qualified Data.Map as Map++import qualified Data.Array.IArray as IAr+import qualified Data.Array.MArray as MAr+import qualified Data.Array.ST as STAr+import qualified Data.Ix as Ix+import qualified Data.Array.Base as BAr++import Control.Monad.ST+import GHC.Arr++#ifdef USE_MPFR+type B = AERN.BAP -- use pure Haskell floats+--type B = AERN.BMPFR -- use MPFR floats+#else+type B = AERN.BAP -- use pure Haskell floats+#endif+type RA = AERN.RA B+type IRA = AERN.IRA B++testMatrixN = 100+incrementGran = (+) 50++-- Hilbert 100x100 matrix:+addOneDiag = False+targetPrec = 167 -- approx 50 decimal digits after the point+initialGran = 2050 -- 100x100+--initialGran = 2388 -- 100x100 Norbert's+--initialGran = 750 -- 50x50+--initialGran = 300 -- 10x10++--targetPrec = 34 -- approx 10 decimal digits after the point+--initialGran = 1350+--initialGran = 50 -- 50x50++-- Hilbert matrix + 1:+--addOneDiag = True+--targetPrec = 167 -- approx 50 decimal digits after the point+--initialGran = 200++--targetPrec = 34 -- approx 10 decimal digits after the point+--initialGran = 50++main =+    do+    AERN.initialiseBaseArithmetic (0 :: RA)+    putStrLn $ +          "Inverting the " ++ show n ++ "x" ++ show n ++ " Hilbert matrix " +          ++ "with target binary precision " ++ show targetPrec ++ "..." +--    putStrLn $ +--        "sorted matrix elements = \n" ++ (unlines $ map show elemsSortedByPrec)+    putStrLn $ +        "sum of all elements in inverted matrix = " ++ show (sum elems)+--    putStrLn $ show (Matrix n n rarr)+    where+    n = testMatrixN+    elems = IAr.elems rarr+    elemsSortedByPrec =+        List.sortBy comparePrec elems+        where+        comparePrec a b =+            compare aPrecLO bPrecLO+            where+            aPrecLO = fst $ AERN.bounds $ aHI - aLO+            (aLO, aHI) = AERN.bounds a+            bPrecLO = fst $ AERN.bounds $ bHI - bLO+            (bLO, bHI) = AERN.bounds b+    rarr =+        STAr.runSTArray $+            do+            mInv@(Matrix _ _ rowsInv) <- +                invert testMatrix+--            m <- testMatrix initialGran+--            mUnit@(Matrix _ _ rowsUnit) <- multM m mInv+            return rowsInv+++testMatrix ::+    Granularity -> +    ST s (STMatrix s IRA)+testMatrix gran =+    do+    marr <- MAr.newArray ((1,1),(n,n)) 0+    mapM (updateCell marr) assocsGran+    return $ Matrix n n marr+    where+    assocsGran = map (mapSnd $ AERN.setMinGranularityOuter gran) assocs+    assocs = +--        assocsMini+        assocsHilbert gran n+    assocsMini = +        [((1,1),1),+         ((1,2),3),+         ((2,1),2),+         ((2,2),0)+        ]+    n = testMatrixN+    updateCell marr (ix, el) =+        do+        unsafeMatrixWrite marr n ix el ++assocsHilbert gran n =+    [((i,j), coeff i j)| i <- [1..n], j <- [1..n]]+    where+    coeff i j +        | addOneDiag && i == j = +            1 + oneOverIplusJ+        | otherwise =+            oneOverIplusJ+        where+        oneOverIplusJ =+            recip $ (AERN.setMinGranularityOuter gran $ iRA + jRA + 1)+        iRA = fromInteger $ toInteger i+        jRA = fromInteger $ toInteger j++    +--invert ::+--    Precision ->+--    () ->+invert getMatrix =+    do+    gaussElim getMatrixI+    where+    n = testMatrixN+    getMatrixI gran =+        do+        m <- getMatrix gran+        mI <- addIdentity m+        return mI++gaussElim getMatrix =+    elimWithMinGran initialGran+    where+    elimWithMinGran workingGran =+        do+        mI@(Matrix colN rowN _) <- getMatrix workingGran+        idPerm <- MAr.newListArray (1,rowN) [1..rowN]+        elimAtRow mI 1 idPerm+        where+        elimAtRow mI@(Matrix colN rowN mIarr) i perm =+            do+            success <- ensureNonZeroDiag -- make sure (i,i) is non-zero by permuting+            case success of+                False -> -- failed - all elements contain 0 -> try larger granularity+                    unsafePrint ("failed to divide at granularity " ++ show workingGran) $+                        elimWithMinGran (incrementGran workingGran)+                True ->+                    do+                    normaliseRow+                    eliminateColumn+                    case i == rowN of+                        True -> +                            do+                            mInv <- permuteRowsDropCols perm testMatrixN mI+                            mPrec <- getMatrixPrecision mInv+                            case mPrec >= targetPrec of+                                False -> -- resulting precision insufficient+                                    unsafePrint +                                    ("insufficient precision " ++ show mPrec ++  +                                     " at granularity " ++ show workingGran) $+                                        elimWithMinGran (incrementGran workingGran)+                                True -> +                                    unsafePrint +                                    ("precision " ++ show mPrec ++ +                                     " succeeded at granularity " ++ show workingGran)+                                    return mInv+                        False -> elimAtRow mI (i+1) perm+            where+            ensureNonZeroDiag =+                do+                maybeNonZeroIx <- findNonZeroRow+                case maybeNonZeroIx of+                    Nothing ->+                        return False+                    Just ii ->+                        do+                        case ii > 0 of+                            True -> swap i (i + ii) perm+                            False -> return ()+                        return True+            findNonZeroRow =+                do+                elems <- mapM getElemPerm [(i,rowIx) | rowIx <- [i..rowN]]+                return $ List.findIndex (\e -> not $ 0 `AERN.refines` e) elems+            getElemPerm (colIx,rowIx) =+                do+                rowIxPerm <- unsafePermRead perm rowIx+                unsafeMatrixRead mIarr rowN (colIx, rowIxPerm)++            normaliseRow =+                do+                rowIxPerm <- unsafePermRead perm i+                e <- unsafeMatrixRead mIarr rowN (i, rowIxPerm)+                unsafeMatrixWrite mIarr rowN (i, rowIxPerm) 1+                mapM (divideCellBy e rowIxPerm) [(i+1)..colN]+            divideCellBy e rowIxPerm colIx =+                do+                e2 <- unsafeMatrixRead mIarr rowN (colIx, rowIxPerm)+                unsafeMatrixWrite mIarr rowN (colIx, rowIxPerm) (e2/e)+                +            eliminateColumn =+                do+                iRowPerm <- unsafePermRead perm i+                mapM (eliminateColumnRow iRowPerm) $ [1..(i-1)] ++ [(i+1)..rowN]+            eliminateColumnRow iRowPerm rowIx =+                do+                rowIxPerm <- unsafePermRead perm rowIx+                c <- unsafeMatrixRead mIarr rowN (i, rowIxPerm) -- remember old element for scaling i'th row+                unsafeMatrixWrite mIarr rowN (i,rowIxPerm) 0 -- at column i we set 0+                mapM (eliminateColumnRowColumn iRowPerm rowIxPerm c) [(i+1)..colN]+            eliminateColumnRowColumn iRowPerm rowIxPerm c colIx =+                do+                ei <- unsafeMatrixRead mIarr rowN (colIx, iRowPerm) -- at i'th row+                er <- unsafeMatrixRead mIarr rowN (colIx, rowIxPerm) -- at current row+                unsafeMatrixWrite mIarr rowN (colIx, rowIxPerm) (er - c * ei) -- eliminate by i'th row+               + +swap ::+    Int ->+    Int ->+    (STAr.STUArray s Int Int) ->+    ST s ()+swap i1 i2 perm =+    do+    a1 <- unsafePermRead perm i1+    a2 <- unsafePermRead perm i2+    unsafePermWrite perm i1 a2+    unsafePermWrite perm i2 a1+            ++unsafePermWrite permArr i e =+    do+    BAr.unsafeWrite permArr (i - 1) e+                +unsafePermRead permArr i =+    do+    BAr.unsafeRead permArr (i - 1)+                ++addIdentity ::+    (STMatrix s IRA) ->+    ST s (STMatrix s IRA)+addIdentity (Matrix colN rowN marr) =+    do+--    (_, (colN,rowN)) <- MAr.getBounds marr+    mElems <- MAr.getElems marr+    mIarr <- MAr.newListArray ((1,1),(colN+rowN,rowN)) $ mElems ++ (idElems rowN)+    return $ Matrix (colN + rowN) rowN mIarr+    where+    idElems m =+        1 : (concat $ replicate (m-1) $ (replicate m 0) ++ [1])+++data Matrix marr el =+    Matrix+    {+        mxRowN :: Int,+        mxColN :: Int,+        mxRows :: marr (ColIx,RowIx) el+    }++type ColIx = Int +type RowIx = Int ++type IMatrix el = +    Matrix Array el+    +type STMatrix s el =+    Matrix (STArray s) el+    +instance +    (IAr.IArray marr el,-- IAr.IArray marr (marr Int el), +     Show el) => +    Show (Matrix marr el)+    where+    show (Matrix colN rowN rows) =+        "\nMatrix:\n" ++ +        (concat $ map showCol [1..colN])+        where+--        (_,(colN,rowN)) = IAr.bounds rows+        showCol colIx =+            unlines $+                map showCell [(colIx, rowIx) | rowIx <- [1..rowN]] +        showCell ix@(colIx, rowIx) =+            (show ix) +++            (replicate colIx '.') ++  +            (show $ (IAr.!) rows ix)+    +getMatrixPrecision (Matrix _ _ marr) =+    do+    elems <- MAr.getElems marr+    return $ foldl1 min $ map AERN.getPrecision elems++unsafeMatrixWrite marr rowN (i,j) e =+    do+    BAr.unsafeWrite marr (rowN*(i-1) + j-1) e+--    MAr.writeArray marr (i,j) e++unsafeMatrixRead marr rowN (i,j) =+    do+    BAr.unsafeRead marr (rowN*(i-1) + j-1)+--    MAr.readArray marr (i,j)+    +permuteRowsDropCols ::+    (STAr.STUArray s Int Int) ->+    Int {-^ drop this many first columns -} ->+    (STMatrix s IRA) ->+    ST s (STMatrix s IRA)+permuteRowsDropCols perm dropN (Matrix colN rowN marr) =+    do+--    (_, (colN,rowN)) <- MAr.getBounds marr+    (_, permN) <- MAr.getBounds perm    +    rarr <- MAr.newArray ((1,1),(colN - dropN, permN)) 0+    mapM (copyElem marr rarr rowN) [(colIx, rowIx) | colIx <- [1..colN - dropN], rowIx <- [1..permN]]+    return (Matrix (colN - dropN) permN rarr)+    where+    copyElem marr rarr rowN (colIx, rowIx) =+        do+        permRowIx <- unsafePermRead perm rowIx+        e <- unsafeMatrixRead marr rowN (colIx + dropN, permRowIx)+        unsafeMatrixWrite rarr rowN (colIx, rowIx) e+        +    +addM m1 m2 +    | mxColN m1 == mxColN m2 && mxRowN m1 == mxRowN m2 =+        do+        marr <- MAr.newArray ((1,1),(colN, rowN)) 0+        mapM (addCell marr) [(c,r) | c <- [1..colN], r <- [1..rowN]]+        return (Matrix colN rowN marr)   +    | otherwise =+        error "Matrix: addM mismatch"+    where+    colN = mxColN m1+    rowN = mxRowN m1+    marr1 = mxRows m1+    marr2 = mxRows m2+    addCell marr (colIx, rowIx) =+        do+        elem1 <- unsafeMatrixRead marr1 rowN (colIx, rowIx)+        elem2 <- unsafeMatrixRead marr2 rowN (colIx, rowIx)+        unsafeMatrixWrite marr rowN (colIx, rowIx) (elem1 + elem2)++multM m1 m2 +    | colN1 == rowN2 =+        do+        marr <- MAr.newArray ((1,1),(colN, rowN)) 0+        mapM (multCell marr) [(c,r) | c <- [1..colN], r <- [1..rowN]]+        return (Matrix colN rowN marr)   +    | otherwise =+        error "Matrix: multM mismatch"+    where+    colN1 = mxColN m1+    rowN1 = mxRowN m1+    colN2 = mxColN m2+    rowN2 = mxRowN m2+    colN = colN2+    rowN = rowN1+    marr1 = mxRows m1+    marr2 = mxRows m2+    multCell marr (colIx, rowIx) =+        do+        elems1 <- mapM (getCell1 rowIx) [1..colN1]+        elems2 <- mapM (getCell2 colIx) [1..rowN2]+        unsafeMatrixWrite marr rowN (colIx, rowIx) (sum $ zipWith (*) elems1 elems2)+    getCell1 rowIx colIx =+        do+        unsafeMatrixRead marr1 rowN1 (colIx, rowIx)+    getCell2 rowIx colIx =+        do+        unsafeMatrixRead marr2 rowN2 (colIx, rowIx)+        

@@ -0,0 +1,43 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE DeriveDataTypeable #-}+module Main++where++import qualified Data.Number.ER.Real as AERN+import Data.Number.ER.Real (ConvergRealSeq(..), convertFuncRA2Seq)+import Data.Number.ER.BasicTypes+import Data.Number.ER.Misc++import Data.Maybe++#ifdef USE_MPFR+--type B = AERN.BMPFR -- use MPFR floats+type B = AERN.BAP -- use pure Haskell floats+#else+type B = AERN.BAP -- use pure Haskell floats+--type B = AERN.BMAP -- use combination of double and pure Haskell floats+#endif+type RA = AERN.RA B+type IRA = AERN.IRA B+++decimalPrec = 1000+binaryPrec =+    fromInteger $ toInteger $+    snd $ AERN.integerBounds $+        (fromInteger decimalPrec :: RA) * (AERN.log 100 10)/(AERN.log 100 2)++main =+    do+    AERN.initialiseBaseArithmetic (0 :: RA)+    putStrLn $ +        show decimalPrec +        ++ " decimal digits of pi = \n" +        ++ (AERN.showConvergRealSeqAuto binaryPrec pi)+    where+    pi :: ConvergRealSeq RA+    pi = ConvergRealSeq AERN.pi+

@@ -1,6 +1,6 @@ {-|     Module      :  Data.Number.ER.BasicTypes-    Description :  generic types for exact real number processing +    Description :  auxiliary types for exact real number processing      Copyright   :  (c) Michal Konecny     License     :  BSD3 @@ -8,12 +8,12 @@     Stability   :  experimental     Portability :  portable -    generic types for exact real number processing+    auxiliary types for exact real number processing -} module Data.Number.ER.BasicTypes  where -import qualified Data.Number.ER.ExtendedInteger as EI+import qualified Data.Number.ER.BasicTypes.ExtendedInteger as EI  {-|     Precision represents an upper bound on the measure of 

@@ -0,0 +1,192 @@+{-# LANGUAGE MultiParamTypeClasses  #-}+{-# LANGUAGE FunctionalDependencies  #-}+{-|+    Module      :  Data.Number.ER.BasicTypes.DomainBox+    Description :  portions of many-dimensional domains   +    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mik@konecny.aow.cz+    Stability   :  experimental+    Portability :  portable++    Abstractions of the 'Box' datatype, often used to represent+    sections of multi-dimensional function domains.+    +    To be imported qualified, usually with prefix DBox.+    +    VariableID(..) and DomainBox +    are usually imported separately and not qualified.+-}+module Data.Number.ER.BasicTypes.DomainBox+(+    VariableID(..),+    getNVars,+    DomainBox(..),+    DomainBoxMappable(..),+    DomainIntBox(..)+)+where++import Data.Number.ER.BasicTypes++import qualified Data.Set as Set+import qualified Data.Map as Map++import Prelude hiding (lookup)+++{-| +    A class abstracting a type of variable identifiers +    for axes in function domains, polynomials etc.+-}+class (Ord varid) => VariableID varid+    where+    newVarID :: Set.Set varid -> varid+    defaultVar :: varid+    defaultVar = newVarID Set.empty+    showVar :: varid -> String++getNVars :: (VariableID varid) => Int -> [varid]+getNVars n =+    aux (Set.empty) n+    where+    aux prevVars n +        | n > 0 = +            aux (Set.insert (newVarID prevVars) prevVars) (n - 1)+        | n == 0 =+            Set.toAscList $ prevVars ++{-|+    A class abstracting a type of many-dimensional points, intervals+    or anything indexed by a subset of dimensions.+    +    More generally, this class abstracts most of 'Data.Map.Map'.+-}+class (VariableID varid) => DomainBox box varid val+    | box -> varid val, varid val -> box+    where+    noinfo :: box+    isNoinfo :: box -> Bool+    size :: box -> Int+    {-| constructor using 'defaultVar' -}+    unary :: val -> box+    singleton :: varid -> val -> box+    toList :: box -> [(varid, val)]+    fromList :: [(varid, val)] -> box+    toAscList :: box -> [(varid, val)]+    fromAscList :: [(varid, val)] -> box+--    toMap :: box -> Map.Map varid val+--    fromMap :: Map.Map varid val -> box+    compare :: (val -> val -> Ordering) -> box -> box -> Ordering+    adjust :: (val -> val) -> varid -> box -> box+    insert :: varid -> val -> box -> box+    insertWith :: (val -> val -> val) -> varid -> val -> box -> box+    delete :: varid -> box -> box+    member :: varid -> box -> Bool+    notMember :: varid -> box -> Bool+    union :: box -> box -> box+    unionWith :: (val -> val -> val) -> box -> box -> box+    keys :: box -> [varid]+    elems :: box -> [val]+    filter :: (val -> Bool) -> box -> box+    fold :: (val -> a -> a) -> a -> box -> a+    foldWithKey :: (varid -> val -> a -> a) -> a -> box -> a+    {-| +        for all variables that appear in both boxes,+        apply the function and add the result to the list +     -}+    zipWith :: (val -> val -> a) -> box -> box -> [(varid, a)] +    {-| +        For all variables that appear in either of the two boxes,+        apply the function and add the result to the list.+        +        Supply the default value when the variable is missing from either box. +     -}+    zipWithDefault :: val -> (val -> val -> a) -> box -> box -> [(varid, a)] +    {-| +        For all variables that appear in the first box,+        apply the function and add the result to the list.+        +        Supply the default value when the variable is missing from the second box. +     -}+    zipWithDefaultSecond :: val -> (val -> val -> a) -> box -> box -> [(varid, a)] +    findWithDefault :: val -> varid -> box -> val+    {-|+        Pick the extents of a single variable in a domain box.+        If there is no information for this variable, assume the+        variable ranges over the whole real line.+    -}+    lookup ::     +        String {-^ identification of caller location to use in error messages -} ->+        varid ->+        box ->+        val+        +{-|+    A class linking two domain box types that share the+    index type so that boxes of the two types can be+    converted etc.+-}+class (DomainBox box1 varid val1, DomainBox box2 varid val2) => +    DomainBoxMappable box1 box2 varid val1 val2+    where+    map :: (val1 -> val2) -> box1 -> box2+    mapWithKey :: (varid -> val1 -> val2) -> box1 -> box2+    intersectionWith :: (val1 -> val2 -> val1) -> box1 -> box2 -> box1+    difference :: box1 -> box2 -> box1 ++{-|+    A class abstracting a type of many-dimensional intervals.+-}+class (DomainBox box varid ira) => DomainIntBox box varid ira+    | box -> varid ira, varid ira -> box+    where+    {-|+        Check whether the two domains specify the same+        interval for each variable that they share.+    -}+    compatible ::+        box ->+        box ->+        Bool+    {-|+        Assuming that two domains are compatible, take the+        most information from both of the domains about the+        ranges of variables.+    -}+    unify ::+        String {-^ identification of caller location to use in error messages -} ->+        box ->+        box ->+        box+    {-|+        Find the variable with the largest interval+        and return it together with the default splitting point+        in its domain.+    -}+    bestSplit ::+        box  {-^ box considered for splitting -} ->+        (varid, (ira, ira))+        {-^ variable with widest domain, its domain and default split point -}+    split ::+        box {-^ box to split -} ->+        varid {-^ direction to split in -} ->+        Maybe ira  {-^ point to split the domain of variable @varid@ at, if absent use default -} ->+        (box, box)+    classifyPosition ::+        box {-^ domain @d1@ -} ->+        box {-^ domain @d2@ -} ->+        (Bool, Bool, Bool, Bool) +            {-^ +                Answers to these (mutually exclusive) questions:+                +                * is @d1@ outside and /not/ touching @d2@?+            +                * is @d1@ outside and touching @d2@?+            +                * is @d1@ intersecting and not inside @d2@?+            +                * is @d1@ inside @d2@?+            -}+            

@@ -0,0 +1,207 @@+{-# LANGUAGE MultiParamTypeClasses  #-}+{-# LANGUAGE FlexibleInstances   #-}+{-# LANGUAGE TypeSynonymInstances   #-}+{-|+    Module      :  Data.Number.ER.BasicTypes.DomainBox.IntMap+    Description :  implementation of DomainBox based on Data.IntMap   +    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mik@konecny.aow.cz+    Stability   :  experimental+    Portability :  portable++    A simple implementation of the 'VariableID' and 'DomainBox' classes.+-}+module Data.Number.ER.BasicTypes.DomainBox.IntMap +(+    VarID, Box+)+where++import qualified Data.Number.ER.Real.Approx as RA+import qualified Data.Number.ER.BasicTypes.DomainBox as DBox+import Data.Number.ER.BasicTypes.DomainBox (VariableID(..), DomainBox, DomainBoxMappable, DomainIntBox)++import Data.Number.ER.Misc++import qualified Data.IntMap as IMap+import qualified Data.Set as Set++type VarID = Int+type Box ira = IMap.IntMap ira++instance VariableID VarID+    where+    newVarID prevVars +        | Set.null prevVars = 0+        | otherwise =+            1 + (Set.findMax prevVars)+    showVar v+        | v == 0 = "x"+        | otherwise = "x" ++ show v++instance (Show val) => (DomainBox (Box val) VarID val)+    where+    noinfo = IMap.empty+    isNoinfo = IMap.null+    size = IMap.size+    unary r = IMap.singleton defaultVar r+    singleton = IMap.singleton+    toList = IMap.toList+    fromList = IMap.fromList+    toAscList = IMap.toAscList+    fromAscList = IMap.fromAscList+--    toMap = id+--    fromMap = id+    compare compareVals b1 b2 =+        compareListsWith comparePairs (IMap.toList b1) (IMap.toList b2)+        where+        comparePairs (k1,v1) (k2,v2) =+            compareComposeMany+                [+                    compare k1 k2,+                    compareVals v1 v2+                ]+             +    adjust = IMap.adjust+    insert = IMap.insert+    insertWith = IMap.insertWith+    delete = IMap.delete+    member = IMap.member +    notMember = IMap.notMember+    union = IMap.union +    unionWith = IMap.unionWith +    elems = IMap.elems+    keys = IMap.keys+    filter = IMap.filter+    fold = IMap.fold+    foldWithKey = IMap.foldWithKey+    zipWith f b1 b2 = +        applyF (IMap.toAscList b1) (IMap.toAscList b2)+        where+        applyF [] _ = []+        applyF _ [] = []+        applyF bl1@((k1,v1):rest1) bl2@((k2,v2):rest2) +            | k1 == k2 = +                (k1, f v1 v2) : (applyF rest1 rest2)+            | k1 < k2 = applyF rest1 bl2+            | otherwise = applyF bl1 rest2 +    zipWithDefault defaultValue f b1 b2 = +        applyF (IMap.toAscList b1) (IMap.toAscList b2)+        where+        applyF [] [] = []+        applyF bl1@((k1,v1):rest1) [] =+            (k1, f v1 defaultValue) : (applyF rest1 [])+        applyF [] bl2@((k2,v2):rest2) =+            (k2, f defaultValue v2) : (applyF [] rest2)+        applyF bl1@((k1,v1):rest1) bl2@((k2,v2):rest2) +            | k1 == k2 = +                (k1, f v1 v2) : (applyF rest1 rest2)+            | k1 < k2 = +                (k1, f v1 defaultValue) : (applyF rest1 bl2)+            | otherwise =  +                (k2, f defaultValue v2) : (applyF bl1 rest2)+    zipWithDefaultSecond defaultValue f b1 b2 = +        applyF (IMap.toAscList b1) (IMap.toAscList b2)+        where+        applyF [] _ = []+        applyF bl1@((k1,v1):rest1) [] =+            (k1, f v1 defaultValue) : (applyF rest1 [])+        applyF bl1@((k1,v1):rest1) bl2@((k2,v2):rest2) +            | k1 == k2 = +                (k1, f v1 v2) : (applyF rest1 rest2)+            | k1 < k2 = +                (k1, f v1 defaultValue) : (applyF rest1 bl2)+            | otherwise =  +                applyF bl1 rest2+    findWithDefault = IMap.findWithDefault+    lookup locspec var dom =+        IMap.findWithDefault err var dom+        where+        err =+            error $+                locspec ++ "DomainBox.IntMap lookup: domain box " ++ show dom +                ++ " ignores variable " ++ show var++instance (Show val1, Show val2) => +    (DomainBoxMappable (Box val1) (Box val2) VarID val1 val2)+    where+    map = IMap.map+    mapWithKey = IMap.mapWithKey+    intersectionWith = IMap.intersectionWith+    difference = IMap.difference++instance (RA.ERIntApprox ira) => DomainIntBox (Box ira) VarID ira+    where+    compatible dom1 dom2 =+        foldl (&&) True $ map snd $+            DBox.zipWith RA.equalIntervals dom1 dom2+    unify locspec dom1 dom2+        | DBox.compatible dom1 dom2 =+            IMap.union dom1 dom2+        | otherwise =+            error $+                locspec ++ "incompatible domains " ++ show dom1 ++ " and " ++ show dom2+    bestSplit domB =+        (var, (varDom, pt))+        where+        pt = +            RA.defaultBisectPt varDom+        (_, (varDom, var)) = +            foldl findWidestVar (0, err) $ IMap.toList domB+        err =+            error $ "DomainBox: bestSplit: failed to find a split for " ++ show domB +        findWidestVar (prevWidth, prevRes) (v, d)+            | currWidth `RA.leqSingletons` prevWidth = (prevWidth, prevRes)+            | otherwise = (currWidth, (d, v))+            where+            currWidth = snd $ RA.bounds $ domHI - domLO+            (domLO, domHI) = RA.bounds d+    split domB var maybePt = +        (IMap.insert var varDomL domB, +         IMap.insert var varDomR domB)+        where+        varDomL = varDomLO RA.\/ pt+        varDomR = pt RA.\/ varDomHI+        pt = +            case maybePt of+                Nothing -> varDomMid+                Just pt | pt `RA.refines` varDom -> pt+                Just pt -> +                    error $  +                        "ER.DomainBox.IntMap: split given an invalid split point " +                        ++ show pt ++ " for the domain box " ++ show domB +                        ++ " and split variable " ++ show var +        (varDomLO, varDomMid, varDomHI, _) = RA.exactMiddle varDom+        varDom = DBox.lookup "DomainBox.IntMap: split: " var domB+    classifyPosition dom sdom =    +        (away, touch, intersect, inside)+            where+            (away, touch, inside, intersect) =+                foldl addDimension (True, True, True, False) awayTouchInsides+            addDimension +                    (prevAway, prevTouch, prevInside, prevIntersect) +                    (thisAway, thisTouch, thisInside, thisIntersect) =+                (prevAway && thisAway, +                 (prevTouch || prevAway) && (thisTouch || thisAway) && (prevTouch || thisTouch),+                 prevInside && thisInside,+                 prevIntersect || thisIntersect)+            awayTouchInsides =+                map snd $+                    DBox.zipWith classifyRA dom sdom+            classifyRA d sd =+                (outsideNoTouch, outsideTouch, inside,+                 not (outsideNoTouch || outsideTouch || inside))+                 where+                 outsideNoTouch = sdR < dL || dR < sdL+                 outsideTouch = sdR == dL || dR == sdL+                 inside = sdL =< dL && dR =< sdR+                 (==) = RA.eqSingletons+                 (<) = RA.ltSingletons+                 (=<) = RA.leqSingletons+                 (dL, dR) = RA.bounds d +                 (sdL, sdR) = RA.bounds sd +        ++    

@@ -0,0 +1,125 @@+{-|+    Module      :  Data.Number.ER.BasicTypes.ExtendedInteger+    Description :  integer with infinities +    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    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.BasicTypes.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"+        

@@ -0,0 +1,47 @@+{-# LANGUAGE DeriveDataTypeable   #-}+{-|+    Module      :  Data.Number.ER.BasicTypes.PlusMinus+    Description :  mini sign datatype+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    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.BasicTypes.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

@@ -0,0 +1,92 @@+{-|+    Module      :  Data.Number.ER.BasicTypes.Tests.Generate+    Description :  (testing) generating values for tests+    Copyright   :  (c) 2007-2008 Michal Konecny+    License     :  BSD3++    Maintainer  :  mik@konecny.aow.cz+    Stability   :  experimental+    Portability :  portable+    +    Instances of Arbitrary so that values+    can be randomly generated for use in QuickCheck tests.+-}++module Data.Number.ER.BasicTypes.Tests.Generate +where++import Data.Number.ER.BasicTypes+import Data.Number.ER.BasicTypes.ExtendedInteger+import Data.Number.ER.BasicTypes.PlusMinus++import Test.QuickCheck hiding (two, three)++data Nat10 = Nat10 Int deriving (Show)+data Nat100 = Nat100 Int deriving (Show)++data Ix10 = Ix10 EffortIndex deriving (Show)+data Ix20 = Ix20 EffortIndex deriving (Show)++data Gran100 = Gran100 Granularity deriving (Show)+data Gran1000 = Gran1000 Granularity deriving (Show)++data SmallRatio = SmallRatio Int Int deriving (Show)++instance (Arbitrary Nat10)+    where+    arbitrary =+        do+        ix <- choose (0,10)+        return $ Nat10 ix+    coarbitrary (Nat10 ix) =+        error "ER.BasicTypes.Tests.Generate: coarbitrary not implemented for Nat10"++instance (Arbitrary Nat100)+    where+    arbitrary =+        do+        ix <- choose (0,100)+        return $ Nat100 ix+    coarbitrary (Nat100 ix) =+        error "ER.BasicTypes.Tests.Generate: coarbitrary not implemented for Nat100"++instance (Arbitrary Ix20)+    where+    arbitrary =+        do+        ix <- choose (2,20)+        return $ Ix20 ix+    coarbitrary (Ix20 ix) =+        error "ER.BasicTypes.Tests.Generate: coarbitrary not implemented for Ix20"++instance (Arbitrary Ix10)+    where+    arbitrary =+        do+        ix <- choose (1,10)+        return $ Ix10 ix+    coarbitrary (Ix10 ix) =+        error "ER.BasicTypes.Tests.Generate: coarbitrary not implemented for Ix10"++instance (Arbitrary PlusMinus)+    where+    arbitrary = +        do+        isPlus <- arbitrary+        case isPlus of+            True -> return Plus+            False -> return Minus+    coarbitrary pm =+        error "ER.BasicTypes.Tests.Generate: coarbitrary not implemented for PlusMinus"+    +instance (Arbitrary SmallRatio)+    where+    arbitrary =+        do+        num <- choose (-1000000,1000000)+        denom <- choose (1,1000000)+        return $ SmallRatio num denom+    coarbitrary (SmallRatio num denom) =+        error "ERChebPoly: Generate: Arbitrary: coarbitrary not implemented for SmallRatio"+        +    

@@ -1,125 +0,0 @@-{-|-    Module      :  Data.Number.ER.ExtendedInteger-    Description :  integer with infinities -    Copyright   :  (c) Michal Konecny-    License     :  BSD3--    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"-        

@@ -0,0 +1,42 @@+{-|+    Module      :  Data.Number.ER.Misc.STM+    Description :  some STM extras +    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mik@konecny.aow.cz+    Stability   :  experimental+    Portability :  portable+    +    Miscelaneous utilities related to concurrency.+-}+module Data.Number.ER.Misc.STM where++import Control.Concurrent as Concurrent+import Control.Concurrent.STM as STM++modifyTVar tv update =+    do+    value <- readTVar tv+    let newValue = update value+    writeTVar tv newValue+    return newValue++modifyTVarGetOldVal tv update =+    do+    value <- readTVar tv+    writeTVar tv $ update value+    return value++modifyTVarHasChanged tv update =+    do+    value <- readTVar tv+    let newValue = update value+    if value == newValue+        then return False+        else +            do+            writeTVar tv $ update value+            return True+    +    

@@ -0,0 +1,54 @@+{-|+    Module      :  Data.Number.ER.Misc.Tests+    Description :  some QuickCheck extras +    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mik@konecny.aow.cz+    Stability   :  experimental+    Portability :  portable+    +    Miscelaneous utilities related to testing.+-}+module Data.Number.ER.Misc.Tests ++where++import Data.Number.ER.Misc++import Test.QuickCheck+import Test.QuickCheck.Batch++import System.IO++erRunTests testsetName options initialise tests =+    do+    mapM (mkRunTest $ length tests) $ zip [1..] tests+    return ()+    where+    mkRunTest testCount (n, (testName, test)) =+        do+        initialise+        putStr testDescr+        result <- test options+        putStrLn $ "  result: " ++ show result+--        runTests testDescr options [test]+        hFlush stdout+        where+        testDescr = +            "(" ++ show n ++ "/" ++ show testCount ++ ") " ++ testsetName ++ ": " ++ testName ++ "\n" ++instance Show TestResult+    where+    show result =+        case result of+            TestOk msg ntest stamps ->+                msg ++ " " ++ show ntest ++ " " -- ++ show stamps+            TestExausted msg ntest stamps ->+                msg ++ " " ++ show ntest ++ " " -- ++ show stamps+            TestAborted exception ->+                "aborted: " ++ show exception+            TestFailed args ntest ->+                "failed after " ++ show ntest ++ " tests" +                ++ "\n args = " ++ show args+                    

@@ -1,31 +0,0 @@-{-|-    Utilities related to concurrency.--}-module Data.Number.ER.MiscSTM where--import Control.Concurrent as Concurrent-import Control.Concurrent.STM as STM--modifyTVar tv update =-    do-    value <- readTVar tv-    writeTVar tv $ update value--modifyTVarGetOldVal tv update =-    do-    value <- readTVar tv-    writeTVar tv $ update value-    return value--modifyTVarHasChanged tv update =-    do-    value <- readTVar tv-    let newValue = update value-    if value == newValue-        then return False-        else -            do-            writeTVar tv $ update value-            return True-    -    

@@ -1,47 +0,0 @@-{-# LANGUAGE DeriveDataTypeable   #-}-{-|-    Module      :  Data.Number.ER.PlusMinus-    Description :  mini sign datatype-    Copyright   :  (c) Michal Konecny-    License     :  BSD3--    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

@@ -16,47 +16,70 @@     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 +            Approximating a real number by a *set* of real numbers+            that includes the approximated number.            +            Precision is measured using some fixed measure on the sets.+            Operations are "safe" wrt inclusion.+            The sets can sometimes be "anti-consistent" - being smaller than+            the empty set in the inclusion order.+                    +        * 'ERInnerOuterApprox' = +            Like 'ERApprox' with the addition of operations that are "inner rounded"+            in the sense that each element of the rounded result set can +            be obtained by the same operation performed on some elements of the arument set(s).++        * 'ERIntApprox' =+            Like ERApprox but assuming that the sets are +            *intervals* of real numbers with finitely+            representable endpoints.          To be imported qualified, usually with the synonym RA. -} module Data.Number.ER.Real.Approx (     ERApprox(..),-    ERIntApprox(..),-    effIx2ra,-    splitIRA,---    checkShrinking,     eqSingletons,     leqSingletons,     ltSingletons,+    effIx2ra,+    ERInnerOuterApprox(..),+    ERIntApprox(..),+    splitIRA,     equalIntervals,     exactMiddle,     maxExtensionR2R,-    maxExtensionInnerR2R+    maxExtensionInnerR2R,+    ERApproxApprox(..) ) where  import Data.Number.ER.BasicTypes-import qualified Data.Number.ER.ExtendedInteger as EI+import qualified Data.Number.ER.BasicTypes.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.+    A type whose elements represent sets that can be used+    to approximate a single extended real number with arbitrary precision.++    Operations are "safe" with respect to inclusion, which means that+    for any numbers admitted by the operand approximations the result+    of the operation is admitted by the result approximation.+   +    The sets can sometimes be "anti-consistent" - being smaller than+    the empty set in the inclusion order.  +    This can be understood as indicating that not only there is no correct real number+    approximated here, but some numbers (ie those in interior of the set)+    are excluded more strongly than the others.+    Prime examples of such sets are directed "inverted" intervals such as [2,1].  +    Such sets arise naturally from "inner rounded" operations - see 'ERInnerOuterApprox'. -} class (Fractional ra) => ERApprox ra  	where     initialiseBaseArithmetic :: ra -> IO ()     getPrecision :: ra -> Precision      {-^ -            Precision is a measure of the set size.+            Precision is a measure of the set size.  It can be infinite.                          The default interpretation:             @@ -65,53 +88,69 @@     -}     getGranularity :: ra -> Granularity     -- ^ the lower the granularity the bigger the rounding errors-    setGranularity :: Granularity -> ra -> ra+    setGranularityOuter :: Granularity -> ra -> ra     -- ^ increase or safely decrease granularity-    setMinGranularity :: Granularity -> ra -> ra+    setMinGranularityOuter :: 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+    -- ^ true if this approximation holds no information, ie it admits any real number +    bottomApprox :: ra +    -- ^ the bottom approximation - it admits any real number     isExact :: ra -> Bool -    -- ^ true if this is a singleton+    -- ^ true if this approximation admits only one real number+    isConsistent :: ra -> Bool+    {- ^ true iff this approximation admits at least one real number -}+    isAnticonsistent :: ra -> Bool+    {- ^ true if this approximation is anti-consistent, which is a computational error +         unless we used inner rounded operations -}+    toggleConsistency :: ra -> ra+    {- ^ +        Toggle consistency - anti-consistency of the approximation. +        Top is toggled with bottom.  +        Exact approximations are the only fixed points for this operation.+    -} +    isTop :: ra -> Bool+    -- ^ true if this approximation is the most anti-consistent one+    topApprox :: ra +    -- ^ the top approximation - strongly rejects all real numbers     isDisjoint :: ra -> ra -> Bool-    isDisjoint a b = isEmpty $ a /\ b+    isDisjoint a b = not $ isConsistent $ a /\ b     isInteriorDisjoint :: ra -> ra -> Bool-    isInteriorDisjoint a b = -        isEmpty isect || isExact isect -        where-        isect = a /\ b-    isBounded :: ra -> Bool -    -- ^ true if the approximation excludes infinity-    bottomApprox :: ra -    -- ^ the bottom element - any number-    emptyApprox :: ra -    -- ^ the top element - error+    isInteriorDisjoint a b = isAnticonsistent $ a /\ b+    isBounded :: ra -> Bool+    {- ^ +        True iff the approximation excludes infinity+        and, if anti-consistent, does not strongly exclude infinity.+    -}+    plusInfinity :: ra+    -- ^ an exact approximation admitting only the positive infinity     refines :: ra -> ra -> Bool      -- ^ first arg is a subset of the second arg+    maybeRefines :: ra -> ra -> Maybe Bool +    -- ^ like 'refines' but usable for types where 'refines' is only partially decidable     (/\) :: ra -> ra -> ra -    -- ^ join; combining two approximations of the same number+    -- ^ join; combining the information in two approximations of the same number     intersectMeasureImprovement ::         EffortIndex -> ra -> ra -> (ra, ra)     {-^ -            Like intersection but the second component:+            First component of result is the intersection and the second component:             -             * measures improvement of the intersection relative to the first of the two approximations+             * measures precision improvement of the intersection relative to the first argument                            * is a positive number: 1 means no improvement, 2 means doubled precision, etc.      -}-    equalReals :: ra -> ra -> Maybe Bool -    -- ^ nothing if overlapping and not singletons+    equalReals :: ra -> ra -> Maybe Bool+    -- ^ semantic semi-decidable equality test     compareReals :: ra -> ra -> Maybe Ordering-    -- ^ nothing if overlapping and not singletons+    -- ^ semantic semi-decidable comparison     leqReals :: ra -> ra -> Maybe Bool-    -- ^ nothing if overlapping on interior or by a wrong endpoint+    -- ^ semantic semi-decidable less-than-or-equal comparison     equalApprox :: ra -> ra -> Bool-    -- ^ syntactic comparison+    -- ^ syntactic equality test     compareApprox :: ra -> ra -> Ordering     -- ^ syntactic linear ordering     double2ra :: Double -> ra+    -- ^ safe approximate conversion     showApprox ::          Int {-^ number of relevant decimals to show -} ->         Bool {-^ should show granularity -} ->@@ -120,26 +159,99 @@         String      {-|-    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)+    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 +        +{-|    +    This function converts+    an effort index to a real number approximation.+    +    Useful when an effort index is used in a formula+    mixed with real approximations.  +-}+effIx2ra :: +    (ERApprox ra) =>+    EffortIndex -> ra+effIx2ra = fromInteger . toInteger++{-|+    A type whose elements represent some kind of nominal sets of real numbers+    over which one can perform two kinds of arithmetic:+   +    * "outer rounded": arithmetic that approximates maximal extensions from outside (ie the 'ERApprox' arithmetic)+   +    * "inner rounded": arithmetic that approximates maximal extensions from inside, potentially leading to+      anti-consistent set specifications (eg intervals whose endpoints are not in the usual order)++    Another explanation of the difference:++    * `outer': the approximation contains all the number(s) of interest+    * `inner': all numbers eligible for the approximation are numbers of interest++    Ie inner rounded operations have the property that each real number admitted by the result can+    be obtained as the exact result of the same operation performed on some real numbers admitted+    by the operand approximations.+    +    While in "outer rounded" operations it is desirable to make the result set as small as+    possible in order to reduce the amount of bogus result numbers, +    in "inner rounded" operations it is desirable to make the result set as large as possible+    to lose less of the genuinely feasible result numbers.+     +    Inner rounded arithmetic is useful eg for proving/disproving inclusions "f(x) subset g(x)"+    where f and g are expressions using arithmetic extended to sets.+    For proving the inclusion, we need an inner rounded approximation of g(x)+    and for disproving the inclusion we need an inner rounded approximation of f(x).+   +    This is an abstraction of Kaucher's extended interval arithmetic    +    [Kaucher, E.: Interval Analysis in the Extended Interval Space IR, +     Computing, Suppl. 2, 1980, pp. 33-49].+-}+class (ERApprox xra) => ERInnerOuterApprox xra +    where+    (+:) :: xra -> xra -> xra+    -- ^ inner rounded addition+    (-:) :: xra -> xra -> xra+    -- ^ inner rounded subtraction+    a -: b = a +: (negate b)+    (*:) :: xra -> xra -> xra+    -- ^ inner rounded multiplication+    (/:) :: xra -> xra -> xra+    -- ^ inner rounded division+    setGranularityInner :: Granularity -> xra -> xra+    -- ^ increase or safely decrease granularity+    setMinGranularityInner :: Granularity -> xra -> xra+    -- ^ ensure granularity is not below the first arg++{-|    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+--class (ERApprox sra) => ERSetApprox sra where --    (\/) :: sra -> sra -> sra -- ^ union; either approximation could be correct  {-|@@ -174,34 +286,6 @@     (\/) :: ira -> ira -> ira      {-|-    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 -        -{-|     Return true if and only if the two intervals have equal endpoints. -} equalIntervals ::@@ -214,18 +298,6 @@     (d2L, d2U) = bounds d2  -{-|    -    This function converts-    an effort index to a real number approximation.-    -    Useful when an effort index is used in a formula-    mixed with real approximations.  --}-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.@@ -268,8 +340,8 @@     gran = max (getGranularity domL) (getGranularity domR)     domM = (domL + domR) / 2     higherGran = gran + 1-    domLhg = setMinGranularity higherGran domL-    domRhg = setMinGranularity higherGran domR+    domLhg = setMinGranularityOuter higherGran domL+    domRhg = setMinGranularityOuter higherGran domR     domMhg = (domLhg + domRhg) / 2               @@ -278,16 +350,21 @@     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 its 'EffortIndex' tending to infinity.+    For finite effor indices the function may add *outer* rounding+    but it should be reasonably small. -} maxExtensionR2R ::     (ERIntApprox ira) =>     (EffortIndex -> ira -> [ira]) -        {-^ returns a safe approximation of all extrema within the interval -} ->+        {-^ returns an *outer* approximation of all extrema within the interval -} ->     (EffortIndex -> ira -> ira) -        {-^ a function behaving well on sequences that intersect to a point -} ->+        {-^ an *outer* rounding 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 -}+        {- ^ an outer rounding function behaving well on sequences that intersect to a non-empty interval -} maxExtensionR2R getExtremes f ix x+    | not $ isConsistent x =+        toggleConsistency $+            maxExtensionInnerR2R getExtremes f ix $ toggleConsistency x      | getPrecision x < effIx2prec ix =         foldl1 (\/) $ [f ix xL, f ix xR] ++ (getExtremes ix x)     -- x is thin enough (?), don't bother evaluating by endpoints and extrema:@@ -301,47 +378,44 @@     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 its 'EffortIndex' tending to infinity.+    For finite effor indices the function may include *inner* rounding+    but it should be reasonably small. -} maxExtensionInnerR2R ::     (ERIntApprox ira) =>-    (EffortIndex -> ira -> ([ira], (Maybe Bool, Maybe Bool)))-        {-^ returns a safe approximation of all extrema within the interval-            and an indication whether the function is increasing or decreasing -            at the endpoints of the queried real approximation -} ->+    (EffortIndex -> ira -> [ira]) +        {-^ returns an *outer* 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 -}-maxExtensionInnerR2R getExtremesAndDirections f ix x =-    case (isIncreasing, isDecreasing, compareReals leftVal rightVal) of---        (True, _, Just GT) -> emptyApprox---        (True, _, Nothing) -> emptyApprox---        (_, True, Just LT) -> emptyApprox---        (_, True, Nothing) -> emptyApprox-        (True, _, _) -> fromBounds (leftVal, rightVal)-        (_, True, _) -> fromBounds (rightVal, leftVal)-        _ -> -            (/\) ((-1) \/ 1) $ -                foldl1 (\/) $ [leftVal, rightVal] ++ extremes+        {-^ an *outer* rounding function behaving well on sequences that intersect to a point -} ->+    (EffortIndex -> ira -> ira) +        {- ^ an inner rounding function behaving well on sequences that intersect to a non-empty interval -}+maxExtensionInnerR2R getExtremes f ix x+    | not $ isConsistent x =+        toggleConsistency $+            maxExtensionR2R getExtremes f ix $ toggleConsistency x+    | otherwise =+        foldl1 (\/) $ map toggleConsistency $ [f ix xL, f ix xR] ++ (getExtremes ix x)     where-    (extremes, (maybeLowIncreasing, maybeHighIncreasing)) =-        getExtremesAndDirections ix x -    (Just lowIsIncreasing) = maybeLowIncreasing -    (isIncreasing, isDecreasing)-        | null extremes = (lowIsIncreasing, not lowIsIncreasing)-        | otherwise = (False, False)-    leftVal =-        case maybeLowIncreasing of-            Just True -> snd $ bounds $ f ix xL-            Just False -> fst $ bounds $ f ix xL-            Nothing -> emptyApprox-    rightVal =-        case maybeHighIncreasing of-            Just True -> fst $ bounds $ f ix xR-            Just False -> snd $ bounds $ f ix xR-            Nothing -> emptyApprox     (xL, xR) = bounds x         -        -        -        +{-|+   A type whose elements are thought of as sets of approximations of real numbers.+   +   Eg intervals of intervals, eg [[0,3],[1,2]] containing all intervals+   whose left endpoint is between 0 and 1 and the right endpoint is between 2 and 3.+   The upper bound interval can sometimes be anti-consistent,+   eg [[0,3],[2,1]] containing all intervals (consistent as well as anti-consistent) +   with a left endpoint between [0,2] and the right endpoint between [1,3].+-}+class ERApproxApprox xra +    where+    safeIncludes :: xra -> xra -> Bool+    -- ^ safe inclusion of approximations+    safeNotIncludes :: xra -> xra -> Bool+    -- ^ safe negation of inclusion of approximations+    includes :: xra -> xra -> Maybe Bool+    -- ^ like 'safeIncludes' but usable for types where 'safeIncludes' is only partially decidable+    includes aa1 aa2 +        | safeIncludes aa1 aa2 = Just True+        | safeNotIncludes aa1 aa2 = Just False+        | otherwise = Nothing

@@ -12,13 +12,15 @@ -} module Data.Number.ER.Real.Approx.Elementary  (-    ERApproxElementary(..)+    ERApproxElementary(..),+    ERInnerOuterApproxElementary(..) ) where  import Prelude hiding (exp, log, sin, cos)  import qualified Data.Number.ER.Real.Approx as RA +import Data.Number.ER.Real.Approx ((+:),(-:),(*:),(/:))  import Data.Number.ER.BasicTypes  import Data.Number.ER.Real.Arithmetic.Elementary@@ -38,6 +40,8 @@     min ix = Prelude.min     max :: EffortIndex -> ra -> ra -> ra     max ix = Prelude.max+    sqrt :: EffortIndex -> ra -> ra+    sqrt = erSqrt_IR     exp :: EffortIndex -> ra -> ra     exp = erExp_IR     log :: EffortIndex -> ra -> ra@@ -54,6 +58,38 @@     tan ix r = (sin ix r) / (cos ix r)      atan :: EffortIndex -> ra -> ra     atan = erATan_IR+    +{-|+    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, RA.ERInnerOuterApprox ra, Ord ra) => (ERInnerOuterApproxElementary ra) +    where+    absInner :: EffortIndex -> ra -> ra+    absInner ix = Prelude.abs+    minInner :: EffortIndex -> ra -> ra -> ra+    minInner ix = Prelude.min+    maxInner :: EffortIndex -> ra -> ra -> ra+    maxInner ix = Prelude.max+    sqrtInner :: EffortIndex -> ra -> ra+    sqrtInner = erSqrt_IR_Inner+    expInner :: EffortIndex -> ra -> ra+    expInner = erExp_IR_Inner+    logInner :: EffortIndex -> ra -> ra+    logInner = erLog_IR_Inner+    (**:) :: EffortIndex -> ra -> ra -> ra+    (**:) ix b e = expInner ix $ e *: (logInner ix b)+    sinInner :: EffortIndex -> ra -> ra+    sinInner = erSine_IR_Inner+    cosInner :: EffortIndex -> ra -> ra+    cosInner = erCosine_IR_Inner+    tanInner :: EffortIndex -> ra -> ra+    tanInner ix r = (sinInner ix r) /: (cosInner ix r) +    atanInner :: EffortIndex -> ra -> ra+    atanInner = erATan_IR_Inner               

@@ -1,5 +1,6 @@ {-# LANGUAGE DeriveDataTypeable   #-} {-# LANGUAGE ScopedTypeVariables  #-}+{-# LANGUAGE FlexibleInstances  #-} {-|     Module      :  Data.Number.ER.Real.Approx.Interval     Description :  safe interval arithmetic@@ -16,17 +17,16 @@ module Data.Number.ER.Real.Approx.Interval  (     ERInterval(..),-    normaliseERInterval,-    intervalTimesInner,-    intervalPlusInner,-    intervalDivideInner+    normaliseERIntervalOuter,+    normaliseERIntervalInner ) where  import qualified Data.Number.ER.Real.Approx as RA+import Data.Number.ER.Real.Approx ((+:),(-:),(*:),(/:)) 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 qualified Data.Number.ER.BasicTypes.ExtendedInteger as EI  import Data.Number.ER.BasicTypes import Data.Number.ER.Misc@@ -44,26 +44,22 @@     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+--    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+        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+  put (ERInterval a b) = putWord8 0 >> 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)+      0 -> get >>= \a -> get >>= \b -> return (ERInterval a b)       _ -> fail "no parse" {- the above has been generated by BinaryDerive -}     @@ -73,52 +69,66 @@          * no NaNs as endpoints     -    * @l <= r@-    -    * no (-Infty, +Infty)+    Note that inverted intervals are fully supported using Warmus-Kaucher arithmetic.+    This version interprets NaN's as bottomApprox.  -}-normaliseERInterval :: +normaliseERIntervalOuter ::      (B.ERRealBase b) =>      ERInterval b -> ERInterval b-normaliseERInterval (ERInterval minusInfty plusInfty) -    | B.isPlusInfinity plusInfty && B.isPlusInfinity (- minusInfty) = -        ERIntervalAny-normaliseERInterval (ERInterval nan1 nan2) +normaliseERIntervalOuter (ERInterval nan1 nan2)      | B.isERNaN nan1 && B.isERNaN nan2 =-        ERIntervalAny-normaliseERInterval (ERInterval nan r) +        RA.bottomApprox+normaliseERIntervalOuter (ERInterval nan r)      | B.isERNaN nan =          ERInterval (- B.plusInfinity) r-normaliseERInterval (ERInterval l nan) +normaliseERIntervalOuter (ERInterval l nan)      | B.isERNaN nan =          ERInterval l (B.plusInfinity)-normaliseERInterval (ERInterval l r)-    | l > r = ERIntervalEmpty-normaliseERInterval i = i+normaliseERIntervalOuter i = i  {-|+    convert to a normal form, ie:+    +    * no NaNs as endpoints+    +    Note that inverted intervals are fully supported using Warmus-Kaucher arithmetic.+    This version interprets NaN's as topApprox. +-}+normaliseERIntervalInner :: +    (B.ERRealBase b) => +    ERInterval b -> ERInterval b+normaliseERIntervalInner (ERInterval nan1 nan2) +    | B.isERNaN nan1 && B.isERNaN nan2 =+        RA.topApprox+normaliseERIntervalInner (ERInterval nan r) +    | B.isERNaN nan = +        ERInterval (B.plusInfinity) r+normaliseERIntervalInner (ERInterval l nan) +    | B.isERNaN nan = +        ERInterval l (- B.plusInfinity)+normaliseERIntervalInner 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.+    Notice that this is +Infty for singleton and anti-consistent intervals+    and -Infty for unbounded intervals. -}     erintvPrecision ::      (B.ERRealBase b) =>      ERInterval b -> EI.ExtendedInteger-erintvPrecision (ERInterval l r) =-    -1 - (B.getApproxBinaryLog $ (r - l)) -- /(1 + abs r + abs l))-erintvPrecision ERIntervalEmpty = EI.PlusInfinity-erintvPrecision ERIntervalAny = EI.MinusInfinity+erintvPrecision i@(ERInterval l r)+    | not $ RA.isConsistent i = EI.PlusInfinity+    | not $ RA.isBounded i = EI.MinusInfinity+    | otherwise = +        -1 - (B.getApproxBinaryLog $ (r - l)) -- /(1 + abs r + abs l))  erintvGranularity ::      (B.ERRealBase b) =>      ERInterval b -> Int-erintvGranularity ERIntervalAny = 0-erintvGranularity ERIntervalEmpty = 0 erintvGranularity (ERInterval l r) =     min (B.getGranularity l) (B.getGranularity r) @@ -132,9 +142,6 @@     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@@ -142,12 +149,6 @@ 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@@ -165,10 +166,6 @@     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@@ -184,10 +181,6 @@     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@@ -204,10 +197,6 @@     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@@ -237,34 +226,35 @@     Granularity ->      (ERInterval b) ->     (ERInterval b)-erintvDefaultBisectPt gran ERIntervalAny = 0-erintvDefaultBisectPt gran ERIntervalEmpty = ERIntervalEmpty-erintvDefaultBisectPt gran (ERInterval l r) =+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+    m = +        case (B.isMinusInfinity l, B.isPlusInfinity r, B.isPlusInfinity l, B.isMinusInfinity r) of+            (True, True, _, _) -> 0 -- [-oo,+oo] +            (True, _,_,True) -> B.minusInfinity -- [-oo,-oo]+            (_, True,True,_) -> B.plusInfinity -- [+oo,+oo]+            (True, _,_,_) | r > 0 -> 0 +            (True, _,_,_) -> 2 * (B.setMinGranularity gran r) - 1+            (_,True,_,_) | l < 0 -> 0 +            (_,True,_,_) -> 2 * (B.setMinGranularity gran l) + 1  +            (_,_,True,_) | r < 0 -> 0 +            (_,_,True,_) -> 2 * (B.setMinGranularity gran r) + 1+            (_,_,_,True) | l > 0 -> 0 +            (_,_,_,True) -> 2 * (B.setMinGranularity gran l) - 1  +            _ -> ((B.setMinGranularity gran l) + r)/2 -- no infinities       erintvBisect ::-    (B.ERRealBase b, RealFrac b) => +    (B.ERRealBase b) =>      Granularity ->      (Maybe (ERInterval b)) ->     (ERInterval b) ->     (ERInterval b, ERInterval b)-erintvBisect gran maybePt i =-    (l RA.\/ m, m RA.\/ r)+erintvBisect gran maybePt i@(ERInterval l r) =+    (ERInterval l mR, ERInterval mL r)     where-    (l,r) = RA.bounds i+    ERInterval mL mR = m     m =         case maybePt of             Just m -> m@@ -294,20 +284,10 @@         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+        ERInterval (max l1 l2) (max r1 r2)     {- 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+        ERInterval (min l1 l2) (min r1 r2)          instance (B.ERRealBase b) => Show (ERInterval b)      where@@ -316,10 +296,12 @@ erintvShow numDigits showGran showComponents interval =     showERI interval     where-    showERI ERIntervalEmpty = "[NONE]"-    showERI ERIntervalAny = "[ANY]"-    showERI (ERInterval l r) +    showERI (ERInterval l r)+        | (B.isMinusInfinity r) && (B.isPlusInfinity r) =+            "[ANY]"          | l == r = "<" ++ showBase l ++ ">"+        | l > r =+            "[!" ++ showBase l ++ "," ++ showBase r ++ "!]"         | otherwise =              "[" ++ showBase l ++ "," ++ showBase r ++ "]"     showBase = B.showDiGrCmp numDigits showGran showComponents@@ -331,195 +313,201 @@             H.toHtml $ show l         | otherwise =             H.simpleTable [] [] [[H.toHtml l],[H.toHtml r]]-    toHtml i = H.toHtml $ show i   instance (B.ERRealBase b) => Num (ERInterval b) where     fromInteger n =-        normaliseERInterval $ ERInterval (fromInteger n) (fromInteger n)+        ERInterval (B.fromIntegerDown n) (B.fromIntegerUp n)     {- abs -}     abs (ERInterval l r)-        | l < 0 && r > 0 = ERInterval 0 (max (-l) r)+        | l <= 0 && r >= 0 = ERInterval 0 (max (-l) r)+        | l >= 0 && r <= 0 = ERInterval (max l (-r)) 0         | 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+    signum i@(ERInterval l r) =+        error "ER.Real.Approx.Interval: signum not implemented for ERInterval"+--        | 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     {- negate -}     negate (ERInterval l r) = (ERInterval (-r) (-l))-    negate ERIntervalEmpty = ERIntervalEmpty-    negate ERIntervalAny = ERIntervalAny     {- addition -}-    i1 + i2 = fst $ intervalPlusOuterInner i1 i2+    i1@(ERInterval l1 r1) + i2@(ERInterval l2 r2) = +        normaliseERIntervalOuter $+            ERInterval (l1 `plusDown` l2) (r1 `plusUp` r2)     {- multiplication -}-    i1 * i2 = fst $ intervalTimesOuterInner i1 i2 --{-|-    Add two real approximations, assuming the approximations are `inner'-    as opposed to `outer':-    -    * `outer': the approximation contains all the number(s) of interest-    * `inner': all numbers eligible for the approximation are numbers of interest--}-intervalPlusInner ::-    (B.ERRealBase b) =>-    (ERInterval b) -> -    (ERInterval b) -> -    (ERInterval b) -intervalPlusInner i1 i2 = snd $ intervalPlusOuterInner i1 i2--{-|-    Multiply two real approximations, assuming the approximations are `inner'-    as opposed to `outer':-    -    * `outer': the approximation contains all the number(s) of interest-    * `inner': all numbers eligible for the approximation are numbers of interest--}-intervalTimesInner ::-    (B.ERRealBase b) =>-    (ERInterval b) -> -    (ERInterval b) -> -    (ERInterval b) -intervalTimesInner i1 i2 = snd $ intervalTimesOuterInner i1 i2--intervalPlusOuterInner (ERInterval l1 r1) (ERInterval l2 r2) =-        (normaliseERInterval $-         ERInterval (l1 `plusDown` l2) (r1 `plusUp` r2),-         ERInterval (l1 `plusUp` l2) (r1 `plusDown` r2))-intervalPlusOuterInner ERIntervalAny i2 = (ERIntervalAny, ERIntervalAny)-intervalPlusOuterInner l1 ERIntervalAny = (ERIntervalAny, ERIntervalAny)-intervalPlusOuterInner ERIntervalEmpty i2 = (ERIntervalEmpty, ERIntervalEmpty)-intervalPlusOuterInner l1 ERIntervalEmpty = (ERIntervalEmpty, ERIntervalEmpty)--intervalTimesOuterInner (ERInterval l1 r1) (ERInterval l2 r2)-    | haveNan = (ERIntervalAny, ERIntervalAny)-    | otherwise =-        (normaliseERInterval $-         ERInterval minProdOuter maxProdOuter, -         ERInterval minProdInner maxProdInner)-    where-    haveNan = or $ map B.isERNaN (prodsUp ++ prodsDown)-    minProdOuter = foldl1 min prodsDown-    maxProdOuter = foldl1 max prodsUp-    minProdInner = foldl1 min prodsUp-    maxProdInner = foldl1 max prodsDown-    prodsDown = [l1 `timesDown` l2, l1 `timesDown` r2, r1 `timesDown` l2, r1 `timesDown` r2]-    prodsUp = [l1 `timesUp` l2, l1 `timesUp` r2, r1 `timesUp` l2, r1 `timesUp` r2]-intervalTimesOuterInner ERIntervalAny i2 = (ERIntervalAny, ERIntervalAny)-intervalTimesOuterInner l1 ERIntervalAny = (ERIntervalAny, ERIntervalAny)-intervalTimesOuterInner ERIntervalEmpty i2 = (ERIntervalEmpty, ERIntervalEmpty)-intervalTimesOuterInner l1 ERIntervalEmpty = (ERIntervalEmpty, ERIntervalEmpty)+    i1@(ERInterval l1 r1) * i2@(ERInterval l2 r2) = +        normaliseERIntervalOuter $+             intervalTimes timesDown timesUp i1 i2  instance (B.ERRealBase b) => Fractional (ERInterval b) where     fromRational rat =         (fromInteger $ numerator rat)         / (fromInteger $ denominator rat)     {- division -}-    i1 / i2 =-        fst $ intervalDivideOuterInner i1 i2-        -intervalDivideInner ::-    (B.ERRealBase b) =>-    (ERInterval b) -> -    (ERInterval b) -> -    (ERInterval b) -intervalDivideInner i1 i2 = snd $ intervalDivideOuterInner i1 i2+    recip i@(ERInterval l r)+        | not $ RA.isConsistent i = +            RA.toggleConsistency $ +                1 /: (RA.toggleConsistency i)+        | 0 < l || r < 0 =+            normaliseERIntervalOuter $+                ERInterval (1 `divideDown` r) (1 `divideUp` l)+        | otherwise =+            RA.bottomApprox -intervalDivideOuterInner (ERInterval l1 r1) (ERInterval l2 r2)-        | l2 < 0 && r2 > 0 = (ERIntervalAny, ERIntervalAny)-        | haveNan = ---            unsafePrint "ERInterval: /: haveNan" $ -            (ERIntervalAny, ERIntervalAny)-        | l2 == 0 && r2 > 0 && 1/l2 < 0 = -- minus 0-            intervalDivideOuterInner (ERInterval l1 r1) (ERInterval (-l2) r2) -- correct it to +0-        | r2 == 0 && l2 < 0 && 1/r2 > 0 = -- plus 0-            intervalDivideOuterInner (ERInterval l1 r1) (ERInterval l2 (-r2)) -- correct it to -0++instance (B.ERRealBase b) => RA.ERInnerOuterApprox (ERInterval b)+    where+    {- addition -}+    i1@(ERInterval l1 r1) +: i2@(ERInterval l2 r2) = +        normaliseERIntervalInner $+            ERInterval (l1 `plusUp` l2) (r1 `plusDown` r2)+    {- multiplication -}+    i1@(ERInterval l1 r1) *: i2@(ERInterval l2 r2) = +        normaliseERIntervalInner $+             intervalTimes timesUp timesDown i1 i2+    {- division -}+    i1@(ERInterval l1 r1) /: i2@(ERInterval l2 r2) +        | not $ RA.isConsistent i2 = +            (*:) i1 $+                RA.toggleConsistency $ +                    1 / (RA.toggleConsistency i2)+        | 0 < l2 || r2 < 0 = +            (*:) i1 $+                normaliseERIntervalInner $+                    ERInterval (1 `divideDown` r2) (1 `divideUp` l2)         | otherwise =---            unsafePrintReturn---            (---                "intervalDivideOuterInner: "---                ++ "\n divsUp = " ++ show divsUp---                ++ "\n divsDown = " ++ show divsDown---                ++ "\n result = "---            )-            (-             normaliseERInterval $-             ERInterval minDivOuter maxDivOuter-            ,-             ERInterval minDivInner maxDivInner-            )-        where-        haveNan = or $ map B.isERNaN (divsUp ++ divsDown)-        minDivOuter = foldl1 min divsDown-        maxDivOuter = foldl1 max divsUp-        minDivInner = foldl1 min divsUp-        maxDivInner = foldl1 max divsDown-        divsDown = [l1 `divideDown` l2, l1 `divideDown` r2, r1 `divideDown` l2, r1 `divideDown` r2]-        divsUp = [l1 `divideUp` l2, l1 `divideUp` r2, r1 `divideUp` l2, r1 `divideUp` r2]-intervalDivideOuterInner ERIntervalAny i2 = (ERIntervalAny, ERIntervalAny)-intervalDivideOuterInner i1 ERIntervalAny = (ERIntervalAny, ERIntervalAny)-intervalDivideOuterInner ERIntervalEmpty i2 = (ERIntervalEmpty, ERIntervalEmpty)-intervalDivideOuterInner i1 ERIntervalEmpty = (ERIntervalEmpty, ERIntervalEmpty)+            RA.bottomApprox+    {- setMinGranularityInner -}+    setMinGranularityInner gr (ERInterval l r) =+        normaliseERIntervalInner $+        (ERInterval (B.setMinGranularity gr l) (negate $ B.setMinGranularity gr (-r)))+    {- setGranularityInner -}+    setGranularityInner gr (ERInterval l r) =+        normaliseERIntervalInner $+        (ERInterval (B.setGranularity gr l) (negate $ B.setGranularity gr (- r))) +intervalTimes timesL timesR i1@(ERInterval l1 r1) i2@(ERInterval l2 r2) =+    ERInterval l r+    where+    (l,r) = +        case (compare l1 0, compare r1 0, l1 <= r1, compare l2 0, compare r2 0, l2 <= r2) of+            -- i1 negative, i2 positive+            (LT, LT, _, GT, GT, _) -> (l1 `timesL` r2, r1 `timesR` l2)+            -- i1 negative, i2 negative+            (LT, LT, _, LT, LT, _) -> (r1 `timesL` r2, l1 `timesR` l2)+            -- i1 negative, i2 consistent and containing zero+            (LT, LT, _, _, _, True) -> (l1 `timesL` r2, l1 `timesR` l2)+            -- i1 negative, i2 inconsistent and anti-containing zero+            (LT, LT, _, _, _, False) -> (r1 `timesL` r2, r1 `timesR` l2)+            +            -- i1 positive, i2 positive+            (GT, GT, _, GT, GT, _) -> (l1 `timesL` l2, r1 `timesR` r2)+            -- i1 positive, i2 negative+            (GT, GT, _, LT, LT, _) -> (r1 `timesL` l2, l1 `timesR` r2)+            -- i1 positive, i2 consistent and containing zero+            (GT, GT, _, _, _, True) -> (r1 `timesL` l2, r1 `timesR` r2)+            -- i1 positive, i2 inconsistent and anti-containing zero+            (GT, GT, _, _, _, False) -> (l1 `timesL` l2, l1 `timesR` r2)++            -- i1 consistent and containing zero, i2 positive+            (_, _, True, GT, GT, _) -> (l1 `timesL` r2, r1 `timesR` r2)+            -- i1 consistent and containing zero, i2 negative+            (_, _, True, LT, LT, _) -> (r1 `timesL` l2, l1 `timesR` l2)+            -- i1 consistent and containing zero, i2 consistent and containing zero+            (_, _, True, _, _, True) -> +                (l,r)+                where+                l | B.isERNaN l1r2 || B.isERNaN r1l2 = B.minusInfinity+                  | otherwise = min l1r2 r1l2+                  where+                  l1r2 = l1 `timesL` r2+                  r1l2 = r1 `timesL` l2+                r | B.isERNaN l1l2 || B.isERNaN r1r2 = B.plusInfinity+                  | otherwise = max l1l2 r1r2+                  where+                  l1l2 = l1 `timesR` l2+                  r1r2 = r1 `timesR` r2+            -- i1 consistent and containing zero, i2 inconsistent and anti-containing zero+            (_, _, True, _, _, False) -> (0, 0)++            -- i1 inconsistent and anti-containing zero, i2 positive +            (_, _, False, GT, GT, _) -> (l1 `timesL` l2, r1 `timesR` l2)+            -- i1 inconsistent and anti-containing zero, i2 negative +            (_, _, False, LT, LT, _) -> (r1 `timesL` r2, l1 `timesR` r2)+            -- i1 inconsistent and anti-containing zero, i2 consistent and containing zero +            (_, _, False, _, _, True) -> (0, 0) +            -- i1 inconsistent and anti-containing zero, i2 the same +            (_, _, False, _, _, False) ->+                (l,r)+                where+                l | B.isERNaN l1l2 || B.isERNaN r1r2 = B.plusInfinity+                  | otherwise = max l1l2 r1r2+                  where+                  l1l2 = l1 `timesL` l2+                  r1r2 = r1 `timesL` r2+                r | B.isERNaN l1r2 || B.isERNaN r1l2 = B.minusInfinity+                  | otherwise = min l1r2 r1l2+                  where+                  l1r2 = l1 `timesR` r2+                  r1l2 = r1 `timesR` l2++                   -instance (B.ERRealBase b, RealFrac b) => RA.ERApprox (ERInterval b) where+instance (B.ERRealBase b) => RA.ERApprox (ERInterval b) where     initialiseBaseArithmetic _ =         B.initialiseBaseArithmetic (0 :: b)     getPrecision i = erintvPrecision i     getGranularity i = erintvGranularity i     {- setMinGranularity -}-    setMinGranularity gr (ERInterval l r) =-        normaliseERInterval $+    setMinGranularityOuter gr (ERInterval l r) =+        normaliseERIntervalOuter $         (ERInterval (- (B.setMinGranularity gr (-l))) (B.setMinGranularity gr r))-    setMinGranularity _ i = i     {- setGranularity -}-    setGranularity gr (ERInterval l r) =-        normaliseERInterval $+    setGranularityOuter gr (ERInterval l r) =+        normaliseERIntervalOuter $         (ERInterval (- (B.setGranularity gr (-l))) (B.setGranularity gr r))-    setGranularity _ i = i-    {- 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+        B.isMinusInfinity l && B.isPlusInfinity r+    {- bottomApprox -}+    bottomApprox = +        ERInterval B.minusInfinity B.plusInfinity     {- isExact -}-    isExact ERIntervalEmpty = False-    isExact ERIntervalAny = False     isExact (ERInterval l r) = l == r+    {- isConsistent -}+    isConsistent (ERInterval l r) = l <= r+    {- isAnticonsistent -}+    isAnticonsistent (ERInterval l r) = l >= r+    {- toggleConsistency -}+    toggleConsistency (ERInterval l r) = (ERInterval r l)+    {- isTop -}+    isTop (ERInterval l r) =+        B.isPlusInfinity l && B.isMinusInfinity r+    {- topApprox -}+    topApprox =+        ERInterval B.plusInfinity B.minusInfinity     {- isBounded -}-    isBounded ERIntervalEmpty = True-    isBounded ERIntervalAny = False     isBounded (ERInterval l r) = -        (- B.plusInfinity) < l && r < B.plusInfinity+        (- B.plusInfinity) < l && l < B.plusInfinity+        &&+        (- B.plusInfinity) < r && r < B.plusInfinity+    {- plusInfinity -}+    plusInfinity = ERInterval B.plusInfinity B.plusInfinity  +    {- refines -}+    refines (ERInterval l1 r1) (ERInterval l2 r2) =+        l2 <= l1 && r1 <= r2+    {- maybeRefines -}+    maybeRefines i1 i2 = Just $ RA.refines i1 i2+                  {- 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@@ -529,18 +517,10 @@             | 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 (ERInterval l r) -        | B.isPlusInfinity r && B.isPlusInfinity (-l) = True-    refines ERIntervalAny _ = False-    refines _ ERIntervalEmpty = False-    refines (ERInterval l1 r1) (ERInterval l2 r2) =-        l2 <= l1 && r1 <= r2+        (isecL, isecH) = RA.bounds $ RA.setMinGranularityOuter gran isec  +        (i1L, i1H) = RA.bounds $ RA.setMinGranularityOuter gran i1+        gran = effIx2gran ix+               {- semantic comparisons -}     equalReals = erintvEqualReals     compareReals = erintvCompareReals@@ -556,50 +536,27 @@     {- formatting -}     showApprox = erintvShow -instance (B.ERRealBase b, RealFrac b) => RA.ERIntApprox (ERInterval b)+instance (B.ERRealBase b) => RA.ERIntApprox (ERInterval b)     where-    doubleBounds ERIntervalAny = (- infinity, infinity)-        where-        infinity = 1/0-    doubleBounds ERIntervalEmpty = -        error "ERInterval: doubleBounds: empty interval"     doubleBounds (ERInterval l r) =-        (B.toDouble l, B.toDouble r) -    floatBounds ERIntervalAny = (- infinity, infinity)-        where-        infinity = 1/0-    floatBounds ERIntervalEmpty = -        error "ERInterval: floatBounds: empty interval"+        (negate $ B.toDouble (-l), B.toDouble r)      floatBounds (ERInterval l r) =-        (B.toFloat l, B.toFloat r) -    integerBounds ERIntervalAny = (- infinity, infinity)-        where-        infinity = EI.PlusInfinity-    integerBounds ERIntervalEmpty = -        error "ERInterval: integerBounds: empty interval"+        (negate $ B.toFloat (-l), B.toFloat r)      integerBounds (ERInterval l r) = -        (- (mkEI (- l)), mkEI r)+        (negate $ mkEI (- l), mkEI r)         where         mkEI f              | B.isPlusInfinity f = EI.PlusInfinity-            | B.isPlusInfinity (-f) = EI.MinusInfinity+            | B.isMinusInfinity f = EI.MinusInfinity             | otherwise = ceiling f-    defaultBisectPt dom = erintvDefaultBisectPt  (RA.getGranularity dom + 1) dom+    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)     {- RA.fromBounds -}@@ -607,9 +564,11 @@         | l1 == r1 && l2 == r2 = ERInterval l1 l2     fromBounds i1i2 =         error $-            "ER.Real.Approx.Interval: fromBounds: bounds not thin: "+            "ER.Real.Approx.Interval: fromBounds: bounds not exact: "             ++ show i1i2 -instance (B.ERRealBase b, RealFrac b) => RAEL.ERApproxElementary (ERInterval b)+instance (B.ERRealBase b) => RAEL.ERApproxElementary (ERInterval b)+instance (B.ERRealBase b) => RAEL.ERInnerOuterApproxElementary (ERInterval b) -- all operations here have appropriate default implementations-+    +    

@@ -0,0 +1,56 @@+{-# OPTIONS_GHC -fno-warn-missing-methods #-}+{-|+    Module      :  Data.Number.ER.Real.Approx.OI+    Description :  outer and inner approximations of approximations  +    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mik@konecny.aow.cz+    Stability   :  experimental+    Portability :  portable++    This module offers a transformation of a safely rounded real approximation type into+    a type that approximates these approximations from outside as well as *inside*. +-}+module Data.Number.ER.Real.Approx.OI where++import qualified Data.Number.ER.Real.Approx as RA ++{-|+    A pair of approximations that form an "interval" in the lattice of+    approximations. +    +    Eg outer = [1,4] inner = [3,2] can be thought of as the set of all+    generalised intervals where the left endpoint is between 1 and 3+    and the right endpoint is between 2 and 4 (eg [1,4], [3,4],+    [3,2], [3,3]).+-}+data ERApproxOI ra = +    ERApproxOI+    {+        eroiOuter :: ra,+        eroiInner :: ra+    }+    deriving (Eq, Ord)++instance (RA.ERApprox ra) => (Show (ERApproxOI ra))+    where+    show (ERApproxOI oi ii) =+        "{ outer = " ++ show oi ++ "; inner = " ++ show ii ++ "}" ++instance (RA.ERApprox ra) => RA.ERApproxApprox (ERApproxOI ra)+    where+    safeIncludes (ERApproxOI oi1 ii1) (ERApproxOI oi2 ii2) =+        oi2 `RA.refines` ii1+    safeNotIncludes (ERApproxOI oi1 ii1) (ERApproxOI oi2 ii2) =+        not $ ii2 `RA.refines` oi1++{- TODO when required: -}    +instance (RA.ERApprox ra) => (Num (ERApproxOI ra))+instance (RA.ERApprox ra) => (Fractional (ERApproxOI ra))+instance (RA.ERApprox ra) => RA.ERApprox (ERApproxOI ra)+    where+    (ERApproxOI oi1 ii1) `leqReals` (ERApproxOI oi2 ii2) =+        oi1 `RA.leqReals` oi2++    

@@ -174,7 +174,7 @@     fromInteger n = ConvergRealSeq sq         where         sq ix =-            RA.setMinGranularity (effIx2gran ix) $ fromInteger n+            RA.setMinGranularityOuter (effIx2gran ix) $ fromInteger n     abs = pointwiseConvergRealSeq1 $ abs     signum = pointwiseConvergRealSeq1 $ signum     negate = pointwiseConvergRealSeq1 $ negate@@ -188,7 +188,7 @@     fromRational q = ConvergRealSeq sq         where         sq ix =-            (RA.setMinGranularity (effIx2gran ix) num) / denom+            (RA.setMinGranularityOuter (effIx2gran ix) num) / denom         num = fromInteger $ numerator q         denom = fromInteger $ denominator q     recip = pointwiseConvergRealSeq1 $ recip

@@ -0,0 +1,177 @@+{-|+    Module      :  Data.Number.ER.Real.Approx.Tests.Generate+    Description :  (testing) generating real approximations+    Copyright   :  (c) 2009 Michal Konecny+    License     :  BSD3++    Maintainer  :  mik@konecny.aow.cz+    Stability   :  experimental+    Portability :  portable+    +    Generic instances of 'Arbitrary' class for generating (almost) random instances. +-}++module Data.Number.ER.Real.Approx.Tests.Generate where++import qualified Data.Number.ER.Real.Approx as RA+import Data.Number.ER.BasicTypes++import Test.QuickCheck++import qualified Data.List as List++newtype RAThin ira = RAThin ira deriving (Show)+newtype RAConsistent ira = RAConsistent ira deriving (Show)+newtype RADirected ira = RADirected ira deriving (Show)+++instance (RA.ERIntApprox ira) => Arbitrary (RAThin ira)+    where+    arbitrary = +        sized arbitrarySized+        where+        arbitrarySized n +            | n < 28 =+                do+                gran <- choose (8,20)+                (f1,f2,f3) <- arbitrary+                isInfty <- choose (-inftyChance,inftyChance)+                pow <- choose (-10,10)+                return $ RAThin $ constructThinRA isInfty gran (f1,f2,f3) pow+            | n <= 68 =+                do+                gran <- choose (30,100)+                (f1,f2,f3) <- arbitrary+                isInfty <- choose (-inftyChance,inftyChance)+                pow <- choose (-100,100)+                return $ RAThin $ constructThinRA isInfty gran (f1,f2,f3) pow+            | otherwise =+                do+                gran <- choose (400,1000)+                (f1,f2,f3) <- arbitrary+                isInfty <- choose (-inftyChance,inftyChance)+                pow <- choose (-10000,10000)+                return $ RAThin $ constructThinRA isInfty gran (f1,f2,f3) pow+    coarbitrary _ =+        error "ER.Real.Approx: Tests: coarbitrary not implemented"++inftyChance = 15+                +constructThinRA ::+    (RA.ERIntApprox ra) =>+    Granularity ->+    Int ->+    (Double, Double, Double) ->+    Int ->+    ra+constructThinRA gran isInfty (f1,f2,f3) pow +    | isInfty == inftyChance =+        RA.setGranularityOuter gran $ RA.plusInfinity+    | isInfty == - inftyChance =+        RA.setGranularityOuter gran $ negate $ RA.plusInfinity+    | isInfty == 0 =+        RA.setGranularityOuter gran 0+    | otherwise =+        fst $ RA.bounds $ -- ensure thinness+            (\ (Just a) -> a) $ List.find RA.isBounded results+    where+    results = [result1, result2, result3, result4, result5, 0]+    result1 = (b1/b2) ^^ pow + b3+    result2 = b1 * b2 + b3+    result3 = b1 ^^ pow - b2+    result4 = b1 - b2+    result5 = b1+    [b1,b2,b3] = map cvt [f1,f2,f3]+    cvt f = RA.setGranularityOuter gran $ RA.double2ra f++instance (RA.ERIntApprox ira) => Arbitrary (RAConsistent ira)+    where+    arbitrary = +        sized arbitrarySized+        where+        arbitrarySized n +            | n < 28 =+                do+                gran <- choose (8,20)+                (f11,f12,f13) <- arbitrary+                isInfty1 <- choose (-inftyChance,inftyChance)+                pow1 <- choose (-10,10)+                (f21,f22,f23) <- arbitrary+                isInfty2 <- choose (-inftyChance,inftyChance)+                pow2 <- choose (-10,10)+                let t1 = constructThinRA isInfty1 gran (f11,f12,f13) pow1+                let t2 = constructThinRA isInfty2 gran (f21,f22,f23) pow2+                return $ RAConsistent $ t1 RA.\/ t2+            | n <= 68 =+                do+                gran <- choose (30,100)+                (f11,f12,f13) <- arbitrary+                isInfty1 <- choose (-inftyChance,inftyChance)+                pow1 <- choose (-100,100)+                (f21,f22,f23) <- arbitrary+                isInfty2 <- choose (-inftyChance,inftyChance)+                pow2 <- choose (-100,100)+                let t1 = constructThinRA isInfty1 gran (f11,f12,f13) pow1+                let t2 = constructThinRA isInfty2 gran (f21,f22,f23) pow2+                return $ RAConsistent $ t1 RA.\/ t2+            | otherwise =+                do+                gran <- choose (400,1000)+                (f11,f12,f13) <- arbitrary+                isInfty1 <- choose (-inftyChance,inftyChance)+                pow1 <- choose (-10000,10000)+                (f21,f22,f23) <- arbitrary+                isInfty2 <- choose (-inftyChance,inftyChance)+                pow2 <- choose (-10000,10000)+                let t1 = constructThinRA isInfty1 gran (f11,f12,f13) pow1+                let t2 = constructThinRA isInfty2 gran (f21,f22,f23) pow2+                return $ RAConsistent $ t1 RA.\/ t2+    coarbitrary _ =+        error "ER.Real.Approx: Tests: coarbitrary not implemented"++instance (RA.ERIntApprox ira) => Arbitrary (RADirected ira)+    where+    arbitrary = +        sized arbitrarySized+        where+        arbitrarySized n +            | n < 28 =+                do+                gran <- choose (8,20)+                (f11,f12,f13) <- arbitrary+                isInfty1 <- choose (-inftyChance,inftyChance)+                pow1 <- choose (-10,10)+                (f21,f22,f23) <- arbitrary+                isInfty2 <- choose (-inftyChance,inftyChance)+                pow2 <- choose (-10,10)+                let t1 = constructThinRA isInfty1 gran (f11,f12,f13) pow1+                let t2 = constructThinRA isInfty2 gran (f21,f22,f23) pow2+                return $ RADirected $ RA.fromBounds (t1, t2)+            | n <= 68 =+                do+                gran <- choose (30,100)+                (f11,f12,f13) <- arbitrary+                isInfty1 <- choose (-inftyChance,inftyChance)+                pow1 <- choose (-100,100)+                (f21,f22,f23) <- arbitrary+                isInfty2 <- choose (-inftyChance,inftyChance)+                pow2 <- choose (-100,100)+                let t1 = constructThinRA isInfty1 gran (f11,f12,f13) pow1+                let t2 = constructThinRA isInfty2 gran (f21,f22,f23) pow2+                return $ RADirected $ RA.fromBounds (t1, t2)+            | otherwise =+                do+                gran <- choose (400,1000)+                (f11,f12,f13) <- arbitrary+                isInfty1 <- choose (-inftyChance,inftyChance)+                pow1 <- choose (-10000,10000)+                (f21,f22,f23) <- arbitrary+                isInfty2 <- choose (-inftyChance,inftyChance)+                pow2 <- choose (-10000,10000)+                let t1 = constructThinRA isInfty1 gran (f11,f12,f13) pow1+                let t2 = constructThinRA isInfty2 gran (f21,f22,f23) pow2+                return $ RADirected $ RA.fromBounds (t1, t2)+    coarbitrary _ =+        error "ER.Real.Approx: Tests: coarbitrary not implemented"++        

@@ -0,0 +1,266 @@+{-|+    Module      :  Data.Number.ER.Real.Base.Tests.Properties+    Description :  (testing) properties to check for real approximations+    Copyright   :  (c) 2009 Michal Konecny+    License     :  BSD3++    Maintainer  :  mik@konecny.aow.cz+    Stability   :  experimental+    Portability :  portable+    +    Properties of real approximations we want to check in tests. +-}++module Data.Number.ER.Real.Approx.Tests.Properties +where++import Data.Number.ER.Real.Approx.Tests.Generate+import Data.Number.ER.Real.Approx.Tests.Reporting+import Data.Number.ER.BasicTypes.Tests.Generate++import qualified Data.Number.ER.Real.Approx as RA+import Data.Number.ER.Real.Approx ((+:),(-:),(*:),(/:))+import qualified Data.Number.ER.Real.Approx.Elementary as RAEL++import Data.Number.ER.BasicTypes++import Data.Number.ER.Misc++import Test.QuickCheck++type RAPropTupleUnary ira =+    ira ->+    String ->+    ((Ix20, RAThin ira) -> Bool, +     (Ix20, RAConsistent ira) -> Bool,+     (Ix20, RAThin ira) -> Bool, +     (Ix20, RAConsistent ira) -> Bool, +     (Ix20, RADirected ira) -> Bool)++props_ra_AMinusA_eq_oi ::+    (RA.ERIntApprox ira, RA.ERInnerOuterApprox ira) => RAPropTupleUnary ira+props_ra_AMinusA_eq_oi =+    props_ra_eq_oi_unary 0 (\ix a -> a - a) (\ix a -> a -: a)++props_ra_ADivA_eq_oi ::+    (RA.ERIntApprox ira, RA.ERInnerOuterApprox ira) => RAPropTupleUnary ira+props_ra_ADivA_eq_oi =+    props_ra_eq_oi_unary 1 (\ix a -> a / a) (\ix a -> a /: a)+    +props_ra_AddCommut_eq_oi ::+    (RA.ERIntApprox ira, RA.ERInnerOuterApprox ira) => RAPropTupleUnary ira+props_ra_AddCommut_eq_oi =+    props_ra_eq_oi_unary 0 commutDiff commutDiffInner+    where+    commutDiff ix a =+        (a + b) - (b + a)+        where+        b = 1 / (a + 1)+    commutDiffInner ix a =+        (a +: b) -: (b +: a)+        where+        b = 1 / (a + 1)++props_ra_MultCommut_eq_oi ::+    (RA.ERIntApprox ira, RA.ERInnerOuterApprox ira) => RAPropTupleUnary ira+props_ra_MultCommut_eq_oi =+    props_ra_eq_oi_unary 0 commutDiff commutDiffInner+    where+    commutDiff ix a =+        (a * b) - (b * a)+        where+        b = 1 / (a + 1)+    commutDiffInner ix a =+        (a *: b) -: (b *: a)+        where+        b = 1 / (a + 1)++props_ra_AddAssoc_eq_oi ::+    (RA.ERIntApprox ira, RA.ERInnerOuterApprox ira) => RAPropTupleUnary ira+props_ra_AddAssoc_eq_oi =+    props_ra_eq_oi_unary 0 assocDiff assocDiffInner+    where+    assocDiff ix a =+        ((a + b) + c) - (a + (b + c))+        where+        b = 1 / (a + 1)+        c = (a - 1)+    assocDiffInner ix a =+        ((a +: b) +: c) -: (a +: (b +: c))+        where+        b = 1 / (a + 1)+        c = (a - 1)++props_ra_MultAssoc_eq_oi ::+    (RA.ERIntApprox ira, RA.ERInnerOuterApprox ira) => RAPropTupleUnary ira+props_ra_MultAssoc_eq_oi =+    props_ra_eq_oi_unary 0 assocDiff assocDiffInner+    where+    assocDiff ix a =+        ((a * b) * c) - (a * (b * c))+        where+        b = 1 / (a + 1)+        c = (a - 1)+    assocDiffInner ix a =+        ((a *: b) *: c) -: (a *: (b *: c))+        where+        b = 1 / (a + 1)+        c = (a - 1)++props_ra_Distr_eq_oi ::+    (RA.ERIntApprox ira, RA.ERInnerOuterApprox ira) => RAPropTupleUnary ira+props_ra_Distr_eq_oi =+    props_ra_eq_oi_unary 0 distrDiff distrDiffInner+    where+    distrDiff ix a =+        (a * (b + c)) - (a * b + a * c)+        where+        b = 1 / (a + 1)+        c = (a - 1)+    distrDiffInner ix a =+        (a *: (b +: c)) -: ((a *: b) +: (a *: c))+        where+        b = 1 / (a + 1)+        c = (a - 1)++props_ra_SinCos_eq_oi ::+    (RAEL.ERInnerOuterApproxElementary ira, RAEL.ERApproxElementary ira) => RAPropTupleUnary ira+props_ra_SinCos_eq_oi =+    props_ra_eq_oi_unary 1 sincos sincosInner+    where+    sincos ix a =+        (RAEL.sin ix a)^2 + (RAEL.cos ix a)^2+    sincosInner ix a =+        (sina *: sina) +: (cosa *: cosa)+        where+        sina = RAEL.sinInner ix a +        cosa = RAEL.cosInner ix a ++props_ra_TanATan_eq_oi ::+    (RAEL.ERInnerOuterApproxElementary ira, RAEL.ERApproxElementary ira) => RAPropTupleUnary ira+props_ra_TanATan_eq_oi =+    props_ra_eq_oi_unary 0 tanAtan tanAtanInner+    where+    tanAtan ixP a =+--        unsafePrint +--        (+--            "tanAtan: "+--            ++ "\n ix = " ++ show ix +--            ++ "\n a = " ++ show a+--            ++ "\n atan ix a = " ++ show tana+--            ++ "\n tan ix (atan ix a) = " ++ show tanatana+--        ) $+        tanatana - a+        where+        tanatana = RAEL.tan ix tana+        tana = RAEL.atan ix a+        ix = min 10 ixP+    tanAtanInner ixP a =+        (RAEL.tanInner ix $ RAEL.atanInner ix a) -: a+        where+        ix = min 10 ixP++props_ra_LogExp_eq_oi ::+    (RAEL.ERInnerOuterApproxElementary ira, RAEL.ERApproxElementary ira) => RAPropTupleUnary ira+props_ra_LogExp_eq_oi =+    props_ra_eq_oi_unary 0 logExp logExpInner+    where+    logExp ixP a =+--        unsafePrint +--        (+--            "logExp: "+--            ++ "\n ix = " ++ show ix +--            ++ "\n a = " ++ show a+--            ++ "\n exp ix a = " ++ show expa+--            ++ "\n log ix (exp ix a) = " ++ show logexpa+--        ) $+        logexpa - a+        where+        logexpa = RAEL.log ix expa +        expa = RAEL.exp ix a+        ix = min 10 ixP+    logExpInner ixP a =+        logexpa -: a+        where+        logexpa = RAEL.logInner ix expa+        expa = RAEL.expInner ix a+        ix = min 10 ixP+++{------------------  auxiliary functions ------------------------}++props_ra_eq_oi_unary constRes opOuter opInner sampleRA reportFileName =+    (prop_Eq_Thin, prop_Eq_Consistent, +     prop_OI_Thin, prop_OI_Consistent, prop_OI_Directed)+    where+    prop_Eq_Thin (Ix20 ix, RAThin a) =+        raConsistentWithPrecise sampleRA (reportFileName  ++ "_Eq_Thin") (ix,aId) 0 constRes resOuter+        where+        resOuter = opOuter ix a+        aId = RA.showApprox 10 True True a+    prop_Eq_Consistent (Ix20 ix, RAConsistent a) =+        raConsistentWithPrecise sampleRA (reportFileName ++ "_Eq_Consistent") (ix,aId) 0 constRes resOuter+        where+        resOuter = opOuter ix a+        aId = RA.showApprox 10 True True a+    prop_OI_Thin (Ix20 ix, RAThin a) =+        raIncludedIn sampleRA (reportFileName ++ "_OI_Thin") (ix, aId) 0 resInner resOuter+        where+        resOuter = opOuter ix a+        resInner = opInner ix a+        aId = RA.showApprox 10 True True a+    prop_OI_Consistent (Ix20 ix, RAConsistent a) =+        raIncludedIn sampleRA (reportFileName ++ "_OI_Consistent") (ix,aId) 0 resInner resOuter+        where+        resOuter = opOuter ix a+        resInner = opInner ix a+        aId = RA.showApprox 10 True True a+    prop_OI_Directed (Ix20 ix, RADirected a) =+        raIncludedIn sampleRA (reportFileName ++ "_OI_Directed") (ix, aId) 0 resInner resOuter+        where+        resOuter = opOuter ix a+        resInner = opInner ix a+        aId = RA.showApprox 10 True True a+    +raConsistentWithPrecise sampleRA reportFileName caseId subId preciseVal approxVal +    | result =+        unsafeERTestReport reportFileName+            (caseId, subId, preciseVal, approxVal) $+        result+    | otherwise = +        unsafePrint+        (+            "raAntiIncludes failed"+            ++ "\n caseId = " ++ show caseId+            ++ "\n subId = " ++ show subId+            ++ "\n preciseVal = " ++ show preciseVal+            ++ "\n approxVal = " ++ show approxVal+        ) $+        result+    where+    result = +        (approxVal `RA.refines` preciseVal)+        || +        (preciseVal `RA.refines` approxVal)+    _ = [sampleRA, approxVal]++raIncludedIn sampleRA reportFileName caseId subId innerVal outerVal +    | result =+        unsafeERTestReport reportFileName+            (caseId, subId, innerVal, outerVal) $+        result+    | otherwise = +        unsafePrint+        (+            "raIncludes failed"+            ++ "\n caseId = " ++ show caseId+            ++ "\n subId = " ++ show subId+            ++ "\n innerVal = " ++ show innerVal+            ++ "\n outerVal = " ++ show outerVal+        ) $+        result+    where+    result = innerVal `RA.refines` outerVal +    _ = [sampleRA, innerVal]++    

@@ -0,0 +1,167 @@++module Data.Number.ER.Real.Approx.Tests.Reporting ++where++import qualified Data.Number.ER.Real.Approx as RA++import Data.Number.ER.Misc++import qualified Data.List as List+import Text.Regex+import System.IO+++unsafeERTestReport ::+    (Show tId, Show sId, RA.ERIntApprox ira) =>+    String ->+    (tId, sId, ira, ira) ->+    a -> a+unsafeERTestReport reportFileName (testId, subId, almostPreciseVal, approxVal) =+    unsafeReport reportFileName $ +        stdRepLine (testId, subId) (overestimation, detail)+    where+    overestimation = fst $ getOverestimation almostPreciseVal approxVal+    detail = (almostPreciseVal, approxVal)++stdRepLine (testId, subId) (overestimation, detail) =+    "case=" ++ show testId+    ++ ";pt=" ++ show subId+    ++ ";ovest=" ++ show overestimation+    ++ ";detail=" ++ show detail++getOverestimation ::+    (RA.ERIntApprox ira) =>+    ira -> ira -> (Double, (ira, ira))+getOverestimation model res =+    ((abs $ wMod - wRes) / (1 + (max 0 (wMod))), (model, res))+    where+    wMod = hMod - lMod+    wRes = hRes - lRes+    (lMod, hMod) = RA.doubleBounds model+    (lRes, hRes) = RA.doubleBounds res+    +produceSummary :: String -> IO ()+produceSummary filepath =+    do+    casesInfo <- parseReport filepath+    writeFile summaryFilepath $ formatSummary casesInfo+    return ()+    where+    summaryFilepath = filepath ++ "-summary"+    formatSummary casesInfo =+        "all " ++ show casesCount ++ " cases:"+        ++ "\n  approx. average time per case: " ++ show timeInSeconds ++ " seconds"+        ++ "\n  approx. average per case average overestimation: " ++ show avgOverestimation+        ++ "\n  approx. average per case maximum overestimation: " ++ show maxOverestimation+        ++ "\n\n removing the worst 5% of the cases (for each measure separately):"+        ++ "\n  approx. average time per case: " ++ show timeInSeconds95 ++ " seconds"+        ++ "\n  approx. average per case average overestimation: " ++ show avgOverestimation95+        ++ "\n  approx. average per case maximum overestimation: " ++ show maxOverestimation95+        ++ "\n\n considering only the worst 50% but not the worst 5% of the cases (for each measure separately):"+        ++ "\n  approx. average time per case: " ++ show timeInSeconds45 ++ " seconds"+        ++ "\n  approx. average per case average overestimation: " ++ show avgOverestimation45+        ++ "\n  approx. average per case maximum overestimation: " ++ show maxOverestimation45+        ++ "\n\n considering only the best 50% of the cases (for each measure separately):"+        ++ "\n  approx. average time per case: " ++ show timeInSeconds50 ++ " seconds"+        ++ "\n  approx. average per case average overestimation: " ++ show avgOverestimation50+        ++ "\n  approx. average per case maximum overestimation: " ++ show maxOverestimation50+        ++ "\n\n" ++ (unlines $ map formatSummaryCase casesInfo)+        where+        (allTimes, (allAvgOvers, allMaxOvers)) =+            mapSnd unzip $ unzip $ snd $ unzip casesInfo+        timeInSeconds = (sum allTimes) / casesCountF+        avgOverestimation = (sum allAvgOvers) / casesCountF+        maxOverestimation = (sum allMaxOvers) / casesCountF+        casesCount = length casesInfo+        casesCountF :: Double+        casesCountF = fromInteger $ toInteger casesCount+        +        timeInSeconds95 = (sum allTimes95) / casesCount95F+        avgOverestimation95 = (sum allAvgOvers95) / casesCount95F+        maxOverestimation95 = (sum allMaxOvers95) / casesCount95F+        allTimes95 = drop fivePerCent $ reverse $ List.sort allTimes +        allAvgOvers95 = drop fivePerCent $ reverse $ List.sort allAvgOvers +        allMaxOvers95 = drop fivePerCent $ reverse $ List.sort allMaxOvers+        casesCount95F = fromInteger $ toInteger $ casesCount - fivePerCent+        fivePerCent = max 1 $ (5 * casesCount) `div` 100+        +        timeInSeconds50 = (sum allTimes50) / casesCount50F+        avgOverestimation50 = (sum allAvgOvers50) / casesCount50F+        maxOverestimation50 = (sum allMaxOvers50) / casesCount50F+        allTimes50 = drop fiftyPerCent $ reverse $ List.sort allTimes +        allAvgOvers50 = drop fiftyPerCent $ reverse $ List.sort allAvgOvers +        allMaxOvers50 = drop fiftyPerCent $ reverse $ List.sort allMaxOvers+        casesCount50F = fromInteger $ toInteger $ casesCount - fiftyPerCent+        fiftyPerCent = casesCount `div` 2+        +        timeInSeconds45 = (sum allTimes45) / casesCount45F+        avgOverestimation45 = (sum allAvgOvers45) / casesCount45F+        maxOverestimation45 = (sum allMaxOvers45) / casesCount45F+        allTimes45 =  drop fivePerCent $ reverse $ drop fiftyPerCent $ List.sort allTimes +        allAvgOvers45 = drop fivePerCent $ reverse $ drop fiftyPerCent $ List.sort allAvgOvers +        allMaxOvers45 = drop fivePerCent $ reverse $ drop fiftyPerCent $ List.sort allMaxOvers+        casesCount45F = fromInteger $ toInteger $ casesCount - fiftyPerCent - fivePerCent+    formatSummaryCase (caseId, (timeInSeconds, (avgOverestimation, maxOverestimation))) =+        "case " ++ caseId ++ ":"+        ++ "\n  approximate time = " ++ show timeInSeconds ++ " seconds"+        ++ "\n  average sampled overestimation = " ++ show avgOverestimation +        ++ "\n  maximal sampled overestimation = " ++ show maxOverestimation+    parseReport :: String -> IO [(String, (Double, (Double, Double)))]+    parseReport filepath =+        withFile filepath ReadMode readFirstAndOtherLines+        where+        readFirstAndOtherLines h =+            do+            startLine <- hGetLine h+            firstLine <- hGetLine h+            readCases (firstLine, (getTime firstLine) - (getTime startLine)) h+        readCases (currentLine, caseCompTime) h =+            do+            (caseOverestimations, maybeNextLineAndTime) <- readCase [] 0 currentLine+            let caseInfo = (caseId, (caseCompTime, avgAndMax caseOverestimations))+            case maybeNextLineAndTime of+                Nothing -> return [caseInfo]+                Just (nextLine, nextCaseTime) ->+                    do+                    otherCases <- readCases (nextLine, nextCaseTime) h+                    return $ caseInfo : otherCases+            where     +            avgAndMax ns =+                (sum ns / (fromInteger $ toInteger $ length ns), foldl1 max ns)+            caseId = getCaseId currentLine+            readCase overestimationsSoFar currentTimeStep currentLine+                | currentCaseId /= caseId =+                    return (overestimationsSoFar, Just (currentLine, currentTimeStep))+                | otherwise =+                    do+                    finished <- hIsEOF h+                    case finished of+                        True -> return (currentOverestimations, Nothing)+                        False ->+                            do+                            nextLine <- hGetLine h+                            let nextTimeStep = (getTime nextLine) - (getTime currentLine)+                            readCase currentOverestimations nextTimeStep nextLine+                where+                currentCaseId = getCaseId currentLine+                currentOverestimations = +                    currentOverestimation : overestimationsSoFar+                currentOverestimation = getOverestimation currentLine+        getTime :: String -> Double+        getTime line = +            case reads line of+                [(time,'s':_)] -> time+        getCaseId :: String -> String+        getCaseId line =+            case matchRegex idRegex line of+                Just [caseId] -> caseId+            where+            idRegex = mkRegex "case=([^;]*);"+        getOverestimation :: String -> Double+        getOverestimation line =+            case matchRegex ovestRegex line of+                Just [ovestS] -> read ovestS+            where+            ovestRegex = mkRegex "ovest=([^;]*);"+

@@ -0,0 +1,100 @@+{-|+    Module      :  Data.Number.ER.Real.Approx.Tests.Run+    Description :  (testing) running all function enclosure base tests in a batch+    Copyright   :  (c) 2007-2008 Michal Konecny+    License     :  BSD3++    Maintainer  :  mik@konecny.aow.cz+    Stability   :  experimental+    Portability :  portable+    +    Support for running all function enclosure base tests in a batch.+-}+module Data.Number.ER.Real.Approx.Tests.Run+where++import Data.Number.ER.Real.Approx.Tests.Generate+import Data.Number.ER.Real.Approx.Tests.Properties+import Data.Number.ER.Real.Approx.Tests.Reporting++import qualified Data.Number.ER.Real.Approx as RA+import qualified Data.Number.ER.Real.Approx.Elementary as RAEL++import Data.Number.ER.BasicTypes.DomainBox (VariableID(..), DomainBox, DomainBoxMappable, DomainIntBox)++import Data.Number.ER.Misc.Tests+import Data.Number.ER.Misc++import Test.QuickCheck+import Test.QuickCheck.Batch++import System.Directory+import qualified System.FilePath as FP+import Data.Time.Clock+import Data.Time.Calendar++runRATests :: +    (RAEL.ERApproxElementary ra, +     RAEL.ERInnerOuterApproxElementary ra,+     Ord ra) =>+    String -> ra -> IO () -> IO ()+runRATests title sampleRA initialise =+    do+    (UTCTime (ModifiedJulianDay days) secs) <- getCurrentTime+    let folder = "tests-" ++ title ++ "-" ++ (show days) ++ "-" ++ (show $ floor $ toRational secs)+    createDirectory folder+    erRunTests (title ++ " real approx tests") raTestOptions initialise (raTests sampleRA folder)++raTestOptions = +    TestOptions+      { +--        no_of_tests = 10+--        no_of_tests = 50+--        no_of_tests = 200+        no_of_tests = 500+      , +        length_of_tests = 240 * 3600 -- ie 4h time limit+      ,+        debug_tests = False +      }++raTests sampleRA folder =+    (propTuple "a-a=0" "AMinusAIsZero" props_ra_AMinusA_eq_oi)+    +++    (propTuple "a/a=1" "ADivAIsOne" props_ra_ADivA_eq_oi)+    +++    (propTuple "a+b=b+a" "AddCommut" props_ra_AddCommut_eq_oi)+    +++    (propTuple "a*b=b*a" "MultCommut" props_ra_MultCommut_eq_oi)+    +++    (propTuple "(a+b)+c=a+(b+c)" "AddAssoc" props_ra_AddAssoc_eq_oi)+    +++    (propTuple "(a*b)*c=a*(b*c)" "MultAssoc" props_ra_MultAssoc_eq_oi)+    +++    (propTuple "a*(b+c)=a*b+a*c" "Distr" props_ra_Distr_eq_oi)+    +++    (propTuple "log(exp(a))=a" "LogExp" props_ra_LogExp_eq_oi)+    +++    (propTuple "(sin a)^2+(cos a)^2=1" "SinCos" props_ra_SinCos_eq_oi)+    +++    (propTuple "tan(atan(a))=a" "TanATan" props_ra_TanATan_eq_oi)+    where+    propTuple testName testFileName propGen =+        [+            (testName  ++ ", equality, thin intervals", runR prop_eq_Thin $ filepath ++ "_Eq_Thin"),+            (testName  ++ ", equality, consistent intervals", runR prop_eq_Consistent $ filepath ++ "_Eq_Consistent"),+            (testName  ++ ", inner in outer, thin intervals", runR prop_oi_Thin $ filepath ++ "_OI_Thin"),+            (testName  ++ ", inner in outer, consistent intervals", runR prop_oi_Consistent $ filepath ++ "_OI_Consistent"),+            (testName  ++ ", inner in outer, directed intervals", runR prop_oi_Directed $ filepath ++ "_OI_Directed")+        ]+        where+        (prop_eq_Thin, prop_eq_Consistent, prop_oi_Thin, prop_oi_Consistent, prop_oi_Directed) = +            propGen sampleRA filepath+        filepath = FP.combine folder testFileName+    runR test filepath opts =+        unsafeReport filepath "started" $+            do +            result <- run test opts+            produceSummary filepath+            return result+

@@ -47,6 +47,7 @@  import qualified Data.Number.ER.Real.Approx as RA import Data.Number.ER.BasicTypes+import qualified Data.Number.ER.BasicTypes.ExtendedInteger as EI  import Data.Number.ER.Real.Arithmetic.Taylor -- import Data.Number.ER.Real.Arithmetic.Newton@@ -74,13 +75,10 @@     (RA.ERIntApprox ira, Ord ira) =>     EffortIndex ->      ira -> ira-erSqr_R ix a-    | RA.isEmpty a =-        RA.emptyApprox-    | otherwise = -        max 0 $ a' * a'+erSqr_R ix a =+    max 0 $ a' * a'     where-    a' = RA.setMinGranularity gran a+    a' = RA.setMinGranularityOuter gran a     gran = effIx2gran ix      {-@@ -109,8 +107,6 @@     Integer ->     ira -> ira erPow_R ix p a-    | RA.isEmpty a =-        RA.emptyApprox     | p < 0 =         1 / erPow_R ix (-p) a     | p == 0 = @@ -120,7 +116,7 @@     | otherwise =         a' * (erPow_R ix (div (p - 1) 2) (erSqr_R ix a'))     where-    a' = RA.setMinGranularity gran a+    a' = RA.setMinGranularityOuter gran a     gran = effIx2gran ix  {-@@ -145,7 +141,7 @@     EffortIndex -> ira -> ira erSqrt_IR_Inner =     RA.maxExtensionInnerR2R -        sqrtExtremaAndDirections+        sqrtExtrema         (\ ix x -> erSqrt_R ix x)  sqrtExtrema ix x = fst $ sqrtExtremaAndDirections ix x@@ -163,11 +159,11 @@     EffortIndex -> ira -> ira erSqrtContFr_R ix a     | aR == 0 = 0-    | aL == 1/0 = 1/0-    | aR `RA.ltSingletons` 0 = RA.emptyApprox+    | aL == RA.plusInfinity = RA.plusInfinity+    | aR `RA.ltSingletons` 0 = RA.topApprox     | otherwise =         contFrIter (ix + 3) $-            RA.setMinGranularity gran $ 0 RA.\/ aR -- assuming aR >= 0 +            RA.setMinGranularityOuter gran $ 0 RA.\/ aR -- assuming aR >= 0      where     gran = effIx2gran ix     (aL, aR) = RA.bounds a@@ -185,10 +181,9 @@     (RA.ERIntApprox ira) =>      EffortIndex -> ira -> ira erSqrtNewton_R ix a-    | RA.isEmpty a = RA.emptyApprox     | aR == 0 = 0-    | aL == 1/0 = 1/0-    | aR `RA.ltSingletons` 0 = RA.emptyApprox+    | aL == RA.plusInfinity = RA.plusInfinity+    | aR `RA.ltSingletons` 0 = RA.topApprox     | otherwise =         x_i RA.\/ (a/x_i)     where@@ -198,7 +193,7 @@          x_i =          newtonIter ((ix `div` 10) + 5) $-                RA.setMinGranularity gran aR -- assuming aR >= 0 +                RA.setMinGranularityOuter gran aR -- assuming aR >= 0      newtonIter i x_i         | i == 0 = x_i         | otherwise =@@ -230,7 +225,7 @@     EffortIndex -> Integer -> ira -> ira erRoot_IR_Inner ix p =     RA.maxExtensionInnerR2R -        (rootExtremaAndDirections p)+        (rootExtrema p)         (\ ix x -> erRoot_R ix p x) $             ix rootExtrema p ix x = fst $ rootExtremaAndDirections p ix x@@ -247,10 +242,9 @@     (RA.ERIntApprox ira, Ord 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+    | aL == RA.plusInfinity = RA.plusInfinity+    | aR < 0 && even p = RA.topApprox     | aR < 0 = - erRootNewton_R ix p (-a)     | p > 0 =         x_i RA.\/ (a/x_i_pow_p_minus_1)@@ -265,7 +259,7 @@          (x_i, x_i_pow_p_minus_1) =          newtonIter (ix + 5) $-                RA.setMinGranularity gran $ max 0 aR+                RA.setMinGranularityOuter gran $ max 0 aR       newtonIter i x_0         | i == 0 = @@ -294,10 +288,19 @@      erExp_R ix x      | RA.isBounded x =+--        unsafePrintReturn+--        (+--            "erExp_R: "+--            ++ "\n x = " ++ show x+--            ++ "\n xNear0 = " ++ show xNear0+--            ++ "\n n = " ++ show n+--            ++ "\n erExp_Tay_Opt_R ix xNear0 = " ++ (show $ erExp_Tay_Opt_R ix xNear0)+--            ++ "\n result = "+--        ) $         erPow_IR ix n $          erExp_Tay_Opt_R ix xNear0-    | x `RA.refines` (-1/0) = 0-    | (-1/0) `RA.refines` x =+    | x `RA.refines` (-RA.plusInfinity) = 0+    | (-RA.plusInfinity) `RA.refines` x =         0 RA.\/ (erExp_R ix (snd $ RA.bounds x))     | otherwise = RA.bottomApprox     where@@ -316,7 +319,7 @@      erExp_IR =     RA.maxExtensionR2R-        (\ ix x -> [])+        noExtrema         erExp_R  erExp_IR_Inner :: @@ -324,9 +327,11 @@     EffortIndex -> ira -> ira erExp_IR_Inner =     RA.maxExtensionInnerR2R-        (\ ix x -> ([], (Just True, Just True)))+        noExtrema         erExp_R +noExtrema ix x = []+ {- Log using Newton -}  erLog_R :: @@ -352,7 +357,7 @@      erLog_IR_Inner =     RA.maxExtensionInnerR2R-        logExtremaAndDirections+        logExtrema         (\ ix x -> logDivSeries_R ix x)          logExtrema ix x = fst $ logExtremaAndDirections ix x@@ -361,24 +366,40 @@     case RA.compareReals 0 x of         Just LT -> ([], (Just True, Just True))         Just GT -> ([], (Nothing, Nothing))-        _ -> ([-1/0], (Nothing, Just True))+        _ -> ([-RA.plusInfinity], (Nothing, Just True))          {-| log using a fast converging series, designed to be used with singletons -} logDivSeries_R ::     (RA.ERIntApprox ira, Ord 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+    | posx `RA.refines` 0 = -RA.plusInfinity+    | 0 `RA.refines` posx = RA.bottomApprox+    | posx `RA.refines` (RA.plusInfinity) = RA.plusInfinity     | otherwise =         case RA.compareReals posx 1 of             Just LT ->-                - (logDivSeries_R ix (recip posx))+--                unsafePrint +--                (+--                    "logDivSeries_R: recursion via recip" +--                ) $+                negate $+                    (logDivSeries_R ix posxRecipL) +                    RA.\/ +                    (logDivSeries_R ix posxRecipR)             _ ->-                nearLogx + 2 * t * (series ix (RA.setMinGranularity gran 1))+--                unsafePrint +--                (+--                    "logDivSeries_R: using series"+--                    ++ "\n posx = " ++ show posx +--                    ++ "\n nearLogx = " ++ show nearLogx +--                    ++ "\n remNearLogx = " ++ show remNearLogx +--                    ++ "\n t = " ++ show t +--                ) $+                nearLogx + 2 * t * (series ix (RA.setMinGranularityOuter gran 1))     where     gran = effIx2gran ix-    posx = (RA.setMinGranularity gran x) RA./\ (0 RA.\/ (1/0))+    posx = (RA.setMinGranularityOuter gran x) RA./\ (0 RA.\/ (RA.plusInfinity))+    (posxRecipL, posxRecipR) = RA.bounds $ recip posx     nearLogx =         0.69314718055994530941 * (fromInteger $ intLogUp 2 $ xCeiling)     remNearLogx =@@ -411,14 +432,14 @@ --        _ ->     --            erNewton_FullArgs  --                ( \ i y -> (erExp_RA i y) - posx, erExp_RA) ---                (RA.setMinGranularity gran nearLogx) ---                (RA.setMinGranularity gran 1) +--                (RA.setMinGranularityOuter gran nearLogx) +--                (RA.setMinGranularityOuter gran 1)  --                (fromInteger $ toInteger i) --                i --    where --    gran = effIx2gran i --    posx = ---        RA.setMinGranularity gran x /\ (ira2ra $ 0 RA.\/ (1/0))+--        RA.setMinGranularityOuter gran x /\ (ira2ra $ 0 RA.\/ (RA.plusInfinity)) --    nearLogx =                     --        0.69314718055994530941 * (fromInteger $ intLog 2 $ xCeiling) --    xCeiling @@ -490,7 +511,7 @@     EffortIndex -> ira -> ira  erSine_Tay_R ix x-    | (1/0) `RA.refines` x || (-1/0) `RA.refines` x = +    | (RA.plusInfinity) `RA.refines` x || (-RA.plusInfinity) `RA.refines` x =          (-1) RA.\/ 1       | otherwise =         erTaylor_R ix sine_coefSeq sine_error 0 x@@ -528,14 +549,14 @@     EffortIndex -> ira -> ira       erSine_IR_Inner = -    RA.maxExtensionInnerR2R sineExtremesAndDirections erSine_R+    RA.maxExtensionInnerR2R sineExtremes erSine_R      erCosine_IR_Inner ::     (RA.ERIntApprox ira) =>     EffortIndex -> ira -> ira       erCosine_IR_Inner = -    RA.maxExtensionInnerR2R cosineExtremesAndDirections erCosine_R+    RA.maxExtensionInnerR2R cosineExtremes erCosine_R          sineExtremes ix x = fst $ sineExtremesAndDirections ix x cosineExtremes ix x = fst $ cosineExtremesAndDirections ix x@@ -592,18 +613,14 @@     EffortIndex -> ira -> ira       erATan_IR =-    RA.maxExtensionR2R atanExtremes erATan_R--atanExtremes ix x = []+    RA.maxExtensionR2R noExtrema erATan_R  erATan_IR_Inner ::     (RA.ERIntApprox ira) =>     EffortIndex -> ira -> ira   erATan_IR_Inner =-    RA.maxExtensionInnerR2R atanExtremesAndDirections erATan_R--atanExtremesAndDirections ix x = ([], (Just True, Just True))+    RA.maxExtensionInnerR2R noExtrema erATan_R  {- 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 + ...)))))@@ -619,11 +636,13 @@     EffortIndex -> ira -> ira  atanEuler_R ix x -    | RA.isEmpty x = RA.emptyApprox+    | x `RA.refines` RA.plusInfinity = RA.plusInfinity  +    | x `RA.refines` (- RA.plusInfinity) = - RA.plusInfinity  +    | not $ RA.isBounded x = RA.bottomApprox      | x `RA.refines` ((-1.5) RA.\/ 1.5) =-        (x / xSquarePlus1) * (series ix (RA.setMinGranularity gran 2))+        (x / xSquarePlus1) * (series ix (RA.setMinGranularityOuter gran 2))     | otherwise = -- too far from 0, needs atan(x) = 2*atan(x/(1+sqrt(1+x^2)))-        2 * (atanEuler_R ix $ x / (1 + (erSqrt_R ix $ 1 + x * x)))+        2 * (atanEuler_R ix $ x / (1 + sqrtXQuarePlus1))     where     gran = effIx2gran ix     series termsCount coeffBase @@ -636,6 +655,16 @@     xSquare = abs $ x * x     xSquarePlus1 = xSquare + 1     xSquareOverXSquarePlus1 = xSquare / xSquarePlus1+    sqrtXQuarePlus1 =+        iterateIx 10 EI.MinusInfinity+        where+        iterateIx ix prevPrec +            | prevPrec == currentPrec = result+            | otherwise =+                iterateIx (ix * 2) currentPrec+            where +            result = erSqrt_R ix xSquarePlus1+            currentPrec = RA.getPrecision result       --{- atan using Newton -} --@@ -646,8 +675,8 @@ --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) +--        (RA.setMinGranularityOuter (effIx2gran i) (x))+--        (RA.setMinGranularityOuter (effIx2gran i) 1)  --        (fromInteger $ toInteger i) --        i @@ -718,9 +747,9 @@     -}      where     gran = max 0 (effIx2gran ix) + 10-    r1over64 = (RA.setMinGranularity gran 1) / 64-    r1over1024 = (RA.setMinGranularity gran 1) / 1024-    z = RA.setMinGranularity gran 0+    r1over64 = (RA.setMinGranularityOuter gran 1) / 64+    r1over1024 = (RA.setMinGranularityOuter gran 1) / 1024+    z = RA.setMinGranularityOuter gran 0     bellardTerms n nMax (mult, r4n, r10n)         | n >= nMax = []         | otherwise =

@@ -19,8 +19,8 @@ where  import qualified Data.Number.ER.Real.Approx as RA -import qualified Data.Number.ER.Real.DomainBox as DBox-import Data.Number.ER.Real.DomainBox (VariableID(..), DomainBox, DomainBoxMappable, DomainIntBox)+import qualified Data.Number.ER.BasicTypes.DomainBox as DBox+import Data.Number.ER.BasicTypes.DomainBox (VariableID(..), DomainBox, DomainBoxMappable, DomainIntBox) import Data.Number.ER.BasicTypes  import Data.List@@ -104,10 +104,10 @@     i = DBox.lookup "ER: LinearSolver: splitBox: " k box     k = widestVar box     incrementGranularity x =-        RA.setMinGranularity (RA.getGranularity x + 1) x+        RA.setMinGranularityOuter (RA.getGranularity x + 1) x  widestVar box =-    fst $ DBox.bestSplit box    +    fst $ DBox.bestSplit box  width i =     snd $ RA.bounds (iR-iL)

@@ -13,8 +13,9 @@ 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 qualified Data.Number.ER.BasicTypes.ExtendedInteger as EI import Data.Number.ER.BasicTypes+import Data.Number.ER.Misc   erTaylor_R@@ -41,7 +42,7 @@     -> ira      -> ira erTaylor_R_FullArgs coefSeq derivBounds n a gran x = -    rec_apTaylor (RA.setMinGranularity gran 0) 0+    rec_apTaylor (RA.setMinGranularityOuter gran 0) 0     where     rec_apTaylor i j         | n > j = (coefSeq(j)) + @@ -59,9 +60,24 @@     -> 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))+    | x `RA.refines` (-RA.plusInfinity) = +--        unsafePrintReturn+--        (+--            "erExp_Tay_Opt_R (x `RA.refines` (-RA.plusInfinity)): "+--            ++ "\n x = " ++ show x+--            ++ "\n ix = " ++ show ix+--            ++ "\n result = "+--        ) $+        0 -- -infty is not handled well by the Taylor formula+    | otherwise = +--        unsafePrintReturn+--        (+--            "erExp_Tay_Opt_R: "+--            ++ "\n x = " ++ show x+--            ++ "\n ix = " ++ show ix+--            ++ "\n result = "+--        ) $+        1 + (te ix x (RA.setMinGranularityOuter gran 1))     where     gran = effIx2gran ix     te steps x i@@ -87,13 +103,13 @@                     | xFloor == EI.MinusInfinity =                         0                     | otherwise = -                        ((RA.setMinGranularity gran 26)/10) ^^ xFloor +                        ((RA.setMinGranularityOuter gran 26)/10) ^^ xFloor                              -- lower estimate of e^x                 pow28xCeiling                      | xCeiling == EI.PlusInfinity =-                        (1/0)+                        (RA.plusInfinity)                     | otherwise = -                        ((RA.setMinGranularity gran 28)/10) ^^ xCeiling +                        ((RA.setMinGranularityOuter gran 28)/10) ^^ xCeiling                              -- upper estimate of e^x  {-@@ -109,7 +125,7 @@     -> ira     -> ira erSine_Tay_Opt_R ix x = -    taylor_seg ix x (RA.setMinGranularity gran 1)+    taylor_seg ix x (RA.setMinGranularityOuter gran 1)     where     gran = effIx2gran ix     taylor_seg i x n -- 'i' for iterator@@ -127,7 +143,7 @@     => EffortIndex      -> ira     -> ira-erCosine_Tay_Opt_R ix x = taylor_seg ix x (RA.setMinGranularity gran 1)+erCosine_Tay_Opt_R ix x = taylor_seg ix x (RA.setMinGranularityOuter gran 1) 	where     gran = effIx2gran ix     taylor_seg i x n -- 'i' for iterator

@@ -22,15 +22,15 @@ where  import Data.Number.ER.BasicTypes-import qualified Data.Number.ER.ExtendedInteger as EI+import qualified Data.Number.ER.BasicTypes.ExtendedInteger as EI  import Data.Typeable  {-|     This class is an abstraction of a subset of real numbers-    with upwards rounded operations. +    with *upwards rounded* operations.  -}-class (Fractional rb, Ord rb) => ERRealBase rb +class (RealFrac rb, Ord rb) => ERRealBase rb      where     typeName :: rb -> String     initialiseBaseArithmetic :: rb -> IO ()@@ -49,9 +49,14 @@     isERNaN :: rb -> Bool     erNaN :: rb     isPlusInfinity :: rb -> Bool+    isMinusInfinity :: rb -> Bool+    isMinusInfinity = isPlusInfinity . negate     plusInfinity :: rb     minusInfinity :: rb     minusInfinity = - plusInfinity+    fromIntegerUp :: Integer -> rb+    fromIntegerDown :: Integer -> rb+    fromIntegerDown i = negate $ fromIntegerUp $ - i     fromDouble :: Double -> rb     toDouble :: rb -> Double     fromFloat :: Float -> rb

@@ -23,7 +23,7 @@ where  import qualified Data.Number.ER.Real.Base as B-import qualified Data.Number.ER.ExtendedInteger as EI+import qualified Data.Number.ER.BasicTypes.ExtendedInteger as EI import Data.Number.ER.Real.Base.MachineDouble import Data.Number.ER.Real.Base.Float import Data.Number.ER.BasicTypes@@ -232,6 +232,7 @@     isPlusInfinity =          op1ERMachineAP (== 1/0) B.isPlusInfinity     plusInfinity = B.fromDouble $ 1/0+    fromIntegerUp = fromInteger     fromDouble d =          ERMachineAPMachineDouble (B.defaultGranularity (0 :: b)) d     toDouble = op1ERMachineAP id B.toDouble

@@ -19,8 +19,8 @@ ) where -import qualified Data.Number.ER.ExtendedInteger as EI-import Data.Number.ER.PlusMinus+import qualified Data.Number.ER.BasicTypes.ExtendedInteger as EI+import Data.Number.ER.BasicTypes.PlusMinus import Data.Number.ER.Misc import Data.Number.ER.BasicTypes import qualified Data.Number.ER.Real.Base as B@@ -494,7 +494,8 @@     erNaN = ERFloatNaN (B.defaultGranularity zero)     isPlusInfinity (ERFloatInfty _ Plus) = True     isPlusInfinity _ = False-    plusInfinity = ERFloatInfty (B.defaultGranularity zero) Plus    +    plusInfinity = ERFloatInfty (B.defaultGranularity zero) Plus+    fromIntegerUp i = fromInteger i         fromDouble d         | isNaN d = ERFloatNaN (B.defaultGranularity zero)         | otherwise = (fromRational . toRational) d

@@ -24,7 +24,7 @@ where  import qualified Data.Number.ER.Real.Base as B-import qualified Data.Number.ER.ExtendedInteger as EI+import qualified Data.Number.ER.BasicTypes.ExtendedInteger as EI import Data.Number.ER.Misc  import Data.Binary@@ -65,6 +65,7 @@     erNaN = 0/0     isPlusInfinity x = M.isInfinite x && x > 0     plusInfinity = 1/0+    fromIntegerUp = fromInteger     fromDouble = M.fromDouble M.Up 53     toDouble = M.toDouble M.Up     fromFloat = B.fromDouble . fromRational . toRational  

@@ -12,14 +12,14 @@      Make 'Double' an instance of 'B.ERRealBase' as much as possible.     -}-module Data.Number.ER.Real.Base.MachineDouble +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 qualified Data.Number.ER.BasicTypes.ExtendedInteger as EI import Data.Number.ER.Misc  import Foreign.C@@ -78,6 +78,8 @@             fromInteger $ intLogUp 2 $ ceiling d         | d < 1 =             negate $ fromInteger $ intLogUp 2 $ ceiling $ recip d+        | otherwise = +            error $ "ER.Real.Base.MachineDouble: getApproxBinaryLog: illegal argument " ++ show d      getGranularity _ = 53     setMinGranularity _ = id     setGranularity _ = id@@ -86,6 +88,14 @@     erNaN = 0/0     isPlusInfinity f = isInfinite f && f > 0     plusInfinity = 1/0+    fromIntegerUp i+        | i <= floor nearest = nearest+        | otherwise = nearestIncreased+        where+        nearestCeil = ceiling nearest+        nearest = fromInteger i+        nearestIncreased = encodeFloat (s+1) e+        (s,e) = decodeFloat nearest     fromDouble = fromRational . toRational     toDouble = fromRational . toRational     fromFloat = fromRational . toRational

@@ -24,8 +24,8 @@ 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 qualified Data.Number.ER.BasicTypes.ExtendedInteger as EI+import Data.Number.ER.BasicTypes.PlusMinus import Data.Number.ER.Misc  import Data.Ratio@@ -226,6 +226,7 @@     isPlusInfinity (Infinity Plus) = True     isPlusInfinity _ = False     plusInfinity = Infinity Plus+    fromIntegerUp = fromInteger     fromDouble = fromRational . toRational     toDouble (Infinity Plus) = 1/0      toDouble (Infinity Minus) = -1/0 

@@ -0,0 +1,90 @@+{-|+    Module      :  Data.Number.ER.Real.Base.Tests.Generate+    Description :  (testing) generating base real numbers+    Copyright   :  (c) 2009 Michal Konecny+    License     :  BSD3++    Maintainer  :  mik@konecny.aow.cz+    Stability   :  experimental+    Portability :  portable+    +    Generic instances of 'Arbitrary' class for generating (almost) random instances+    according to different distributions. +-}+module Data.Number.ER.Real.Base.Tests.Generate where++import qualified Data.Number.ER.Real.Base as B+import Data.Number.ER.BasicTypes++import Test.QuickCheck++newtype BGran20 b = BGran20 b deriving Show+newtype BGran100 b = BGran100 b deriving Show+newtype BGran1000 b = BGran1000 b deriving Show++instance (B.ERRealBase b) => Arbitrary (BGran20 b)+    where+    arbitrary =+        do+        gran <- choose (8,20)+        (f1,f2,f3) <- arbitrary+        pow <- choose (-10,10)+        return $ BGran20 $ constructB gran (f1,f2,f3) pow+    coarbitrary _ =+        error "ER.Real.Base: Tests: coarbitrary not implemented"++constructB ::+    (B.ERRealBase b) =>+    Granularity ->+    (Double, Double, Double) ->+    Int ->+    b+constructB gran (f1,f2,f3) pow =+    (b1/b2) ^^ pow + b3 +    where+    [b1,b2,b3] = map cvt [f1,f2,f3]+    cvt f = B.setGranularity gran $ B.fromDouble f++instance (B.ERRealBase b) => Arbitrary (BGran100 b)+    where+    arbitrary = +        sized arbitrarySized+        where+        arbitrarySized n +            | n <= 28 =+                do+                (BGran20 b) <- arbitrary+                return (BGran100 b)+            | otherwise =+                do+                gran <- choose (30,100)+                (f1,f2,f3) <- arbitrary+                pow <- choose (-100,100)+                return $ BGran100 $ constructB gran (f1,f2,f3) pow+    coarbitrary _ =+        error "ER.Real.Base: Tests: coarbitrary not implemented"++instance (B.ERRealBase b) => Arbitrary (BGran1000 b)+    where+    arbitrary = +        sized arbitrarySized+        where+        arbitrarySized n +            | n <= 28 =+                do+                (BGran20 b) <- arbitrary+                return (BGran1000 b)+            | n <= 68 =+                do+                (BGran100 b) <- arbitrary+                return (BGran1000 b)+            | otherwise =+                do+                gran <- choose (400,1000)+                (f1,f2,f3) <- arbitrary+                pow <- choose (-10000,10000)+                return $ BGran1000 $ constructB gran (f1,f2,f3) pow+    coarbitrary _ =+        error "ER.Real.Base: Tests: coarbitrary not implemented"+            +            

@@ -1,181 +0,0 @@-{-# LANGUAGE MultiParamTypeClasses  #-}-{-# LANGUAGE FunctionalDependencies  #-}-{-|-    Module      :  Data.Number.ER.Real.DomainBox-    Description :  portions of many-dimensional domains   -    Copyright   :  (c) Michal Konecny-    License     :  BSD3--    Maintainer  :  mik@konecny.aow.cz-    Stability   :  experimental-    Portability :  portable--    Abstractions of the 'Box' datatype, often used to represent-    sections of multi-dimensional function domains.-    -    To be imported qualified, usually with prefix DBox.-    -    VariableID(..) and DomainBox -    are usually imported separately and not qualified.--}-module Data.Number.ER.Real.DomainBox -(-    VariableID(..),-    DomainBox(..),-    DomainBoxMappable(..),-    DomainIntBox(..)-)-where--import qualified Data.Number.ER.Real.Approx as RA--import Data.Number.ER.BasicTypes--import qualified Data.Set as Set-import qualified Data.Map as Map--import Prelude hiding (lookup)---{-| -    A class abstracting a type of variable identifiers -    for axes in function domains, polynomials etc.--}-class (Ord varid) => VariableID varid-    where-    newVarID :: Set.Set varid -> varid-    defaultVar :: varid-    defaultVar = newVarID Set.empty-    showVar :: varid -> String--{-|-    A class abstracting a type of many-dimensional points, intervals-    or anything indexed by a subset of dimensions.-    -    More generally, this class abstracts most of 'Data.Map.Map'.--}-class (VariableID varid) => DomainBox box varid val-    | box -> varid val, varid val -> box-    where-    noinfo :: box-    isNoinfo :: box -> Bool-    size :: box -> Int-    {-| constructor using 'defaultVar' -}-    unary :: val -> box-    singleton :: varid -> val -> box-    toList :: box -> [(varid, val)]-    fromList :: [(varid, val)] -> box-    toAscList :: box -> [(varid, val)]-    fromAscList :: [(varid, val)] -> box---    toMap :: box -> Map.Map varid val---    fromMap :: Map.Map varid val -> box-    compare :: (val -> val -> Ordering) -> box -> box -> Ordering-    insert :: varid -> val -> box -> box-    insertWith :: (val -> val -> val) -> varid -> val -> box -> box-    delete :: varid -> box -> box-    member :: varid -> box -> Bool-    notMember :: varid -> box -> Bool-    union :: box -> box -> box-    unionWith :: (val -> val -> val) -> box -> box -> box-    keys :: box -> [varid]-    elems :: box -> [val]-    filter :: (val -> Bool) -> box -> box-    fold :: (val -> a -> a) -> a -> box -> a-    foldWithKey :: (varid -> val -> a -> a) -> a -> box -> a-    {-| -        for all variables that appear in both boxes,-        apply the function and add the result to the list -     -}-    zipWith :: (val -> val -> a) -> box -> box -> [(varid, a)] -    {-| -        For all variables that appear in either of the two boxes,-        apply the function and add the result to the list.-        -        Supply the default value when the variable is missing from either box. -     -}-    zipWithDefault :: val -> (val -> val -> a) -> box -> box -> [(varid, a)] -    {-| -        For all variables that appear in the first box,-        apply the function and add the result to the list.-        -        Supply the default value when the variable is missing from the second box. -     -}-    zipWithDefaultSecond :: val -> (val -> val -> a) -> box -> box -> [(varid, a)] -    findWithDefault :: val -> varid -> box -> val-    {-|-        Pick the extents of a single variable in a domain box.-        If there is no information for this variable, assume the-        variable ranges over the whole real line.-    -}-    lookup ::     -        String {-^ identification of caller location to use in error messages -} ->-        varid ->-        box ->-        val-        -{-|-    A class linking two domain box types that share the-    index type so that boxes of the two types can be-    converted etc.--}-class (DomainBox box1 varid val1, DomainBox box2 varid val2) => -    DomainBoxMappable box1 box2 varid val1 val2-    where-    map :: (val1 -> val2) -> box1 -> box2-    mapWithKey :: (varid -> val1 -> val2) -> box1 -> box2-    intersectionWith :: (val1 -> val2 -> val1) -> box1 -> box2 -> box1-    difference :: box1 -> box2 -> box1 --{-|-    A class abstracting a type of many-dimensional intervals.--}-class (DomainBox box varid ira) => DomainIntBox box varid ira-    | box -> varid ira, varid ira -> box-    where-    {-|-        Check whether the two domains specify the same-        interval for each variable that they share.-    -}-    compatible ::-        box ->-        box ->-        Bool-    {-|-        Assuming that two domains are compatible, take the-        most information from both of the domains about the-        ranges of variables.-    -}-    unify ::-        String {-^ identification of caller location to use in error messages -} ->-        box ->-        box ->-        box-    {-|-        Find the variable with the largest interval-        and return it together with the default splitting point-        in its domain.-    -}-    bestSplit ::-        box ->-        (varid, ira)-    split ::-        box ->-        varid ->-        ira ->-        (box, box)-    classifyPosition ::-        box {-^ domain @d1@ -} ->-        box {-^ domain @d2@ -} ->-        (Bool, Bool, Bool, Bool) -            {-^ -                Answers to these (mutually exclusive) questions:-                -                * is @d1@ outside and /not/ touching @d2@?-            -                * is @d1@ outside and touching @d2@?-            -                * is @d1@ intersecting and not inside @d2@?-            -                * is @d1@ inside @d2@?-            -}-            

@@ -1,196 +0,0 @@-{-# LANGUAGE MultiParamTypeClasses  #-}-{-# LANGUAGE FlexibleInstances   #-}-{-# LANGUAGE TypeSynonymInstances   #-}-{-|-    Module      :  Data.Number.ER.Real.DomainBox.IntMap-    Description :  implementation of DomainBox based on Data.IntMap   -    Copyright   :  (c) Michal Konecny-    License     :  BSD3--    Maintainer  :  mik@konecny.aow.cz-    Stability   :  experimental-    Portability :  portable--    A simple implementation of the 'VariableID' and 'DomainBox' classes.--}-module Data.Number.ER.Real.DomainBox.IntMap -(-    VarID, Box-)-where--import qualified Data.Number.ER.Real.Approx as RA-import qualified Data.Number.ER.Real.DomainBox as DBox-import Data.Number.ER.Real.DomainBox (VariableID(..), DomainBox, DomainBoxMappable, DomainIntBox)--import Data.Number.ER.Misc--import qualified Data.IntMap as IMap-import qualified Data.Set as Set--type VarID = Int-type Box ira = IMap.IntMap ira--instance VariableID VarID-    where-    newVarID prevVars -        | Set.null prevVars = 0-        | otherwise =-            1 + (Set.findMax prevVars)-    showVar v-        | v == 0 = "x"-        | otherwise = "x" ++ show v--instance (Show val) => (DomainBox (Box val) VarID val)-    where-    noinfo = IMap.empty-    isNoinfo = IMap.null-    size = IMap.size-    unary r = IMap.singleton defaultVar r-    singleton = IMap.singleton-    toList = IMap.toList-    fromList = IMap.fromList-    toAscList = IMap.toAscList-    fromAscList = IMap.fromAscList---    toMap = id---    fromMap = id-    compare compareVals b1 b2 =-        compareListsWith comparePairs (IMap.toList b1) (IMap.toList b2)-        where-        comparePairs (k1,v1) (k2,v2) =-            compareComposeMany-                [-                    compare k1 k2,-                    compareVals v1 v2-                ]-             -    insert = IMap.insert-    insertWith = IMap.insertWith-    delete = IMap.delete-    member = IMap.member -    notMember = IMap.notMember-    union = IMap.union -    unionWith = IMap.unionWith -    elems = IMap.elems-    keys = IMap.keys-    filter = IMap.filter-    fold = IMap.fold-    foldWithKey = IMap.foldWithKey-    zipWith f b1 b2 = -        applyF (IMap.toAscList b1) (IMap.toAscList b2)-        where-        applyF [] _ = []-        applyF _ [] = []-        applyF bl1@((k1,v1):rest1) bl2@((k2,v2):rest2) -            | k1 == k2 = -                (k1, f v1 v2) : (applyF rest1 rest2)-            | k1 < k2 = applyF rest1 bl2-            | otherwise = applyF bl1 rest2 -    zipWithDefault defaultValue f b1 b2 = -        applyF (IMap.toAscList b1) (IMap.toAscList b2)-        where-        applyF [] [] = []-        applyF bl1@((k1,v1):rest1) [] =-            (k1, f v1 defaultValue) : (applyF rest1 [])-        applyF [] bl2@((k2,v2):rest2) =-            (k2, f defaultValue v2) : (applyF [] rest2)-        applyF bl1@((k1,v1):rest1) bl2@((k2,v2):rest2) -            | k1 == k2 = -                (k1, f v1 v2) : (applyF rest1 rest2)-            | k1 < k2 = -                (k1, f v1 defaultValue) : (applyF rest1 bl2)-            | otherwise =  -                (k2, f defaultValue v2) : (applyF bl1 rest2)-    zipWithDefaultSecond defaultValue f b1 b2 = -        applyF (IMap.toAscList b1) (IMap.toAscList b2)-        where-        applyF [] _ = []-        applyF bl1@((k1,v1):rest1) [] =-            (k1, f v1 defaultValue) : (applyF rest1 [])-        applyF bl1@((k1,v1):rest1) bl2@((k2,v2):rest2) -            | k1 == k2 = -                (k1, f v1 v2) : (applyF rest1 rest2)-            | k1 < k2 = -                (k1, f v1 defaultValue) : (applyF rest1 bl2)-            | otherwise =  -                applyF bl1 rest2-    findWithDefault = IMap.findWithDefault-    lookup locspec var dom =-        IMap.findWithDefault err var dom-        where-        err =-            error $-                locspec ++ "DomainBox.IntMap lookup: domain box " ++ show dom -                ++ " ignores variable " ++ show var--instance (Show val1, Show val2) => -    (DomainBoxMappable (Box val1) (Box val2) VarID val1 val2)-    where-    map = IMap.map-    mapWithKey = IMap.mapWithKey-    intersectionWith = IMap.intersectionWith-    difference = IMap.difference--instance (RA.ERIntApprox ira) => DomainIntBox (Box ira) VarID ira-    where-    compatible dom1 dom2 =-        foldl (&&) True $ map snd $-            DBox.zipWith RA.equalIntervals dom1 dom2-    unify locspec dom1 dom2-        | DBox.compatible dom1 dom2 =-            IMap.union dom1 dom2-        | otherwise =-            error $-                locspec ++ "incompatible domains " ++ show dom1 ++ " and " ++ show dom2-    bestSplit dom =-        (var, pt)-        where-        pt = -            RA.defaultBisectPt varDom-        (_, (varDom, var)) = -            foldl findWidestVar (0, err) $ IMap.toList dom-        err =-            error $ "DomainBox: bestSplit: failed to find a split for " ++ show dom -        findWidestVar (prevWidth, prevRes) (v, d)-            | currWidth `RA.leqSingletons` prevWidth = (prevWidth, prevRes)-            | otherwise = (currWidth, (d, v))-            where-            currWidth = snd $ RA.bounds $ domHI - domLO-            (domLO, domHI) = RA.bounds d-    split dom var pt =-        (IMap.insert var varDomL dom, IMap.insert var varDomR dom)-        where-        varDomL = varDomLO RA.\/ varDomMid-        varDomR = varDomMid RA.\/ varDomHI-        (varDomLO, varDomMid, varDomHI, _) = RA.exactMiddle varDom-        varDom = DBox.lookup "DomainBox.IntMap: split: " var dom-    classifyPosition dom sdom =    -        (away, touch, intersect, inside)-            where-            (away, touch, inside, intersect) =-                foldl addDimension (True, True, True, False) awayTouchInsides-            addDimension -                    (prevAway, prevTouch, prevInside, prevIntersect) -                    (thisAway, thisTouch, thisInside, thisIntersect) =-                (prevAway && thisAway, -                 (prevTouch || prevAway) && (thisTouch || thisAway) && (prevTouch || thisTouch),-                 prevInside && thisInside,-                 prevIntersect || thisIntersect)-            awayTouchInsides =-                map snd $-                    DBox.zipWith classifyRA dom sdom-            classifyRA d sd =-                (outsideNoTouch, outsideTouch, inside,-                 not (outsideNoTouch || outsideTouch || inside))-                 where-                 outsideNoTouch = sdR < dL || dR < sdL-                 outsideTouch = sdR == dL || dR == sdL-                 inside = sdL =< dL && dR =< sdR-                 (==) = RA.eqSingletons-                 (<) = RA.ltSingletons-                 (=<) = RA.leqSingletons-                 (dL, dR) = RA.bounds d -                 (sdL, sdR) = RA.bounds sd -        --    

@@ -1,149 +0,0 @@-{-# LANGUAGE CPP #-}-{-| -    Module      :  Main-    Description :  simple examples of using AERN-Real-    Copyright   :  (c) Michal Konecny-    License     :  BSD3--    Maintainer  :  mik@konecny.aow.cz-    Stability   :  experimental-    Portability :  portable--    Simple examples of using AERN-Real--}-module Main where--import qualified Data.Number.ER.Real as AERN-import Data.Number.ER.Real (ConvergRealSeq(..), convertFuncRA2Seq)--#ifdef USE_MPFR---type B = AERN.BAP -- use pure Haskell floats-type B = AERN.BMAP -- use combination of double and pure Haskell floats---type B = AERN.BMPFR -- use MPFR floats-#else---type B = AERN.BAP -- use pure Haskell floats-type B = AERN.BMAP -- use combination of double and pure Haskell floats-#endif-type RA = AERN.RA B-type IRA = AERN.IRA B-type R = ConvergRealSeq IRA--one :: R-one = 1--two :: R-two = 2--piSeq :: R-piSeq = ConvergRealSeq $ AERN.pi--seqExp = convertFuncRA2Seq $ AERN.exp-seqSine = convertFuncRA2Seq $ AERN.sin-seqCosine = convertFuncRA2Seq $ AERN.cos--main = -    do-    AERN.initialiseBaseArithmetic (0 :: RA)-    putStrLn "****************************"-    putStrLn "Testing interval arithmetic:"-    putStrLn "****************************"-    putStrLn "**** Fractions:"-    putStrLn $-        "(default granularity, show internals) 1/3 =\n  " ++ -        AERN.showApprox 30 True True (1/3 :: RA) -    putStrLn $-        "(granularity 50, show internals) 1/3 =\n  " ++ -        AERN.showApprox 30 True True ((AERN.setGranularity 50 1/3) :: RA) -    putStrLn $-        "(granularity 100, show internals) 1/3 =\n  " ++ -        AERN.showApprox 40 True True ((AERN.setGranularity 100 1/3) :: RA) -    putStrLn $-        "(granularity 100, do not show internals) 1/3 =\n  " ++ -        AERN.showApprox 40 True False ((AERN.setGranularity 100 1/3) :: RA) -    putStrLn $-        "(granularity 100, default show) 1/3 =\n  " ++ -        show ((AERN.setGranularity 100 1/3) :: RA) -    putStrLn "**** Exp:"-    putStrLn $ -        "(effort 5, granularity 50) exp 1 =\n  " ++ -        (show $ AERN.exp 5 (AERN.setGranularity 50 (1::RA)))-    putStrLn $ -        "(effort 10, granularity 50) exp 1 =\n  " ++ -        (show $ AERN.exp 10 (AERN.setGranularity 50 (1::RA)))-    putStrLn $-        "(effort 10, granularity 100) exp 1 =\n  " ++ -        (show $ AERN.exp 10 (AERN.setGranularity 100 (1::RA)))-    putStrLn $ -        "(effort 20, granularity 50) exp 1 =\n  " ++ -        (show $ AERN.exp 20 (AERN.setGranularity 50 (1::RA)))-    putStrLn $-        "(effort 20, granularity 100) exp 1 =\n  " ++ -        (show $ AERN.exp 20 (AERN.setGranularity 100 (1::RA)))-    putStrLn "**** Pi:"-    putStrLn $ -        "(effort 10) pi =\n  " ++ -        (show $ (AERN.pi 10 :: RA))-    putStrLn $ -        "(effort 50) pi =\n  " ++ -        (AERN.showApprox 20 True False $ (AERN.pi 50 :: RA))-    putStrLn $ -        "(effort 100) pi =\n  " ++ -        (AERN.showApprox 35 True False $ (AERN.pi 100 :: RA))-    putStrLn $ -        "(effort 200) pi =\n  " ++ -        (AERN.showApprox 65 True False $ (AERN.pi 200 :: RA))-    putStrLn $ -        "(effort 400) pi =\n  " ++ -        (AERN.showApprox 125 True False $ (AERN.pi 400 :: RA))-    putStrLn "**** Sine:"-    putStrLn $-        "(effort 10, granularity 50) sin 1 =\n  " ++ -        (show $ AERN.sin 10 (AERN.setGranularity 50 (1::RA)))-    putStrLn $-        "(effort 10, granularity 100) sin 1 =\n  " ++ -        (show $ AERN.sin 10 (AERN.setGranularity 100 (1::RA)))-    putStrLn "**** Integration:"-    putStrLn $ -        "(effort 10, granularity 50) integrate exp 0 1 =\n  " ++ -        (show $ AERN.integrateContAdapt_R AERN.exp 10 0 (AERN.setGranularity 50 (1::RA)))-    putStrLn $ -        "(effort 20, granularity 50) integrate exp 0 1 =\n  " ++ -        (show $ AERN.integrateContAdapt_R AERN.exp 20 0 (AERN.setGranularity 50 (1::RA)))---    putStrLn $ ---        "(effort 30, granularity 50) integrate exp 0 1 =\n  " ++ ---        (show $ AERN.integrateContAdapt_R AERN.exp 30 0 (AERN.setGranularity 50 (1::RA)))-    putStrLn "*****************************"-    putStrLn "Testing convergent sequences:"-    putStrLn "*****************************"---    putStrLn $ "1 =\n  " ++ show one---    putStrLn $ "1 + 2 =\n  " ++ (show $ one + two)-    putStrLn "**** Fractions:"-    putStrLn $ -        "(precision 20) 1/3 =\n  " ++ -        (AERN.showConvergRealSeqAuto 20 $ one / 3)-    putStrLn $ -        "(precision 20) 100000000001/300000000000 =\n  " ++ -        (AERN.showConvergRealSeqAuto 20 $ (one + 100000000000)/300000000000 )-    putStrLn $ -        "100000000001/300000000000 =? 1/3:\n  " ++ -        (show $ one/3 == 100000000001/300000000000)---    putStrLn $ "abs -1 = " ++ (show $ abs (- one))---    putStrLn $ "neg 2 = " ++ (show $ negate two)---    putStrLn $ "1 + 2 = " ++ (show $ one + 2)-    putStrLn "**** Elementary:"-    putStrLn $ -        "(precision 30) exp 1 =\n  " ++ -        (AERN.showConvergRealSeqAuto 30 $ seqExp one)-    putStrLn $ -        "(precision 500) pi =\n  " ++ -        (AERN.showConvergRealSeqAuto 500 $ piSeq)-    putStrLn $ -        "(precision 30) cosine(1) =\n  " ++ -        (AERN.showConvergRealSeqAuto 30 $ seqCosine one)    -    putStrLn $-        "(precision 30) sine(1) =\n  " ++ -        (AERN.showConvergRealSeqAuto 30 $ seqSine one)-    putStrLn "**** Integration:"-    putStrLn $ -- very slow for precision > 4-        "(precision 3) integrate exp 0 1 =\n  " ++ -        (AERN.showConvergRealSeqAuto 3 $ AERN.integrateCont AERN.exp 0 one)

@@ -1,384 +0,0 @@-{-# LANGUAGE CPP #-}-{-# LANGUAGE UndecidableInstances #-}-{-# LANGUAGE TypeOperators #-}-{-# LANGUAGE DeriveDataTypeable #-}-module Main--where--import qualified Data.Number.ER.Real as AERN-import Data.Number.ER.BasicTypes-import Data.Number.ER.Misc--import Data.Maybe-import qualified Data.List as List-import qualified Data.Map as Map--import qualified Data.Array.IArray as IAr-import qualified Data.Array.MArray as MAr-import qualified Data.Array.ST as STAr-import qualified Data.Ix as Ix-import qualified Data.Array.Base as BAr--import Control.Monad.ST-import GHC.Arr--#ifdef USE_MPFR-type B = AERN.BMPFR -- use MPFR floats-#else-type B = AERN.BAP -- use pure Haskell floats-#endif-type RA = AERN.RA B-type IRA = AERN.IRA B--testMatrixN = 100-incrementGran = (+) 50---- Hilbert 100x100 matrix:-addOneDiag = False-targetPrec = 167 -- approx 50 decimal digits after the point-initialGran = 2050 -- 100x100---initialGran = 2388 -- 100x100 Norbert's---initialGran = 750 -- 50x50---initialGran = 300 -- 10x10----targetPrec = 34 -- approx 10 decimal digits after the point---initialGran = 1350---initialGran = 50 -- 50x50---- Hilbert matrix + 1:---addOneDiag = True---targetPrec = 167 -- approx 50 decimal digits after the point---initialGran = 200----targetPrec = 34 -- approx 10 decimal digits after the point---initialGran = 50--main =-    do-    AERN.initialiseBaseArithmetic (0 :: RA)-    putStrLn $ -          "Inverting the " ++ show n ++ "x" ++ show n ++ " Hilbert matrix " -          ++ "with target binary precision " ++ show targetPrec ++ "..." ---    putStrLn $ ---        "sorted matrix elements = \n" ++ (unlines $ map show elemsSortedByPrec)-    putStrLn $ -        "sum of all elements in inverted matrix = " ++ show (sum elems)---    putStrLn $ show (Matrix n n rarr)-    where-    n = testMatrixN-    elems = IAr.elems rarr-    elemsSortedByPrec =-        List.sortBy comparePrec elems-        where-        comparePrec a b =-            compare aPrecLO bPrecLO-            where-            aPrecLO = fst $ AERN.bounds $ aHI - aLO-            (aLO, aHI) = AERN.bounds a-            bPrecLO = fst $ AERN.bounds $ bHI - bLO-            (bLO, bHI) = AERN.bounds b-    rarr =-        STAr.runSTArray $-            do-            mInv@(Matrix _ _ rowsInv) <- -                invert testMatrix---            m <- testMatrix initialGran---            mUnit@(Matrix _ _ rowsUnit) <- multM m mInv-            return rowsInv---testMatrix ::-    Granularity -> -    ST s (STMatrix s IRA)-testMatrix gran =-    do-    marr <- MAr.newArray ((1,1),(n,n)) 0-    mapM (updateCell marr) assocsGran-    return $ Matrix n n marr-    where-    assocsGran = map (mapSnd $ AERN.setMinGranularity gran) assocs-    assocs = ---        assocsMini-        assocsHilbert gran n-    assocsMini = -        [((1,1),1),-         ((1,2),3),-         ((2,1),2),-         ((2,2),0)-        ]-    n = testMatrixN-    updateCell marr (ix, el) =-        do-        unsafeMatrixWrite marr n ix el --assocsHilbert gran n =-    [((i,j), coeff i j)| i <- [1..n], j <- [1..n]]-    where-    coeff i j -        | addOneDiag && i == j = -            1 + oneOverIplusJ-        | otherwise =-            oneOverIplusJ-        where-        oneOverIplusJ =-            recip $ (AERN.setMinGranularity gran $ iRA + jRA + 1)-        iRA = fromInteger $ toInteger i-        jRA = fromInteger $ toInteger j--    ---invert ::---    Precision ->---    () ->-invert getMatrix =-    do-    gaussElim getMatrixI-    where-    n = testMatrixN-    getMatrixI gran =-        do-        m <- getMatrix gran-        mI <- addIdentity m-        return mI--gaussElim getMatrix =-    elimWithMinGran initialGran-    where-    elimWithMinGran workingGran =-        do-        mI@(Matrix colN rowN _) <- getMatrix workingGran-        idPerm <- MAr.newListArray (1,rowN) [1..rowN]-        elimAtRow mI 1 idPerm-        where-        elimAtRow mI@(Matrix colN rowN mIarr) i perm =-            do-            success <- ensureNonZeroDiag -- make sure (i,i) is non-zero by permuting-            case success of-                False -> -- failed - all elements contain 0 -> try larger granularity-                    unsafePrint ("failed to divide at granularity " ++ show workingGran) $-                        elimWithMinGran (incrementGran workingGran)-                True ->-                    do-                    normaliseRow-                    eliminateColumn-                    case i == rowN of-                        True -> -                            do-                            mInv <- permuteRowsDropCols perm testMatrixN mI-                            mPrec <- getMatrixPrecision mInv-                            case mPrec >= targetPrec of-                                False -> -- resulting precision insufficient-                                    unsafePrint -                                    ("insufficient precision " ++ show mPrec ++  -                                     " at granularity " ++ show workingGran) $-                                        elimWithMinGran (incrementGran workingGran)-                                True -> -                                    unsafePrint -                                    ("precision " ++ show mPrec ++ -                                     " succeeded at granularity " ++ show workingGran)-                                    return mInv-                        False -> elimAtRow mI (i+1) perm-            where-            ensureNonZeroDiag =-                do-                maybeNonZeroIx <- findNonZeroRow-                case maybeNonZeroIx of-                    Nothing ->-                        return False-                    Just ii ->-                        do-                        case ii > 0 of-                            True -> swap i (i + ii) perm-                            False -> return ()-                        return True-            findNonZeroRow =-                do-                elems <- mapM getElemPerm [(i,rowIx) | rowIx <- [i..rowN]]-                return $ List.findIndex (\e -> not $ 0 `AERN.refines` e) elems-            getElemPerm (colIx,rowIx) =-                do-                rowIxPerm <- unsafePermRead perm rowIx-                unsafeMatrixRead mIarr rowN (colIx, rowIxPerm)--            normaliseRow =-                do-                rowIxPerm <- unsafePermRead perm i-                e <- unsafeMatrixRead mIarr rowN (i, rowIxPerm)-                unsafeMatrixWrite mIarr rowN (i, rowIxPerm) 1-                mapM (divideCellBy e rowIxPerm) [(i+1)..colN]-            divideCellBy e rowIxPerm colIx =-                do-                e2 <- unsafeMatrixRead mIarr rowN (colIx, rowIxPerm)-                unsafeMatrixWrite mIarr rowN (colIx, rowIxPerm) (e2/e)-                -            eliminateColumn =-                do-                iRowPerm <- unsafePermRead perm i-                mapM (eliminateColumnRow iRowPerm) $ [1..(i-1)] ++ [(i+1)..rowN]-            eliminateColumnRow iRowPerm rowIx =-                do-                rowIxPerm <- unsafePermRead perm rowIx-                c <- unsafeMatrixRead mIarr rowN (i, rowIxPerm) -- remember old element for scaling i'th row-                unsafeMatrixWrite mIarr rowN (i,rowIxPerm) 0 -- at column i we set 0-                mapM (eliminateColumnRowColumn iRowPerm rowIxPerm c) [(i+1)..colN]-            eliminateColumnRowColumn iRowPerm rowIxPerm c colIx =-                do-                ei <- unsafeMatrixRead mIarr rowN (colIx, iRowPerm) -- at i'th row-                er <- unsafeMatrixRead mIarr rowN (colIx, rowIxPerm) -- at current row-                unsafeMatrixWrite mIarr rowN (colIx, rowIxPerm) (er - c * ei) -- eliminate by i'th row-               - -swap ::-    Int ->-    Int ->-    (STAr.STUArray s Int Int) ->-    ST s ()-swap i1 i2 perm =-    do-    a1 <- unsafePermRead perm i1-    a2 <- unsafePermRead perm i2-    unsafePermWrite perm i1 a2-    unsafePermWrite perm i2 a1-            --unsafePermWrite permArr i e =-    do-    BAr.unsafeWrite permArr (i - 1) e-                -unsafePermRead permArr i =-    do-    BAr.unsafeRead permArr (i - 1)-                --addIdentity ::-    (STMatrix s IRA) ->-    ST s (STMatrix s IRA)-addIdentity (Matrix colN rowN marr) =-    do---    (_, (colN,rowN)) <- MAr.getBounds marr-    mElems <- MAr.getElems marr-    mIarr <- MAr.newListArray ((1,1),(colN+rowN,rowN)) $ mElems ++ (idElems rowN)-    return $ Matrix (colN + rowN) rowN mIarr-    where-    idElems m =-        1 : (concat $ replicate (m-1) $ (replicate m 0) ++ [1])---data Matrix marr el =-    Matrix-    {-        mxRowN :: Int,-        mxColN :: Int,-        mxRows :: marr (ColIx,RowIx) el-    }--type ColIx = Int -type RowIx = Int --type IMatrix el = -    Matrix Array el-    -type STMatrix s el =-    Matrix (STArray s) el-    -instance -    (IAr.IArray marr el,-- IAr.IArray marr (marr Int el), -     Show el) => -    Show (Matrix marr el)-    where-    show (Matrix colN rowN rows) =-        "\nMatrix:\n" ++ -        (concat $ map showCol [1..colN])-        where---        (_,(colN,rowN)) = IAr.bounds rows-        showCol colIx =-            unlines $-                map showCell [(colIx, rowIx) | rowIx <- [1..rowN]] -        showCell ix@(colIx, rowIx) =-            (show ix) ++-            (replicate colIx '.') ++  -            (show $ (IAr.!) rows ix)-    -getMatrixPrecision (Matrix _ _ marr) =-    do-    elems <- MAr.getElems marr-    return $ foldl1 min $ map AERN.getPrecision elems--unsafeMatrixWrite marr rowN (i,j) e =-    do-    BAr.unsafeWrite marr (rowN*(i-1) + j-1) e---    MAr.writeArray marr (i,j) e--unsafeMatrixRead marr rowN (i,j) =-    do-    BAr.unsafeRead marr (rowN*(i-1) + j-1)---    MAr.readArray marr (i,j)-    -permuteRowsDropCols ::-    (STAr.STUArray s Int Int) ->-    Int {-^ drop this many first columns -} ->-    (STMatrix s IRA) ->-    ST s (STMatrix s IRA)-permuteRowsDropCols perm dropN (Matrix colN rowN marr) =-    do---    (_, (colN,rowN)) <- MAr.getBounds marr-    (_, permN) <- MAr.getBounds perm    -    rarr <- MAr.newArray ((1,1),(colN - dropN, permN)) 0-    mapM (copyElem marr rarr rowN) [(colIx, rowIx) | colIx <- [1..colN - dropN], rowIx <- [1..permN]]-    return (Matrix (colN - dropN) permN rarr)-    where-    copyElem marr rarr rowN (colIx, rowIx) =-        do-        permRowIx <- unsafePermRead perm rowIx-        e <- unsafeMatrixRead marr rowN (colIx + dropN, permRowIx)-        unsafeMatrixWrite rarr rowN (colIx, rowIx) e-        -    -addM m1 m2 -    | mxColN m1 == mxColN m2 && mxRowN m1 == mxRowN m2 =-        do-        marr <- MAr.newArray ((1,1),(colN, rowN)) 0-        mapM (addCell marr) [(c,r) | c <- [1..colN], r <- [1..rowN]]-        return (Matrix colN rowN marr)   -    | otherwise =-        error "Matrix: addM mismatch"-    where-    colN = mxColN m1-    rowN = mxRowN m1-    marr1 = mxRows m1-    marr2 = mxRows m2-    addCell marr (colIx, rowIx) =-        do-        elem1 <- unsafeMatrixRead marr1 rowN (colIx, rowIx)-        elem2 <- unsafeMatrixRead marr2 rowN (colIx, rowIx)-        unsafeMatrixWrite marr rowN (colIx, rowIx) (elem1 + elem2)--multM m1 m2 -    | colN1 == rowN2 =-        do-        marr <- MAr.newArray ((1,1),(colN, rowN)) 0-        mapM (multCell marr) [(c,r) | c <- [1..colN], r <- [1..rowN]]-        return (Matrix colN rowN marr)   -    | otherwise =-        error "Matrix: multM mismatch"-    where-    colN1 = mxColN m1-    rowN1 = mxRowN m1-    colN2 = mxColN m2-    rowN2 = mxRowN m2-    colN = colN2-    rowN = rowN1-    marr1 = mxRows m1-    marr2 = mxRows m2-    multCell marr (colIx, rowIx) =-        do-        elems1 <- mapM (getCell1 rowIx) [1..colN1]-        elems2 <- mapM (getCell2 colIx) [1..rowN2]-        unsafeMatrixWrite marr rowN (colIx, rowIx) (sum $ zipWith (*) elems1 elems2)-    getCell1 rowIx colIx =-        do-        unsafeMatrixRead marr1 rowN1 (colIx, rowIx)-    getCell2 rowIx colIx =-        do-        unsafeMatrixRead marr2 rowN2 (colIx, rowIx)-        

@@ -1,43 +0,0 @@-{-# LANGUAGE CPP #-}-{-# LANGUAGE UndecidableInstances #-}-{-# LANGUAGE TypeOperators #-}-{-# LANGUAGE DeriveDataTypeable #-}-module Main--where--import qualified Data.Number.ER.Real as AERN-import Data.Number.ER.Real (ConvergRealSeq(..), convertFuncRA2Seq)-import Data.Number.ER.BasicTypes-import Data.Number.ER.Misc--import Data.Maybe--#ifdef USE_MPFR-type B = AERN.BMPFR -- use MPFR floats-#else-type B = AERN.BAP -- use pure Haskell floats---type B = AERN.BMAP -- use combination of double and pure Haskell floats-#endif-type RA = AERN.RA B-type IRA = AERN.IRA B---decimalPrec = 100---decimalPrec = 1000-binaryPrec =-    fromInteger $ toInteger $-    snd $ AERN.integerBounds $-        (fromInteger decimalPrec :: RA) * (AERN.log 100 10)/(AERN.log 100 2)--main =-    do-    AERN.initialiseBaseArithmetic (0 :: RA)-    putStrLn $ -        show decimalPrec -        ++ " decimal digits of pi = \n" -        ++ (AERN.showConvergRealSeqAuto binaryPrec pi)-    where-    pi :: ConvergRealSeq RA-    pi = ConvergRealSeq AERN.pi-

@@ -0,0 +1,43 @@+{-# LANGUAGE CPP #-}+{-| +    Module      :  Main+    Description :  laucher for approximated exact real arithmetic tests+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mik@konecny.aow.cz+    Stability   :  experimental+    Portability :  portable++    An executable for easy automated launch of tests +    of approximated exact real arithmetic.+-}+module Main where++import qualified Data.Number.ER.Real.Approx as RA+import Data.Number.ER.Real.Approx.Tests.Run +import Data.Number.ER.Real.DefaultRepr++--import Data.Number.ER.Real.Approx.Tests.Properties+--import Data.Number.ER.Real.Approx.Tests.Generate++main =+    do+    runRATests "interval-double" sampleRABM (RA.initialiseBaseArithmetic sampleRABM)+    runRATests "interval-haskell" sampleRABM (RA.initialiseBaseArithmetic sampleRABAP)+    runRATests "interval-haskell-double" sampleRABM (RA.initialiseBaseArithmetic sampleRABMAP)+#ifdef USE_MPFR+    runRATests "interval-mpfr" sampleRABM (RA.initialiseBaseArithmetic sampleRABMPFR)+#endif++sampleRABM :: RA BM+sampleRABAP :: RA BAP+sampleRABMAP :: RA BMAP+sampleRABM = 0+sampleRABAP = 0+sampleRABMAP = 0++#ifdef USE_MPFR+sampleRABMPFR :: RA BMPFR+sampleRABMPFR = 0+#endif