AERN-Real 0.9.6.1 → 0.9.7

raw patch · 17 files changed

+649/−23 lines, 17 filesdep +hmpfrPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: hmpfr

API changes (from Hackage documentation)

Files

• AERN-Real.cabal +22/−8
• ChangeLog +4/−0
• src/Data/Number/ER/Misc.hs +3/−0
• src/Data/Number/ER/Real.hs +2/−0
• src/Data/Number/ER/Real/Approx.hs +3/−1
• src/Data/Number/ER/Real/Approx/Interval.hs +4/−1
• src/Data/Number/ER/Real/Base.hs +5/−0
• src/Data/Number/ER/Real/Base/CombinedMachineAP.hs +5/−0
• src/Data/Number/ER/Real/Base/Float.hs +1/−0
• src/Data/Number/ER/Real/Base/MPFR.hs +78/−0
• src/Data/Number/ER/Real/Base/MachineDouble.hs +15/−7
• src/Data/Number/ER/Real/Base/Rational.hs +1/−0
• src/Data/Number/ER/Real/DefaultRepr.hs +23/−4
• tests/Demo.hs +13/−2
• tests/ISin3.hs +47/−0
• tests/Matrix.hs +381/−0
• tests/Pi.hs +42/−0

@@ -1,5 +1,5 @@ Name:           AERN-Real-Version:        0.9.6.1+Version:        0.9.7 Cabal-Version:  >= 1.2 Build-Type:     Simple License:        BSD3@@ -13,7 +13,7 @@ Tested-with:    GHC ==6.8.2 Description:     Datatypes and abstractions for approximating exact real numbers-    and a basic arithmetic over such approximations.  +    and a basic arithmetic over such approximations.     The main datatype is interval with arbitrary precision endpoints     supported by safely rounding field and elementary operations     whose precision can be increased arbitrarily, so that they@@ -27,19 +27,32 @@     .     Simple examples of usage can be found in module @Demo.hs@ in folder @tests@. Extra-source-files:-    ChangeLog tests/Demo.hs+    ChangeLog tests/Demo.hs tests/Pi.hs tests/Matrix.hs tests/ISin3.hs  Flag containers-in-base     Default: False +Flag use-hmpfr+    Default: False+ Library   hs-source-dirs:  src   if flag(containers-in-base)-    Build-Depends:-      base < 3, binary >= 0.4+    if flag(use-hmpfr)+      Build-Depends:+        base < 3, binary >= 0.4, hmpfr == 0.1.1+      cpp-options: -DUSE_MPFR+    else+      Build-Depends:+        base < 3, binary >= 0.4   else-    Build-Depends:-      base >= 3, containers, binary >= 0.4+    if flag(use-hmpfr)+      Build-Depends:+        base >= 3, containers, binary >= 0.4, hmpfr == 0.1.1+      cpp-options: -DUSE_MPFR+    else+      Build-Depends:+        base >= 3, containers, binary >= 0.4   Exposed-modules:     Data.Number.ER,     Data.Number.ER.Real,@@ -48,6 +61,7 @@     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,@@ -64,5 +78,5 @@     Data.Number.ER.BasicTypes,     Data.Number.ER.Misc,     Data.Number.ER.ExtendedInteger-  Extensions: DeriveDataTypeable, ForeignFunctionInterface, ScopedTypeVariables, MultiParamTypeClasses, FunctionalDependencies, TypeSynonymInstances, FlexibleInstances+  Extensions: CPP, DeriveDataTypeable, ForeignFunctionInterface, ScopedTypeVariables, MultiParamTypeClasses, FunctionalDependencies, TypeSynonymInstances, FlexibleInstances   

@@ -1,3 +1,7 @@+0.9.7: 30 September 2008+    * made it easier to switch among various bases (double, mpfr, pure haskell floats...)+    * added MPFR backend via hmpfr (cabal install -f "use-hmpfr")+ 0.9.6.1: 7 August 2008     * revamped package description to make it much shorter and linked it       to the main module

@@ -21,6 +21,9 @@         putStrLn $ "unsafe: " ++ msg         return val +unsafePrintReturn msg a =+    unsafePrint (msg ++ show a) a+ {-|     Compose as when defining the lexicographical ordering. -}

@@ -35,6 +35,8 @@     precisions.  Rational numbers can be used as one of the alternatives.     Augustsson's Data.Number.BigFloat can be easily wrapped as an instance     of ERRealBase except that it uses a different method to control precision.+    Optionally, one can also have MPFR floating point numbers via package+    hmpfr if compiled with USE_MPFR.          ERIntApprox is implemented via outwards-rounded arbitrary precision interval arithmetic.       Any instance of ERRealBase can be used for the endpoints of the intervals.

