AERN-Real 0.10.0.1 → 0.10.0.2
raw patch · 8 files changed
+581/−580 lines, 8 filesPVP ok
version bump matches the API change (PVP)
API changes (from Hackage documentation)
Files
- AERN-Real.cabal +3/−3
- ChangeLog +1/−0
- demos/Demo.hs +0/−149
- demos/Matrix.hs +0/−385
- demos/Pi.hs +0/−43
- examples/Demo.hs +149/−0
- examples/Matrix.hs +385/−0
- examples/Pi.hs +43/−0
AERN-Real.cabal view
@@ -1,5 +1,5 @@ Name: AERN-Real-Version: 0.10.0.1+Version: 0.10.0.2 Cabal-Version: >= 1.2 Build-Type: Simple License: BSD3@@ -26,13 +26,13 @@ . For an architectural overview, see module "Data.Number.ER.Real". .- Simple examples of usage can be found in folder @demos@.+ Simple examples of usage can be found in folder @examples@. . 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+ examples/Demo.hs examples/Pi.hs examples/Matrix.hs tests/RunERIntervalTests.hs ChangeLog
ChangeLog view
@@ -1,3 +1,4 @@+0.10.0.2: 29 July 2009: renamed "demos" folder to "examples" 0.10.0.1: 28 July 2009: fixed a few bugs in meta-data 0.10.0: 28 July 2009 * switching to beta status
− demos/Demo.hs
@@ -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.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)
− demos/Matrix.hs
@@ -1,385 +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.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)-
− demos/Pi.hs
@@ -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-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-
+ examples/Demo.hs view
@@ -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)
+ examples/Matrix.hs view
@@ -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)+
+ examples/Pi.hs view
@@ -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+