@@ -51,7 +51,9 @@    A type whose elements represent sets that can be used    to approximate a single extended real number with arbitrary precision. -}-class (Fractional ra) => ERApprox ra where+class (Fractional ra) => ERApprox ra +	where+    initialiseBaseArithmetic :: ra -> IO ()     getPrecision :: ra -> Precision      {-^              Precision is a measure of the set size.

@@ -1,4 +1,5 @@ {-# LANGUAGE DeriveDataTypeable   #-}+{-# LANGUAGE ScopedTypeVariables  #-} {-|     Module      :  Data.Number.ER.Real.Approx.Interval     Description :  safe interval arithmetic@@ -104,7 +105,7 @@     (B.ERRealBase b) =>      ERInterval b -> EI.ExtendedInteger erintvPrecision (ERInterval l r) =-    - (B.getApproxBinaryLog $ (r - l)/(1 + abs(r) + abs(l)))+    -1 - (B.getApproxBinaryLog $ (r - l)) -- /(1 + abs r + abs l)) erintvPrecision ERIntervalEmpty = EI.PlusInfinity erintvPrecision ERIntervalAny = EI.MinusInfinity @@ -398,6 +399,8 @@     i1 / ERIntervalEmpty = ERIntervalEmpty              instance (B.ERRealBase b, RealFrac b) => RA.ERApprox (ERInterval b) where+    initialiseBaseArithmetic _ =+        B.initialiseBaseArithmetic (0 :: b)     getPrecision i = erintvPrecision i     getGranularity i = erintvGranularity i     {- setMinGranularity -}

@@ -32,6 +32,11 @@ -} class (Fractional rb, Ord rb) => ERRealBase rb      where+    typeName :: rb -> String+    initialiseBaseArithmetic :: rb -> IO ()+    initialiseBaseArithmetic x = +		do+		putStrLn $ "Base arithmetic: " ++ typeName x     defaultGranularity :: rb -> Granularity     getApproxBinaryLog :: rb -> EI.ExtendedInteger     getGranularity :: rb -> Granularity

@@ -190,6 +190,11 @@          instance (B.ERRealBase b) => B.ERRealBase (ERMachineAP b)     where+    typeName _ = "auto switching double and " ++ (B.typeName (0::b))+    initialiseBaseArithmetic x = +		do+		putStr $ "Base arithmetic:" ++ B.typeName x ++ "; "+		initMachineDouble     defaultGranularity _ = (B.defaultGranularity (0 :: b))     getApproxBinaryLog =          op1ERMachineAP doubleBinaryLog B.getApproxBinaryLog

@@ -482,6 +482,7 @@      instance B.ERRealBase ERFloat     where+    typeName _ = "ERFloat (pure Haskell implementation)"     defaultGranularity _ = 10     getApproxBinaryLog = apFloatExponent     getGranularity = apfltGran

@@ -0,0 +1,78 @@+{-# LANGUAGE CPP #-}+-- #define USE_MPFR+{-|+    Module      :  Data.Number.ER.Real.Base.MPFR+    Description :  enabling MPFR dyadics as interval endpoints+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mik@konecny.aow.cz+    Stability   :  experimental+    Portability :  non-portable (requires fenv.h)++    Make Ales Bizjak's Haskell interface to MPFR an instance of +    'B.ERRealBase'.+     +    If compiled without USE_MPFR, this module is empty.+-}+module Data.Number.ER.Real.Base.MPFR+(+#ifdef USE_MPFR+    MPFR+#endif+)+where++import qualified Data.Number.ER.Real.Base as B+import qualified Data.Number.ER.ExtendedInteger as EI+import Data.Number.ER.Misc++import Data.Binary++#ifdef USE_MPFR+import qualified Data.Number.MPFR.Up as M++type MPFR = M.MPFR++instance Binary M.MPFR+    where+    get = error "Data.Number.Dyadic: Binary not implemented yet"+    put = error "Data.Number.Dyadic: Binary not implemented yet"+++instance B.ERRealBase M.MPFR+    where+    typeName _ = "MPFR float"+    defaultGranularity _ = 30+    getApproxBinaryLog d +        | d < 0 =+            error $ "ER.Real.Base.MPFR: getApproxBinaryLog: negative argument " ++ show d +        | d == 0 = EI.MinusInfinity +        | d >= 1 =+            fromInteger $ intLogUp 2 $ ceiling d+        | d < 1 =+            negate $ fromInteger $ intLogUp 2 $ ceiling $ recip d+    getGranularity = mPrec2gran . M.getPrec+    setMinGranularity g x +        | g > xGran = B.setGranularity g x+        | otherwise = x+        where+        xGran = B.getGranularity x  +    setGranularity g = M.set M.Up (gran2mPrec g)+    getMaxRounding _ = +        error "ER.Real.Base.MPFR: getMaxRounding undefined"+    isERNaN = M.isNaN+    erNaN = 0/0+    isPlusInfinity x = M.isInfinite x && x > 0+    plusInfinity = 1/0+    fromDouble = M.fromDouble M.Up 53+    toDouble = fromRational . toRational+    fromFloat = B.fromDouble . fromRational . toRational  +    toFloat = fromRational . toRational . B.toDouble+    showDiGrCmp numDigits _showGran _showComponents f = +        M.toStringExp (int2word numDigits) f+#endif+    +mPrec2gran = fromInteger . toInteger+gran2mPrec = fromInteger . toInteger+int2word = fromInteger . toInteger

@@ -56,20 +56,28 @@     currentRndMode <- fegetround     case currentRndMode == eFE_UPWARD of         True -> -            putStrLn "initMachineDouble: already rounding upwards" +            putStrLn "already rounding upwards"          False ->             do             fesetround eFE_UPWARD-            putStrLn "initMachineDouble: switched to upwards rounding" +            putStrLn "switching to upwards rounding"   instance B.ERRealBase Double     where+    typeName _ = "double"+    initialiseBaseArithmetic x = +		do+		putStr $ "Base arithmetic:" ++ B.typeName x ++ "; "+		initMachineDouble     defaultGranularity _ = 53-    getApproxBinaryLog f -        | f == 0 =-            EI.MinusInfinity-        | otherwise =-            intLogUp 2 (abs $ ceiling f)+    getApproxBinaryLog d +        | d < 0 =+            error $ "ER.Real.Base.MachineDouble: getApproxBinaryLog: negative argument " ++ show d +        | d == 0 = EI.MinusInfinity +        | d >= 1 =+            fromInteger $ intLogUp 2 $ ceiling d+        | d < 1 =+            negate $ fromInteger $ intLogUp 2 $ ceiling $ recip d     getGranularity _ = 53     setMinGranularity _ = id     setGranularity _ = id

@@ -206,6 +206,7 @@          instance B.ERRealBase ExtendedRational     where+    typeName _ = "extended rationals"     defaultGranularity _ = 10     getApproxBinaryLog (Finite r)         | r == 0 =

@@ -1,3 +1,5 @@+{-# LANGUAGE CPP #-}+-- #define USE_MPFR {-|     Module      :  Data.Number.ER.Real.DefaultRepr     Description :  concise names for default real representations@@ -18,8 +20,10 @@ -} module Data.Number.ER.Real.DefaultRepr (-    initMachineDouble,-    B, BM, BAP, BMAP, BR, +    B, BM, BAP, BMAP, BR,+#ifdef USE_MPFR+	BMPFR,+#endif          RA, IRA ) where@@ -35,6 +39,8 @@ import Data.Number.ER.Real.Base.MachineDouble import Data.Number.ER.Real.Base.CombinedMachineAP +import Data.Number.ER.Real.Base.MPFR+ type BAP = ERFloat  {-| @@ -45,6 +51,10 @@ -} type BM = Double +#ifdef USE_MPFR+type BMPFR = MPFR+#endif+ {-|         Use machine 'Double' while the granularity is up to its significant bit length         and when the granularity grows beyond that, use 'ERFloat'.@@ -61,18 +71,27 @@ {-|      the default base type -}++#ifdef USE_MPFR+type B = BMPFR+--type B = BMAP --type B = BAP --type B = BM+--type B = BR+#else type B = BMAP+--type B = BAP+--type B = BM --type B = BR+#endif  {-|      the default instance of 'Data.Number.ER.Real.Approx.ERApprox'  -}-type RA = ERInterval B+type RA b = ERInterval b  {-|      the default instance of 'Data.Number.ER.Real.Approx.ERIntApprox'  -}-type IRA = ERInterval B+type IRA b = ERInterval b 

@@ -1,3 +1,4 @@+{-# LANGUAGE CPP #-} {-|      Module      :  Main     Description :  simple examples of using AERN-Real@@ -13,8 +14,18 @@ module Main where  import qualified Data.Number.ER.Real as AERN-import Data.Number.ER.Real (RA, IRA, ConvergRealSeq(..), convertFuncRA2Seq)+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@@ -32,7 +43,7 @@  main =      do-    AERN.initMachineDouble+    AERN.initialiseBaseArithmetic (0 :: RA)     putStrLn "****************************"     putStrLn "Testing interval arithmetic:"     putStrLn "****************************"

@@ -0,0 +1,47 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE DeriveDataTypeable #-}+module Main++where++import qualified Data.Number.ER.Real as AERN+import qualified Data.Number.ER.RnToRm as AERNFunc+import Data.Number.ER.BasicTypes+import Data.Number.ER.Misc+import Data.Number.ER.RnToRm.TestingDefs++import Data.Maybe+import qualified Data.List as List+import qualified Data.Map as Map++#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++main =+    do+    AERN.initialiseBaseArithmetic (0 :: RA)+    putStrLn $ "integ(sin(sin(sin(x)))dx = " ++ show result+    putStrLn $ " precision = " ++ show (AERN.getPrecision result)+    where+    result = +        head $ +            AERNFunc.eval (AERNFunc.unary 1) $ +                AERNFunc.integrateUnary 0 +                    (sin3 120 0 42 160) +                    0 (0 AERN.\/ 1) [0]+    sin3 ix depth deg gran = +        AERN.sin ix $ +            AERN.sin ix $ +                AERN.sin ix $ +                    AERNFunc.bisectUnbisectDepth depth $ +                        AERN.setMinGranularity gran $ +                            AERNFunc.setMaxDegree deg fapwUPX0+

@@ -0,0 +1,381 @@+{-# 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 $ +--        "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)+        

@@ -0,0 +1,42 @@+{-# 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 = 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+